1、非线性有限元 第6章 求解方法和稳定性,计算固体力学,第6章 求解方法和稳定性,引言 显式方法 平衡解答和隐式时间积分 线性化 稳定性和连续方法 数值稳定性 材料稳定性,1 引言,描述非线性有限元的求解过程，瞬态问题的显式和隐式求解方法，以及平衡问题的解决方法，并且检验它们的编程和性质。展示了计算结果的稳定性、数值过程的稳定性和材料的稳定性。显式时间积分的中心差分方法，编程方法，相关技术如质量缩放、子循环和动态松弛。,以Newmark -方法为模型描述了隐式方法，静态问题的平衡求解。应用Newton方法求解离散方程，包括收敛性检验和线性搜索技术。在隐式系统解答中的临界步长和平衡问题是控制方程的
2、线性化。作为平衡方程的特殊情况，描述了运动方程的线性化过程。比较了显式和隐式方法以及评价了它们的相对优越性。描述了平衡方法的特殊技术，如连续性方法(参数法和弧长法)。,双曲线型偏微分方程，典型问题是波的传播,在双曲线型系统中，信息以有限的速度传播，波速为cx/t的直线斜率。 一个力在 t0 时刻施加在杆的左端，在右侧 x处的观察者直到波传播到该点时才有感觉。,2 显式方法,应力波传播的描述,2 显式方法,应力波传播的描述从概念上理解应用显式动力学算法时应力是如何在模型中传播的。在这个例子中，考虑应力波沿着一个由三个单元构成的杆件模型传播的过程，随着时间增量的变化，将研究杆件的各个状态。,自由端
3、作用有集中力P 的杆件的初始状态,自由端作用有集中力的杆件在第一时间增量段结束时的状态,2 显式方法,应力波传播的描述,杆件在第二个时间增量段开始时的状态,杆件在第三个时间增量段开始时的状态,高速动力学（high-speed dynamic）事件最初发展显式动力学方法是为了分析那些用隐式方法分析起来可能极端费时的高速动力学事件。例如，分析一块钢板在瞬时爆炸载荷下的响应。在迅速施加的巨大载荷作用下，结构的响应变化得非常快。对于捕获动力响应，精确地跟踪板内的应力波是非常重要的，由于应力波与系统的最高阶频率相关联，因此为了得到精确的解答需要许多小的时间增量。 复杂的接触（contact）问题应用显式
4、方法建立接触条件的公式要比应用隐式方法容易得多，因此能够比较容易地分析包括许多独立物体相互作用的复杂接触问题，特别适合于分析受冲击载荷并随后在结构内部发生复杂相互接触作用的结构的瞬间动态响应问题。,2 显式方法,适用求解的问题,复杂的后屈曲（postbuckling）问题不稳定的后屈曲问题。随着载荷的施加，结构的刚度会发生剧烈的变化。在后屈曲响应中常常包括接触相互作用的影响。 高度非线性的准静态（quasi-static）问题能够有效地解决某些在本质上是静态的问题。准静态过程模拟包括复杂接触的问题；如锻造、滚压和薄板成型等过程一般都属于这类问题。 材料退化（degradation）和失效（fa
5、ilure）在隐式程序中，材料的退化和失效常常导致严重的收敛困难，但是显式方法能够很好地模拟这类材料。混凝土开裂的模型是一个材料退化的例子，其拉伸裂缝导致了材料的刚度成为负值。金属的延性失效模型是一个材料失效的例子，其材料刚度能够退化并且一直降低到零，在这段时间中，单元从模型中被完全除掉。这些类型分析的问题都有可能包含温度和热传导的影响。,适用求解的问题,2 显式方法,中心差分方法,定义时间增量,关于速度的中心差分公式为,第n时间步的时间和位移,和,差分公式转化为数值积分公式,采用变化时间步的算法。在大多数实际计算中是必要的，当网格变形和由于应力而波速变化时，稳定时间步随之而改变。,2 显式方
6、法,中心差分方法,时间显式算法的图谱,2 显式方法,中心差分方法,在时间间隔的中点定义速度，称之为半步长或中点步长，,在等时间步长的情况下，上述公式简化为,是已知的关于函数二阶导数的中心差分公式。,考虑半离散运动方程的时间积分，在第n 时间步给出为,位移边界条件和在模型中的其它约束,I=1到nC,内部和外部节点力是节点位移和时间的函数。外部载荷通常描述成为时间的函数，也可能取决于结构的构形，是节点位移的函数，比如当压力施加在大变形的表面上。内部节点力对位移的依赖性是显而易见的：内部节点力取决于应力，而应力通过本构方程依赖于应变和应变率，而它们依次取决于位移和它们的导数。内部节点力也可以直接取决
7、于时间，例如：当温度设定为时间的函数时，则应力和相应的节点力也直接是时间的函数。,2 显式方法,中心差分方法,当质量矩阵M为对角阵时，实现节点速度和位移的更新不用求解任何方程。这是显式方法的一个突出特征：对离散动量方程的时间积分不需要求解任何方程，关键在于应用了对角质量矩阵。,2 显式方法,编程,在这个流程中的主要相关变量是速度和Cauchy应力。对于速度、Cauchy应力和所有材料状态参数，必须给出初始条件。一般假设初始位移为零（除了应力取决于变形历史的超弹性材料外，初始位移对于非线性分析是没有意义的，如果初始条件为位移，编程非常复杂，需要罚函数或Lagrange乘子法）。,1. 初始条件和
8、初始化设定v0，0和其它材料状态参数的初始值；d0=0， n=0， t=0；计算质量M， 2. 给出作用力,3. 计算加速度(阻尼不同时发生),2 显式方法,编程,4. 时间更新,5. 第一次部分更新节点速度,6. 强迫速度边界条件，在边界上节点I 处,7. 更新节点位移,8 给出作用力,9. 再次计算加速度,2 显式方法,编程,10. 第二次部分更新节点速度,11. 在第n+1时间步检查能量平衡,12. 更新增量步骤数目：,13. 输出；如果模拟没有完成，返回4, 更新时间。,与5. 第一次部分更新节点速度时的能量比较，检查能量平衡,2 显式方法,编程,子程序给出作用力,集合单元节点位移和速
9、度,ii. 当前步的内力初始值,iii. 对高斯积分点Q循环计算,1. 如果n=0，转到4,2. 计算变形的度量：,3. 通过本构方程计算应力,4. 在节点处集合内力,结束积分点循环,iv 计算单元外部节点力,v 合成单元节点力,vi 计算,如果,则,vii 离散,到整体,结束单元循环,计算总体外部节点力,2 显式方法,编程,流程的主要部分是计算节点力，子程序中的主要步骤；通过集合过程，从整体数组中提取单元的节点位移和速度。在单元的每个积分点上计算应变度量，通过本构方程计算应力。在整个单元域上，对B 矩阵和Cauchy应力的乘积进行积分，计算内部节点力，向节点汇集节点内力。离散单元的节点力进入
10、整体数组。在第一时间步中，没有计算应变度量和应力，直接从初始应力计算内部节点力。,2 显式方法,编程,显式方法的优点：方法的简单性和避免了方程的求解。容易编程，因此，显式时间积分是非常强健的，很少因为数值运算的失败而终止。显式积分的缺点：条件稳定性，如果时间步长超过了一个临界值,计算结果将会发散,线性系统的最大频率,单元e的特征长度,当前波速,非线性不稳定影响的折减系数,这种简单的稳定极限定义提供了某些直觉上的理解。稳定极限是当膨胀波通过由单元特征长度定义的距离时所需要的时间。如果已知单元最小尺寸和材料波速，就能够估算稳定极限。例如，最小单元尺寸是5 mm，膨胀波速是5000 m/s，稳定的时
11、间增量就在110-6 s的量级上。,2 显式方法,材料模型通过对膨胀波波速的限制作用来影响稳定极限。在线性材料中，波速是常数；所以，在分析过程中稳定极限的唯一变化来自于最小单元尺寸的变化。在非线性材料中，例如产生塑性的金属材料，当材料屈服和材料的刚度变化时波速发生变化。在整个分析过程中，应用在每个单元中的当前材料状态估算稳定性。在屈服之后刚度下降，减小了波速并因而相应地增加了稳定极限。,网格和材料对稳定时间增量的影响,下图给出了具有不同应力-应变曲线的平行联结的两根杆件，加载时保持相同的应变。当应变达到杆B的屈服点时（应变轴上的a点），杆件A已经经历了一定量的塑性应变。,若移去外载荷F，杆B的
12、应变将按照原加载线下降；但是杆A则循着一条平行于初始弹性线的轨迹下降。当A的应力下降到零时（b点），B仍然受到拉伸（c点）。为使二杆内力自相平衡，A和B二杆必须分别承受压应力和拉应力（两杆截面积相同），如e点和d点所示。,残余应力已留在杆件中，由于二杆并联为整体结构，故具有相同的拉应变。本例由于考虑了塑性变形，才使得两平行组合杆件，尽管应变相同，却具有大小相等方向相反的应力。,2 显式方法,临界时间步长随着网格的细划和材料刚度的增加而减小。显式模拟的计算成本是独立于频域的，它仅取决于模型的尺寸和时间步的数目。对于弹塑性材料，在塑性响应中减慢波速是否能够增加时间步长。基于经验，答案是否定的。弹塑
13、性材料可以在任何时刻卸载，由于数值的干扰在计算中常常发生卸载。在弹性卸载的过程中，临界时间步长取决于弹性波速，一个超过临界时间步长的时间步将导致不稳定。从单元时间步长得到网格时间步长。对于每一个单元，计算单元时间步长，并选择最小的步长作为网格时间步长。,网格和材料对稳定时间增量的影响,材料对稳定时间增量的影响,同样的波扩展分析在不同的材料中进行将需要不同的CPU时间，这取决于材料的波速。例如将材料从钢改为铝，波速将从5.15103 m/s变成,2 显式方法,由于刚度和密度几乎改变了相同的量级，所以从铝到钢对稳定时间增量只有微小的影响。在铅的情况下，差别则变得非常大，其波速减小为,这个值大约是钢
14、中波速的五分之一。铅棒的稳定时间增量是钢棒的5倍。,2 显式方法,能量平衡,从对线性运动方程积分的稳定性分析中，得到了上述稳定性条件。现在，还没有稳定性原理可以涵盖在工程问题中遇到的全部非线性现象，例如接触冲击，撕裂等等。即使满足临界时间步长，不稳定也可能进一步发展。对比线性问题，这里的不稳定将导致结果的指数性增长，因此不可轻视，非线性问题的不稳定解答有时不容易被识别。例如，Belytschko(1983)描述了一个数值现象，称为抑制失稳：由非线性诱发了不稳定，例如几何硬化，而材料是弹性的。这种不稳定引起了结果的局部指数增长，从而依次导致了塑性行为。而塑性反应软化了结构，并降低了波速，从而使积
15、分又得到稳定。,2 显式方法,能量平衡,在低阶方法中，如中心差分方法，一般是通过类似的阶数对时间积分得到能量，如梯形法则。内能和外能积分如下：,通过检验能量平衡可以察觉它们。在伪能量生成中，任何不稳定的结果将导致能量守恒的破坏。因此，在非线性计算中，通过检验能量平衡可以知道是否保持着稳定性。,动能给出为,2 显式方法,能量平衡,能量守恒要求,在单元或积分点水平，也可以计算内能,这里 是一个很小的容许极限，一般阶数为10-2,如果系统是非常庞大，在105节点量级或者更多，能量平衡必须在模型的子区域中进行，相邻子区域的内力可以视为每一个子区域的外力。,精确性,中心差分方法在时间上是2阶的，即在位移
16、上的截断误差具有 阶。,在L2范数的位移上，其空间误差具有h2 阶，这里h是单元尺寸。尽管在误差的两种度量之间存在一定的技术差别，其结果是相似的。由于时间步长和单元尺寸必须是相同阶数才能满足稳定条件，对于中心差分时间积分，时间积分误差和空间误差具有相同的阶数。然而，对于迅速改变硬度的材料，如粘塑性材料，中心差分方法的精确性有时是不够的。在这种情况下，建议采用Runge-Kutta法(4阶精度)。不需要所有的方程都应用Runge-Kutta法：它可以应用于本构方程，而通过中心差分方法对运动方程进行积分。,2 显式方法,2 显式方法,质量缩放、子循环和动态松弛,当模型中包含几个非常小或者比较硬的单
17、元时，显式积分的效率受到严重损害，因为整个网格的时间步长是根据这些小或硬的单元决定的。有几种技术可以避开这一困难：质量缩放：增加较硬单元的质量，使得时间步长不会因为这些单元而减小。质量缩放必须应用于高频效果不重要的问题。例如在金属薄板成型中，是一个静态过程，不会发生困难。 子循环：对于较硬的单元采用较小的时间步长。,在显式程序中，经常采用动态松弛以获得静态的解答。基本思路是：非常缓慢地施加外力，并应用足够的阻尼求解动态方程，从而使振荡最小。,质量缩放、子循环和动态松弛,由于质量密度影响稳定极限，在某些情况下，缩放质量密度能够潜在地提高分析的效率。例如，许多模型需要复杂的离散，因此有些区域包含着
18、控制稳定极限的非常小或者形状极差的单元。这些控制单元常常数量很少并且可能只存在于局部区域。通过仅增加这些控制单元的质量，就可以显著地增加稳定极限，而对模型的整体动力学行为的影响是可以忽略的。程序中的自动质量缩放功能，可以阻止这些有缺陷的单元对稳定极限的影响。质量缩放可以采用两种基本方法：直接定义一个缩放因子或者给那些质量需要缩放的单元逐个定义所需要的稳定时间增量。然而要非常小心，因为模型质量的显著变化可能会改变问题的物理模型。,2 显式方法,小结 应用中心差分方法对时间进行动力学显式积分，需要许多小的时间增量。因为不必同时求解联立方程，每个增量计算成本很低。随着模型尺寸的增加，显式方法比隐式方
19、法能够节省大量的计算成本。 稳定极限是能够用来前推动力学状态并仍保持精度的最大时间增量。随着材料刚度增加，稳定极限降低；随着材料密度的增加，稳定极限提高。对于单一材料的网格，稳定极限大致与最小单元的尺寸成比例。,2 显式方法,3 平衡解答和隐式时间积分,桁架问题的离散化模型,每个节点的分离图,对于隐式方法，这些平衡方程需要同时进行求解以获取每个节点的位移，最好采用矩阵算法求解。,显式与隐式相反，并不需要同时求解一套方程组或整体刚度阵。求解是通过动态方法从一个增量步前推到下一个增量步得到的。,3 平衡解答和隐式时间积分,离散动量方程,开关设置为,残数,对于运动方程的隐式更新和平衡方程，该离散方程
20、是节点位移,的非线性代数方程 。,在平衡问题中，残数对应于非平衡力；加速度可以忽略的问题为静态问题。上述的结果称为平衡点，并且结果的连续轨迹称为平衡分支或者平衡路径。,3 平衡解答和隐式时间积分,Newmark -方程,一个普遍应用的时间积分器Newmark -方法，更新的位移和速度为,应用参数 控制人为的粘性，即由于数值方法引进的阻尼，抑制结果的振荡,积分器无附加阻尼,积分器附加了,的人工阻尼比例,无条件稳定,离散动量方程,条件稳定,临界阻尼的分数,3 平衡解答和隐式时间积分,时间隐式算法的图谱,时间显式算法的图谱,3 平衡解答和隐式时间积分,Newton-Raphson方法,公式的求解是一
21、个迭代过程。迭代的次数由希腊字母的下角标表示,在时间步n+1上迭代 次的位移,开始迭代过程，必须选择未知量的初始值；通常选择上一步的结果,在动态解答中应用Newmark-方法，一个较好的初始值是,关于节点位移当前值的残数的Taylor展开，并设定计算的残数等于零，给出,略去比线性高阶的项，则公式给出一个关于,的线性方程,称为非线性方程的线性模型或线性化模型,3 平衡解答和隐式时间积分,Newton-Raphson方法,线性模型是非线性残差函数的正切；获得模型的过程称为线性化 。,对于位移增量，求解这个线性模型给出,在Newton过程中，通过迭代求解一系列线性模型，获得非线性方程的解答。在迭代的
22、每一步中，重写公式为,获得未知数的更新值，如图持续这一过程直到达到理想的精确度水平为止。,Jacobian矩阵,考虑平衡方程,定义残差,通过Newton迭代求解和线性化,内力Jacobian切线刚度,外力Jacobian载荷刚度,整体Jacobian为二者之差,ABAQUS的UMAT就是求内力的Jacobian矩阵切线刚度,3 平衡解答和隐式时间积分,3 平衡解答和隐式时间积分,Newton-Raphson方法,在计算力学中，称Jacobian矩阵为等效切线刚度矩阵，惯性力、内部和外部节点力的贡献被线性地分开。Newmark 积分器的Jacobian矩阵为,内部节点力的Jacobian矩阵称为
23、切线刚度矩阵 Kint,外部节点力的Jacobian矩阵称为载荷刚度矩阵，给出为,载荷刚度的名称相当奇怪，因为载荷不可能有刚度，名称是沿袭下来的,Jacobian矩阵记为,应用于将节点力的微分与节点位移的微分,3 平衡解答和隐式时间积分,切线刚度矩阵,内部节点力的Jacobian矩阵称为切线刚度矩阵 Kint ，分解为材料刚度和几何刚度。,不能应用于本构方程，因为它不具有客观性。,回顾,，因此有,则,材料,几何,3 平衡解答和隐式时间积分,材料刚度,几何刚度,切线刚度矩阵,材料,几何,TL格式的切线刚度矩阵组成为(证明见6.4.2和6.4.3节)：,对于UL格式的切线刚度矩阵，注意代换:,小结
24、：三种切线刚度矩阵，材料刚度，几何刚度，载荷刚度,3 平衡解答和隐式时间积分,Newton方法的编程,完全Newton算法，在每一次迭代过程中，Jacobian矩阵被赋值并求逆。修正Newton算法，Jacobian矩阵被组合和三角化，只在每一步的开始或者在中间赋值。例如，仅当迭代过程似乎不能很好地收敛时，或者在一个时间步开始迭代时，Jacobian矩阵可能被三角化。这些修正的算法更快捷，但是缺乏强健性。,隐式方法从施加初始条件开始，初始条件可以按照显式方法同样处理。通过迭代过程得到位移dn+1。开始迭代过程，需要d的初始值；通常应用前一时间步的结果。为了这个初始值则计算残数。在平衡解答中，残
25、数仅依赖于内部和外部的节点力。在瞬态隐式解答中，残数也依赖于加速度。,隐式程序框图,1. 初始条件和初始化设定v0，0 和其它材料状态参数的初始值；d0=0， n=0， t=0；计算质量M， 2. 得到,3. 计算初始加速度,3 平衡解答和隐式时间积分,4. 估计下一步解答：,或,5. 第n+1步的Newton迭代：,6. 更新位移，计数器和时间：,7. 检验能量平衡,8. 输出；如果模拟没有完成，返回到第4步。,a. 给出作用力计算,3 平衡解答和隐式时间积分,5. 第n+1步的Newton迭代：,隐式程序框图,见框6.1,b. 更新速度和加速度,c. 计算残值,d. 计算Jacobian
26、A(d),e. 对于基本边界条件，修正A(d),f. 解线性方程,g.,h. 检验收敛准则；如果没有满足，返回到第a步。,a. 给出作用力单元循环此部分与显式计算相同,集合单元节点位移和速度,ii. 当前步的内力初始值,iii. 对高斯积分点Q循环计算,1. 如果n=0，转到4,2. 计算变形的度量：,3. 通过本构方程计算应力,4. 在节点处集合内力,结束积分点循环,iv 计算单元外部节点力,v 合成单元节点力,vi 计算,如果,则,vii 离散,到整体,结束单元循环,隐式程序框图,3 平衡解答和隐式时间积分,保守问题,通过设势能的导数等于零，得到平衡解答,当势能为局部最小值时，一个平衡点是
27、稳定的。因此，通过势能W最小化，可以找到稳定平衡解答。通过数学程序和优化技术，以这种形式得到平衡解答，如最速下降法。,如果希望得到稳定和非稳定解答，必须找到W(d)的驻值点，则离散问题成为,找到 d,根据,3 平衡解答和隐式时间积分,约束,罚函数法 Lagrange乘子法 Lagrangian增广法 Lagrangian摄动法,位移边界条件和其它约束是节点位移的线性或非线性代数方程：,将约束附加到采用Lagrange乘子的目标函数上。保守问题中的目标函数是势能，即势能函数最小化。根据约束对一个函数最小化，可以形成一个Lagrange乘子问题：函数的最小化对应于函数的驻值点和由Lagrange乘
28、子权重的约束，解答是与寻找驻值点等价的,Lagrange乘子,Lagrange乘子法,3 平衡解答和隐式时间积分,Lagrange乘子法,WL表示Lagrange乘子修正的势能，在平衡点上有,注意到稳定平衡点对应于d(W)的最小值和对应于 的最大值，即鞍点。对应于d(W) 和 ，在驻值点处的导数为零，因此：,I=1 to nc,上面共有,个代数方程组成的系统；方程重写为,约束引入了附加力是Lagrange乘子的线性组合。如果约束是线性的，附加力将独立于节点位移。,3 平衡解答和隐式时间积分,Lagrange乘子法,为了得到的线性模型，取其Taylor展开，并将结果设置为零给出,注意到对重复指标
29、求和。为了将此写成矩阵标记，定义,线性模型成为,可以看到，由于约束，线性模型具有nc个附加方程。即使当矩阵A是正定时，方程的扩展系统将可能是非正定的，因为在矩阵右下角的对角线上有零元素。,3 平衡解答和隐式时间积分,例6.1,考虑一段长为,并铰接在原点O上的非线性弹性杆OA。应用超弹性，Kirchhoff材料模型。该杆受常数外力fext（保守荷载）作用，并从A点拉伸到A点。杆的右端被限定搁置在圆形槽中。,(a) 单轴向应变状态写成为,并计算E11,(b) 建立非约束问题的势能及它的一阶导数(残数)和二阶导数(刚度)。(c) 建立Lagrange乘子法公式。,(a) Green应变,(b) 总势
30、能给出为,W对应于ux1和uy1的偏微分给出残数r,3 平衡解答和隐式时间积分,例6.1,对于保守荷载,(c) 借助于Lagrange乘子施加约束。给出修正势能为,g=0是约束，要求节点1保持在圆心(2a,0)处，以节点坐标的形式，约束为,3 平衡解答和隐式时间积分,例6.1,残数对应于ux1和uy1的偏微分给出切线刚度,由于初始条件 Y1=0和X1=a，所以，约束的梯度是,由公式(6.3.40)给出线性化方程，其中矩阵A和G的定义同上，而H给出为,Lagrange乘子 对应于约束g=0。,3 平衡解答和隐式时间积分,例6.1,3 平衡解答和隐式时间积分,收敛准则,对于隐式求解，在Newton
31、算法中终止迭代是由收敛准则决定的。这些准则适用于迭代求解,应用三种l2范数收敛准则终止迭代： 1根据残数r的量级的准则。,3能量误差准则(残数r是力)。,2根据位移增量d的量级的准则。,4 线性化,传统Newton方法应用于固体力学问题的困难：,1 对于弹塑性材料，节点力不是节点位移的连续可微函数。,2 对于路径相关材料，迭代过程的中间结果不是真实加载路径。,3 对于大变形和转动，线性化增量引入误差。,4 导数和积分需要离散，离散形式的切线取决于增量步长。,对Newton方法做出修正：,1 应用方向导数建立切线刚度。,2 应用正割方法，取代切线刚度，用最后收敛解答作为迭代点。,3 采用依赖于增
32、量大小的公式，联系力与位移的增量。,目的是寻求本构积分中增量的线性刚度。,承受载荷的二杆桁架，应力均为相等压力，达到屈服，考虑材料非线性。当施加载荷增量，切线刚度矩阵取决于增量位移，内部节点力的改变依赖于位移增量，其残数不是增量节点位移的连续可微函数。在Jacobian中有4条不连续的线段。,原因是如果位移增量导致了拉伸应变增量，则压杆将弹性卸载，切线模量从塑性模量Hp(连续加载)变化到弹性模量E(卸载)。,4 线性化,切线刚度矩阵,当解答不平滑时，必须应用算法模量，称为算法的一致切线刚度，适用于在本构方程中缺少平顺性和分区线性。,象限1,象限2,象限3,象限4,位移增量的切线刚度：图中代表了
33、不同节点力行为的区域，显然，用一个标准的导数(切线刚度)无法计算，四个象限中关于位移增量的切线刚度分别为：,4 线性化,隐式方法的精度和稳定性 隐式优于显式是对于线性瞬态问题，隐式积分器是无条件稳定的。尽管对于特殊情况的理论结果表明，至少对于一些非线性系统，无条件稳定是适用的，经验表明隐式积分器的时间步长远远大于显式积分器的时间步长，并没有导致失稳。在隐式方法中，对于时间步长的主要限制来自于对精度的要求，以及当时间步长增加时降低了Newton过程的强健性。Newmark方法是二阶精确的，与中心差分方法是同阶的。较大的时间步长也削弱了Newton方法的收敛性，特别是在一些有非常粗糙响应的问题中，
34、如接触冲击。应用较大的时间步长，初始迭代可能远远偏离解答，因此，增加了Newton方法收敛失败的可能性。较小步长则改善了Newton算法的强健性。作为增强稳定性的报答，隐式计算付出高昂的代价：需要在每一个时间步求解非线性代数方程，存储这些方程需要大量的内存。,隐式小结,5 稳定性和连续方法,在非线性问题中，需要考虑解答的稳定性。稳定性是一个概念，取决于观察者和他的目的，如数值稳定性问题：,首先给出稳定性定义并探索它的内涵。考虑一个由演化方程控制的过程，诸如运动方程或热传导。对于初始条件,设结果用,表示。现在考虑对于初始条件,的结果,在某些范数下,如果对于所有初始条件满足公式，解答是稳定的，则结
35、果满足,初始条件位于围绕,的球上，如结果 位于包围结果,的球上，则结果是稳定的，否则是不稳定的。,稳定性,5 稳定性和连续方法,收敛性是用范数度量的。L2范数：度量位移的精度,范数,能量误差：度量应变和应力的精度,能量的误差一般低于L2范数的误差一个数量级。,5 稳定性和连续方法,稳定性,在下面的例子中将解释这一定义与稳定性直观概念的关系。考虑一根梁轴向加载的过程，如图所示。通过一个 或者2 的距离，摄动加载位置，则平衡路径如图(c)所示。可以看到当荷载低于屈曲荷载时，对于不同初始条件的路径保持着与AC接近。然而，当荷载超过屈曲荷载时，结果将分岔开，任何过程都是不稳定的。尽管没有在图中显示，不
36、稳定分支的方向依赖于初始缺陷的迹象；图中只显示了一个方向。,5 稳定性和连续方法,稳定性,当梁十分笔直的时候，数值计算结果一般保持在路径AC上，如图(b)所示。即使当荷载超出屈曲荷载时横向位移仍然是零。在有限元模拟中，在增量静态求解或者动态求解中，无论应用显式还是隐式算法，一根笔直的梁一般将不发生屈曲。仅当舍入误差引起的数值缺陷或者当数据中引入缺陷，将导致笔直梁在模拟中屈曲，即需要有缺陷以破坏对称性，否则“压不垮”。,5 稳定性和连续方法,稳定性,一个系统的状态稳定与过程稳定多少有所不同。通过检验应用于平衡状态的摄动增长，确定一个平衡状态的稳定性。其结果不如过程稳定的概念直观，例如图中的梁，平
37、衡状态的任何扰动在分支 AB和 BD上不增长，即这些分支上的任意平衡状态是稳定的(在例6.4中说明)。另一方面，在分支BC上平衡结果的任何扰动可以看到增长，即它是不稳定的。在B点上，分支AC从稳定变为不稳定，称为临界点，或者是屈曲点。对应于Euler梁的屈曲荷载音叉。,5 稳定性和连续方法,稳定性,这种方法也应用于研究在管中液体流动的稳定性。当流速低于临界Reynolds数时，流动是稳定的，流动的扰动导致小的变化。当流速超过临界Reynolds数时，由于流动从层流变为湍流，小的扰动将导致大的改变。因此，高于临界雷诺数的流动是不稳定的。在动力学中，经常给出稳定和不稳定的例子，如图所示。状态 A
38、是稳定的，因为球的位置的小扰动将不能显著地改变系统的状态。状态 B 是不稳定的，因为任何一个小扰动将导致很大的变化：球可能向左或右滚动。在教材中，状态 C 通常称为中性稳定。由范数定义，状态 C 是不稳定的，因为速度的小变化将导致位置在长时间内的大变化，尽管不像状态 B 那样大。,5 稳定性和连续方法,平衡解答的分支,为了更好地理解一个系统，必须确定它的平衡路径或者分支，以及它们的稳定性。在结构力学中，通过直接获得动态解答可以回避与不稳定行为相关的困难，如关于位移摄动的系数行列式为零的非平凡解答。当结构加载超过它的临界点或在动态模拟中的分支点时，结构动态地通过了最接近的稳定分支。然而，不稳定是
39、不容易出现的，并且缺陷敏感性的可能性也是不清晰的。因此，为了理解结构的行为或者整个过程，必须小心地检查它的平衡行为。,5 稳定性和连续方法,平衡解答的分支,一个具有反对称分支的系统对于缺陷是非常敏感的，在图(e)中，最大荷载随着缺陷变化。在实际结构中缺陷是不可避免的，一个完好结构的理论分岔点不是强度的真正度量；实际结构可以在一个远远低于理论值的荷载下屈曲。,许多结构行为的奇怪表现可能被动态模拟掩盖了。例如在承受轴向荷载的圆柱壳中，交叉分支是反对称的，如图(e)所示。,5 稳定性和连续方法,平衡解答的分支,为了达到确定非线性系统性质的目的，通常将平衡解答分成若干分支，当一个参数变化时，它是连续地
40、线性描述系统的响应。称这些分支为平衡分支。追踪模型的平衡分支作为参数 的函数。,寻找,使得残数为零（即系统是处于平衡的）,非线性系统显示了三种分岔行为： 1 转折点：结构分析中的临界点，在该点上分支斜率改变符号。 2 静止分岔：简单分岔，在这些点上两个平衡分支交叉。 3 Hopf分岔：在这些点上平衡分支与一个周期运动的分支交叉。,5 稳定性和连续方法,平衡解答的分支,图示浅桁架行为展示了一个转折点（或者临界点）；A和B是转折点。在转折点之后，一个分支可以是稳定或者是不稳定的。在这种情况下，在例6.4会看到，第一个转折点A 后面的分支是不稳定的，而第二个转折点B 后面的分支是稳定的。(例6.4中
41、字母有区别),5 稳定性和连续方法,平衡解答的分支,在图a中显示的梁是分岔的典型例题。两个分支的交叉点B是分岔点。在分岔点后，基本分支的继续部分BC成为不稳定的。分岔点B对应于Euler梁的屈曲荷载。这种类型的分支通常称为音叉。Hopf分岔在率无关结构中不常见，发生在率相关材料或者主动控制的问题。在Hopf分岔中，在一个分支的终点不可能得到稳定平衡的结果。代之是这个平衡分支分成为两个分支，其结果对时间是周期性变化的。,5 稳定性和连续方法,连续方法和弧长法,跟踪平衡分支的方法称为连续方法。平衡分支的跟踪是相当困难的；对于连续性，仍然没有实现强健和自动化的过程。下面描述基于参数化的连续方法-弧长
42、法。,如果采用载荷增量求解，希望在A附近收敛，有可能跳到另一个替代平衡点A；同样情况也适用于位移增量求解，有可能从B跳到B。如何跟踪真正的平衡轨迹，采用弧长法。,5 稳定性和连续方法,连续方法和弧长法,在跟踪分支时，荷载参数通常从零开始并逐渐增加。对于参数 的每一个增量，计算平衡方程解答，即找到关于方程的解。,或者,在弧长法中，除了逐渐增加荷载参数 之外，在位移荷载参数空间中的弧长度量也需要逐渐增加，这是通过在平衡方程中增加参数化方程实现的。,弧长增量,载荷增量,首先讨论如何应用荷载参数跟踪分支。,位置变化,5 稳定性和连续方法,连续方法和弧长法,在步骤中弧长被横切的近似值,载荷增量比例因数,
43、一个可能的比例因数与材料刚度矩阵的对角线单元有关，即,产生许多其它类型的参数化方程。方程的整个系统包括平衡方程和,弧长方程,参数化方程,(1),(2),5 稳定性和连续方法,连续方法和弧长法,这样，有了一个附加的方程和一个附加的未知数,使得弧长s(圆的半径)是逐渐增加的。弧长法既不是力加载，也不是位移加载，而是载荷位移混合加载。,假设节点外部荷载的分布不随模型的变形而改变；当荷载因子是从属荷载的系数时，如压力，根据荷载几何的影响，必须修改方法以考虑到荷载分布的变化。上述不包括惯性项，因为连续方法仅适用于平衡问题。在结构力学中，弧长法经常称为Riks方法（Riks,1972）,对于一个自由度的问
44、题，这一过程是最容易解释的。图示浅桁架，假设在点n得到平衡解答。在荷载位移平面中，弧长方程(1)是一个围绕点n的圆。参数方程(2)的解答是与该圆交叉的平衡分支。在点n增量荷载参数是无效的，因为它又带回到分支上。 在弧长法中，以沿分支的弧长形式重新表示问题，因此一个分支可能跟随一个下降的荷载。在步骤中不需要增加荷载，而仅需要增加弧长参数，这是跟踪分支的完美的自然方式。在两个自由度的问题中，弧长方程可以围绕平衡点定义一个球面或球，以及分支与球交叉的一个解答。,5 稳定性和连续方法,连续方法和弧长法,5 稳定性和连续方法,连续方法和弧长法,关于对称桁架的参数方程可以设置如下：找到如下方程一个解答,根
45、据,在上面方程中，d 是竖向位移，r 是对应残数；假设水平位移为零。另外，可以将上述方程以位移和荷载参数增量的形式写出：找到一个解答根据 r = 0 得到,因此，对于含一个未知数的一组原始离散方程可以通过一个方程扩展，并增加一个未知数 。方程的解答可以通过标准Newton方法得到，设,正确的结果通常是,的一个最大值，即相同方向上前一步的解答,5 稳定性和连续方法,例6.4,具有稳定和不稳定路径问题的一个简单例子是铰接浅桁架，如图所示。单元的初始横截面面积为A0，两个单元的初始长度为l0，,竖向荷载 p 施加在节点1上。由于这是唯一荷载，设 p为荷载参数。材料服从Kirchhoff法则，,因为系
46、统是保守的，通过寻找势能的驻值点，将确定平衡路径。通过线性稳定性分析，检验分支的稳定性。,5 稳定性和连续方法,例 6.4,通过中心节点的当前竖向坐标 y1，描述桁架的变形。与以位移形式的方程相比，它导致了简单的方程。势能成为,两个单元的Green应变给出,两个单元的内能给出为,式中,组合上式和外力势能得到整体势能为,5 稳定性和连续方法,例 6.4,应用势能驻值原理得到平衡方程。令上式的导数为零得到,节点力是中心点竖向位置的三次函数。平衡分支有两个转折点B和C，三个分支分别由AB、BC和CD表示。通过求解关于状态 y1的摄动和线性化的运动方程，检查分支的稳定性。通常应用切线刚度矩阵获得线性化
47、方程，在本例中，只是简单地对节点力求导数。令节点1的位置为y1；按照公式以 y1的形式，给出荷载,5 稳定性和连续方法,例 6.4,摄动解答为,式中 是一个小参数。运动方程给出为,替换了,展开上述方程，并舍弃在,中高于线性的项，得到摄动的运动方程,由于y1是平衡状态，荷载p抵消了上式括号中的第一项，运动方程成为,5 稳定性和连续方法,例 6.4,将公式代入上式，服从,对于,在式 中的一个参数 是实数且为正，于是摄动解答将增长。因此由,所定义的分支BC是不稳定的,对于,的任何其它值，系数 是虚数，因此，摄动解答是具有常数幅值 的调和函数，并且平衡点是线性稳定的。,5 稳定性和连续方法,例 6.4
48、,上述稳定性分析的结果可以直接由势能函数的二阶导数得到，从公式给出,检验上述Jacobian是否为正定,当,其它情况,因此，对于其它通过摄动分析得到的结果，稳定条件是一致的。总之，平衡分支AB和CD是稳定的，而分支BC是不稳定的。,RIKS算法对于塑性屈曲和后屈曲问题，需要通过跟踪加载过程，由荷载-变形曲线识别其屈曲荷载. 对这类问题采用有限元进行分析时，常用的荷载增量法不适用，如前所述ABC段，增量可以不断增加，而没有侧向变形迹象。采用在有限元软件中应用的RIKS法。在采用 RIKS法进行屈曲或后屈曲分析时，需给结构引入适当的初始扰动，如结构的几何或材料的初始缺陷，可以采用实际结构的缺陷，如
49、制造、安装误差等。在结构设计阶段，这难以做到，常用的方法是采用弹性屈曲模态的线性组合，作为假想的初始缺陷，见欧洲钢壳结构稳定性设计标准。,5 稳定性和连续方法,应用实例,对于内压作用下的压力容器，由弹性屈曲分析的结果显示在内压作用下结构不失稳，故不能引入弹性屈曲模态的线性组合作为假想的初始缺陷。这里采用过渡段周向受压区局部引入与静力位移相关量作为初始几何缺陷，分析荷载 变形历程，识别屈曲荷载及后屈曲。,5 稳定性和连续方法,应用实例,如图所示蝶形封头压力容器，由球壳(R1=402), 柱壳(R2=201), 过渡环壳(母线半径R3=30), 厚度(t=0.8), 内压2.0MPa。弹性和弹塑性分析结果表明，在内压作用下，在过渡区局部向内瘪进，证明该区域有局部压应力，其余大部分区域向外鼓胀。因此，不会因为内压引起大部分区域失稳，只有外压作用才会失稳。,5 稳定性和连续方法,应用实例,在局部区域周向受压，出现周向稳定性问题，在内压达到一定值后过渡区可能形成周向皱褶。引入与静力位移相关的几何初始缺陷, 预制局部突起，根据弹性静力计算所得到变形量确定径向突起量,
