材料非线性有限元法(共27页)

1、精选优质文档-倾情为你奉上第四章 材料非线性有限元法 以上三章分别研究了线性弹性有限元法，材料非线性本构方程和非线性方程组解法，本章就可以研究材料非线性有限元法了。 在材料非线性基本方程中，除第二章所述的本构方程外，与线性弹性一样，而非线性有限元法又归结为一系列线性弹性问题。因此，只要在第一章中改用第二章的本构方程，就可建立材料非线性有限元法的基本内容。§4-1 非线性弹性有限元法 第二章提到，非线性弹性本构方程与形变理论弹塑性本构方程在形式上相同，所以与第二章一样，这里也按塑性力学形变理论，研究非线性弹性有限元法，以便把二者统一起来。 1 非线性弹性基本方程 为了便于以后直接引用，

2、这里列出全量形式的非线性弹性（或形变理论弹塑性）基本方程，并用矩阵表示。 几何方程： (1.14) 本构方程：=D (2.13) 平衡方程：(在内) (1.20) 边界条件： （在A上） （1.22） （在A上） (1.23) 虚功方程： (1.28) 位能变分方程：=0 (1.31)其中 (1.32) (4.1) 2 非线性方程组的建立 由于虚功方程本身不涉及材料性质，所以第一章由虚功方程得到的单元平衡方程（1.48）式和总体平衡方程（1.109）式完全适用于非线性弹性（或形变理论弹塑性）问题。可见，只要把非线性弹性（或形变理论弹塑性）本构方程代入单元或总体的平衡方程，就可以建立非线性方程组

3、。 （1）割线刚度方程 仿照线性弹性有限元法，把（1.36）式代入（2.13）式后，再把（2.13）代入（1.48）式便得单元割线刚度方程，即 (4.2)其中单元割线刚度矩阵 (4.3)而割线本构矩阵，如（2.14）式所示。 仿照（1.113）式的推导，同样可得总体割线刚度方程即 (4.4)其中总体割线刚度矩阵 (4.5)而总体节点载荷P仍如（1.110）式所示。 由（4.5）式可知，总体割线刚度矩阵K取决于各单元的等效应变 ；又由（2.5）式可知，等效应变 是由应变计算出来的；再由（1.36）和（1.106）式可知，应变与总体节点位移U有关。可见，总体割线刚度矩阵K是总体节点位移U的函数，所

4、以总体割线刚度方程（4.4）式是一个非线性方程组。 必须指出，建立非线性方程组（4.4）式，只是为了说明非线性弹性（或形变理论弹塑性）有限元方程的非线性性质。实际求解时并不用（4.4）式。因为求解（4.4）式要用直接迭代法，而正如 3-2指出，直接迭代法不但计算量太大，而且常常不收敛。 （2）切线刚度矩阵 由 3-2-3-6可知，在求解非线性方程组时，除上述直接迭代法外，都要用到切线刚度矩阵（至少要用到初始切线刚度矩阵k和K）。为此，这里讨论一下建立非线性弹性（或形变理论弹塑性）有限元方程中的切线刚度矩阵问题。 由（1.109）和（3.11）式可知 (4.6)于是由（3.10）和（4.6）式可

5、得 (4.7)由于由（2.16），（1.36）和（1.106）式，并考虑到符号d和d分别是d和d，有 (4.8) (4.9) (4.10)所以把（4.8）-（4.10）式代入式（4.7）式便得总体割线刚度矩阵，即 (4.11)其中单元切线刚度矩阵 (4.12) （3）具有初应变理论或初应力的刚度方程 仿照线性弹性有限元法，把形式上相同的（3.101）式代入（2.13）式，并令=0或=0，再把（2.13）式代入（1.48）式便得单元刚度方程，即 (4.13)或 (4.14)其中单元刚度矩阵和初应变，初应力节点载荷，仍分别如（1.50）和（1.53）、（1.54）式所示。但要强调，这里k的含义是单

