第7章能量泛函的转换形式及其应用(16K)

1、第 7 章 能量泛函的转换形式及其应用§ 7.1 总位能泛函转换形式及其应用由§ 5.1 节中的（ 5-16）式，定义了总位能泛函，即P A(ij ) Fi ui dVTi ui dS（ 5-16）VS该泛函为单变量变分原理，其自变量要求满足位移应变关系及位移边界条件，即ij1 (ui , ju j ,i ) ，在 V 内2uiui 0 ，在 S 上所以，这种变分原理是有条件的，并可以进一步证明总位能原理是极小值原理，解的收敛性得到保证。 这种原理是目前广为流行的绝大部分有限元素模型的基础，比较理想的情形是 “保续元”的建立，而放松某些边界协调条件则构成了有限元素法中的“非

2、保续元”。【例 7-1】 弯曲梁的总位能泛函及其变换。图 7-1 所示的一维梁，承受横向分布载荷p(x) ，简支端（ x L ）作用一集中力矩 M ，梁的另一端为固持。显然，其边界条件为x0： w(0)w (0) 0xL ： w(L)0及M(L) M（ 7-1）总位能泛函根据定义可写为pU V（ 7-2）其中图 7-1一维弯曲梁U1L2 dx （应变能）（ 7-3）EJ (w )20VLM w ( L) （外力位能）（ 7-4）pwdx0上面各式中， w 表示挠度，它是坐标x 的函数，而 w 与 w 分别代表 dw 及 d2 w。dx dx 2现在对总位能取一阶变分，UVLLMw (L) pE

3、Jw w dxpwdx（ 7-5）00当弯曲刚度 EJ 沿长度不变时，可将它放在积分号之前，再利用Green 公式，可得LEJw w 0LEJw w 0LLEJw w dxEJw(4 ) wdx（ 7-6）00将（ 7-6）式代入（ 7-5）式中，利用条件（7-1）式，整理后可得 pL EJw(4 )pwdx EJw x LM w（ 7-7）0现令（ 7-7）式的 p0 ，利用变分法中的预备定理，可得到-100-EJw (4)p0（ 7-8）EJw x LM0（ 7-9）（ 7-8）式即为平衡方程，与材料力学所导出的公式完全一致， （ 7-9）式为力的边界条件，即相当于（ 7-1）式中的最后一

4、个公式。以上的分析再次验证了总位能泛函的驻值条件是等价于平衡方程的。应当指出，方程（7-8）对自变量即挠度w 要求它具有四阶可微，而泛函（7-2）中最高可微阶次为两次。显然，定义泛函p 的自变量的因次可能满足不了平衡方程（7-8）的要求，从这一点来说，直接利用泛函（7-2）来导出的离散型式有限元素法模型，对自变量阶次的要求可能要低得多，这对选择自变量的函数形式带来方便。在连续体力学中所求寻的解一般都具有高阶可微性，且满足微分方程及所有的边界条件。有限元素法情形却不一样，它的解是用有限个自由度来表示的，且是分片光滑函数，这些函数的可微性一般均低于微分方程式中导数的最高阶数。【例 7-2】 图 7

5、-2为一维梁元素，节点位移分别为w1 , 1, w2 , 2，下标 1 代表节点 1 的，下标 2 代表节点2 的，节点位移列阵为 w11w22 T（ 7-10）因为节点位移有四个，我们以3 次多项式表达挠度 w ，即w a1a2 x a3 x2a4 x 3（ 7-11）或w x（ 7-12）式中： x1x x2x3 ，图 7-2一维梁元素 a1a2a3a4 T显然，（ 7-11）式的阶次并不满足平衡方程式（7-8）。利用节点位移（7-10）式，可得 c（ 7-13）则（ 7-12）式化为w x cN（ 7-14）式（ 7-14）中的矩阵 N 为位移插值函数，其物理涵意在一般有限元书中均有说明

6、。下面由式（ 7-14）式导出几何矩阵 B ，梁的弯曲应变为2 wz B（ 7-15）z（ 7-15）式中的 B 阵为x22 B2 N（ 7-16）x将（ 7-14）式中的 N 代入（ 7-16 ）式，可求出几何矩阵B 为 BL 62( 2x1)2 ( 3x2)62( 2x 1)2 ( 3x 1)LLL LLLLL最后，利用（5-16）式求出梁的总位能泛函为 T N T p dx T F e （ 7-17）p1 TLBT D Bdx L200式中 DEJ 为梁的抗弯模量，Jz2 dydz 为梁横截面关于y 轴的惯性矩。A-101-由泛函p 的驻值条件，即 p 0，可得 eK （ 7-18） F

