非线性有限元读书报告

非线性有限元读书报告1 对张量分析及其在结构非线性力学和有限元中应用的认识。

张量分析将张量中的每一项的计算过程（这些计算过程可能是相似的）用张量之间的计算方法代替，并统一，提供了一种更加方便的计算方法。

张量标记中，指标不出现，因此张量标记的表示是独立于坐标系统的，并且可以用于笛卡尔坐标系，柱坐标系，曲线坐标等。

如同矩阵一般，或者说张量就是矩阵的更一般的形式，而矩阵是二维的张量。

此外，在张量标记中的方程非常容易记忆。

2 对各种应变、应力张量和材料本构关系的认识。

本段回答除本课程讲义外大部分参考文献[1] 2.1 应变张量为了使刚体运动时，特别是刚体转动时，应变度量为0，工程应变张量被抛弃了。

根据采用连续介质力学中的欧拉（Eulerian ）方法或拉格朗日方法（Lagrangian ），得到的应变分别为阿尔曼西（Almansi ）应变（或称欧拉应变）张量（E ij ε）与格林（Green ）应变（或称拉格朗日应变）张量（G ij ε），它们与工程应变（或称柯西小应变张量）（此处记为ij ε）均不相同。

三者的几何方程表达式为：阿尔曼西（Almansi ）应变（或欧拉应变）张量：()11,,,,11()22E ij ij ki kj i j j i k i k j F F u u u u εδ--=-=+-； 格林（Green ）应变（或拉格朗日应变）张量：(),,,,11()22G ij ki kj ij i j j i k i k j F F u u u u εδ=-=++； 工程应变（柯西小应变）张量：,,1()2ij i j j i u u ε=+。

其中u 为位移，ij F 为变形梯度矩阵，即i ij j x F X ∂⎡⎤=⎢⎥∂⎣⎦另外，阿尔曼西（Almansi ）应变（或欧拉应变）张量（Eij ε）和格林（Green ）应变（或拉格朗日应变）张量（Gij ε）之间的转换关系为：1122G Ekl ik jl ijF F εε--= 其中i ij j x F X ∂⎡⎤=⎢⎥∂⎣⎦，x 为欧拉坐标，X 为拉格朗日坐标。

2.2 应力张量当平衡方程建立在当前构型时，由柯西假定得到的应力张量为柯西应力（Cauchy ）张量（ij σ，在本课程讲义中称为欧拉应力张量）。

这种应力张量是真实的应力张量，存在 T ij ij σσ=。

但是在大变形几何非线性问题中，当前构型的边界尚未确定。

因此在初始构型上定义另外两种应力以避免这种困难。

将柯西应力（Cauchy ）张量（在本课程讲义中称欧拉应力张量）中平衡方程的面积和法线的形式直接变为参考构型，则得到了拉格朗日（Lagrange ）应力张量（N ij σ，或称名义应力张量）。

它不是对称的，转置被称为第一Piola-Kirchhoff应力。

将拉格朗日（Lagrange ）应力张量中的外力项乘以1ij F -便得到了第二Piola-Kirchhoff 应力张量（K ij σ，在本课程讲义中称为克希荷夫（Kirchhoff ）应力张量），这种变换使得第二Piola-Kirchhoff （PK2）应力张量变为对称的，并且，它与格林（Green ）应变张量共轭。

三者的平衡方程表达式分别为：柯西（Cauchy ）应力张量（ij σ，本课程讲义中称为欧拉应力张量，也是真实应力）：i i ij dT tdA n dA σ==；拉格朗日（Lagrange ）应力张量（Nij σ，或称名义应力张量）： 0000Ni i ij dT t dA n dA σ==；第二Piola-Kirchhoff 应力张量（PK2），在本课程讲义中称为克希荷夫（Kirchhoff ）应力张量（N ij σ）： 1100000K i i ij F dT F t dA n dA σ--==；其中i dT 为外力的微分，dA 为面积的微分，i n 为方向向量，ij F 为变形梯度矩阵，下标中带0的为建立在参考构型上的方程，不带0的为建立在当前构型上的方程。