6、元初始切线刚度矩阵；中的初应变或中的初应力随迭代过程而变。 仿照线性弹性有限元法，同样可得总体刚度方程，即 (4.15)或 (4.16)其中总体刚度矩阵和总体初应变、初应力节点载荷、在形式上均与线性弹性有限元法相同。 3 等效应力、等效应变关系 由（4.11）-（4.16）式可知，要建立并求解非线性弹性（或形变理论弹塑性）有限元方程，关键是要具体知道材料的本构矩阵。而由（2.14）和（2.18）式可知，只要（2.15）和（2.19）式中的函数关系是已知的，那么本构矩阵就是显式的。 根据单一曲线假设，和的关系与单向拉伸时相同，即 (4.17)再考虑体积不可压缩条件（），则 (4.18)其中取决于

7、所采用的简化模型。 理想塑性（见图4-1）： (4.19)线性强化塑性（见图4-2）： (4.20)幂次强化塑性（见图4-3）： ， （4.21）4 迭代公式的具体化 由于非线性弹性（或形变理论弹塑性）有限元方程一般都写成全量形式，所以这里只相应的列出几种迭代类型解法的具体迭代公式。（1）Newton-Raphson法 由（1.36）、（1.106）和（2.5）式以及（3.17）和（3.18）式，有 （4.22） （4.23） （4.24） (4.25) （4.26） （4.27）（2）初应变迭代法 由（2.10）、（2.13）和（3.99）、（3.101）式可知 （4.28） （4.29）

8、（4.30）所以仿照Newton-Raphson法，并考虑到（3.109）和（3.110）式，有 （4.31） （4.32） （4.33） （4.34） （4.35）（3）初应力迭代法 由（2.10）、（2.13）和（3.118）、（3.120）式可知 （4.36） （4.37） （4.38）所以仿照Newton-Raphson法，并考虑到（3.126）和（3.127）式，有 （4.39） (4.40) (4.41) (4.42) (4.43)§4-2 非线性弹性手算例题为了熟悉非线性弹性(或形变理论弹塑性)有限元法及其非线性方程组的求解过程,这里以图4-4(a)所示的弹塑性拉压超静

9、定问题为例，用Newton-Raphson法、初应变迭代法、初应力迭代法进行手算。其中用Newton-Raphson法的求解作较详细的叙述，以便了解非线性弹性(或形变理论弹塑性)有限元分析的全过程，而用其他方法的求解只给出主要计算过程和计算结果。在该拉压超静定杆中划分的节点和单元如图4-4（a）中所示。单元和分别由线性强化材料和线性弹性材料制成，如图4-4（b）和（c）所示。两个单元的截面积均为，长度均为，弹性模量均为，单元的强化模量为。节点2所受集中载荷，其中为单元的屈服极限。 1 2 3 图4-4（a）1 用Newton-Raphson法求解（1）非线性有限元方程的形成 首先去掉两固定端约

10、束，用其约束反力和代替，所以 （a）由于 （b） 所以 （c） （d）设单元内任意一点位移为而 所以其中几何矩阵 （e）单元节点位移向量这里的单元都是单向应力状态，即 （f） （g）所以由（4.24）、（4.25）和（4.26）式可得 （h） （i）把（c）（g）式代入（h）和（i）式，并考虑到，有 （j） （h）把（a）、（j）和（k）式代入（3.17）式，并考虑到，有 （）最后由边界条件和，用第一章所述的消行降阶法便得约束处理后的非线性有限元方程，即 （m）（2）切线刚度迭代公式的建立 由于 所以无量纲迭代公式为 （n） （o） （p） （q） （r） （s）（3）节点位移和单元应力的求解

11、 按迭代公式（n）（s）式，求节点位移和单元应力的迭代过程和计算结果如表4-1所示： 表4-1n101.5001.5001.250-1.50022.7500.1661.6661.333-1.66633.00002 用初应变迭代法求解 若令（j）式中，则 （A）仿照用Newton-Raphson法的求解，由（3.103）式有 （B）约束处理后可得 （C）因为只有单元具有塑性，所以初应变节点载荷只有由单元提供。由于 所以 （D）由图4-4（b）可知因此 （E）于是由（n）、（C）（E）和（3.109）式可得无量纲迭代公式，即 （F） （G） (H) （I） （J） 按迭代公式（F）（J）式，求节点

