DSGZ模型

最新推荐文章于 2023-07-17 16:19:03 发布

battlestar

最新推荐文章于 2023-07-17 16:19:03 发布

阅读量1.7k

点赞数

分类专栏： Python学习 CAE学习

本文链接：https://blog.csdn.net/battlestar/article/details/105200661

版权

Python学习同时被 2 个专栏收录

22 篇文章 5 订阅

订阅专栏

CAE学习

9 篇文章 1 订阅

订阅专栏

[学习]DSGZ模型 ¹

一种基于唯象理论的材料本构模型，包含了四种模型（Johnson-Cook模型，G’Shell-Jonas模型，Matsuoka模型和Brooks模型）。该模型是关于应变，应变率和温度的函数，可以描述玻璃化和半结晶化聚合物的形变特性，尤其是内在应变软化和后续的硬化阶段。该模型通过五个应力应变点获得本构模型的八个参数.
三个应力-应变点在同一曲线上，分别为上屈服点，下屈服点（软化点），硬化点；第四个点为不同应变率下，同一温度下的上屈服点；第五个点为同一应变率下，不同温度下的上屈服点。

基本模型

$\sigma(\varepsilon,\dot {\varepsilon},T) = K \{ f(\varepsilon)+[\frac {\varepsilon\times e^{(1-\frac {\varepsilon} {C_3 \times h(\dot \varepsilon,T)})}} {C_3\times h(\dot \varepsilon,T)}-f(\varepsilon)] \times e^{[ln(g(\dot \varepsilon,T)-C_4)] \times \varepsilon} \} \times h(\dot \varepsilon,T)$

这里，
$f(\varepsilon)=(e^{-C_1 \times \varepsilon}+\epsilon^{C_2})(1-e^{-\alpha \times \varepsilon})$
$h(\dot \varepsilon,T)=(\dot \varepsilon)^m e^{\frac a T}$

参数说明

$C_1$ 和 $C_2$
$C_1$ 和 $C_2$ 确定基准应力应变曲线的走势，由方程求解获得：

