1.概述
固体力学中有3组基本方程，即本构方程、几何运动方程和平衡方程。经典线性理论基于3个基本假定，即材料的应力、应变关系满足广义虎克定律；位移是微小的；约束是理想约束。这些假定使得3组基本方程成为线性。只要基本假定中任何一个失效，问题就转化为非线性。表1-1给出了非线性问题的分类及基本特点。
表1-1非线性问题的分类及基本特点
几何非线性理论将平衡方程建立在结构变形后位置上。以图1-5所示结构为例，按线性理论求解就无法找到平衡位置。按几何非线性分析方法处理，在外力P作用下，B点产生竖向位移，当位移达到一定值δ时，AB、BC两杆件中轴力的竖向分力与P平衡，δ即为B点位移的解。可见，受力状态因变形而发生明显改变时，就必须用几何非线性方法进行分析。
图1-5 受集中力的二力杆
凡是在本构关系中放弃材料线性关系假定的理论，均属材料非线性范畴。根据不同的材料性态，又可以分成表1-2给出的几种不同的材料非线性问题。
桥梁结构以钢和混凝土作为主要建材，因此，涉及的材料非线性主要是非线性弹塑性问题和混凝土徐变问题。
表1-2几种材料非线性问题
2.几何非线性分析
在整个分析过程中，以t=0时的构形作为参考，且参考位形保持不变，这种列式称为总体拉格朗日列式（T.L列式）。
以杆系结构为例，对于任意应力-应变关系与几何运动方程，杆单元的平衡方程可由虚功原理推导得到：
式（1-70）中，（σ）——单元的应力向量；
（f）——单元杆端力向量；
V——单元体积分域，对T.L列式V是变形前的单元体积域；
（B）——应变矩阵，是单元应变与节点位移的关系矩阵。即：
d（ε）=（B）d（δ） （1-71）
（δ）——杆端位移向量。
在有限位移情况下，（B）是位移（δ）的函数矩阵，可分解为与杆端位移无关的部分[B0]和与杆端位移有关的部分（BL）两部分，即：
（B）=（B0）+（BL） （1-72）
采用增量列式法将式（1-70）写成微分形式：
或
根据式（1-72），式（1-74）左边第一项可写成：
当材料满足线弹性时，有：
（σ）=（D）（（ε）−（ε0））+（σ0） （1-76）
式（1-76）中：（ε0）——单元的初应变向量；
（σ0）——单元的初应力向量；
（D）——弹性矩阵。
于是，单元的应力、应变增量关系可表示成：
d（σ）=（D）d（ε） （1-77）
将式（1-71）和式（1-72）代入式（1-77）得：
d（σ）=（D）（（B0）+（BL））d（δ） （1-78）
于是，式（1-74）左边第二项可表示为：
记：
式（1-74）最后可表达为：
（0（k）0+（k）L+0（k）σ）d（δ）=0（K）Td（δ）=d（f） （1-82）
式（1-82）就是增量形式T.L列式的单元平衡方程。
式中0（K）T是3个刚度矩阵之和，称为单元切线刚度矩阵，它表示载荷增量与位移增量之间的关系，也可理解为单元在特定应力、变形下的瞬时刚度。0（k）0与单元节点位移无关，是单元弹性刚度矩阵；0（k）L称为单元初位移刚度矩阵或单元大位移刚度矩阵，是由大位移引起的结构刚度变化，是d（δ）的函数；0（k）σ称为初应力刚度矩阵，它表示初应力对结构刚度的影响，当应力为压应力时，单元切线刚度减小，反之单元切线刚度增大。
将各单元切线刚度方程按节点力平衡条件组集成结构增量刚度方程，即有：
0（K）Td（Δ）=d（P） （1-83）
式（1-83）中：0（K）T为结构切线刚度矩阵，可以由单元切线刚度矩阵按常规方法进行组集形成；d（P）为载荷增量。由于载荷增量一般取为有限值而不可能取成微分形式，结构在求得的位移状态下，抗力与总外载荷之间有一差量，即失衡力，结构必须产生相应位移以改变结构的抗力来消除这个失衡力。在计算中，一般通过迭代法来求解。
在建立t+Δt时刻物体平衡方程时，如果选择的参照构形不是未变形状态t=0时的构形，而是最后一个已知平衡状态，即以本增量步的起始时刻t的构形作为参照构形，这种列式法称为更新的拉格朗日列式法（U.L列式）。
由于采用了U.L列式，平衡方程式（1-74）中的积分须在t时刻单元体积内进行，且t（k）L的积分式是t（k）0的一阶或二阶小量，因此，代表[k]L的积分式可以略去。这是U.L列式与T.L列式的一个重要区别。最后，增量形式的U.L列式平衡方程可写成：
（t（k）0+t（k）σ）d（Δ）=d（P） （1-84）
3.材料非线性分析
桥梁结构材料非线性主要是非线性弹塑性问题和混凝土徐变问题，本节介绍非线性弹塑性问题的分析方法，混凝土徐变问题读者可参阅相关资料。
根据实验结果，单轴应力下材料的应力、应变关系可归结为如下几点。
1）应力在达到比例极限前，材料为线弹性；应力在比例极限和弹性极限之间，材料为非线性弹性。
2）应力超过屈服点，材料应变中出现不可恢复的塑性应变：
ε=εe+εp （1-85）
应力和应变间为非线性关系：
σ=ϕ（ε） （1-86）
3）应力在某一应力σ0（σ0＞σs，σs为材料的屈服点）下卸载，则应力增量与应变增量之间存在线性关系，即：
dσ=Edε （1-87）
为了判断是加载还是卸载，用如下加载准则：
●当σdσ.0时为加载，满足式（1-86）；
●当σdσ＜0时为卸载，满足式（1-87）。
4）在卸载后某应力σ下重新加载，则：
σ＜σ0时，dσ=Edε （1-88）
σ0为卸载前材料曾经受到过的最大应力值，称后屈服应力，若：
●σ0=σs，则材料称为理想塑性的；
●σ0＞σs，则材料称为硬化的。
5）从卸载转入反向力加载，应力、应变关系继续采用式（1-87）或式（1-88），一直到反向屈服。在复杂应力状态下，判断材料是否屈服，可以用应力的某种函数表示：
F（σij，K）=0 （1-89）
式（1-89）中：σij——应力状态；K——硬化函数。
若以σij为坐标轴建立一坐标空间，则式（1-89）的几何意义为空间超曲面。任一应力状态在此空间中代表一个点，当此点落在屈服面之内时，F（σij，K）＜0，材料呈弹性状态；F（σij，K）=0时，材料开始进入塑性。
常用的屈服条件如下。
1）屈雷斯卡（Tresca）屈服条件：假定最大剪应力达到某一极限值时，材料开始屈服，相当于材料力学中的第三强度理论。
2）密赛斯（Von Mises）屈服条件：假定偏应力张量的第二不变量达到某一极限时，材料开始屈服，相当于材料力学中的第四强度理论。
此外还有Drucker-Prager屈服准则、Zienkiewicz-Pande屈服准则等。
在弹塑性增量理论中，讨论仍限于小变形情况。于是，其应变-位移几何运动方程和平衡方程与线性问题相同，不需要进行任何变动。需要改变的只是在塑性区范围内用塑性材料的本构关系矩阵（Dep）代替原来的弹性系数矩阵（De）。因此，可直接得到弹塑性分析有限元平
衡方程：
（tKT）（Δtu）=（ΔtR） （1-90）
式（1-90）中：
（ΔtR）=（ΔtF）+（ΔtT）+（ΔtFc）−（ΔtFI） （1-92）
其中，（ΔtF）和（ΔtT）分别表示与结构面载荷f及体载荷t对应的等效节点力增量；（ΔtFc）为节点集中外载荷增量；（ΔtFI）为初应力或初应变增量引起的外载荷增量，它们在t-Δt至t时间的增量为：
对于初应力问题：(www.xing528.com)
对于初应变问题：
式（1-90）～式（1-96）给出了小变形弹塑性分析的有限元方程，式中[tKT]代表了载荷与位移增量的切线刚度，随不同加载历程而变化。求解这一问题的关键是计算单元的切线刚度矩阵和应力，由于本构关系[Dep]是当前应力的函数，即当前位移的隐函数，因此计算时要引入一个材料模型的子程序来处理塑性问题。这个子程序的主要计算内容与步骤如下。
1）由前边迭代的位移结果计算应变增量：
Δtε=Δtε（tu，t−Δtu） （1-97）
式（1-97）中：tu、t−Δtu——t与t−Δt时刻结构的位移。
2）暂假定Δε是弹性的，计算
Δtσe=DeΔtε （1-98）
3）由此推出新的应力状态为
tσ=t−Δtσ+Δtσe=t−Δtσ+DeΔtε （1-99）
4）核对在第二步中的假设是否符合事实。将式（1-99）代入加载函数中，计算当前的加载函数值：
tF=F（tσ，K） （1-100）
5）若tF-0，说明Δtε确实是弹性的，第二、三步中的计算正确，此子程序的执行可以结束。
6）若tF＞0，说明Δtε中包括了（或甚至全部是）塑性变形，则执行以下计算。
7）若本次迭代开始时的应力是弹性的，则本次迭代的应力增量中有一部分是弹性的，而另一部分是弹塑性的。将弹性部分记为：
mΔtσt=mDeΔtε （1-101）
显然，m＜1，将式（1-101）代入到加载函数中可解出m。
F（t−Δtσ+mΔtσe，K）=0 （1-102）
8）计算塑性部分应变增量及当前应力
Δtεp=（1−m）Δtε （1-103）
9）计算应变增量之塑性部分Δtεp所引起的应力。由于材料刚度矩阵是非线性的，因此这一计算应是积分过程。作为数值计算，可改为逐段线性化求和。为此，将Δtεp再细分为M个小的增量：
Δ（Δtεp）=Δtεp/M（ 1-105）
10）在每一个小的子增量Δ（Δtεp）（i）中，先根据子增量起始时的应力计算[Dep]（i），而
Δ（Δtσ）=[Dep]（i）Δ（Δtεp）（i） （1-106）
于是新的应力状态为：
tσ（i）=tσ（i−1）+Δ（Δtσ）（i） （1-107）
由tσ（i）可计算下一个子增量时的[Dep]（i+1），并重复以上步骤，结果
由此可形成最终状态的[Dep]。
以上方法将平衡迭代与本构迭代分开，主步进行平衡迭代，子步进行本构迭代，故称之为子增量法。
4.非线性方程组的求解
结构非线性有限元分析最终归结为一组非线性代数方程的求解。非线性代数方程组的求解方法很多，要根据问题的非线性程度、计算结果等因素来选择恰当的方法。下面介绍几种常用的求解方法。
（1）直接求解法
直接求解法是基于全量列式的求解过程，应用最多的是直接迭代法，由虚功原理建立的非线性有限元平衡方程为：
（K（{δ}））（δ）=（P） （1-109）
当设定位移向量（δ）的初值（δ0）后，改进的近似解可由下式得到：
（δ1）=（K0（（δ0）））−1（P） （1-110）
整个迭代过程可用下式表示：
（δn）=（Kn−1（（δn−1）））−1（P） （1-111）
当迭代结果满足预定的收敛准则时，就得到了所要求的节点位移向量。图1-6a为取（δ0）=（0）时单自由度问题的迭代过程取得收敛的示意图。
直接迭代法应用简单，运算速度一般也较快，可应用于具有轻微非线性的问题。这一求解过程的成功与否很大程度上取决于对初值位移{δ0}的正确估计。图1-6b表示的是直接迭代法迭代过程发散时的情形。为改善收敛性和收敛速度，可以对载荷进行分级。
图1-6 直接迭代法收敛和发散过程
（2）增量法
增量形式的有限元列式方法具有一个共同的特点：将整个载荷变形过程划分为一连串增量段，每一增量段中结构的载荷反应被近似地线性化。简单增量法将每一级增量载荷下直接求得的状态变量视作结构平衡状态，计算相应的切线刚度矩阵，进而进行下一级载荷计算，并不断累加其位移增量。图1-7描述了简单增量法的求解过程。
几何非线性问题的有限元分析最初多采用简单增量法进行，虽然这种求解方法对每一级载荷作用时的计算速度较快，但由于每一级载荷作用前结构并未精确地到达平衡位置，所求得的解答会随着增量过程的继续而越来越偏离真实的载荷-变形过程。为了保证计算精度，常常将增量区间划分得相当小。
此外，为了评价解的精度，一般要对同一问题在进一步细分增量区间后再次求解，通过两次解的比较判定是否收敛。这样就需要消耗大量的计算时间。
作为对这一方法的改进，可将不平衡力作为一种修正载荷并入下一级载荷增量。这就是有一阶自校正的增量法。一阶自校正增量法求解过程的示意图如图1-8所示。一阶自校正增量法具有较快的求解速度，同时也比简单增量法的计算精度高。这一方法在求解非线性问题特别是求解塑性问题时得到广泛的应用。
图1-7简单增量法的求解过程
图1-8 一阶自校正增量法的收敛过程
（3）Newton-Raphson法
对于多自由度体系，同样导出相应的迭代公式：
Δn为失衡力，式（1-112）即为Newton-Raphson法（N⋅R法）求解结构非线性问题的基本形式，其收敛过程如图1-9a所示。
由式（1-112）可见，N⋅R法在每次迭代后都要重新形成（K）T，对于大跨度桥梁结构进行这一过程很费计算机时间。为了减少形成总刚及其三角化分解的次数，有时用（K（δ0））T代替（K（δn））T，这样，仅进行一次切线刚度矩阵和三角化分解计算，后面的迭代只是线性方程组的回代，这种方法称为修正的N⋅R法（M⋅N⋅R法）。
图1-9b给出了该方法的迭代过程。M⋅N⋅R法在每次迭代中均用同一斜率，收敛较N⋅R法差。图1-10给出了N⋅R法和M⋅N⋅R法求解非线性方程组的流程，编程时可将这两种方法结合使用。
图1-9 N⋅R法的收敛过程
（4）收敛准则
在迭代计算中，为了中止迭代过程，必须确定一个收敛标准。在实际应用中，可以从结构的不平衡力向量和位移增量向量两方面来判断迭代计算的敛散性。
数的大小可以用其绝对值来衡量，而对于一个结构，无论其节点力还是节点位移都是向量，其大小一般用该向量的范数来表示。
设列向量（v）=（v1，v2，v3，…，vn）T，该向量的范数可以定义为：
图1-10 N⋅R法和M⋅N⋅R法迭代流程图
●各元素绝对值之和
●元素中绝对值最大者
这3个范数记为
（P=1，2，∞），应用中可任选其中的一种。
有了列向量的范数，无论是节点力向量还是节点位移向量，其“大小”均可按其范数的大小来判断。所谓足够小就是指其范数已小于预先指定的某个很小的数。
取位移增量为衡量收敛标准的准则称为位移准则，若满足下列条件，就认为迭代收敛：
式（1-116）中：αd——位移收敛容差
Δui+1——位移增量向量的某种范数
实践证明，对有些问题，前后两次迭代所得到的位移向量范数之比会出现剧烈跳动，以导致收敛不可靠。
取不平衡节点力为衡量收敛标准的准则称为平衡力准则，若满足下列条件，就认为迭代收敛：
式（1-117）中：P为外载荷向量；ΔPi为不平衡力向量；αp为不平衡力收敛容差。
上面公式中取哪一种范数，理论上可以任选。有学者认为，在用平衡力准则时，取
比较好；在用位移准则时，取
更为方便。在非线性比较严重的问题中，用位移准则更合适。有的学者还用能量（P）T（Δu）作为收敛标准，综合了力与位移两个方面，但要增加更多的计算量。
免责声明：以上内容源自网络，版权归原作者所有，如有侵犯您的原创版权请告知，我们将尽快删除相关内容。
