非线性方程组的解法

第四章非线性方程组的解法4・1非线性方程组的一般形式从上面两章中，我们研究了离散化结构中任一单元在t-1+At的时间增量步内，由材料非线性引起的单元切线刚度阵是线性的，(如第三章得出的增量平衡方程kAq=Ap⑺1(假定t时刻的状态已知))，由此集合而成结构的增量平衡方程也是线性的K村=Ap，这就为求解整个的非线性过程准备了条件。即只要确定每一步的切线刚度，通过求解一系列的线性方程组，累加起来就得到了解的全过程。结构总的平衡方程是非线性的：K(q)q=P (1)i.eq=K-1P。令R=K(q)q(1)—F=P-R(q)=0 (1)’分段线性化是求解非线性问题中一个普遍有效的技术，但作为具体的解法还有许多种，主要的有：1、 增量法一纯增量法2、 迭代法一直接迭代法(刚线刚度法)、Newton-Raphson迭代法(切线刚度法)3、 .混合法一增量/迭代型方法4・2载荷增量法(纯增量法)1、基本思想方法：若将外载荷分成N个增量步，A为载荷系数i方法：若将外载荷分成N个增量步，A为载荷系数i而每个增量载荷为Ap=A\%，q1 q2 q3(或称载荷因子)，则总载荷p=xP；0其中：M=U牌p为基准载荷.i=1上面的结构平衡方程为F(q)=P-R(q)=0 (1)'i.e F(q)=XP-R(q)=0 (2)

上式两边对人微分得i.e上式两边对人微分得i.edF(4)dq=K-1(q以。P-K)qdPdq=K-1(q以。P-K)qdP..Ae将(5)式写成增量形式便有以下求解格式Aq=[q(q.1)]-iAP