三个应力张量之间转换关系如下：11NKij ijik kj ijik kl jl F F F F F σσσ--==；1N Kkj ij ik kj ik jk F F F σσσ-==； 111K N kl ij ik kl jl kj jl F F F F σσσ---==；2.3 本构关系本构关系即为物理方程，对于各项同性材料，在线弹性假定中，得到了名义应力张量（拉格朗日（Lagrange ）应力张量，N ij σ）和名义应变张量（柯西小应变张量，ij ε）之间的三个关系，它们是等价的，即：N ij ijkl kl E σε=1,2N ij ijkl ij kl ij U U E σεεε∂==∂ N ijij ijklE ttσε∂∂=∂∂式中，ijkl E 即为广义胡克定律规定得到的本构关系。

讨论有限应变时的本构关系也是从这三个形式向外推广的。

2.3.1对于弹性介质弹性介质是对过去的历史没有记忆的材料，仅考虑等温过程。

对于Cauchy（Euler ）应力张量（ij σ）与Almansi （Euler ）应变张量（Eij ε），存在以下本构关系：Eij ijkl klC σε= 仅当小应变时，基于当前构型的本构关系ijkl C 才退化为广义胡克定律。

根据Cauchy （Euler ）应力张量（ij σ）和第二Piola-Kirchhoff 应力张量（Kij σ）之间的关系，以及Almansi （Euler ）应变张量（Eij ε）和Green （Lagrange ）应变张量之间的关系，得到第二Piola-Kirchhoff 应力张量（N ij σ）和Green （Lagrange ）应变张量（G kl ε）之间的本构关系：11110K G Gij ij im jn mnpq kp lq kl ijkl klF F F C F F C σεε----== 2.3.2对于超弹性（hyperelastic ）介质单位质量应变能的变化率等于应力做功率，这种由能量关系来定义本构的关系称为超弹性介质，通过求应变能的时间变化率，由定义：011E Gij ijK ij ijd d dU dt dt dtεεσσρρ== 故，采用欧拉方法时：1ij EijUσρε∂=∂ 采用拉格朗日方法时：01K ij GijUσρε∂=∂ ρ为瞬时密度，0ρ为参考构型的密度，U 为应变能。

可以证明，ij σ与E ij ε之积等于Kij σ与G kl ε之积，均为单位质量应变能变化率。

因此前两者共轭，后两者共轭。

2.3.3对于次弹性（亚弹性）（hypoelastic ）介质应力变化率是应变变化率的齐次线函数，这种定义本构的关系称为次弹性介质，有三种频繁应用在当前有限元软件中的Cauchy 应力的客观变化率，这些客观变化率存在的原因很多是源于消除刚体转动的影响：(1)Jaumann 变化率ij JTij ik kj ik kj D W W Dtσσσσ∇=--(2)Truesdell 变化率ij Tk i i ij ij kj ik k k kD v v vDtx x x σσσσσ∇∂∂∂=+--∂∂∂(3)Green-Naghdi 变化率ij G T ijik kj ik kj D Dtσσσσ∇=-Ω-Ω（以上具体符号的意思尚未理解）而对于kirchhoff 应力Kij σ，它的变化率率与Cauchy 应力的变化率直接的关系为：11K ijpq q p k ip jq pq pk qk k k k D D v v v F F F Dt Dt x x x σσσσσ--∂∂⎛⎫∂=+-- ⎪∂∂∂⎝⎭一般来说，前面所述的几种本构关系并不等价，实际材料适合哪种定义，必须通过试验来确定。

3 对TL 法、UL 法有限元列式以及切线刚度矩阵的认识。