12、位移和单元应力的迭代过程和计算结果如表4-2所示： 表4-2n1001.5001.500-1.50020.5000.2501.7501.250-1.75030.2500.1251.6251.375-1.62540.3750.1881.6881.313-1.68750.3130.1561.6561.343-1.65660.3430.1721.6721.329-1.67270.3290.1641.6641.335-1.66480.3350.1681.6681.333-1.66890.3330.1661.6661.333-1.666100.3333 用初应力迭代法求解 仿照用初应变迭代法求解，有 （

13、K）约束处理后可得 （L）由于只有单元具有塑性，所以初应力节点载荷只由单元提供。由于所以 （M）于是由（n）、（L）、（M）和（3.126）式可得无量纲迭代公式，即(N) （O） （P） （Q）按迭代公式（N）（Q）式，求节点位移和单元应力的迭代过程和计算结果如表4-3所示： 表4-3n101.5001.250-1.50020.1251.6251.312-1.62530.1561.6561.328-1.65640.1641.6641.332-1.66450.1601.6671.333-1.66760.166§4-3 弹塑性有限元法在§4-1中讨论了基于形变理论的弹塑性有限元

14、法。这里将讨论基于流动理论的弹塑性有限元法。为了跟踪加载历史求出位移、应变和应力的全量，基于流动理论弹塑性有限元方程只能取增量形式。同非线性弹性（或形变理论弹塑性）问题一样，第一章由虚功方程得到的平衡方程也完全适用于流动理论的弹塑性问题。所以由（1.109）、（3.11）和（3.59）式便得总体增量刚度方程，即 （4.44）其中总体切线刚度矩阵 （参考（4.11）式） （4.45）单元切线刚度矩阵 （参考（4.12）式） （4.46）（4.46）式中的切线弹塑性矩阵，按（2.63）式（正则屈服面）或（2.72）和（2.88）式（非正则屈服面）计算。几种常见材料模型的切线弹塑性矩阵已在§

15、;2-3 §2-5中给出显示形式。由于流动理论弹塑性有限元方程只能取增量形式，所以在第三章所述的解法中，只有增量类型的解法才能用来对它进行求解，其中常用的是Euler-Cauchy法、Euler一次迭代法、初应变增量法和初应力增量法。1 逐步求解公式的具体化 这里列出一求得的第I增量步的具体公式。（1）Euler-Cauchy法 对比（3.59）、（3.61）式和（4.44）（4.46）式，并仿照非线性弹塑性有限元法，Euler-Cauchy法逐步求解公式为： （4.47） （4.48） （4.49） （4.50）Euler一次迭代法：同样，仿照非线性弹塑性有限元法，Euler一次迭

16、代法逐步求解公式为： （4.51） （4.52） （4.53） （4.54） （4.55）（2）初应变增量法 仿照§4-1中的初应变迭代法，并考虑到（3.102）和（3.114）（3.116）式，有 （4.56） （4.57） （4.58） （4.59）对于Mises模型，若把塑性应变增量直接看作是初应变增量，并写成有限增量形式，则由（2.126）、（2.135）和（2.136）式可得： （4.60）初应变增量法：仿照§4-1中的初应力迭代法，并考虑到（3.121）和（3.135）（3.137）式，有： （4.61） （4.62） （4.63） （4.64）2 Mises模

17、型和Drucher-Prager模型 （1）单元弹、塑性状态的确定 由于同一增量步的不同单元可能处于不同的状态弹性或塑性，而不同状态下单元的应力增量应按不同公式计算，所以要计算（4.47）（4.64）式中的单元应力增量，必须参照§2-2中的加卸载准则，确定第I增量步单元的弹、塑性状态。为了确定第I增量步单元的状态，可先假定该单元在此增量步内不产生新的塑性变形（弹性、卸载或中性变载），于是有： （4.65）据此可以计算第i步终了时的加载函数值或屈服函数值。为了计算方便，可把加载函数写成另一种形式。由（2.102）、（2.148）和（2.165）式可知，Mises等向强化加载函数值为：

