1、第六章第六章 非线性有限元法（几何非线性）非线性有限元法（几何非线性）1、变形体的运动描述x3x1x2P0t0=0tn+1=tn+tn tn PnPn+1A0An+1 An 变形体上的质点的运动状态变形体上的质点的运动状态可以随不同的坐标选取以下几可以随不同的坐标选取以下几种描述方法：种描述方法：1 1、全拉格朗日列式法全拉格朗日列式法（T.LT.L列式列式法法Total Lagrangian Formulation)Total Lagrangian Formulation)：选取选取t t0 0=0=0时刻未变形物体的构时刻未变形物体的构形形A A0 0作为参照构形进行分析。作为参照构形进行
2、分析。2 2、修正拉格朗日列式法修正拉格朗日列式法（U.LU.L列式法列式法Updated Lagrangian FormulationUpdated Lagrangian Formulation）：选取选取t tn n时刻的物体构形时刻的物体构形A An n作为参照构形。由于作为参照构形。由于A An n随计算而变化，因随计算而变化，因此其构形和坐标值也是变化的，即与此其构形和坐标值也是变化的，即与t t有关。有关。t tn n为非线性增量求解时增为非线性增量求解时增量步的开始时刻。量步的开始时刻。3 3、欧拉描述法欧拉描述法(Eulerian Formulation)(Eulerian F
3、ormulation)：独立变量是质点当前时刻的位置独立变量是质点当前时刻的位置x xn+1n+1与时间与时间t tn+1n+1。几何非线性的有限元方程一几何非线性的有限元方程一般采用般采用T.LT.L或或U.LU.L列式法建立！列式法建立！、变形梯度张量x3x1x2PPP P初始初始/未变形未变形 变形后变形后 位移位移u u x x x x 1 1、首先采用、首先采用LagrangianLagrangian方法，方法，将一个物体的加载过程划分为将一个物体的加载过程划分为一系列平衡状态。一系列平衡状态。iiiuxx位移方程位移方程 初始状态与变形后状态之间坐初始状态与变形后状态之间坐标关系为
4、：标关系为：2 2、然后，考虑材料方向矢量，这个矢量、然后，考虑材料方向矢量，这个矢量描述物体内一段无限小的单元。描述物体内一段无限小的单元。jijjjiixdFxdxxdxx3x1x2 ixPiidxxPiixdxQ ixP式中，式中，F Fij ij称为变形梯度张量。称为变形梯度张量。jiijxxF初始状态与变形后状态之间材料方向矢量初始状态与变形后状态之间材料方向矢量的关系：的关系：、变形梯度张量jijijiijxuxxxxFjjiiijijFFFFI212由位移方程，得：由位移方程，得：jiijijxuFijijiiFFI1由二阶张量特性，变形梯度张量由二阶张量特性，变形梯度张量的三个
5、不变量为：的三个不变量为：VJdVdFdVij detJFIij det3 由于由于F Fij ij表示从初始状态到变表示从初始状态到变形后状态的一个映射，其逆映射形后状态的一个映射，其逆映射Fij-1一定存在，即：一定存在，即：AdFJNdAnijij1或写为：或写为：体积映射体积映射:面积映射：面积映射：变形前面积变形前面积dA dA Ni(初始面积法向矢量初始面积法向矢量)变形后面积变形后面积dAdAni(变形后面积法向矢量变形后面积法向矢量)映射映射F Fij ij逆映射逆映射F F-1-1ij ijF Fij ij是一个二阶张量。是一个二阶张量。jiijjiijxuxxF1、应变与变
6、形测度 由于变形梯度张量由于变形梯度张量F Fij ij中包含了刚体运动，因此不能直接用于定中包含了刚体运动，因此不能直接用于定义应变测度。而材料方向矢量则不包含刚体运动，因此它的平方值义应变测度。而材料方向矢量则不包含刚体运动，因此它的平方值可以作为衡量从某一状态到变形后状态的一个测度，定义为：可以作为衡量从某一状态到变形后状态的一个测度，定义为：iixdxdsd2初始状态初始状态:一个应变测度应该能反映出材料一段一个应变测度应该能反映出材料一段长度发生的改变。因此，应变张量可以由长度发生的改变。因此，应变张量可以由下式定义：下式定义：iiiixdxddxdxsdds22x3x1x2 ixP
7、iidxxPiixdxQ ixPiidxdxds2变形后状态：变形后状态：提醒：提醒：由于由于GreenGreen应变张量表达式中的变形梯度张量对应于初始状应变张量表达式中的变形梯度张量对应于初始状 态，因此该应变张量也应在初始状态下计算。态，因此该应变张量也应在初始状态下计算。、应变与变形测度、AlmanshiAlmanshi应变张量应变张量1 1、Green Green 应变张量应变张量GreenGreen应变张量采用应变张量采用LagrangianLagrangian运运动描述方法，即按初始状态下的动描述方法，即按初始状态下的构形定义应变张量。构形定义应变张量。iiijiiijkjkii
8、ijkjkiiiiiixdxdexdxdFFxdxdxdFFxdxdxddxdxsdds222ijkjkiijFFe21式中，式中，e eij ij称为称为GreenGreen应变张量应变张量或或Green-LagrangianGreen-Lagrangian应变张量应变张量。AlmanshiAlmanshi应变张量采用应变张量采用EularEular运动运动描述方法，即按当前状态下的构描述方法，即按当前状态下的构形定义应变张量。形定义应变张量。iiijiikjkiijjkjkiiiiiiiidxdxEdxdxFFdxFFdxdxdxxdxddxdxsdds21111221121kjkiiji
9、jFFE式中，式中，E Eij ij称为称为AlmanshiAlmanshi应变张量应变张量或或Almanshi Almanshi EularEular应变张量应变张量。由于大变形问题有由于大变形问题有限元方程主要采用限元方程主要采用T.LT.L列式法列式法或或U.LU.L列式列式法法建立，因此应在初建立，因此应在初始状态下定义应变张始状态下定义应变张量，即采用量，即采用GreenGreen应应变张量。变张量。可以证明可以证明GreenGreen应变张量和应变张量和AlmanshiAlmanshi应变张量都是二阶对称张量。应变张量都是二阶对称张量。、应变与变形测度2 2、Green Green
10、 Lagrangian Lagrangian应变张量应变张量e eij ij与小应变张量与小应变张量ij ij的关系的关系 将变形梯度张量表达式代入到将变形梯度张量表达式代入到GreenGreen应变张量公式中，得：应变张量公式中，得：ijijjkikijjiijjkikijjiijijjkkjikkiijxuxuxuxuxuxuxuxuxuxue21212121ijjiijxuxu21式中：式中：为小变形应变张量；为小变形应变张量；kjkiijFFCjkikijxuxu212 2、GreenGreen变形张量也可写为：变形张量也可写为：为非线性二次项为非线性二次项1 1、GreenGreen
11、应变张量应变张量为小应变张量与一个非线性二为小应变张量与一个非线性二次项之和，这意味所有大变形次项之和，这意味所有大变形分析都是非线性的。分析都是非线性的。ijijijCe21ijijije式中，式中，C Cij ij是是CauchyCauchy变形张量变形张量由于由于CauchyCauchy变形张量是正定对称变形张量是正定对称阵，因此该张量有三个实特征值；阵，因此该张量有三个实特征值；这些特征值的平方根记为材料的这些特征值的平方根记为材料的主轴拉伸。主轴拉伸。、大变形的应力测度1 1、柯西应力张量、柯西应力张量(Cauchys stress(Cauchys stress tensor)ten
12、sor)取三维空间笛卡儿坐标系，在取三维空间笛卡儿坐标系，在t t时刻时刻的现时构形中截取一个四面体素，斜面的现时构形中截取一个四面体素，斜面的法线为的法线为n n，另外三个面元与所取坐标，另外三个面元与所取坐标面平行。由四面体素的平衡条件得出其面平行。由四面体素的平衡条件得出其上的应力为：上的应力为：jijinn niin3x2x1x这里这里ij=ji便是便是柯西应力张量柯西应力张量，它是二阶对称张量。，它是二阶对称张量。、柯西、柯西(Cauchy)(Cauchy)应力张量是一种采用欧拉描述法应力张量是一种采用欧拉描述法(是以质点的瞬时是以质点的瞬时坐标坐标x xk k和时间和时间t t作为
13、自变量描述作为自变量描述)定义在定义在t t时刻的现时构形上的应力张时刻的现时构形上的应力张量量ij ij，又称，又称欧拉应力张量欧拉应力张量。、在大变形、在大变形(有限变形有限变形)情况下，由于变形前的初始构形和变形后情况下，由于变形前的初始构形和变形后的现时构形差别较大，柯西的现时构形差别较大，柯西(Cauchy)(Cauchy)应力张量难于适应。应力张量难于适应。柯西应力是定义在现柯西应力是定义在现时构形（变形后状态时构形（变形后状态下）的单位面积上的下）的单位面积上的力，是与变形相关的力，是与变形相关的真实应力。真实应力。3、大变形的应力测度2 2、一阶、一阶Piola-Kirchof
14、fPiola-Kirchoff应力张量应力张量 一阶一阶Piola-KirchoffPiola-Kirchoff应力张量的定义是建立在总应力张量的定义是建立在总力相等的基础上。即：在参考状态下该应力张量力相等的基础上。即：在参考状态下该应力张量能给出与变形后状态下柯西应力张量相同的力。能给出与变形后状态下柯西应力张量相同的力。变形后状态下：变形后状态下：dAndPjiji称为一阶称为一阶Piola-KirchoffPiola-Kirchoff应力张量应力张量或或名义应力名义应力参考后状态下：参考后状态下：AdNTdPjiji变形前面积变形前面积dA dA Ni(参考面积法向矢量参考面积法向矢量
15、)变形后面积变形后面积dAdAni(变形后面积法向矢量变形后面积法向矢量)AdFJNdAnijij1将面积映射关系：将面积映射关系：代入上式，得：代入上式，得：1jkikijJFTAdNTdAnjijjijAdNTAdFJNjijkjkij1同样，柯西应力张量也可以由一同样，柯西应力张量也可以由一阶阶Piola-KirchoffPiola-Kirchoff应力张量表示：应力张量表示：ikjkijTFJ1 从该式可以看出，一阶从该式可以看出，一阶Piola-KirchoffPiola-Kirchoff应力应力张量张量提供了以参考状态表示实际力的形式。但提供了以参考状态表示实际力的形式。但是，直接
16、应用一阶是，直接应用一阶Piola-KirchoffPiola-Kirchoff应力张量应力张量可能可能存在以下两个困难：存在以下两个困难：1 1、从能量角度上，、从能量角度上，T Tij ij不适合与不适合与GreenGreen应变张量应变张量共同使用。因为共同使用。因为T Tij ij乘以乘以GreenGreen应变张量不会产应变张量不会产生与生与CauchyCauchy应力张量与小应变张量相同的能量应力张量与小应变张量相同的能量密度。密度。2 2、T Tij ij不对称，因而较难应用到有限元分析中。不对称，因而较难应用到有限元分析中。、大变形的应力测度3 3、二阶、二阶Piola-Kir
17、choffPiola-Kirchoff应力张量应力张量 如不采用变形后状态如不采用变形后状态dPdP推导应力张量，而推导应力张量，而是将作用在变形后状态下的是将作用在变形后状态下的dPdP映射到未变形映射到未变形状态上（映射是采用逆变形梯度张量），即：状态上（映射是采用逆变形梯度张量），即：jijidPFPd1AdNSPdjiji这样可以定义另一个应力张量这样可以定义另一个应力张量S S，它给出了，它给出了未变形状态下作用在未变形面积上的总力：未变形状态下作用在未变形面积上的总力：现在，变换柯西应力张量，使：现在，变换柯西应力张量，使：dAnFdPFPdkjkijjiji11AdNJFFPdr
18、rkjkiji11AdFJNdAnijij1将面积映射关系将面积映射关系 代入上式：代入上式：(1)(1)(2)(2)(3)(3)(4)(4)对比对比(2)(2)、(4)(4)式可得：式可得：11rkjkijijJFFSrkijijjkFJSF1S Sij ij称为称为二阶二阶Piola-KirchoffPiola-Kirchoff应力张量应力张量或或伪应力伪应力同样，由上式可得：同样，由上式可得：二阶二阶Piola-KirchoffPiola-Kirchoff应力张量应力张量S Sij ij的性质：的性质：S Sij ij是对称阵；是对称阵；S Sij ij在能量角度下与在能量角度下与Gre
19、enGreen应变张应变张量协调，即：量协调，即：该表达式的优点在于等式右边是该表达式的优点在于等式右边是在参考状态下计算的。在参考状态下计算的。1.1.S Sij ij与与T Tij ij有以下关系：有以下关系：ijijijijeSrjijriTSF二阶二阶Piola-KirchoffPiola-Kirchoff应力张量应力张量的物理意义是明确的：真实的物理意义是明确的：真实的力元可以看成是由的力元可以看成是由S Sij ij定义定义的力元经与变形相同的方式的力元经与变形相同的方式被被“拉长和转动拉长和转动”后得到的。后得到的。、大变形的应力测度4 4、三个应力张量的比较、三个应力张量的比较
20、 张量张量 作用力作用力 作用面积作用面积柯西应力张量柯西应力张量ij ij 变形后状态下的力变形后状态下的力 变形后状态下的面积变形后状态下的面积一阶一阶P-KP-K应力张量变形后状态下的力未变形状态下的面积应力张量变形后状态下的力未变形状态下的面积二阶二阶P-KP-K应力张量未变形状态下的力未变形状态下的面积应力张量未变形状态下的力未变形状态下的面积因此，虽然二阶因此，虽然二阶P-KP-K应力张量有其应用上的优点，但其本身的物理应力张量有其应用上的优点，但其本身的物理意义很难理解。它主要是起到求解大变形问题的桥梁作用，通过意义很难理解。它主要是起到求解大变形问题的桥梁作用，通过它计算出柯西
21、应力张量。它计算出柯西应力张量。、几何非线性有限元方程的建立如前所述，几何非线性的有限元方程一般采用如前所述，几何非线性的有限元方程一般采用T.LT.L或或U.LU.L列式法建立：列式法建立：1 1、全拉格朗日列式法全拉格朗日列式法（T.LT.L列式法列式法)：选取选取t t0 0=0=0时刻未变形物体的构形时刻未变形物体的构形A A0 0作为参照构形进行分析。作为参照构形进行分析。2 2、修正拉格朗日列式法修正拉格朗日列式法（U.LU.L列式法）列式法）：选取选取t tn n时刻的物体构形时刻的物体构形A An n作为参照构形。由于作为参照构形。由于A An n随计算而变化，因随计算而变化，
22、因此其构形和坐标值也是变化的，即与此其构形和坐标值也是变化的，即与t t有关。有关。t tn n为非线性增量求解时增为非线性增量求解时增量步的开始时刻。即增量分析。量步的开始时刻。即增量分析。x3x1x2P0t0=0tn+1=tn+tn tn PnPn+1A0An+1 An 图示物体同时作用有体积力图示物体同时作用有体积力fib和和面力面力fiS，在时刻，在时刻t tn+1n+1=t=tn n+t tn n的平的平衡方程可以按虚功原理建立：衡方程可以按虚功原理建立：S*iSiV*ibiVijijdSufdVufdV*提醒：提醒：该方程此时不可解，因为应力该方程此时不可解，因为应力和应变在变形后
23、状态下表示未知。和应变在变形后状态下表示未知。、几何非线性有限元方程的建立2 2、在外力作用点和方向都不改、在外力作用点和方向都不改变的条件下，也可以将体积力变的条件下，也可以将体积力fib和面力和面力fiS定义到初始状态下：定义到初始状态下：dSfSdfdSfSdfbibiSiSi S*iSiV*ibiVijijdSufdVufdV*S*iSiV*ibiVijijSdufVdufVdeS*提醒：提醒：上式给出的虚功方程是从上式给出的虚功方程是从变形后状态下变形后状态下的虚功方程转换的虚功方程转换而来，因此是准确的，但是已经完全定义在初始状态下了。而来，因此是准确的，但是已经完全定义在初始状态
24、下了。为了求解，需将以上变形后状态下表示的虚功方程转换到为了求解，需将以上变形后状态下表示的虚功方程转换到初始状态下表达。初始状态下表达。VijijVijijVdeSdV*1 1、采用二阶、采用二阶PiolaPiola应力张量和应力张量和GreenGreen应变张量将虚应变能转换应变张量将虚应变能转换到初始状态下表示：到初始状态下表示：将以上关系代入到虚功方程中：将以上关系代入到虚功方程中：得：得：(a)(a)、几何非线性有限元方程的建立ijR21表示该张量对应的时刻：表示该张量对应的时刻：1 1代表初始代表初始状态时刻，状态时刻，2 2为变形后状态时刻；为变形后状态时刻；如该标识缺省，则表示
25、从初始状态如该标识缺省，则表示从初始状态变化到变形后状态该张量的增量。变化到变形后状态该张量的增量。代表定义该张量所对应的构形：代表定义该张量所对应的构形：1 1为初始状态构形，为初始状态构形，2 2为变形后为变形后状态构形；如该标识缺省，则状态构形；如该标识缺省，则为初始状态构形。为初始状态构形。在利用增量法（在利用增量法（修正拉格朗日列式法）修正拉格朗日列式法）求解时，为了分析的方求解时，为了分析的方便，在张量符号的左侧引入上下标，分别该张量对应时刻以及定义便，在张量符号的左侧引入上下标，分别该张量对应时刻以及定义该张量的构形：该张量的构形：S*iSiV*ibiVijijSdufVdufV
26、deS2121*2121当引入以上表示后，当引入以上表示后，按按t t1 1+t t时刻构形建立的虚功方程可以写为：时刻构形建立的虚功方程可以写为：或写为：或写为：QVdeSVijij21*2121式中，式中，表示外力所做的虚功。表示外力所做的虚功。S*iSiV*ibiSdufVdufQ212121、几何非线性有限元方程的建立引入此前引入此前GreenGreen应变张量表达式，可得：应变张量表达式，可得：ijijijijijijeeijjiijxuxu21线线性性项项为为：ijijijSSS11121jkikjkikijxuxuxuxu21非非线线性性项项为为：再将变形后状态下再将变形后状态下
27、KirchoffKirchoff应力张量表示为未变应力张量表示为未变形状态的形状态的KirchoffKirchoff应力张量加上一个应力增量：应力张量加上一个应力增量：(a)(a)(b)(b)ijS11ijijS111注意，式注意，式 (b)(b)中中 为作用在未变形构形上并以未变形状态下表示为作用在未变形构形上并以未变形状态下表示的的KirchoffKirchoff应力张量，实际上就是柯西应力张量：应力张量，实际上就是柯西应力张量：。虚功方程：虚功方程：QVdeSVijij21*2121ijijijSS1121(c)(c)、几何非线性有限元方程的建立S*iSiV*ibiSdufVduf111
28、1VdQV*ijij1111为为t tn n时刻初始构形上时刻初始构形上外力所做虚功。外力所做虚功。QVdeSVdQVdSVdSVdVdVdeSVijijVijijVijijVijijVijijVijijVijijijij21*11*1111*11*11*11*11*1*111将以上将以上 (a)(a)、(c)(c)两式代入到虚功方程中，可得：两式代入到虚功方程中，可得：QQVdVdeSVijijVijij1121*11*11即变形后状态下的虚功方程为：即变形后状态下的虚功方程为：式中：式中：为为t tn n+t t时刻初始构形上外力所做的虚功。时刻初始构形上外力所做的虚功。S*iSiV*ib
29、iSdufVdufQ212121这里，虚功方程中由于包含了非线性二这里，虚功方程中由于包含了非线性二次项，因此方程是非线性方程。这个方次项，因此方程是非线性方程。这个方程还不能直接求解。为了求解这个方程，程还不能直接求解。为了求解这个方程，需要将方程线性化。需要将方程线性化。6、非线性平衡增量方程的线性化 通常，可以假定应变增量和应力增通常，可以假定应变增量和应力增量之间以下线性本构关系：量之间以下线性本构关系：ijijklijeCS111ijS20ije2012ijS20ije10ijS10ije20ijSd20ijed20ijef10斜斜率率ijijijijklijeCS11111将以上关
30、系代入到虚功方程中，得：将以上关系代入到虚功方程中，得：QQVdVdeeCVijijVijijijkl1121*11*111然而上式依然包含有非线性二次项，不可直接求解。一般需要引入然而上式依然包含有非线性二次项，不可直接求解。一般需要引入以下线性化近似：以下线性化近似：则有以下线性化的非线性平衡增量方程：则有以下线性化的非线性平衡增量方程：QQVdVdCVijijVijijijkl1121*11*111从以上线性化过程可以看出，这种线性从以上线性化过程可以看出，这种线性化处理是有局限性的。在分析非线性大化处理是有局限性的。在分析非线性大应变时会造成较大误差。可以采用其他应变时会造成较大误差。
31、可以采用其他有限元格式，如摄动法有限元。有限元格式，如摄动法有限元。7、几何非线性问题的有限元基本方程 kkuSuBCC kVNLijNLVijij1NLuVdBBVdK 1111111 QQuKKkNLL1121221111 kVijNLVijijuVdBVdQ 11111 kVTVijklijklLuVdBCBVdCK 111111111kikikikixNxuNu单元刚度阵的形成单元刚度阵的形成第一步：坐标、位移插值第一步：坐标、位移插值 kuB第二步：应变、应力插值第二步：应变、应力插值将以上插值关系代入到线性化非线性平衡增量方程，得：将以上插值关系代入到线性化非线性平衡增量方程，得：式中：式中：1 1、杆元位移函数、杆元位移函数33221010 xbxbxbbvxaau22yvxuxu注意在以下推导中，杆元的结点位移和结点力实际上为相应杆元的注意在以下推导中，杆元的结点位移和结点力实际上为相应杆元的位移增量和结点力的增量。位移增量和结点力的增量。应变：应变：