本段回答主要参考课件4.A 连续力学基础与几何非线性有限元(自学、参考)；5 小变形几何非线性问题有限元分析 5.1 小变形几何非线性有限元方程2019；5.4.2 平面Euler-Bernoulli 梁-简化TL 列式2019。

在固体力学中，Lagrangian 网格是应用最普遍的，其吸引力在于它们能很容易地处理复杂的边界条件，并能跟踪材料点。

在Lagrangian 有限元的发展中，一般采用两种方法： 3.1 TL 有限元列式以Lagrangian 度量的形式表述应力和应变的公式，导数和积分计算采用相应的Lagrangian 坐标X ，称为完全的Lagrangian 格式。

3.2 UL 有限元列式以Eulerian 度量的形式表述应力和应变的公式，导数和积分计算采用相应的Eulerian 坐标X ，称为更新的Lagrangian 格式。

尽管完全的和更新的Lagrangian 格式表面看起来有很大区别，但这两种格式的本质是相同的，主要区别在于在参考构型还是当前构型上描述变量。

在完全的Lagrangian 格式中习惯采用一个应变的完全度量，而在更新的Lagrangian 格式中常常采用应变的率度量，但是两者可以互换。

这两种格式在理论上是等效的，但在求解的有限元矩阵方程本身和求解步骤上仍有一定的差别。

在有限变形，材料非线性等问题中，刚度矩阵不再是线性的，或者常系数的。

因此需要采用增量分析方法，这就要求获得切线刚度矩阵。

3.3 切线刚度(1) TL 格式中的切线刚度矩阵：采用虚功原理,设一个运动的物体，假定对于时刻m t ，以前的解均已得到，当前的目标是去建立一个方程由方程可解得时刻1m m t t t +=+∆的各未知量．我们知道柯西应力与阿尔曼西应变是一对共扼的应力和应变，因此在时刻1m t +的虚功方程可表示为（采用Euler 应力和Almansi 应变）：11,11m m G m m ij ij V dV W σδεδ++++=⎰，11111m m m m m k k kkV S W F u dV T udS δδδ+++++=+⎰⎰k F ，k T 分别表示物体所受的体积载荷和表面载荷。