7、式中eLT p dx Fe（ 7-19） F N 0为梁元素的等效节点力。利用能量法求近似解的方法较多，其中Rayleigh-Ritz法是一种有效而应用得比较多的一种方法。其主要是选用一系列满足位移边界条件的函数wi(i 1,2,) 来离散实际位移，如nwai wi（ 7-20）i1ai 为待定参数。将上式代入总位能泛函中，得到以ai为独立变量的泛函如pp (a1， a2， ， an )利用泛函驻值条件，np ai p0（ 7-21）i 1ai得到一组代数方程式，p0(i1,2, n)（ 7-22）ai譬如对于图 7-1 所示的一端固持一端简支的梁，（ 7-1）式表示其边界条件。现取w1x2

8、(x L), w2 x 3 ( x L)（ 7-23）显然，（ 7-23）式是满足位移边界条件的两个连续函数。梁的可能挠度w 可取为w a1 x 2 ( x L) a2 x 2 (x L)（ 7-24）这类函数的形式甚多，这里不在列举。【例 7-3】 薄板的总位能泛函及其变换形式。总位能泛函在薄板中也得到广泛应用。 下面我们讨论略去横向剪切效应的 Kirchhoff 板的总位能泛函的形成过程。图 7-3 为板边界的正向边界力的规定， V z , M n , M ns 表示给定的边界力， p 为分布法向载荷。其应变能为图 7-3弯曲板正向边界力U1( M xk xM y k y M xy kxy

9、 ) dxdy（ 7-25）2A-102-式中， kx 、 k y 、 kxy 为板的曲率，2 wk xx2k y2 w（ 7-26）y 2kxy2 wxy2 wM x10x 2M yD102 w（ 7-27）y 2M xy00(1) 222wxyEt 37-26）式和（ 7-27）式代入（ 7-25）式中，式中 D2)，是材料的泊松系数。将（12(1可求得UD(2 w2 w22(1)2 w 2 w2 w2 dxdy2Ax2y2 )x2y2()（ 7-28）x y给定边界上的外力是由以下几部分组成：表面法向载荷p 、法向给定边界力矩M n 及等效给定剪力 VzM ns 组成，于是外力位能V 等

10、于sVpwdxdyw(V zM ns（ 7-29）C1 M n) dsAns式中 C1 表示力的给定边界，而用C 2 表示位移给定边界。总位能泛函为D(2 w2 w22(1)2 w 2 w2 w2 dxdyp22 )22()2Axyxyx ypwdxdyMw(VzM ns（ 7-30）nns)w dxdyAC1（ 7-30）式给出的薄板总位能泛函的一般形式。对于具体薄板（给定位移边界及力边界条件等各种情形），上式应作相应调整。譬如对如图7-4所示的四边简支矩形板，承受横向分布载荷p ，泛函（ 7-30）式只保留前两项积分，即D(2 w2 w22 w 2 w2 w)2 dxdypwdxdyp22

11、 )2(1 )22(2Axyxyx yA（ 7-31）如果利用Rayleigh-Ritz 法求解，可取三角函数来求挠度w ，如下面的形式-103-w（ 7-32）Amn sin m x sin n ym 1 n 1ab图 7-4四边简支矩形板对（ 7-32）式求导，并利用三角函数积分正交性，再代入（7-31）式后，可得422U ab DAmn2mn8m 1 n 1a2b2ab4Amnmn， (m, n为奇数 )2Vpwdydxab m 1 n 1（ 7-33）000， (m, n为偶数 )最后，利用总位能的驻值条件p0（ 7-34）Amn得到一组代数方程，从而求出系数Amn ，并得到挠度值w

12、。总位能原理泛函为位移协调元素模型的建立做出了贡献，其实质是由变分泛函直接形成离散的有限元素模型，而不是通过变分运算得到微分方程。元素刚度矩阵的形成过程在有限元专著中均可查到，这里只简单的回顾一下。【例 7-4】 图 7-5 所示的矩形弯曲薄板。节点位移为 e T（ 7-35）1234其中 wxiyi T(i 1,2,3,4)， w、xi、yi分别表示节点挠度、转角等。通过iii双线性插值函数，完成位移的离散，如w N e（ 7-36）式中N N1N x1N y1N 2N x 2N y2N 3N 3 xN 3 yN 4N 4 xN 4 yN i , N xi , N yi (i1,2,3,4)

13、 都 是 x, y 的 四 次 多 项 式（各式可查阅有限元素法教材）。广义应变与节点位移关系， 是由几何矩阵 B 体现的，如图 7-5矩形弯曲板元 B e（ 7-37）-104-几何矩阵 B 为 3×12矩阵2 N 12 N x12 N y12 N y 4x2x 2x2x2B2 N 12 N x12 N y12 N y 4y2y 2y2y 22 N122 N x12 N y122 N y42yx y2yx yxx现将（ 7-36）式和（ 7-37）式代入总位能泛函，经过整理后，可得1e)TK e( eTe（ 7-38）p( ) F2式中，刚度矩阵 K BT D BdV（ 7-39）

14、V等效节点力 e N TpdAe（ 7-40） F FA（ 7-40）式 F e 为实际作用到节点上的载荷，它组成等效节点载荷的一部分。由驻值条件 p0 ，得到薄板弯曲时的刚度方程为K e e（ 7-41） F【例 7-5】 现在讨论图7-6 所示薄板的屈曲失稳情形。在失稳之前，我们假定在薄板的中面上承受平面应力0x ,0y 和0xy ，这里表示一比例常数。 这些应力可视为初应力， 它们均满足平衡条件和力的边界条件（忽略体力）：ij , j0（在体积 V 内）Tij n j（在 S 或 C1上） （ 7-42）另一种边界为位移边界C 2 （或称 Su ），对于总位能泛函，则需要预先给定，如u

15、vw0及w0（在 C2上）n图 7-6薄板的屈曲失稳（ 7-43）薄板的中面力可以用单位长度上的力N x0 ，N y0和 N xy0表示，这部分力可视为在失稳过程中是不改变大小与方向的常量，由于中面的平面内的变形，这些力所作用的功为V1A N x0 (w ) 2N y0 (w ) 22N xy0ww dxdy（ 7-44）2xyxy注意积分号内的 N x0等以压力为正，所以在积分号内各力在计算时均取正值。将（ 7-44）式代入总位能泛函，则可写出12 w2 w22 w 2 wD(22(1)2 w dxdyp22 )222Axyx yxy-105-1 N x0 (w ) 2N y0 (w )21

16、N xy0ww dxdy（ 7-45）2Axy2xy式中挠度 w 为独立变量，要求w 必须满足给定的边界条件如（7-43）式之 C 2 边界等。经过对（ 7-45）式泛函自变量的离散化，并按有限元素法刚度方程的形成过程，最后可求得一组特征方程，并由此而求出其特征根，确定了失稳临界系数c r。§7.2 总余能泛函转换形式及其应用由§ 5.1 节中的（ 5-22）式定义了总余能泛函为cB(ij )dVTi ui dS（ 5-22）VSu该泛函为单变量泛函，自变量为力或广义力，泛函成立的约束条件是自变量处处满足平衡方程及力的边界条件，这与总位能原理是相对应的，它同总位能原理类似，

17、也是属于两种不同场量的经典变分原理。总余能原理在有限元素法中的应用，不如总位能原理广泛，而它对构造应力杂交模型作出了贡献，为有限元素法开辟了另一领域。为了与总位能泛函有所区分，这里用 U * 和 V * 分别表示余应变能及外力余功，总余能泛函表示为cU *V *（ 7-46）式中U *1ij T D ij dV2VV *T T u dS（ 7-47a,b）Su如果讨论的对象为平面应力板，则应力分量可以表示为 xyxy T（ 7-48）现在引用应力函数，在不考虑体力的情形下，应力函数与应力分量的关系为222xy 2,yx 2 , xy应力函数应满足下面的平衡方程（ 7-49）x yxxyxyy0

18、（ 7-50）xy0,yx现将（ 7-49）式代入（7-48）式，（ 7-48）式又可以表示为yy xx（ 7-51）xy显然，（ 7-51）式中之 表示以二阶导数微分算子前乘应力函数。以上引入应力函数的目的是为了用节点应力函数（包括导数）来离散元素应力，这恰如基于位移法的有限元法中以节点位移来离散位移有相似之处。【例 7-6】 对于图7-7 所示的边长为 a 和 b 的矩形平面元素，节点编号为ij , (i , j 1,2) ，节点应力函数为 eij ,xij ,yij ,xyij T16 1（ 7-52）由（ 7-52）式可以分别表示元素的4 个节点参数。如果11-21 边为应力给定的边，

19、则应有-106-图 7-7矩形平面板元素y|0y yxy|0xy为了保证元素与元素之间的协调，这里采用了Hermitan 插值函数，对自然坐标写出各插值函数为N1 ()13223N 2 ()3223N x1 ()a(223 )N x2 ()a(32 )及N1 ()13 223N 2 ()3223N y1 ()b(223 )N y2 ()b(32 )应力函数可以由下式插值完成形函数 N 展开后，为以下N e16× 1 矩阵形式N N1N1N 1N 2N 1N y1N1N y2N 2N1N 2N 2N 2N y1N 2N y2N x1N 1N x1N 2N x1N y1N x1N y 2

20、N x2N1N x2N 2N x2N y1N x2N y2对应于形函数 N 的节点应力函数参数 e 为( ,)可分别（ 7-53）（ 7-54）（ 7-55）（ 7-56） e 1112y11y122122y21y 22x11x12xy11xy12x21x22xy 21xy22 （ 7-57）-107-式中x , y , xy 均表示对 x 、 y 及 xy 的导数。利用（ 7-51）式和（ 7-55）式，元素应力可表示为N yy N xx e G e（ 7-58）N xy将（ 7-58）式代入余应变能（7-47a）式中，（ 7-47a）式可以转换为U *1 ( e ) T f e（ 7-59

21、）上式中的 f 即为元素的柔度矩阵，为（216× 16）的对称矩阵，且 f G T D 1 GdV（ 7-60）V为了保证元素与元素间的协调，则需要应力函数及其导数沿边界保持连续条件，n这种连续条件可以由双三阶 Hermitan 插值予以保证。这些元素均属于元素相交边界，对于那些应力指定的边界上，这些应力属于已知量或称为边界力，这部分边界力的平衡关系在应力函数中难以满足，因此，对总位能泛函应作松弛处理，或者说以 Lagrange 乘子项对总余能泛函进行修正。对于图 7-7 所示的 11-12 边界上的指定应力为y 及 xy 的情形，2 ey |y 0x2|y 0 N xx （ 7-6

22、1）2 exy |y 0|y 0 N xy （ 7-62）x y如果已知边界上指定的力为Ty（ 7-63）xy从（ 7-61）式及（ 7-62）式中，边界应力x, y 还需要应力函数的二阶偏导yy 与xx ，所以上述的插值对边界来说是不够的，故必须引入满足二阶导数条件的五次插值函数，如图7-8所示。为了与一般元素的应力函数有别，对应力边界元素现在以b 表示其应力函数。插值函数为图 7-8五次锸值函数-108-Q1 (s)110s315s46s5Q 2 (s)10s315s46s5Q3 ( s)s6s38s43s5Q 4 (s)4s37s43s5Q5 ( s)1 (s23s33s4s5 )2Q

23、6 ( s)1 ( s32s4s5 )2且b Q1 ( s)(0)Q 2 ( s)(1) Q 3 (s) s(0)Q 4 ( s)s (1)Q 5 (s)ss (0)Q 6 (s)ss (1)为了更简洁，这里的一般插值函数表示为（ 7-64）（ 7-65）C1 (s)1 3s22s3C 2 ( s)3s22s3（ 7-66）C 3 ( s)s2s2s3C 4 ( s)ss3则应力函数b 为b C 1 (s) (0) C 2 (s)(1)C 3 (s) s (0) C 4 ( s) s (1)（ 7-67）如果图 7-7 所示的元素为边界上的一个元素，为了满足（ 7-51）式的要求，该元素的应力

24、函数可表示为bC1 ( y) Q1 ( x)11Q 2 (x)21Q3 ( x)x114Q ( x)2C ( x)2C ( x)2C ( x)Nx21Q 5 ( x)22 C 3 (x)y21C 3 ( x)y22C 3 ( x)xx11Q 6 ( x)xx21 C 2 ( y) C 1 (x)12x12C 4 ( x)x 22 C 3 ( y) C 1 ( x)y11xy11C 4 ( x)xy 21 C 4 ( y) C 1 (x)y12xy12C 4 (x)xy22 （ 7-68）这里形函数 N 的顺序可以重新按 的对应顺序排列一下及11212212x11x21x 22x12y11y 2

25、1y 22y12xy11xy21xy22xy12xx11xx21xyxyxx 18 1 N C 1 ( y)C 1 ( x)C 4 ( y)C 4 (x)（ 7-70） N 式见下页。现在利用（ 7-68）式，将（ 7-68 ）式代入（ 7-61）式与（ 7-62）式，并合并此两式，如下-109-N xx C （ 7-71）N xy再利用力的边界条件 Tiij n j ，将（ 7-71）式及（ 7-63 ）式代入后，可得C T （ 7-72）（ 7-72）式即为边界力平衡方程。依上类似，如果边界元素两边均承受已给定的应力，如x ( y) 、y ( x) 及xy ，如图 7-9 所示。边界力为x

26、 0，x ( y) 及 xy ( y)y0 ，y ( x) 及 xy (x)对 11-21 及 11-12 两条边均需引用（ 7-64）式的插值函数，为清晰起见，可用下表表示：图 7-9两边承受给定应力的边界元素-110-N i Q 1 (x)Q1 ( y)C 1 ( x)Q 2 ( y)C 1 (x)Q 3 ( y)C 1 ( x)Q 4 ( y)C 1 ( x)Q 5 ( y)C 1 ( x)Q 6 ( y)11Q 2 ( x)C 1 ( y)21Q 3 ( x)C 1 ( y)x11Q 4 ( x)C 1( y)x21Q 5 ( x)C 1 ( y) xx11 Q 6 ( x)C 1 (

27、 y) xx2112C 2 ( x)C 2 ( y)22C 3 ( x)C 2 ( y)x12C 4 (x)C 2( y)x22y11C 2 ( x)C 3 ( y)y21C 3 ( x)C 3 ( y)xy11C 4 (x)C 3 ( y)xy21y12C 2 (x)C 4 ( y)y22C 3 ( x)C 4 ( y)xy12C 4 ( x)C 4 ( y)xy22yy11yy12（ 7-73）-111-（ 7-73）式中的 Q 1 ( x)Q1 ( y) 为 10 次多项式， 这样的形函数为10 次函数。 而且这种排列是十分规律的，对任意扩大的各种力的边界都十分方便，其中包含有Q 5 (

28、s) 与 Q 6 (s) ，所有的边界应力均可得到满足。边界元素的节点参数由（7-73）式不难写出11212212x11x21x22x12y112 y1y22y12xy11xy 21xy22xy12xx11xx12yy11yy11（ 7-74）下一步是按照有限元素法常规过程进行，将各元素的局部自由度 e ，向结构总体自由度 k 过渡，如 eLk （ 7-75）由于自变量应力函数对边界条件（7-72）是松弛的，故（ 5-22）式的总余能泛函必须将部分松弛条件以 Lagrange 乘子相乘计入总余能泛函，参考（7-47a）式，形成如下形式的松弛泛函，如c1 (k ) T f k ( TC k )（

29、 7-76）20 ，自变量为 k 及 取式（ 7-76）的驻值条件即 c ，于是有以下各式c0 ： f k CT 0（ 7-77）kc0 ：T Ck 0（ 7-78）将（ 7-77）式与（ 7-78）式合并，得 f CTk 0（ 7-79） C0 T 由（ 7-79）式的第一式，可得k f 1C T （ 7-80）再取（ 7-79）式第二式，并将（7-80）式代入，可得 C f 1 C T T 取FC f 1 C T ，则上式可转化为F T 及 F 1T（ 7-81）将（ 7-81）式代入（ 7-77）式中，则有 f k （ 7-82）式中 CTF 1T（ 7-83）（ 7-82）式为最终公式， 不妨称为广义位移， f 为广义柔度矩阵，显然，它具有对称带状等特点，因此，对计算带来很大的方便。§7.3 混合泛函变分原理及其