18、（4.66）Mises随动强化（线性）加载函数值为： （4.67）Drucher-Prager理想塑性加载函数（即屈服函数）值为： （4.68）参照§2-2中的加卸载准则，若则应力在加载面内（弹性或卸载）或在加载面上（中性变载），说明上述假定是正确的。我们把这种单元简称为弹性单元。若 则应力在加载面外，这当然是不可能的，说明上述假定是错误的。但这种情况下的单元只有两种可能状态：全塑性或弹塑性，而单元的确切状态又取决于第i-1步终了时（或第i步开始时）的单元状态。若第i-1步终了时单元是塑性的，即则第i步单元是全塑性的（有一种塑性状态进入另一种塑性状态，即完全加载状态），我们把这种单元

19、简称为塑性单元。若第i-1步终了时单元是弹性的，即则第i步单元是弹塑性的（有弹性状态进入塑性状态，即部分加载状态），我们把这种单元简称为弹塑性单元（或过度单元）。（2）单元应力增量的计算 对于弹性单元，显然有： （4.69）对于塑性单元，当单元应变增量较小时，有： （4.70）当单元应变增量较大时，单元应力增量应当用切线模量按增量法求下列积分的近似值： （4.71）对于弹塑性单元，单元应力增量由弹性和塑性两部分组成。若用r表示弹性应力增量与总应力增量之比参看图4-5（a）和相应的一维问题图4-5（b）,则可通过满足加载条件确定这个比例因子。确定出r以后，当单元应变增量较小时，有 （4.72）

20、（4.73）当单元应变增量较大时，单元应力增量应当用切线模量按增量法求下列积分的近似值： （4.74）显然，当r=0时，此式即为（4.71）式。（3）比例因子r的确定 假定用或表示弹性应力偏张量，于是把（4.66）（4.68）式中的或换 2 ri-1 i-1 10 (a) i-1 ri-1 ri-1 i-1(b)图45为或，并令其值为零便得 （4.75）其中A、B和C，Mises模型等向强化材料为（4.76）（4.77）（4.78）Mises模型随动强化（线化）材料为 （4.79）（4.80）（4.81）Drucker-Prager模型为理想塑性材料为（4.82）（4.83）（4.84）（4）

21、积分近似值的计算（4.71）和（4.74）式中的积分很难积得显式。为此，我们可以把积妥对应的塑性变形部分再分为若干子增量，用分段线性的计算结果去逼近积分结果。由于（4.71）式是（4.74）的r=0的特殊情况，所以这里只讨论（4.74）式中积分的近似计算问题。若把，再分为N个相等的子增量（例如图46中N5），则有（4.85）123456(1-r)i-1 图46由于初始子增量步开始时的单元应力为（4.86）第j子增量步开始时的单元应力为（4.87）所以（4.88）（5）第i增量步的单元计算过程在第I增量步，已知的是和，要计算的是和。若用表示判定单元应变增量大小的标准，则计算过程如图47所示。3T

22、resca模型由于Tresca模型的屈服面比较复杂，所以这里只限于平面应力问题的Euler-Cauchy解法。由§25中的和可知，Tresca模型弹塑性矩阵各元素均与主应力方向有关。当单元开始进入塑性状态以后，主方向随载荷的增加而改变，因而每一增量步都要重新对它进行计算。为了使应力满足加载条件，可以按下述插值方法对它进行计算第i增量步的主方向值。（1）按计算主应力和塑性功用（塑性单元）和（弹性单元）以及它们的组合（弹塑性单元）形成，求得位移增量以后，有 （4.89）（4.90）（4.91）（4.92）（4.93）（4.94）（4.95）（2）按计算主应力和塑性功仿照（1）中的计算，有

23、（4.96）（4.97）（4.98）（4.99）（4.100）（4.101）（3）按插值方法计算主方向由、和按一种插值方法求，使得用它求得的应力满足加载条件。对于图26中的边6，有（4.102）§44弹塑性手算例题我们仍以§42中的弹塑性拉压问题为例，用Euler-Cauchy法、 Euler一次迭代法、初应变增量法和初应力增量法进行手算。为了反映应力与加载历史的相关性，除Euler-Cauchy法外，均用叙述方式。（1）用Euler-Cauchy法求解仿照§42中用Newton-Raphson法的求解，有（R）约束处理后可得（S）若以总载荷为参考载荷，并按两个增量步取，则有于是无量纲逐步求解公式为（T） >0返回 =0 &g