GTN模型详尽推导及相关计算流程

1、GTN模型GTN模型的屈服函数为：2**221332()cosh(1())2eqmqqffqff(1)其中，m为静水应力，eq为宏观的等效应力，为微观层面的等效应力。/3mii(2)devmij(3)1/23:2eqdevdev(4)*111()()1CCCCCFFCFfffqfffffffffffqff(5)2、状态变量的演化一般来说，屈服函数涉及应力张量的第一及第二不变量以及损伤参数，可以写为,,0meqiH(6)其中损伤参数由,1,2,,iHin表示。其为一系列的状态变量。张这些状态变量可以是等效塑性应变，空穴体积分数，微裂纹密度等等。接下来简要的探讨本章用到的两个状态变量。2.1、微观等效塑性应变空穴间的材料，即基体材料假定在微观层面服从一般的、各项同性的硬化率plhf，其中pl为微观等效塑性应变。hf定义为：()plf(7)损伤模型的第一个状态变量为微观等效塑性应变p。假设基体材料的等效塑性功率0(1)pf与塑性功率等效:p，即:(1)ppf(8)由上述假设，微观等效塑性应变定义为:(1)ppf(9)2.2、空穴体积分数用于损伤模型的第二个状态变量为空穴体积分数f。空穴体积分数的增长率可以写为growthnucleationfff(10)方程（10）表明变形增加的过程中，空穴的变化率部分是由于已有空穴的增长，部分是由于新的空穴的成核。空穴增长的公式表示为：(1):plgrowthffεI(11)空穴成核的定义常采用下述模型：()plnucleationNmfAB(12)其中A和B为以下模型定义的参数。Chu和Needleman假设对于成核存在平均等效塑性应变n。成核应变对于其平均值服从正态分布。因此，对于应变主导的成核A和B的方程可以写为2nNnnf1A=expandB=022pnss(13)其中ns是分布的标准平均差。此处我们考虑塑性应变导致的成核，因此Eq4.4变为N(1):ppffIA(14)3、弹塑性本构关系3.1、弹塑性本构方程在不含损失的情况下，经典的弹塑性本构模型把小变形情况下的应变率分解为弹性部分和塑性部分：epε=ε+ε(15)其中为总应变率，,ep分别为弹性及塑性应变率。对于各项同性线弹性，eeσ=C:ε(16)其中σ为柯西应力，eC为四阶张量22()3eijklikjlijklCGKG(17)K和G分别为体积模量和弹性剪切模量。通过对大多数金属材料的实验发现，应力率σ受弹性应变率的支配。换而言之，弹性响应不受塑性流动的影响。因此eeσ=C:ε(18)根据关联流动率塑性应变率可以写为pφεσ(19)其中为塑性乘子，为屈服函数。对方程（3）、（4）运用链式法则，变为：13meqφφφInσ(20)其中32deveqn(21)将（21）代入（20），流动法则可以写为1()3pppmeqIn(22)其中，由静水应力引起的体积应变率定义为*21231()3()sinh()2pmmmmqDqqffφ(23)由等效应力引起的等效塑性应变率定义为22()eqpeqeqeqDφ(24)根据式（22），两个标量,meqDD完全定义了塑性应变率p的变化。应力张量同样可以根据静水应力m及等效应力eq。因此23meqIn(25)采用径向返回算法确定,meqDD：在本构方程的积分过程中，应力及状态变量需要在每个应变增量步结束时被计算。通常非线性问题只能靠增量方法进行求解。因此，上述状态变量的积分方法需要同时考虑计算的精度及稳定性和计算的效率。通过（16），时间增量步结束时刻1nt的弹性方程可以写为11::nneeeepttCCt(26)其中试探弹性应力定义为：:()neeetCt(27)径向返回应力可以写为两部分：体积增量及偏增量，1:2neeptemeqCtKtDGtDIn(28)可以写为123ntmeqIn(29)其中emmmKtD(30)3eeqeqeqGtD(31)下标t+1为了简化故省略。3.2、率形式本构方程的数值积分为了计算（30）、（31）中的,meq，需要求解,meqDD。通过（23）、（24）消除塑性乘子得到0meqeqmDD(32)（32）为含有两个未知量,meqDD的相容条件。为了求解两个未知量，上述方程需要和屈服函数联合求解。采用联合流动法则假设的屈服面函数可以由塑性势函数得到12,,,0meqHH(33)因为,,meqiH仅为,meqDD的函数，选自,meqDD为(32)、(33)的基本未知量。可知1(),(),(,,,,)0ppppmmeqeqimeqmeqjmeqeqmfHHDD(34)且GTN屈服函数可以写为212(),(),(,,,,),,,ppppmmeqeqimeqmeqjmeqfHHHH(35)因为(35)中未知量为隐式形式出现，因此采用Newton迭代法进行求解。此方法为将上述非线性方程进行一阶Taylor展开进而线性化：111110kkkkkkmeqmeqffffDDDD(36)122220kkkkkkmeqmeqffffDDDD(37)1kkkmmmDDD(38)1kkkeqeqeqDDD(39)上述方程可以写为：11111221222kmkeqDRKKKKRD(40)两个内变量即状态变量为：1plH(41)2Hf(42)状态变量对时间的微分为：11(1)plplmmeqeqHhtf(43)221(1)(1)plplplmNmNHhfAfAht(44)对时间的积分可以简单的写为：1,11,1nnHHth(45)2,12,2nnHHth(46)假如从试探应力张量探测到有塑性出现，以下的迭代流程得出时间步结束时的应力。a)将迭代因子设置为k=0b)初始化。假如t=0，则0meqDD(47)120HH(48)否则，采用先前时间步结束步的值。c)计算ijK。d)计算iRe)对ijK求逆，并得到kmkeqDD。f)更新变量1kkmmmDDD(49)1kkeqeqeqDDD(50)11kekmmmKtD(51)113kekeqeqeqGtD(52)111111(1)kkkkmmeqeqkpDDHf(53)12(1)()(1)mmeqeqkmeqNDDHffDDAfIn(54)21exp22pNNNNNfAss(55)111,11,1kknnHHtH(56)112,12,2kknnHHtH(57)对于每个迭代步中得到的微观等效塑性应变11,1knH，即pl，采用局部的牛顿迭代法求出相应的微观层面的等效应力1ky，用于下一次迭代计算使用。本文采用的微观等效应力和微观等效塑性应变之间的关系式为00000,,N(58)其中为微观等效应变，00,分别为初始屈服应力和应变，00E。注：这里我不太明白何时从微观层面上讲认为材料在微观上进入了塑性区，也就是这个计算如何起步的问题。式（58）中的是如何得来的？是否认为满足了GTN的屈服函数即默认为微观层面上已经产生了塑性变形，抑或是需要从（58）判断其处于弹性或是弹塑性的状态？我这个问题考虑的不是很清楚，为了简化模型，在计算GTN的屈服函数判断是否产生塑性变形时，我假设微观等效应力为初始屈服极限，求解出微观等效塑性应变后再利用局部牛顿迭代求出新的相对应的微观等效应力替代原来的等效应力，不知道这样的假设是否合理？g)检查收敛情况。假如11,ff同时满足小于810则认为收敛；否则不收敛，k=k+1，转到第三步重新开始计算，反复迭代至收敛为止。假如收敛则计算应力如下：1(2)nemeqtKDIGDn(59)以下内容为迭代中需要用到的系数：222111122221(,)meqimimeqeqmeqimieqimmimfDDHKDtKDHDHDtKHD(60)222112212221(,)33meqimieqmeqeqieqieqimeqmieqfDDHKDGtDHDHDGtHD(61)22211m(,)meqiimimfDDHKtKDHD(62)22221eq(,)3meqiieqieqfDDHKGtDHD(63)11(,)meqmeqeqmRfDDDD(64)2212(,),,meqmeqRfDDHH(65)状态变量对时间的微分为1112(,,,,)meqmeqHhDDHH(66)2212(,,,,)meqmeqHhDDHH(67)由以上两式可以得到111111112112111111212()(()0)eqmmmmmeqmmmmmeqmmddHHhhhhHhHDhDDDHDHDhhhhHhHttKDHDHD(68)222222122112222212212()[()0]eqmmmmmeqmmmmmeqmmddHHhhhhHhHDhDDDHDHDhhhhHhHttKDHDHD(69)写成矩阵形式为1111112222221211mmmmmmHhhhhtKttDDHHthhHhhtttKHHDD(70)111111112112111111212()(0(3))eqmeqeqmeqeqeqeqeqeqmeqmmddHHhhhhHhHDhDDDHDHDhhhhHhHtGtDHDHD(71)222222122212222212212()(0(3))eqmeqeqmeqeqeqeqeqeqmeqmmddHHhhhhHhHDhDDDHDHDhhhhHhHtGtDHDHD(72)写成矩阵形式为111111222222123113eqeqeqeqeqeqHhhhhGtttDDHHthhHhhttGtHHDD