q=q+AqAP=AXP(5)(6)2、求解步骤将载荷分成若十个增量步P=歹备匚，准备位移量累加器［Q］并置零.i=1施加第1个载荷增量AP=AP，计算k(q)=竺线性1 10 t0 6q求解Aq1=［q(q0)］-1气q1=Aq1 并送入位移量累加器［Q］施加第2个增量步AP=AX2P0用q，求K(q)即在q处的切线刚度矩阵1 T1 1求解Aq2=［q(q1)］-1AP2q2=q1+Aq2 在位移量累加器［Q］中完成累加.4)重复3)直至(N)个载荷施加完毕，在位移量累加器［Q］中得到总位移q=寸Aqi=1入OP0P3P△2入OP0P3P△2入2P1△入1q1 q2q3优缺点：优点：了解加载过程，当APT足够小，总能收敛到真实解缺点：实际不可能无限小，因此累积误差，且无法估计，造成极大偏离而失真4・3迭代法1直接迭代法1）基本思想：将载荷一次加上，并假设一个初始解代入方程组求出第一次近似解；将其再代入方程组求解，得出第二次近似解，反复迭代逐次修正解，直至满足方程组（类似于对过渡单元加权平均。任中m的迭代）。2） 步骤①假设近似解q0；代入方程（1）得q⑴二]K。（0），P；依次类推得q（k）=[K（q（T，TP （7）3）几何意义（一维问题） 迭代过程是调整其平行割线刚度的过程。4） 优缺点求解方法简单，对原有弹性分析程序稍加修改即可，适用于非线性弹性。但迭代效率低，对一些问题不收敛。举例

2牛顿一莱夫森迭代法（Newton—Raphson）1）基本思想对非线性方程（1）函数F（q）在某一近似解处qG）作一阶台劳展开四G=四G=四侦))4电］(q-q"—'-qG))+O侦)=0(8)为当前迭代解与准确解之间的差值。由（8）得Aq为当前迭代解与准确解之间的差值。由（8）得Aq(j)=-(10)代入（1）两边求导得(dF(q))I6q设外力（大小、方向）与位移无关(dF(q))I6q设外力（大小、方向）与位移无关Jq=q(/)由(6)得K(q(j))=-TW()q=q(j)(11)（12）式代回（10）、（11）得Aq(j)3®jjq=q(j)(12)=(k(q)A(OT(13)于是得到迭代公式于是得到迭代公式(14)(14)2）N2）N-R法迭代步骤①取初始状态｛q。｝=0故R{q0}=0计算Klq0计算Klq0求解方程（8）代入（9）{q0}=K0即线弹性刚度矩阵(与当前位移无关)TAqJ=K0(q0)-1(p}-R0(q0位K0曷T Tq1技0挫q1}={Aq1},.,线性解由{q1}计算K111）和&2K「'1Md从而得到q2}=K「&1）11板以1（q1）》和£2技1M2}从而得到③重复2③重复2的做法^,}=K。一11《}-&一1）》=K。一11z。一1）T T"顷1》XU0}+EL}j3）几何意义（以一维问题为例）按基本公式（12），（13），（14）凶』为p-q曲线上，在£j》对应点的斜率，也就是切线刚度；{p}为当前状态下的外载向量；b、{。）》=虹［j】［j}=fBTDBdv{qj}{r}称为当前抗力，{p}-{r}={z}称为当前的不平衡£}力；所以，每一次迭代都是在上一次迭代终（当前）的变形，应力状态下，形成新的切线刚度矩阵作为求下一次迭代解的切线刚度，用当前的不平衡力求解。作为对上一次解的修正，不断修正，使不平衡力越来越小，最后达到平衡。4）优缺点及适用范围从上面迭代原理可知：该方法使变刚度法，计算工作量较大，但是收敛速度比较快，而且收敛性比较好，但是有时特例也不一定能保证收敛。为了改进N-R方法，克服计算量大的缺点，产生了等刚度法3修正的N-R迭代法（等刚度法）顾名思义，等刚度就是在迭代过程中取一个（相等得）切线刚度矩阵，如果取初始的K0T或某一个值，取KkTi.e.取KjTKjTKo=K(q)T（j=1，2，。。。）hqj+1}=「Ko（q）］-1（{p}-{R（qj）》）1 T作了这一改进之后，只需在开始第一步迭代时第k步计算各单元的单刚，组装一次总刚后，作为以后每一次迭代的总刚而无需每次重新计算重新组装，只需计算由当前位移

R(qj)引起的结构反力，然后回代求解相应的位移增量。分段等刚度迭代这一改进明显减少了计算工作量，但收敛速度必然会放慢，总的效果有时是好的。以后又有人作了进一步的改进，在等刚度迭代几次以后，改用一次变刚度，形成新的刚度阵再作等刚度迭代。分段等刚度迭代等刚度迭代总结N-R迭代法&。)｝七侦分妇小S)｝LT」(10)(6)(15)qg上q(j顼挫qg(10)(6)(15)纯增量法(\q｝=K(q)L曷｝i Ti-1 i上｝=妃｝+明｝将纯增量法和N-R迭代法结合起来：[LLK。-』1UR。-1)》ji 3｝+£Mi〔 i i-1，'4.4混合法1基本思想结合增量法保证收敛的优点和N-R迭代法不断校正解的优点，可取得好的效果。即首先将总载荷分成若十个增量步，在每个增量步内采用N-R迭代。在一个增量步内达到平衡之后，再进入下一个增量步迭代。下面以第i个增量步内的迭代为例说明混合法的全过程。

2求解步骤1)开始，施加第一个增量载荷，按N-R迭代法求解。设10}=0，求K(/0『，求解方程{\q，}=K(/0HG-R(0»，{q"=(^q"，线性解。以(q"计算K(q，)] 假如迭代了k次收敛，则(q}i={q}0+2L}j=12)在第i个增量步内的迭代(i>2),各个增量步内的迭代过程相同，已知第G-1)步终达到平衡状态(q，1}{.Jb.Ji-1求%I?)R(q?)i-1求得(qQ)}和(qQ)L(q}+(qQ)}TOC\o"1-5"\h\zi i i1 i(1)以(q}为i迭代步的初位移-(qdi(2)解方程dQ)}=K(/0MG-r(/0》i Ti(1)以(q}为i迭代步的初位移-(qdi(2)解方程dQ)}=K(/0MG-r(/0》i Ti-1i i-1(3)以此类推，求出(q(3)}，((q(4)}，……即(q(3)}，(q(4)}，……迭代到k次收敛i i i i最终达到平衡时(q^k)}=qL・}+2kqk}j=1ii4)直至全部载荷施加完毕。3几何意义3)以第i步终达到平衡状态时的(q}{}和(7}作为初始状态，进入第(+1)步的迭代。ii4)直至全部载荷施加完毕。3几何意义4优缺点优点：通常能保证收敛，且收敛速度快。目前应用最广泛。缺点：对一些特殊问题失效，主要是有极值点的问题。4.5迭代收敛准则和增量步的选取1收敛准则(ConvergentCriterion)用迭代法求解非线性方程时，给出一个收敛判据（准则）是非常必要的，否则迭代将无法终止，实践证明，收敛准则将直接影响求解精度和速度，如果准则选的不合适将导致计算失效，目前常用的有下面三种判据：（1）位移准则（见AIAA；10（8）.1982）（a）1节 Aq=——乙 k£(b)i（a）1节 Aq=——乙 k£(b)iNk=14.ref£2冬'k.ref(41)(42)(c)£=max(43)最常用的是位移收敛准则，£取10-2〜10-5,只有当结构出现严重硬化时，位移增量的微小变化会引起不平衡力的很大变化，这一准则不予收敛。（2）不平衡力准则zj={p}―"侦-1）j<£{p}| （44）当物体严重软化时或材料接近于理想塑性时，平衡力的微小变化又将引起位移增量的很大偏差，这时又不宜用不平衡力准则。（3）能量准则能量已同时考虑了位移和不平衡力，能量准则是比较好的判据，它是把每次迭代后的内能增量（不平衡力在位移增量上做功）与初始能量增量做比较。i.e.(45)(46)膈}S-"侦位&1*Gp}4侦）}）i.e.(45)(46)TOC\o"1-5"\h\z，^-j^ ，^r^i ii i采用弧长法时可写作：Dlj{p}^Dqj}1{p》{q}'■'i i-10ai

2增量步长的选择求解非线性方程组时要合理选择增量步长，步长过大使计算结果不会收敛或不可靠，步长过小使计算时间太长，因此要根据具体问题的非线性程度来选择步长，不同结构形式的非线性形态不同，同一物体不同的加载方式、载荷大小不同，非线性程度也不一样。P.G.Bergan等提出的方法设第/步刚度度量S*二"V} (i=1,2,...) (47)^ {贤尸{△〃,}则初始刚度S*={Api}"V} (48)o {△”{△〃]}Ap=p kq=q(49)于是，第i增量步的刚度参数用下面无量纲参数定义(49)S=S*/S*由于{由于{p}=人{p}，0叫=&{p0}c{Ap}T{Ap}{Aq}T{Ap}S= i i 1 1—i{Aq}t{Ap}{Ap}t{Ap}ii11(50)加{p}『{p}{q}m{应［{Mh={Aqk