在拉格朗日体系中，采用第二Piola-Kirchhoff 应力和格林应变这一对共轭的应力应变张量进行分析，并且以0t =时刻的构型为参考构型，此时虚功方程为：1,1,1m K m G m ij ij V W dV δσδε+++=⎰，0111m m m kkk k V S W F udV T u dS δδδ+++=+⎰⎰离散化后，物体由m t 时刻变化到1m t +时刻时，参照构形有：,1,,1K m K m K m ij ij ij σσσ++=+∆，,1,,1G m G m G m ij ij ijεεε++=+∆，11m m m k k k u u u ++=+∆ 由格林应变的定义：()()()()()12,1,1,1111,,,,,,,,,,,,111111,,,,,,,,,1,1112212G m G m G mij ij ij m m m m m m m m m m m m i j i j j i j i k i k i k j k j i j j i k j k j m m m m m m m m i j j i k i k j k i k j k i k j GL m GL m ij ij u u u u u u u u u u u u u u u u u u u u εεεεε++++++++++++++∆=-=+∆++∆++∆+∆-++=∆+∆+∆+∆+∆∆=∆+∆,1,1,1GNL m GL m GNL m ij ij ijεεε++++∆=∆+∆ 其中()1111,,,,,,12GL m m m m m m ij i j j i k i k j k i k j u u u u u u ε++++∆=∆+∆+∆+∆，为线性部分 11,,12GNL m m ijk i k j u u ε++∆=∆∆，为非线性部分 设材料的本构关系有以下增量形式：0K Eij ijkl klC σε∆=∆则虚功方程可写为（注意到,0G mij δε=）： ()()()()01,1,1,,1,,1,,1,,1,10,1,1,10,1,1,m K m G m K m K m G m G m ij ij ijij ij ij V V K m GL m K m GNL m ij ij ijij V V TTGL m m GL m GNL m m GL m klijkl ij kl ijkl ij VV GL m klW dV dVdV dVC dV C dV δσδεσσδεεσδεσδεεδεεδεε+++++++++++++==+∆+∆=∆+∆+∆∆+∆∆+∆⎰⎰⎰⎰⎰⎰()()0010,1,1,10,1,1T Tm GNL m GNL m m GNL m ijkl ij klijkl ij V V C dV C dV δεεδε++++++∆+∆∆⎰⎰在有限元中，位移增量是由节点位移增量与形函数线性表达的，即：[]u N q ∆=∆而格林应变是位移的非线性函数，故（省略下标）：GL L ij B q ε⎡⎤∆=∆⎣⎦，GNL NL ij B q ε⎡⎤∆=∆⎣⎦代入整理后的虚功方程，得：()()0001,1,1,,1,,1,1,,1,,10,1,1,10,1,1m K m G m K m K m G m G m ij ij ijij ij ij V V TTL m K mNL m K mij ij V VT TTL m m L m TNL m m NL m ijklijklVW dV dV B qdV B qdVq B CBqdV q BCBδσδεσσδεεσδσδδ+++++++++++++==+∆+∆⎡⎤⎡⎤=∆+∆⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎡⎤+∆∆+∆⎣⎦⎣⎦⎣⎦⎣⎦⎰⎰⎰⎰⎰000,10,1,1,10,1,1V TTTL m m L m T NL m m NL m ijkl ijkl VVqdVq B C B qdV q B C B qdVδδδ++++++∆⎡⎤⎡⎤⎡⎤⎡⎤+∆∆+∆∆⎣⎦⎣⎦⎣⎦⎣⎦⎰⎰⎰同时，注意到0m k u δ=，故：[][]011111m m m m m kkk k kkk k V S V S k k V S W F udV T u dS F udV T u dSF N qdV T N qdSδδδδδδδ+++++=+=∆+∆=∆+∆⎰⎰⎰⎰⎰⎰则由前两式得到：[][]000,1,,1,,10,1,1,10,1,1,10,1,1TL m K mk k ij V S V TTNL m K m TL m m L m ijijkl V VT TTNL m m NL m TL m m L m ijkl ijklVF N qdV T N qdS B qdV B qdV q B C B qdV q B C B qdV q BCBδδσδσδδδδ+++++++++++⎡⎤∆+∆-∆⎣⎦⎡⎤⎡⎤⎡⎤=∆+∆∆⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎡⎤+∆∆+∆⎣⎦⎣⎦⎣⎦⎣⎦⎰⎰⎰⎰⎰⎰00,10,1,1VTTNL m m NL m ijklVqdVq B CBqdVδ+++∆⎡⎤⎡⎤+∆∆⎣⎦⎣⎦⎰⎰该式对于任意符合边界条件的q δ∆均成立，故：[][]0000,1,,1,,10,1,1,10,1,1,10,1,1,1TL m K mk k ij V S V TTNL m K mTL m m L m ijijkl V VT TTNL m m NL m TL m m L m ijklijklVVT NL m F N dV T N dS B dV BdV q B C B dVq B C BdV q BCBdVq B σσ++++++++++++⎡⎤+-⎣⎦⎡⎤⎡⎤⎡⎤=+∆⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎡⎤+∆+∆⎣⎦⎣⎦⎣⎦⎣⎦⎡⎤+∆⎣⎰⎰⎰⎰⎰⎰⎰00,1,1Tm NL m ijkl VC B dV ++⎡⎤⎦⎣⎦⎰将能够连乘的线性项合并，将[]0k V F N dV ⎰记为{}F ，将[]k S T N d S ⎰记为{}T ，,,1K m L m ij V B dV σ+⎡⎤⎣⎦⎰为应力在节点上的等价合力，记为{}P 。