$\left\{ \begin{array} {lr} \frac {e^{-C_1 \times \varepsilon_1}+\varepsilon^{C_2}} {e^{-C_1 \times \varepsilon_2}+\varepsilon^{C_2}}=\frac {\sigma_1}{\sigma_2} &\\ \frac {e^{-C_1 \times \varepsilon_1}+\varepsilon^{C_2}} {e^{-C_1 \times \varepsilon_3}+\varepsilon^{C_2}}=\frac {\sigma_1}{\sigma_3} \end{array} \right.$

K
比例系数。
m
描述同意温度下，两个不同应变率曲线之间的差异。
$(\frac {\dot \varepsilon_1} {\dot \varepsilon_2})^m=\frac {\sigma_1} {\sigma_2}$
a
描述同一应变率下，两个不同温度曲线之间的差异。
$\frac {ln(\frac {\sigma_1} {\sigma_2})} {\frac {1} {T_1} - \frac {1} {T_2}}$
$C_3$
描述不同应变利率和温度下上屈服点的差异。
注意： $C_3 \times h(\dot \varepsilon,T)=\varepsilon_y$ , $\varepsilon_y$ 和 $\sigma_y$ 为软化前的最大应边和应力，即材料的上屈服点。
$\alpha$ 和 $C_4$
两个参数都是用来拟合应力应变曲线的初始阶段（硬化前）。
$1-e^{-\alpha \times\varepsilon}=0.97$ ， $\varepsilon$ 应选区软化过程的终点（下屈服点）。
$C_4$ 针对不同聚合物材料选取不同的参数。
玻璃化聚合物：
$C_4 = 7+ln(\dot \varepsilon^m \times e^{\frac \alpha T})$
半结晶化聚合物：
$C_4 = 200+ln(\dot \varepsilon^m \times e^{\frac \alpha T})$

参数的影响规律

为了便于对模型个部分的理解：
令
$j=\frac {\varepsilon\times e^{(1-\frac {\varepsilon} {C_3 \times h(\dot \varepsilon,T)})}} {C_3\times h(\dot \varepsilon,T)}$

$l=e^{[ln(g(\dot \varepsilon,T)-C_4)] \times \varepsilon}$
绘制出这两个函数随应变的变化趋势。

在这里插入图片描述
$j$ 和 $l$ 两个函数都是0~1的函数。
$l$ 函数为随 $\epsilon$ 衰减的函数，其影响范围为应力应变曲线的初始阶段。
$j$ 函数为应力应变的前期函数，其峰值表明上屈服点的位置以及软化的速度。

测试参数

测试参数：

strain = [0.1,0.3,1.0];
stress =[112.,95.,160.,105.,80.]
strainrate = [0.001,0.0005]
T = [296.,323]
K = 4.5

	$c_1$	$c_2$	$m$	$a$	$K$	$c_3$	$c_4$	$\alpha$
论文	1.91	1.49	0.064	1191	3.9	0.0029	11	11.7
测试	1.35	2.09	0.093	1191.5	4.5	0.003398	10.38	11.69

测试结果对比I 试验及论文参数的对比

代码（Python）

# -*- coding: utf-8 -*-
"""
Author:Yujin.Wang
Date:2019.04.07
The lib. is used to translate the F-D curve to MAT24 effective plastic strain - true stress.
Usage:
return TrueStrain,TrueStress,[PeakStrain,PeakStress]
################################################################################
Parameters:
    A - float, Specimen Crosssection area
	l0 - float, Gage lengrh
    FDfile - Filename, Force-Disp curve
    strain_rate - list, Strain rate
    curvenum - int, Curve number in .key file
    ratio - float, Display ratio
#################################User Define###################################
    Ymould_user - float, Young modulus
    yieldpoint - float, Yield point
    pointnum - int, the number of output points
    strainturnlist - list, Turn strain 
    scaleturnlist - list,  Scale of the turn strain 
    curvescalelist - list, Curve scale equalent to SFO
    alignstrain - float, align the strain to user defined value(>max(strainturnlist))
    extend - float, Extend strain to a user defined value
    scaleextend - float, Scale of the extend strain 
"""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import interpolate,optimize
import scipy.signal as signal
from sympy import symbols,nsolve

def f_strain(strain,c1,c2,a):
    return (np.power(np.e,-c1*strain) + np.power(strain,c2))*(1-np.power(np.e,-a*strain))

def h_strainrate_T(strainrate,m,T,alpha):
    return np.power(strainrate,m)*np.power(np.e,alpha/T)

def DGSZ(strain,strainrate,T,K,c1,c2,c3,c4,a,m,alpha):
    f = f_strain(strain,c1,c2,a)
    h = h_strainrate_T(strainrate,m,T,alpha)
    j = (strain*np.power(np.e,(1.-strain/(c3*h)))/(c3*h))
    l =  np.power(np.e,(np.log(h)-c4)*strain)
    plt.figure(2)
    plt.plot(strain,f)
    plt.title('f')
    plt.figure(3)
    plt.title('j')
    plt.plot(strain,j)
    plt.figure(4)
    plt.title('j-f')
    plt.plot(strain,j-f)
    plt.figure(5)
    plt.title('l')
    plt.plot(strain,l)
    plt.figure(6)
    plt.title('(j-f)*l')
    plt.plot(strain,((j-f)*l))
    plt.figure(7)
    plt.title('(f+(j-f)*l)')
    plt.plot(strain,(f+(j-f)*l))
    plt.figure(8)
    plt.plot(strain,K*(f+(j-f)*l)*h)
    plt.figure(1)
    value = K* ( f + (j-f)*l )*h
    return value

def DGSZ_c1c2(strain,stress):
    c1,c2 = symbols('c1,c2')
    equations = [(np.power(np.e,-c1*strain[0]) + np.power(strain[0],c2))     \
                 /(np.power(np.e,-c1*strain[1]) + np.power(strain[1],c2)) - stress[0]/stress[1], \
                 (np.power(np.e,-c1*strain[2]) + np.power(strain[2],c2))     \
                 /(np.power(np.e,-c1*strain[1]) + np.power(strain[1],c2)) - stress[2]/stress[1]]
    result = nsolve(equations,[c1,c2],[1,1],verify=True,prec=150)
    return result
    
def DGSZ_m(strainrate,stress):
    return np.log(stress[0]/stress[1]) / np.log(strainrate[0]/strainrate[1])

def DGSZ_alpha(T,stress):
    return np.log(stress[0]/stress[1])/(1./T[0]-1./T[1])

def DGSZ_K(stress,strain,strainrate,m,T,alpha,c1,c2):
    h = h_strainrate_T(strainrate,m,T,alpha)
    return stress/(h*np.power(np.e,-c1*strain)+np.power(strain,-c2))

def DGSZ_c3(strain,strainrate,m,T,alpha):
    h = h_strainrate_T(strainrate,m,T,alpha)
    print (strain)
    return strain/h
    
def DGSZ_c4(strainrate,m,alpha,T):
    return 7 + np.log(np.power(strainrate,m)*np.power(np.e,alpha/T))

def DSGZ_a(strain):
    return -np.log(0.03)/strain

def CowperSymondsFunc(x,c,p):
    rate,ES0 = x
    return ES0*(1+np.power((rate/c),1./p))

def ConditionRemove(data,condition):
    a =[]
    isflag = 1
    string ="if "+ condition +":\n"
    string +="                     for k in range(len(data)):\n   \
                        data[k].pop(i+1)\n   \
                  isflag = 1\n"
    while isflag == 1:
        a.extend(data)
        len0 = len(a[0])
        try:
            for i in range(len(a[0])):
                exec(string)
        except:
            if len(data[0]) == len0:
                break
    return data
    
def MAT89LCSS(fdfile,A,l0):
    '''fdfile -- Force-Displacement Curve
       A -- Area of the Cross section of specimen
       l -- Length of the specimen (gage length or grip length or equalent length)
       Output : TrueStrain,TrueStress,[PeakStrain,PeakStress]
     '''
    data = pd.read_csv(fdfile)
    x = data.Disp
    y = data.Force
    EngStrain = x/l0
    EngStress = y/A
    EngStressFilt =  pd.Series(signal.medfilt(EngStress,201))
    PeakStress = max(EngStressFilt[:800])
    PeakStrain = EngStrain[EngStressFilt[:800].idxmax()]
    print ("PeakStress:%f\tPeakStrain:%f\t Elastic modulus:%f" %(PeakStress,PeakStrain,PeakStress/PeakStrain))
    plt.figure(1)
    plt.grid("on")
    plt.scatter(PeakStrain,PeakStress, marker='o')
    plt.title("Engineering Strain VS. Engineering Stress Curve")
    plt.xlabel("Engineering Strain[-]")
    plt.ylabel("Engineering Stress[GPa]")
    plt.legend(["0.01mm/ms","0.1mm/ms","1mm/ms","10mm/ms"])
    TrueStrain =  np.log(1+EngStrain)
    TrueStress = EngStressFilt*(1+EngStrain)
    plt.figure(2)
    plt.grid("on")
    plt.title("True Strain VS. True Stress Curve")
    plt.xlabel("True Strain[-]")
    plt.ylabel("True Stress[GPa]")
    plt.plot(TrueStrain,TrueStress)
    plt.legend(["0.01mm/ms","0.1mm/ms","1mm/ms","10mm/ms"])
    return TrueStrain,TrueStress,[PeakStrain,PeakStress]

def CurveKey1(x,y,strainrate_list,curvenum):
    # x,y = CowperSymondsCurve(FDfile[0],A,l0,Ymould_user,80,ratelist,x_p,y_p)
    fout = open("Curve.key",'w')
    fout.write('*KEYWORD\n')
    fout.write( "$%10s%10s%10s%10s%10s%10s%10s%10s\n" %('c1','c2','m','a','K','c3','c4','alpha'))
    fout.write( "$%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n" %(c1,c2,m,a,K,c3,c4,alpha))
    fout.write('$ Created: ' + time.strftime("%d.%m.%Y %H:%M:%S") + '\n')
    fout.write("*DEFINE_TABLE\n%d\n" %(curvenum))
    for i in strainrate_list:
        fout.write("%f\n" %(np.log(i)))
    for index,strainrate in enumerate(strainrate_list):
        res = [ x[index], y[index]]
        # condition = "(res[0][i]<res[0][i+1] or res[1][i]>res[1][i+1])"
        # x[index], y[index]= ConditionRemove(res,condition)
        curvenum += 1
        flag = 0 
        fout.write('*DEFINE_CURVE_TITLE\nRate %.5f\n' %(strainrate_list[index]))
        fout.write('$     LCID      SIDR       SFA       SFO      OFFA      OFFO    DATTYP\n')
        fout.write('      %d         0    1.0000    1.0000    0.0000    0.0000\n' %(curvenum))
        plt.figure(5)
#        plt.plot(x[index],y[index])
        plt.grid('on')
        for i in range(len(x[index])):
            if x[index][i] > 0:
                if flag == 0 :
                    fout.write("%f,%f\n" %(0,y[index][i]))
                    flag = 1
                if flag != 0 :
                    fout.write("%f,%f\n" %(x[index][i],y[index][i]))
    fout.write("*END\n")
    fout.close()
    plt.savefig('curve')
    return

if __name__ == '__main__':
    Ymould_user = 1.5 #定义弹性模量
    curvenum = 2350 #key文件中曲线的编号起始编号
    pointnum = 80#
    #############################################################
    strain = [0.1,0.3,1.0];
    stress =[112.,95.,160.,105.,80.]
    
    strainrate = [0.001,0.0005]
    T = [296.,323]
    c1,c2 = DGSZ_c1c2(strain[:3],stress[:3])
    m = DGSZ_m(strainrate,[stress[0],stress[3]])
    alpha = DGSZ_alpha(T,[stress[0],stress[4]])
#    K = DGSZ_K(stress[0],strain[0],strainrate[0],m,T[0],alpha,c1,c2)
    K = 4.5
    c3 = DGSZ_c3(strain[0],strainrate[0],m,T[0],alpha)
    c4 = DGSZ_c4(strainrate[0],m,alpha,T[0])
    a = DSGZ_a(strain[1])
    plt.scatter(strain[:3],stress[:3])
    res_stress = [];res_strain = []
    res_stress_89 = [];res_strain_89 = []
    print ("c1\tc2\tm\ta\tK\tc3\tc4\talpha")
    print ("%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f" %(c1,c2,m,alpha,K,c3,c4,a))
        #
    strain_curve = np.linspace(0.0001,1.,pointnum)
    strainrate_output = [0.0001,0.0005,0.001]
    for index,sr in enumerate(strainrate_output):
        stress = DGSZ(strain_curve,sr,T[0],K,c1,c2,c3,c4,a,m,alpha)
        plt.plot(strain_curve,stress)
        strain_temp_list = [];stress_temp_list = []
        for i in range(pointnum):
            try:
                if stress[i] == stress[i-1]:
                   stress[i] = stress[i]*2.
            except:
                pass
            strain_temp = strain_curve[i] - stress[i]/Ymould_user
            if strain_temp > 0.2:
                if strain_temp_list == []:
                   strain_temp_list.append(0)
                   stress_temp_list.append(stress[i]/3.)
                strain_temp_list.append(strain_curve[i])
                stress_temp_list.append(stress[i])
        res_strain.append(strain_temp_list)
        res_stress.append(stress_temp_list)
        res_strain_89.append(strain_curve)
        res_stress_89.append(stress)   
        plt.grid('on')

特定温度和应变率下的三维图形。

# -*- coding: utf-8 -*-
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import interpolate,optimize
import scipy.signal as signal
from sympy import symbols,nsolve
from mpl_toolkits.mplot3d import Axes3D

def f_strain(strain,c1,c2,a):
    return (np.power(np.e,-c1*strain) + np.power(strain,c2))*(1-np.power(np.e,-a*strain))

def h_strainrate_T(strainrate,m,T,alpha):
    return np.power(strainrate,m)*np.power(np.e,alpha/T)

def DGSZ(strain,strainrate,T,K,c1,c2,c3,c4,a,m,alpha):
    f = f_strain(strain,c1,c2,a)
    h = h_strainrate_T(strainrate,m,T,alpha)
    print ('h:\n',strainrate)
    print (h)
    j = (strain*np.power(np.e,(1.-strain/(c3*h)))/(c3*h))
    l =  np.power(np.e,(np.log(h)-c4)*strain)
    value = K* ( f + (j-f)*l )*h
    return value

def DGSZ_c1c2(strain,stress):
    c1,c2 = symbols('c1,c2')
    equations = [(np.power(np.e,-c1*strain[0]) + np.power(strain[0],c2))     \
                 /(np.power(np.e,-c1*strain[1]) + np.power(strain[1],c2)) - stress[0]/stress[1], \
                 (np.power(np.e,-c1*strain[2]) + np.power(strain[2],c2))     \
                 /(np.power(np.e,-c1*strain[1]) + np.power(strain[1],c2)) - stress[2]/stress[1]]
    result = nsolve(equations,[c1,c2],[1,1],verify=True,prec=150)
    return result
    
def DGSZ_m(strainrate,stress):
    return np.log(stress[0]/stress[1]) / np.log(strainrate[0]/strainrate[1])

def DGSZ_alpha(T,stress):
    return np.log(stress[0]/stress[1])/(1./T[0]-1./T[1])

def DGSZ_K(stress,strain,strainrate,m,T,alpha,c1,c2):
    h = h_strainrate_T(strainrate,m,T,alpha)
    return stress/(h*np.power(np.e,-c1*strain)+np.power(strain,-c2))

def DGSZ_c3(strain,strainrate,m,T,alpha):
    h = h_strainrate_T(strainrate,m,T,alpha)
    return strain/h
    
def DGSZ_c4(strainrate,m,alpha,T):
    return 7 + np.log(np.power(strainrate,m)*np.power(np.e,alpha/T))

def DSGZ_a(strain):
    return -np.log(0.03)/strain

    



if __name__ == '__main__':
    Ymould_user = 1.5 #定义弹性模量
    curvenum = 2350 #key文件中曲线的编号起始编号
    pointnum = 20#
    #############################################################
    strain = [0.1,0.3,1.0];
    stress =[112.,95.,160.,105.,80.]
    
    strainrate = [0.001,0.0005]
    T = [296.,323]
    c1,c2 = DGSZ_c1c2(strain[:3],stress[:3])
    c1 = np.float(c1)
    c2 = np.float(c2)
    m = DGSZ_m(strainrate,[stress[0],stress[3]])
    alpha = DGSZ_alpha(T,[stress[0],stress[4]])
#    K = DGSZ_K(stress[0],strain[0],strainrate[0],m,T[0],alpha,c1,c2)
    K = 4.5
    c3 = DGSZ_c3(strain[0],strainrate[0],m,T[0],alpha)
    c4 = DGSZ_c4(strainrate[0],m,alpha,T[0])
    a = DSGZ_a(strain[1])
    plt.scatter(strain[:3],stress[:3])
    res_stress = [];res_strain = []
    res_stress_89 = [];res_strain_89 = []
    print ("c1\tc2\tm\ta\tK\tc3\tc4\talpha")
    print ("%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f" %(c1,c2,m,alpha,K,c3,c4,a))
        #
    x = np.linspace(0.0001,1.,pointnum)
    y1 = np.linspace(0.1,10.,pointnum)
    y2 = np.linspace(238.,358.,pointnum)
    
    
    X,Y = np.meshgrid(x,y1)
    
    # strainrate,temperature = np.meshgrid(x,y)
    stress = DGSZ(X,Y,273.,K,c1,c2,c3,c4,a,m,alpha)    
    fig = plt.figure(1)
    ax = Axes3D(fig) 
    ax.plot_surface(X,Y, stress, rstride = 1, cstride = 1, cmap = plt.get_cmap('ocean'),alpha=1)
    ax.set_xlabel("Strain")
    ax.set_ylabel("Strainrate")
    ax.set_zlabel("Stress")    
    
    X,Y = np.meshgrid(x,y2)
    strainrate,temperature = np.meshgrid(x,y2)
    stress = DGSZ(X,0.1,Y,K,c1,c2,c3,c4,a,m,alpha)
    fig = plt.figure(2)
    ax = Axes3D(fig) 
    ax.plot_surface(X,Y, stress, rstride = 1, cstride = 1, cmap = plt.get_cmap('ocean'),alpha=1)
    ax.set_xlabel("Strain")
    ax.set_ylabel("Temperature")
    ax.set_zlabel("Stress")    
    plt.show()

Y. Duan, A.Sagal, R. Greif. A uniform Phenomenological Constitutive Model for Glassy and Semicrystalline Polymers. Polymer Engineering and Science,vol.41.2001.08 ↩︎

battlestar

关注

0
点赞
踩
4

收藏

觉得还不错? 一键收藏
0
评论
DSGZ模型

一个基于唯象理论的材料本构模型，包含了四种模型（Johnson-Cook模型，G’Shell-Jonas模型，Matsuoka模型和Brooks模型）。该模型是关于应变，应变率和温度的函数，可以描述玻璃化和半结晶化聚合物的形变特性，尤其是内在应变软化和后续的硬化阶段。该模型通过五个应力应变点获得本构模型的八个参数.三个应力-应变点在同一曲线上，分别为上屈服点，下屈服点（软化点），硬化点；第四个点为不同应变率下，同一温度下的上屈服点；第五个点为同一应变率下，不同温度下的上屈服点。
复制链接

扫一扫