第九章二维梁的有限变形和算法(徐春晖、李明瑞)

第九章 二维梁的有限变形和算法 §1 引论 实际存在的结构都是三维体结构。梁是力学上一种简化的结构模型。当三维体的其 中某一维尺寸远大于其他二维尺寸就可将其简化为一维结构。如果该一维结构的受力变形 形式为拉压。这就是我们在上一章中研究过的二力杆。如果该一维结构的受力变形形式为 弯曲，这种结构构件就是梁。所以要建立梁的模型应符合两个条件：10 几何条件：细而长 的构件。20 受力变形形式：以弯曲变形为主。 在力学简化构件中还一类很特殊的构件称之为索，它的简化模型的条件是：10 几何 条件：细而长的构件。20 只受拉，不受压，不受弯。以著名的美国金门桥的拉索结构为 例，它的断面尺寸为直径36英寸，全长7650英尺。显然，说它不受弯的条件是过分勉强 了。这是由于当时对梁的大变形问题不会分析，不得已而简化成索来处理。一般情况下， 将所谓的索当成大变形的梁是完全正确的，也提高了计算的精度。 建立简化的力学结构模型的好处主要在于将三维体的力学分析简化为一条一维曲线 的力学分析。这条曲线就是梁的中轴线。从三维体的有限变形观点来看，不管哪一种梁的 大变形理论只能是一种近似理论。但是在简化结构模型的观点看却完全可能建立起梁的精 确的大变形理论，这就梁的有限变形理论。就梁在受力变形后的曲线形式看，如果是一条 二维曲线，就称之为二维梁。如果是一条三维曲线，就称之为三维梁。我们知道在连续体 力学模型中，一个空间点的初始位置R由三它的三个坐标决定。该点在变形后的位置r，即 当前位置的三个坐标，则可由该点的初始位置的三个坐标再加上该点的三个线位移u决 定。如图1所示。 1 u 图1 三维连续介质中的任一点在变形前后的位置与位移的关系 r R r = R + u 。 (9.1) 由于一个点可以产生三个线位移，所以也称一个点有三个自由度。 梁作为一种简化的结构模型，与连续体模型的最大不同处在于引进了转动自由度的 概念。这个概念在三维梁的线性理论与有限变形理论中有较大差异。由于本章只讨论二维 梁的问题，所以只对二维梁的转动自由度作详细研究。对于转动自由度的更一般性的研究 则将在三维梁与板壳的有限变形理论中去进行。 梁的变形后的位形由两部分来描述。第一部分是梁的中轴线的变形。这部分的变形 与连续体的变形的描述完全相同。第二部分是梁中不在中轴线上的点在变形后的位置的描 述。如图2 a. 所示为一个梁段的初始位形。取梁段的轴线为x轴。设A0为变形前轴线上任意 一点，线段A0BB0垂直于轴线，实质上它代 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 了一个在初始位形中垂直于轴线的横截面。 x y A0 B0 ϕ ϕ γ θ A1 λ B1 a. t0时刻的位形(初始位形) b. t时刻的位形(当前位形) 图2 Timoshenko梁段变形示意图 2 3 BB0点就是要研究的梁中不在中轴线上的一个任意点。一般而言，二维梁可以发生拉压、弯 曲、剪切的变形形式。由于拉压有限变形已在前面的杆单元的有限变形中作了研究，本章 中将主要研究梁的弯曲、剪切的变形形式。图2 b.表示了梁段在变形后的位形。此时线段 A0B0B 在变形后的映像为线段A1BB1。如果在纯弯曲的情况下线段A1B1B 将依旧垂直于变形后的 轴线，而在一般情况下由于剪切的存在，A1BB1将不再垂直于变形后的轴线。设在A1点作变 形后轴线的切线λ.。令切线与x轴的夹角为ϕ ，线段A1B1B 与y轴的夹角为θ，由于剪切的存 在，如图2有 ϕ = θ + γ 。θ 是线段A1BB1的转角，也就是在初始位形中垂直于轴线的横截面 的转角。γ 是剪切角。为了简化以上的说法，称A0点有线位移A0A1和角位移θ。当轴线上每 一点的线位移和角位移知道后，整个梁段的变形也就完全确定。因此在二维梁的理论中， 我们可以主要研究中轴线的变形，而且认为中轴线上的每一个点都有3个自由度，即2个线 位移与1个角位移。从严格的连续介质力学理论来说，一个点只有位置没有大小，所以只 能有线位移而没有转动。现在对于梁的模型中轴线上的点给与转动自由度，这个转动自由 度实质上表示的是垂直于中轴线的截面的转角，而不是中轴线上的点具有任何特殊性，产 生了转动。所以中轴线上的点与一般的连续介质力学的点的在概念上并无不同，还是只有 线位移自由度。但为了简化说法，我们以后可以约定成俗地说梁的中轴线有线位移与转角 自由度。对这种不严格的简化说法我们应作上述的正确理解就行了。在以上的变形分析中 隐含着如下两条非常重要的基本假设假设： 10 变形前的直线段A0BB0在变形后仍为直线段。这就是梁的平截面假设。我们在有限变形理 论中仍然采用。实际上，由于剪切的存在，截面将发生翘曲。如此，仅靠中轴线的线位移 与角位移就无法确定梁中不在中轴线上点的位移，也就无法确定梁的整体位形，而必须当 作连续体来分析了。但是当梁的形状足够细长时，截面的翘曲仅是局部效应，对梁的整体 变形影响极小。这是简化结构模型仍然采取了平截面假设的缺点同时也是优点所在。 20 在初始位形中垂直于中轴线的任意一条线段在变形过程中只产生刚体平动与刚体转动， 线段长度不变。 以上两条假设中，第一条是主要的。没有这一条假设就不存在梁这种力学上的一维 简化结构模型。第二条仅是为了分析问题方便而引进的。因为在第一条假设的前提下，该 线段还可以有伸长缩短的变形。引起该线段的伸长缩短变形的因素可以有2个。其一：横 向挤压。这个因素由于不符合梁的变形形式的定义，所以不予考虑；其二：由于轴线的拉 压而引起横向的伸缩。这时，其横向伸缩量可由波桑比确定。所以当拉压与弯曲同时存在 时，可以分别考虑，分别处理。在单独考虑弯曲变形时，就完全可以将其当作一条刚性线 段。这就是二维梁模型在弯曲有限变形时的刚性线段假设。 在梁的线性理论中，通常采用的理论有： Bernoulli 理论：忽略剪切因素，即假设剪切角 γ = 0。此时就有 ϕ = θ ，亦即梁的横截面在 变形过程中永远垂直于梁的中轴线。由于对于变形形式作了更多的限制，在有限变形理论 中并不合适。但在线性理论中由于转角很小还可进一步假设有：dv/dx = tanθ ≈ θ ，可以减 少独立变量，这在有限变形理论中更是不可取了。 Timoshenko理论：考虑剪切因素。 以下我们将把梁的线性Timoshenko理论推广到有限变形理论中。对于梁而言，有限 变形的变形特征是：不仅有大位移，而且有大转动和大曲率。 §2 二维梁的有限变形运动学分析 B1 A0 B0 A1 图3 线段AB自时刻t0至t1的运动 uA uB θ 4 5 如图3所示为一条在初始位形中垂直于梁的中轴线的线段A0BB0，在变形过程中经过一 个刚体平动与一个刚体转动后运动到当前位形中位于A1B1B 位置。刚体平动的位移可由位于 中轴线上的A0点的位移uA表示。刚体转动可由转角θ表示。我们所要研究的梁中不在中轴 线上的点BB0的位移为uBB。设A0，BB0，A1，B1B 四点的矢径分别为： rA0 = (x, 0) (9.2a) rB0 = (x, y) (9.2b) rA1 = rA0 + uA (9.2c) rB1 = rB0 + uB = rB A1 + A1B1B (9.2d) 其中 uA = (u, v) 是位于中轴线上的点的位移。uB = (U, V), 是 B0B B 点的位移。由基本假设线段 A0BB0 的长度与A1B1B 的相等，则由图3可以看出 A1BB1 = y(-sinθ , cosθ)T 将 uA, uB 与 AB 1B1B 的表达式代入(2d) 中rB1的表达式，可以得到 U = u - ysinθ (9.3a) V = v - y(1-cosθ) (9.3b) 这个式子说明垂直于梁中轴线的横截面上任意一点的位移可用中轴线的位移（刚体平动） 及横截面的一个转动（刚体转动）来表示。也就是说梁的变形后的位形完全可由中轴线的 位移（线位移+转动）确定。这就是引进基本假设将一个连续体力学问题降阶为一个一维 的力学问题的最大好处所在。同时我们还应注意到一个非常重要的事实，那就是，中轴线 的位移(u, v)与转角θ 独立无关性。那就是说，位移与转动的次序无关，梁的最终位形都由 （3）式决定。所以尽管中轴线在变形过程中可以发生伸长缩短的变形，而使梁的横截面 也发生伸长缩短的变形，但在转动时还是可以把横截面的运动看作刚体转动。这样就使我 们可以在同时有轴向变形与弯曲变形的复合变形时，将两种变形形式分开处理，保证了第 二条基本假设的成立。下面将位移(U,V)分别对x与y取导。按通常规定可用右下角的一撇代 表取导，例如：∂u/∂x = u,x 。可以得到： U,x = u,x - yθ,xcosθ (9.4a) 6 U,y = -sinθ (9.4b) V,x = v,x - yθ,x sinθ (9.4c) V,y = -1 + cosθ (9.4d) 于是可以由此定义 梁上任意一点的Green-Lagrange 应变： Exx = U,x + (U,x2 + V,x2)/2 = u,x + (u,x2 + v,x2)/2 + y 2θ,x2/2 - yθ,x λs (9.5a) Exy = V,x + U,y + U,xU,y + V,xV,y = λn (9.5b) Eyy = V,y + (U,y2 + V,y2)/2 = -1+ cosθ + ( sin2θ+1−2 cosθ+ cos2θ)/2 = 0 (9.5c) 其中 λs = (1+u,x)cosθ + v,xsinθ (9.6a) λn = v,xcosθ - (1+u,x)sinθ (9.6b) 注意，代表剪切的Exy是应变张量的分量E12 的2倍。横向伸长应变Eyy = 0，是因为刚体转动 的必然结果。但由于Syy = 0 ，Exx ≠ 0 ，必然引起横向伸长应变Eyy ≠ 0，此时可由波桑比计 算Eyy。这种处理方法，与平面应力问题处理沿厚度方向的应变类似。显然(5a)中的前3项表 示了中轴线的Lagrange应变： Exx (y=0) = u,x + (u,x2 + v,x2)/2 (5a)中的后2项则是由于截面转动产生的，而且不对称于中轴线。如果略去曲率θ,x的二次 项并令λs=1，则就得到了线性理论中由于截面转动产生的项。对以上表达式还可以作进一 步的几何解释。定义线段的伸长比： λ = tds / 0ds, 其中0ds 是线段在初始时刻 t0时的原长, tds 是线段在当前时刻t 时变形后的长度。于是 λ 表示了线段的相对伸长。设令线段的原长为 0ds=dx，即假设该线段位于初始位形中梁的中轴线上。则有 tds2 =(dx+du)2 + dv2 以及有 λ = ( (1+u,x)2 + v,x2)1/2 (9.7a) 将(7a)与中轴线的Lagrange应变 Exx (y=0)比较，可以发现，有：2Exx (y=0) = λ2 – 1。事实上 Lagrange应变的定义也正是用类似的方法导出的。我们也可求出在当前时刻 t 时位于A1 点 处的切向量 λ = (1+u,x , v,x)T = (λ.x, λ.y)T (9.7b) 显然切向量λ 的长度即为λ.。切向量λ 完全描述了中轴线的变形运动。令 ϕ 表示切向量λ 与 x 轴的夹角，也就是法线与y 轴的夹角，θ 表示变形后的线段A1BB1 与y轴的夹角，即截面转 角，如图2所示。则有 tanϕ = v,x / (1+u,x) (9.8a) 并有 sinϕ = v,x /λ (9.8b) cosϕ = (1+u,x) /λ (9.8c) 于是方程(6a)与(6b)可重写为 λ s = λcosϕcosθ +λsinϕsinθ = λ cos(ϕ−θ) = λ cosγ (9.9a) λ n = λsinϕcosθ - λcosϕsinθ = λ sin(ϕ−θ) = λ sinγ (9.9b) 将(8)(9)可以合写成为矩阵形式： ⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡ ++− ++= ⎥⎦ ⎤⎢⎣ ⎡ ⎥⎦ ⎤⎢⎣ ⎡ −=⎥⎦ ⎤⎢⎣ ⎡ +⎥⎦ ⎤⎢⎣ ⎡ −=⎥⎦ ⎤⎢⎣ ⎡ γλ γλ θθ θθ λ λ θθ θθ θθ θθ λ λ sin cos cos,sin),1( sin,cos),1( cossin sincos , ,1 cossin sincos xx xx y x x x n S vu vu v u (9.9c) 由(9a)(9b) 可以明显的看出λ s与λ n的几何意义。他们分别是切向量λ在变形后的截面 和垂直于截面方向的投影。其中 γ = ϕ −θ 表示变形后的截面与法线的偏差角，也就是剪切 角。在Bernoulli 梁的理论中忽略了剪切的影响。在线性理论中可以简单地表示为 γ = ϕ − θ = 0。而在有限变形理论中则应更精确的表示为 λ n = v,xcosθ - (1+u,x)sinθ = 0 (9.9d) 方程 (9d) 可以称之为梁的有限变形理论中的无剪切的 Bernoulli 条件。线性理论再进一步假 设转角θ 很小，以及 v,x = θ 。这是以下一系列线性化假设的结果： xxx vuv ≈+=≈≈ )1/(tantan ϕθθ 由于有限变形必须考虑大转动问题，以上的近似简化显然是是不可取的。而不利用挠度导 数与转角的关系，也就达不到减少变量的目的，所以采用Bernoulli理论就没有任何好处， 这就是我们在有限变形理论中只采用Timoshenko理论的原因。在Timoshenko理论中的无剪 切的 Bernoulli 条件则相当于此时产生了纯弯曲。也就是说中轴线除了弯曲外没有其他变 7 形。因此应该有中轴线的Lagrange应变Exx (y=0) = 0。也就是中轴线的相对伸长 λ = 1。 于 是由(9)式有 λ s = 1；λ n = 0。所以得到纯弯曲的Lagrange应变为： Exx = y 2θ,x2/2 - yθ,x (9.10) Exy = 0。 还应注意，除了在建立梁的模型时作出简化假设外，在梁的有限变形理论中是不允许额外 再做任何的简化假设的，不允许忽略任何微小的量，包括：假设位移、转角、曲率、应变 很小，忽略高阶小量，以及限制加载步长的大小在内。引进任何假设的理论都将是某种近 似的几何非线性理论，而不是有限变形理论。所以梁的有限变形理论是在模型意义下的精 确几何非线性理论。 §3 虚功方程与有限元列式 一个弹性物体在受力变形后达到平衡状态时有如下的虚功方程成立： δW(i) = δW(e) (9.11a) 其中δW(i) 与 δW(e) 分别表示物体在平衡状态时内力虚功与外力虚功。在应用有限元法时， 可将外力处理为等效节点外力向量P 。令δq 表示节点虚位移向量。于是外力虚功可表示为 δW(e) = δqTP (9.11b) 内力虚功的一般表示应为 δW(i) = dV V ijij∫ δεσ 其中V是变形后的体积，σij是变形后的当前位形中Cauchu应力， δεij = (δui,j+δuj,i)/2是与虚位 移对应的虚应变。取导是对于当前坐标取的。由于Green-Lagrange 应变Eij 与 2nd Piola-Kirchhoff 应力Sij 在初始位形下能量共轭，所以二维梁的内力虚功可以用对初始位形 的体积分表示为： δW(i) = (9.11c) ∫∫ += V xyxyxxxxijV ij VdSESEVdES 00 00 )( δδδ 8 9 不失一般性，我们假定 2nd Kirchhoff 应力 Sxx, Sxy 与Lagrange 应变Exx, Exy 建满足本构关系 （在弹性范围内） Sxx = EExx , 与 Sxy = G Exy (9.11d) 由于本章主要研究梁变形的几何关系与运动学关系。采用(11d)形式的本构关系仅是为了叙 述的方便。当采用其他的本构关系不会对所介绍的有限变形理论有任何影响。 在已介绍过的有限变形算法中，我们分别介绍了TL, UL, TE, TE等各种算法。其区别在于 在不同时刻的参考位形不同，而无实质的不同。为了简单起见，以后我们将只采用TL 算 法，也就是在任何阶段都以初始位形为参考位形。为了导出Lagrange虚应变，先写出(6)式 的变分 δλs = cosθ δu,x + sinθ δv,x + λnδθ (9.12a) δλn = -sinθ δu,x + cosθ δv,x - λsδθ (9.12b) 于是可得Lagrange虚应变为 δExx = Lxxδe (9.13a) δExy = Lxyδe (9.13b) 其中 e = (u,x, v,x, θ,x, θ)T (9.13c) Lxx = (1+u,x - yθ,xcosθ， v,x -yθ,xsinθ， y2θ,x - yλs， -yθ,xλn) (9.13d) Lxy = (-sinθ， cosθ， 0， -λs) (9.13e) 定义量中轴线上的 u,v 与 θ 的插值函数为 u = ∑Nui ui (9.14a) v = ∑Nvivi (9.14b) θ = ∑Nθιθιi (9.14c) u,v 与 θ 的插值函数可以相同，也可以不同。令q = (q 1T, q 2T,…, q nT )T；qi = (ui，vi，θi　)Τ是 轴线上的第i个节点的节点位移。如果取相同的插值函数，则有： δe = Hδq (9.15a) H是由插值函数及其导数构成的矩阵，其形式为为： H = [H1 , H2,..., Hi,..., Hn] (9.15b) 其中 Hi = ( i =1,2, ...,n) (9.15c) ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ i xi xi xi N N N N 00 00 00 00 , , , 其中 n 是两单元的节点数。方程 (13) 可表示为： δExx = LxxHδq = BBxx δq (9.16a) δExy = LxyHδq = BBxyδq (9.16b) 将方程 (16) 代入 (11) 并由节点虚位移在满足约束条件下的任意性，于是可推得对于初始位 形列写的如下形式的平衡方程： F(q) = (9.17) 0)( 0 0 =−+∫ PBB VdSSV xyTxyxxTxx 对于梁单元可以定义如下形式的Kirchhoff内力： N =∫A Sxx dA (9.18a) Q = ∫A Sxy dA (9.18b) M =∫A y Sxx dA (9.18c) T = ∫A y2 Sxx dA (9.18d) 其中，面积A都是指的初始位形时的截面积。将方程 (5) 与(11) 代入 (18) 式就得到内力的表 达式： N = EAu,x + EA(u,x2 + v,x2)/2 + EIθ,x2/2 (9.19a) Q = GAs λ sinγ = GAs λn (9.19b) M = -EIθ,x λcosγ = -EIθ,x λS (9.19c) T = EI[u,x + EA(u,x2 + v,x2)/2] + EI4θ,x2/2 (9.19d) 10 As 一般称之为有效剪切面积，是考虑到剪应力的不均匀分布而引进的。其计算方法可以参 考线性理论中的计算方法。I 是截面惯矩，其定义与线性理论相同。I4 = ∫A y4 dA。是截面 的4次矩。 (19a)式是Kirchhoff轴力，其第一项与线性理论相同。第一项与第二项之和是由中轴线的伸 长有限变形引起的。第三项与曲率有关，则是由弯曲有限变形引起的。 (19b)是Kirchhoff剪力，其值与轴线的伸长比和剪切角的正弦值有关。在线性理论中不考虑 伸长比，或是说令 λ ≈ 1，并将剪切角看作是微小值，取γ ≈ sinγ ≈ 剪应变，则Gγ就是剪应 力，再乘以有效剪切面积就是线性理论中的剪力了。可见将有限变形理论得出的结果，经 过一系列简化可以推出线性理论的值。 (19c)表示的是Kirchhoff弯矩。将其与线性理论的弯矩表达式，可以看出多了λcosγ 项。这 一项也与轴线的伸长比和剪切角的余弦值有关。如果令λ ≈ 1；cosγ ≈ 1，则也可以退化成 线性理论的结果。 (19d)是只有在有限变形理论中才出现的二次矩。 再一次强调有限变形理论的是绝对不可以对任何变量的值的大小作出限制，不可以采用任 何忽略高阶小量的方法来进行理论推导。正因为如此，有限变形得出的理论结果的适用性 基本上也不受任何限制。唯一的限制就是模型本身。所以说梁的有限变形理论是模型意义 下的精确理论。 引入Kirchhoff内力后，梁内力虚功(11c)是可以在对截面进行精确积分，写成为： δW(i) = dxTMQEN xL xSNxx )]2/,(),()0([ 2θδθλδδλδ∫ +++ 其中Exx(0)是中轴线的应变。 §4 牛顿迭代法求解有限元平衡方程 对于用有限元方法建立起来的平衡方程(17)一般均用牛顿迭代法求解。为了简单和统一 起见，如果不作特殊申明，我们一律采用TL法求解。积分区域均对初始位形而言，不再标 11 注左角的 “0”。下面介绍牛顿迭代法求解过程的基本思想与具体细节。设q是一个近似的试 探解，将方程(17)在q附近用Taylor 级数展开如下： F(q+Δq) = F(q) + ΔF +O(||Δq2||) = 0 其中O(Δq2) 是与||Δq2||同阶的高阶小量。设将其忽略，而用迭代法逐步逼近真解，则有： F(q) + ΔF ≈ 0。 (9.20a) 增量函数ΔF 可写成 ΔF = KTΔq (9.20b) 于是方程(17)成为： KTΔq = -F(q) = P - (9.20c) dVSS xyTxyxxV Txx )( BB +∫ KT 称之为切线刚度矩阵。由方程(17) ΔF 可写成： ΔF = (9.20d) dVSSSS xyTxyxyTxyxxTxxxxTxxV )ΔΔΔ(Δ BBBB +++∫ 注意到算子Δ 仅对节点位移有影响，所以由方程可得 ΔBBxxT = H GT xx HΔq (9.21a) ΔBBxyT = HGxy HΔq (9.21b) 其中 Gxx = (9.22a) ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ −− −− Sxynyxyxy yyy λθλθθθθ θθ ,cos,sin, 2sincos 10 1 对称 Gxy = (9.22b) ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ −−− nλθθ 0sincos 000 00 0 对称 在积分过程中将会出现以下4种形式的积分： ∫ ∫ ∫ ∫A A A A dAydAydAyydA 432 ,,, 12 当中轴线经过截面形心时第一个积分将为零。第二个积分为惯性矩。如果梁截面具有对称 平面而且是中性层时，第三个积分也为零。第四个积分称为4次惯性矩。利用(11d)的本构 关系，有 ΔSxx = EΔExx = EBBxxΔq ΔSxy = GΔExy = GBBxyΔq 将其代入ΔF中，且现刚度矩阵可写成 KT = Kxx + Kxy + Kσ + Κτ (9.23a) 其中 Kxx = E (9.23b) dVxxV Txx BB∫ Kxy = G (9.23c) ∫V xyTxy dVBB Kσ = (9.23d) ∫V xxxxT dVS HGH Kτ = (9.23e) dVS xyxyV T HGH∫ (23d) 与 (23e) 通常称之为几何刚度矩阵或称初应力矩阵。在迭代过程中，每当求得一个增 量位移Δq(i) ，就可以作位移更新： q(i+1) = q(i) + Δq(i)。 以上就是标准的牛顿迭代法。对于牛顿迭代法的其他变形如保持切线刚度不变的修正牛顿 迭代法，加一个增量位移因子的位移更新法： q(i+1) = q(i) + c Δq(i)等就不作介绍了。 在阐述用牛顿迭代法求解非线性方程时我们利用了Tayloy级数展开，将非线性方程 (17)略去增量位移的高阶项，近似地写成为切线形式的(20a)式迭代求解。由于这个原因， 增量位移就不可能太大，否则(20a)式就不能成立。也就是说载荷增量步也要受到一定限 制。这是牛顿迭代法本身带来的缺点，与有限变形理论无关。也就是说有限变形理论本身 没有对增量步长有所限制，它的结果的精确性可以不受载荷增量步长的限制，但是由于受 到求解方法的限制，使得在求解过程中不能采用过大的步长。这就可以得出一个结论，如 果在牛顿迭代法允许的载荷增量步长的范围内，一次加载的结果与多次加载的结果应该完 全一样。也就是说在弹性范畴内，如果加载过程不越过分叉点前，计算结果与路径无关。 理论上可证明标准的牛顿迭代法在牛顿吸引区内具有渐近二次收敛率。牛顿吸引区的存在 13 也就是前面讨论过的由于利用Taylor级数展开，抛去位移增量的高阶小量的采用一阶近似 方程(20a)的结果。 以上结果不论对于初始位形时是直梁还是曲梁一律适用。唯一差别在于对于直梁单 元可以对一个单元只选取一个单元局部坐标系，而对于曲梁单元理论上应该对每一个所要 研究的点都建立一个局部坐标系。这个坐标系以量的中轴线的切线为x轴，以垂直于该切 线的直线即在横截面上的直线为y轴。在实际计算时只要在Gauss积分点上建立局部坐标系 即可。 §5 截面上的Cauchy应力 Kirchhoff应力是Lagrange描述下的应力，是一种计算应力，这种应力与Lagrange应变 在初始位形下满足能量共轭条件。由此可以方便地相对于初始位形建立平衡方程。但在 迭代求解过程后得到位移解的同时还求出真应力, 即在当前位形下的Cauchy应力。因为计 算应力可以有多种多样，参考位形更可以任意选择。对于不同作者得到的不同计算结 果，不仅应该以同一标准来进行比较，也因为强度准则与塑性准则也必须以Cauchy应力 为标准。 在二维情况下Cauchy应力张量与Kirchhoff应力张量之间的转换关系是： T x x xxxx zzyzxz yzyyxy xzxyxx SSS SSS SSS J FF ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ zzyzz yzyyy zy1 σσσ σσσ σσσ (9.24) 其中，Kirchhoff应力矩阵的第三行和第三列均为零。所以得到的Cauchy应力张量的第三行 和第三列也均为零。如此得到的Cauchy应力张量是在初始位形坐标系下表达的，还应该 将其转换到变形后的坐标系中。 梁上任意一点的变形梯度和Jacobi行列式为 14 (9.25a) ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ + + + = z yx yx W VV UU ,100 0,1, 0,,1 F ),1](,,),1)(,1[( zyxyx WUVVUJ +−++== F (9.25b) 初始位形上的横截面在梁变形过程中转过了θ角，沿变形后截面法向n及沿平行截面方向s (n-s右手标架)取微元体。正因为Cauchy应力张量只有4个不为零的量，其Cauchy应力张 量在坐标变换的关系可以简化为 Txxxx xxxx T yyxx xxxx ssns nsnn SS SS J )( 000 0 0 1 000 0 0 000 0 0 QFQFQQ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ σσ σσ σσ σσ (9.26a) 其中Q为坐标转动矩阵 ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ −= 100 0cossin 0sincos θθ θθ Q (9.26b) 令 ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ + = zW aa aa ,100 0 0 43 21 QF (9.27a) 其中， θθ sin,cos),1(1 xx VUa ++= θθ sin),1(cos,2 yy VUa ++= (9.27b) θθ cos,sin),1(3 xx VUa ++−= θθ cos),1(sin,4 yy VUa ++−= 将(27)代入(26)，并可缩减为： 15 ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ +++++ +++++= ⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡ ⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡ yyxyxxyyxyxx yyxyxxyyxx yyxy xyxx ssns sn SaSaaSaSaaSaaaaSaa SaaSaaaaSaaSaaaSa J aa aa SS SS aa aa J 2 443 2 342413231 42413231 2 221 2 1 42 31 43 21nn 2)( )(21 1 σσ σσ (9.28a) 写成向量形式 ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ + = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ xy yy xx n ss nn S S S aaaaaaaa aaaa aaaa J 41324231 43 2 4 2 3 21 2 2 2 1 2 2 1 σ σ σ (9.28b) 对二维梁来说，在(27b)式中 ⎥⎦ ⎤⎢⎣ ⎡−⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡ θθ θθ sin, cos, , , , , x x x x x x y y v u V U ⎥⎦ ⎤⎢⎣ ⎡ +− −=⎥⎦ ⎤⎢⎣ ⎡ θ θ cos1 sin , , y y V U 将其代入(27b)式中，得 xsxxx yyvua ,,sin,cos),1(1 θλθθθ −=−++= a 2 0= nxx vua λθθ =++−= cos,sin),1(3 (9.29) a4 1= J = a1(1+W,z) 将(29)代入式(28b)并考虑到Syy=0，则有 ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ + ++=⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ +=⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ xyxx xyxx xx z xy xx z ns ss nn SSa SSa a a Sa W S S aaa aa a Wa 3 3 1 3 1 131 3 2 3 2 1 1 )2( ),1( 10 0 21 00 ),1( 1 σ σ σ (9.30a) 式(30a)还可以写成 16 ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ + += ++−= −+= ),1( )1( )1( ),1)(,( ),( ),1( 2 z xxn ns xx zxs n ss xxxs z nn W EE E Wy E Ey W E λσ θλ λσ θλσ (9.30b) 现在来看 的意义和求法。显然),1( zW+ ),1( zW+ 的几何意义是截面边宽的伸长比。由于Szz = Ezz − ν(Exx + Eyy) = 0。而且有 Eyy = 0；Ezz = W,z + (W,z)2/2。解二次方程： (W,z)2/2 + W,z − νExx = 0 可以得到 1 + W,z = xxνE2 (9.31) 如果对σnn作当前面积的积分，并利用dA = (1+V,y)(1+W,z)d0A,以及1+V,y) = cosθ，可得 Cauchy轴力为 = d ∫∫ ++== A zynnA nnC AdWVdAN 0 0),1)(,1(σσ θθλ cos),(0 xSA xx yS −∫ 0A = NKλScosθ + MKθ,xcosθ (9.32) 类似可定义Cauchy弯矩MC与Cauchy切力QC，并可求得为 MC = − = M∫A nndAyσ KλScosθ +TKθ,xcosθ (9.33) QC = ∫A ns dAσ = λnNK cosθ + QK cosθ (9.34) 这三个式子体现了Cauchy内力与Kirchhoff内力之间的关系。 §6 数值算例 在下面的几个算例中，均设定收敛条件为ru<ε，rf<ε，允许误差ε =1.e-10，其中 r uu = Δ / u 是位移增量与总位移的欧几里得模之比， P/Rrf ΔΔ= 是不平衡力与加载步 载荷增量的欧几里得模之比。算法均用Newton切线法求解，采用Gauss数值积分，积分点 数=单元节点数–1。 17 算例1. 悬臂直梁加一始端弯矩M将其弯成一圆。 悬臂直梁自由端加一弯矩使其弯成半圆或圆，这是人们在梁几何非线性研究中经常引 用的例子。由Bernoulli梁理论可知，当M=2πEI/L时，直梁将弯成一圆。 18 E=720 G=360 A=1 I=5/30 L=120 M=2π 图4 例1之悬臂直梁弯曲成圆 M 分别用2、3、4和5节点单元进行了试验，计算表明， 2、3、4和5节点单元的收敛性 差别不大, 但精度差别较大。4节点单元和5节点单元的精度比较高且基本一样，3节点单元 的精度稍差，2节点单元的精度最差。我们可以大概分析如下。有限元方法总是把作用于 梁单元上的各种载荷统一地处理为节点载荷。因此在线性理论中梁的变形最高也就是一条 3次曲线。因此4节点的梁单元可以得出梁模型的精确解。而有限变形理论是线性理论的自 然推广。所以推荐采用4节点单元或5节点单元，从计算效率看，最好采用4节点单元。 当采用两个4节点单元(总节点数为7) 或采用两个5节点单元(总节点数为9)时，只用4步 加载可以弯曲成圆。用2或3节点单元则须用较多单元才能准确的弯曲成圆。。图41是悬臂 直梁4步弯曲成圆的示意图。 算例2. 悬臂梁受两个横向载荷。 该算例Crivelli曾细加研究[54]。文献[54]用10步加载求得结果。用有限变形理论通过一 步加载即可。有限元理论已经证明了只要不断细化单元，其数值结果将收敛于精确解。这 个结论对于线性理论成立，对于非线性理论也同样成立，只要这是一种精确的非线性理论 就行。有限变形理论就是这样一种唯一的精确的非线性理论。所以在计算中还采用不断细 化的方法求得4位有效数字的精确解。计算时分别选用2、4、6和8个4节点单元，发现选用 多于6个单元后已不再能提高精度，所以可以认为6个单元的4位有效数字结果是精确的。 表1列出了一步加载梁自由端位移的计算结果。从表1可以看出，在单元数很少的情况下， 精确理论也能达到很高的精度。 表1 算例2有限变形理论一步加载梁端位移的计算结果 2个单元 4个单元 6个单元 8个单元 u 30.85 30.80 30.79 30.79 v 67.11 67.03 67.03 67.03 θ 1.043 1.043 1.043 1.043 E=3.0e+7 G=1.1538e+7 A=0.2 I=0.1667 P1=850 P2=1350 L1=52.03 L2=102.75 图5 算例2之悬臂梁受两个横向载荷 P2 P1 19 。 算例3. 一端固支—端铰接拱受垂直载荷作用。 L1 L2 20 该算例R.D.Wood与O.C.Zienkiewicz等曾加以研究[120,121,89]。文献[121]应用平面应力 (paralinear)单元计算该题，求得拱发生跃变的临界载荷为9.24P(该文引用的解析解是 8.97P)，并给出了拱顶端的载荷位移曲线。本文将一端固支—端铰接拱剖分成16个4节点梁 单元进行计算，当载荷增大到一定数值时，拱的变形将发生跃变，下面给出用有限变形理 论计算得到的发生跃变时的临界载荷与拱顶端的位移计算结果。 应用平面应力 (paralinear)单元计算该题，求得拱发生跃变的临界载荷为9.24P(该文引用的解析解是 8.97P)，并给出了拱顶端的载荷位移曲线。本文将一端固支—端铰接拱剖分成16个4节点梁 单元进行计算，当载荷增大到一定数值时，拱的变形将发生跃变，下面给出用有限变形理 论计算得到的发生跃变时的临界载荷与拱顶端的位移计算结果。 Pcr = 8.94P；u = -60.99；v = -113.2；θ = -.05294 P 算例4. 该算例由作者首次提出。一均质直梁，密度为ρ，初始位置垂直向下，在自重下有 初始伸长变形， 然后每次转动 Δα角度，转动数次后回至原位。若刚度一定，密度ρ 较小 时，垂直向上位置是一个稳定平衡位置；密度ρ 太大时，垂直向上位置将是一个不稳定平 衡位置。图7是用有限变形理论计算该例的变形图。作者，采用三个4节点单元，已可求得 临界密度的数值精确值为 算例4. 该算例由作者首次提出。一均质直梁，密度为ρ，初始位置垂直向下，在自重下有 初始伸长变形， 然后每次转动 Δα角度，转动数次后回至原位。若刚度一定，密度ρ 较小 时，垂直向上位置是一个稳定平衡位置；密度ρ 太大时，垂直向上位置将是一个不稳定平 衡位置。图7是用有限变形理论计算该例的变形图。作者，采用三个4节点单元，已可求得 临界密度的数值精确值为 cr = 8.94P；u = -60.99；v = -113.2；θ = -.05294 3804.7 gAL EI cr =ρ (当Ε=1.5e+6，G=.75e+6时，ρcr =1.631)。 Timoshenko用弹性小变形理论求出的下端固定上端自由等截面压杆在均布载荷作用下的临 界载荷[122]为 3cr 8377 gAL EI.=ρ ，比本文结果略大。这是可以理解的，因为Timoshenko的解 来自线性理论，本文结果来自有限变形精确理论，从最小势能原理可知，精确理论的载荷 值最小，位移值最大。 R β β P E=2e+6 G=1e+6 A=1.0 I=0.5 R=100 β=107.5° P=100 图6 例3之初始结构图及其变形图 L=130 A=36 I=108 Δα=π/4 g=9.8 (a). Ε一定，ρ较小时 (b). Ε一定，ρ较大时 图7 例4之变形图 α = π 7，随动载荷的载荷刚度矩阵及计算例题 图6 梁端作用的随动力 θ0 QN QS x 载荷的作用方式可以分做两种，一种是作用方向保持空间不变的载荷，即死载(Dead loading)；一种是作用方向依赖于结构位移的载荷，即随动载荷(Follower loading)。前者对 应于保守系统，后者对应于非保守系统。尽管对于后者，结构的稳定性分析时常需求助于 21 动力学方法，但对于屈曲前的结构的非线性行为，仍可采用静力学方法。我们在前面建立 的有限元平衡方程(17) F(q) = 0)( 0 0 =−+∫ PBB VdSSV xyTxyxxTxx 并用牛顿迭代法由(20a)导出 (20b) ΔF = KTΔq ，此时隐含着载荷P的位置以初始位形为参 考位形时是不变的即假定为死载。当载荷P的位置随时变化时，这种载荷称为随动载荷。 当具有随动载荷时(20b)应推广为： ΔF = KTΔq + KPΔq。 KP可称之为随动载荷切线刚度矩阵，简称为载荷刚度矩阵，它的一般表达式为： KP = 。其具体形式随问题而异。下面举例说明。 qP ∂∂− / 如图6所示，梁之一端作用有互垂的一对力，QS与QN。在梁的变形过程中QN永远沿梁截面 方向，QS在初始位形中沿梁中轴线切线方向，而在变形后永远与梁截面垂直。在一般情 况下，初始位形中的该处的中轴线与总体坐标系的x轴有夹角θ0。于是QS与QN在总体坐标 系中的分力Qx与Qy可表示为： ⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡ − y x N Q Q Q QS 00 00 cossin sincos θθ θθ 设在变形后梁端截面转动了θ角，此时的QS与QN也随之而转动θ角。具有这样性质的力称 之为随动力。其时随动力QS与QN在总体坐标系中的分力Qx与Qy可表示为： ⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡ ++ +−+=⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡ −⎥⎦ ⎤⎢⎣ ⎡ − y x NN Q Q Q Q Q Q S 00 00S 00 00 )cos()sin( )sin()cos( cossin sincos cossin sincos θθθθ θθθθ θθ θθ θθ θθ 随动力形成的随动载荷刚度矩阵 θ∂−∂= /PK P 是应该在局部坐标系中表示的，所以初始 转角不应被考虑进去。就有： ，这个矩阵仅出现在关于 该节点的对角线子块处。显然，当有随动载荷作用时，切线刚度矩阵是不对称的。 ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − −− −= 000 sincos00 cossin00 NS NS P QQ QQ θ θθ K 下面是用本文的有限变形理论处理随动载荷的几个算例，计算使用不对称求解器完成 的。 22 算例5. 悬臂梁在自由端受一垂直随动载荷作用。该算例由J.H.Argyris提出[16]，J.C.Simo也曾 加以研究[89]。文献[16]给出了134步加载(步长P=1000)的结构变形图和载荷—位移曲线。文 献[89]采用Reissner-Simo理论，用5个3节点单元计算，只给出了130步加载(步长P=1000)的 结构变形图和载荷—位移曲线。 本文采用精确理论，也用5个3节点单元进行计算，图5.5是计算所得的197步加载(步长 P=1000)的变形图