变分原理的直接方法

为了书写方便,我们在这里先引入内积空间和线性算子两个概念。1.内积空间H是复线性空间，H上定义一个两元函数x,y:HHC，x,yH满足(A)对称性，x,yy,x(B)双线性，ax1bx2,yax1,ybx2,y，a,bC,x1,x2H,(C)正定性，x,x0，而且只有当x0时等号才成立那么我们称该两元函数,定义了线性空间H上的一个内积。定义了内积的线性空间称为内积空间。nnHR,x,yR,x(x1,x2,...,xn),y(y1,y2,...,yn)，那么下面定义的例7.1:就是内积x,yxiyi

i1

n

例7.2:H是定义在[a,b]上连续函数所组成的线性空间C[a,b]，(x),(x)H，那么下面定义的就是内积b

,(x)(x)dx

a

通常，在连续函数空间中按如下来定义内积,w(x)(x)(x)dxab其中w(x)0是个权函数。2.线性算子A是定义在内积空间H上的一个线性映射A:H,AH；如果存在另一线性映射A，使得*

,H,A,,A**

，则A称为A的伴映射。当H是函数空间*

*

时，A称为（线性）算子，A称为共轭算子；特别当AA，A称为对称算子。如果当对称算子A满足A,0

，并且等号仅当0时成立，则A称为对称正定算子。7.1里兹方法（Ritz）

由线性对称正定算子A及函数f所确定的一个线性泛函为(u)Au,u2f,u该泛函的变分为(7.1.1)2Au,u2f,u2Auf,u泛函极值问题(也就是变分问题)其所对应的Euler方程为Auf0(Auf)次的边界条件:(7.1.2)不失一般性,我们假设泛函的边界条件是齐次的，否则我们总是可以通过函数变换来实现齐uu0u

其中u0非齐次的边界条件,那么u满足齐次的边界条件。现选定一组满足泛函齐次边界条件的函数序列u1,u2,...,un

,那么由该函数序列所张成的子空间为nUspan(u1,u2,...,un)u|uaiui,aiRi1该子空间上的每个函数都满足齐次边界条件.里兹法的核心思想就是用上述函数序列所张成的一个线性空间Uspan(u1,u2,...,un)

u1,u2,...,un不是来近似地替代原泛函的定义域空间，然后在线性空间U中找到一个使得泛函最小的一个函数,该函数就是原问题的一个近似解。显然满足齐次边界条件的函数序列唯一的，如果我们选择了比较合适的函数序列的近似解就能很好地逼近原问题的解。具体地讲，由u1,u2,...,un

，而且该序列的个数足够多时(当然函数序列的个数越多,其张成的子空间就越逼近原来的定义域空间)，那么里兹法所得到u1,u2,...,un

n

线性组合成的一个函数为~auuii

i1

其中系数a1,a2,...,an为待定的常数，对应的泛函为)Au,u2f,u(u

Aaiui,aiui2f,aiui

i1

i1

i1

n

n

n

由于A是线性正定对称算子，那么)aiAui,aiui2f,aiui(ui1i1i1nnnaiajAui,uj2aif,uii1j1nni1nnnAijaiaj2fiaii1j1i1n式中AijAui,uj,这是一个关于fif,uia1,a2,...,an的一个二次型。选择的a1,a2,...,an要使得该函数取到最小值，也就是说

0,as

从而有ns1,2,...,n

Aafisii1s0,s1,2,...na1,a2,...,an，由此得到原泛函极值这是关于一个线性代数方程组。解此代数方程组后得到问题(或者微分方程边值问题)的近似解为n

~auuii

i1

如果用向量的形式来表示u1,u2,...,un,aa1,a2,...,an那么TTaTu

)AaTaT2f,aT(a)(u

aTKa2aTF

其中K[Kij],KijAui,ujFf1,f2,...,fn,(a)

0a

这也就是说Tfif,uiK是nn的矩阵，F是n1向量。要使(a)取到最小值，必须KaF该方程的解为aK1F由于A是对称正定算子,可以证明K是对称正定矩阵,上述解必定存在。所以说，通过里兹法，我们可以把一个泛函的极值问题转化成一个函数的极值问题，求解该函数极值问题所对应的代数方程组，就可以得到原问题的近似解。里兹法的关键在于函数序列的选择,如果选择合适的函数序列是该算法最核心之处。例7.3求变分问题1

J[y](x2y2xy)dx,

0

y(0)y(1)0

的近似解。~~

取y(x)a1x(1x),那么y'(x)a1(12x)

J[a1][a12(x24x34x4)a1(x2x3)]dx01dJ0da令1得到

由此可以求得1

0

[2a1(x24x34x4)(x2x3)]dx0

5

a116

一阶近似解为5

(x)16yx(1x)

更进一步,可以取近似解为(x)a1x(1x)iyi1n例7.4设J[y(x)](y'2y22xy)dx,y(0)y(1)0

0

1

求变分问题的近似解。该变分问题的精确解为y(x)n

sinxxsin1~

n1,如果,也就是说y(x)a1x(1x),上面求法一样,得现取~y(x)a1x(1x)k

i1

5

a118到,也就是说一阶近似解为5

(x)18yx(1x)

JJ

0~aa2

如果n2,y(x)x(1x)(a1a2x),代入泛函表达式,并令1,得到3

10310

31

a120a21213a2105a2

1

20

由此可以求得a171369

,

a2741

也就是说两阶近似解为717

(x)x(1x)(369y41x)

与精确解相比,两阶近似解误差已经非常小。例7.5长度为l，抗弯刚度为EI的简支梁，受均布载荷q的作用。图7.1例7.5图取位移（挠度）的试函数为xxxx~(x)awaaa1234llll

为了满足两端位移的简支边界条件，取234

a1(a2a3a4)

那么xxxx(x)a2a3a4a2a3a4wllll梁内的应变能为234d2w1

U1EI2dx

20dx

l

2

作用在梁上的外力势能为U2qwdx

0

l

把挠度的试函数表达式代入总势能表达式中UU1U2

2

216x12x

U1EI2a22a3a42dx

20lllll

234

lxxxx

U2qa2a3a4a2a3a4dx

0

llll

l

2

由UUU

0a2a3a3

得到a20,a3从而梁挠度的近似解为112ql4,a4EI124ql4EIw

1

24ql4x2x3x434EIlll

这里只是满足了位移边界条件,但是没有满足力的边界条件。练习:另一种解法：取位移（挠度）的试函数为~(x)asinxasin2x...w12

ll

7.2康托罗维奇法（Kantorovich）

康托罗维奇法是里兹法在多元自变函数变分问题中的推广。假设泛函的自变函数是个关于x1,x2,...,xn(n2)的多元函数，在康托罗维奇法中，取试函数为~a(x)u(x,x,...x)uini12n1i1k(7.2.1)也就是说,现在的基函数为定的关于ui(x1,x2,...xn1)

它要满足相应的齐次边界条件，而ai(xn)

是待xn的函数。将该试函数代入到原泛函(u)Au,u2f,ua1(xn),a2(xn),...,ak(xn)的新泛函*(a1(xn),a2(xn),...,ak(xn))a1(xn),a2(xn),...,ak(xn)，使得新泛函能取到极小值。这是关于多个得到一个关于于是问题就变为求函数一元函数的变分问题，相应的Euler方程一般为常微分方程组（而原来变分问题得到的Euler方程一般为偏微分方程），求解该常微分方程的边值问题就得到了原变分问题的近似解。和里兹法相比，康托罗维奇法稍显麻烦，因为里兹法最终得到的是代数方程，而康托罗维奇法最终得到的是常微分方程组。但是由于里兹法中的试函数一般都不满Euler方程，而康托罗维奇法中有一部分函数是通过求Euler方程的边值问题得到，所以康托罗维奇法的精度一般要比里兹法来得高。在实际应用中，一般把变化比较复杂而较难选择试函数的那个变量选作为xn。和里兹法一样，要提高康托罗维奇法的精度，一般要增加试函数中的项数，但是相应的Euler方程个数也会有所增加(微分方程组的维数升高)，在实际计算中有很多麻烦。另外一种较为合理的方法是变量轮换法，也就是说交替轮换算精度。例7.6矩形截面的柱体扭转问题用应力函数表示的应变余能泛函为ak(xn)中的变量来提高计22()4abxy

a

b

dxdy

取一阶近似解为(x,y)(b2y2)u(x)22

v(y)(by)满足yb边界上的边界条件(x,y)0这里新的泛函为*(u)该泛函的Euler方程为16a15a32163b5u'28bu33budxu\"(x)

相应的边界条件为55u(x)2b22b2

u(a)0从而得到一阶近似解为(x,y)(b2y2)u(x)这里u(x)ch1kachkx1,k

也可以直接用里兹法进行求解,取1b52(x,y)(b2y2)(a2x2)Amnxmyn对于矩形截面,取一项(x,y)(b2y2)(a2x2)A00代入到22

()4ab

xy

a

b

dxdy

令d

0dA00

得到相应的解。7.3伽辽金法（Galerkin）

伽辽金法的基础是虚功原理（虚位移原理）：一个平衡力系在任何虚位移中，外力的虚功等于虚应变能TTT

dpudBfud

B



(7.3.1)其中f为体积力，p为表面力，为与外力系所平衡的可能应力。u为虚位移，为与虚位移对应的虚应变。如果在前面的恒等式(4.2.1)中取,uu，得到Td[E(n)]TudB[E()]TudB从而有TT[E()f]ud[pE(n)]udB0B(7.3.2)如果选择的位移试函数不仅满足位移边界条件，同时还满足力的边界条件，即E(n)p，B2上B1上u0，那么就有T[E()f]ud0(7.3.3)和里兹法一样，我们取位移的试算函数为aiuiu

i1

n

那么虚位移为aiuiu

i1

n

因为ai是的变分，因而有[E()f]ud0

i



T

(7.3.4)这里是用位移来表示的应力。更一般地讲，如果微分方程为L(u)f

取试算函数为~auuii

i1

n

其中ui满足所有的边界条件。那么伽辽金法的积分为)fudx0,L(u

i

i1,2,...,n

与里兹法不同的是，在伽辽金法中位移试算函数除了要满足位移边界条件外，一般还需要满足应力边界条件(也就是说要满足所有边界条件)。如果试算函数不满足应力边界条件，那么近似计算的结果可能不是很理想，甚至可能完全是错误的。有时我们也可以根据问题的需要，在伽辽金积分式中添加一个权函数W(x)0以保证收敛速度)fWudx0,L(uii1,2,...,n该方法也称为加权残值法。例7.7y\"y2x,y(0)y(1)0取试函数为yk(x)(1x)xk,

两阶近似为k1,2,

~y(x)a1(1x)xa2(1x)x2

那么~y\"(x)~y(x)2x(2xx2)a1(26xx2x3)a22x代入到伽辽金方程组中(2xx(2xx

1010

2)a1(26xx2x3)a22xx(1x)dx0)a1(26xx2x3)a2

2

2

2xx

(1x)dx0

积分并求解得到14

a1142369,a241于是近似解为214(x)142y369(1x)x41(1x)x例7.8长度为l，抗弯刚度为EI的悬臂梁，受均布载荷q的作用。图7.2例7.8图梁的平衡方程为2

EIw\"12q(lx)0

或者d4w

EI4q0dx

取试算函数为x

w\"(x)a1sin

2l

那么该试函数满足右端的边界条件(M0,Q0)

EIw\"(l)0,EIw'''(l)0

把位移的表达式积分,再根据左边的边界条件(w0,w'0)得到122l2x2l

w(x)a2xsin

2l

21EIw\"q(lx)20w(x)dx0从中可以求得a，从而得到相应的近似解。lx

这样的试函数可以满足所有的力和位移边界条件。把它代入到伽辽金积分式中7.4有限元法

前面讲到，基于最小势能原理的里兹法要求试算函数在整个区域内满足位移边界条件，这对于一些形状比较复杂或者边界条件比较复杂的问题就很难处理，有限元法正好能弥补里兹法的这种缺陷。下面我们用最小势能原理所对应的位移有限元方法加以说明。假如我们在整个区域内取一些点(称为节点)，用这些节点把整个区域划分成一个个子区域（称为单元），单元和单元之间通过节点连接。图7.3单元和节点假设这些节点的位移已经知道，我们记为d，它是所有节点位移分量所组成的一个向量。每个单元（子区域内）包含其中某几个节点，单元内节点的位移向量为d，显然它是总体位ee

wwd移向量的一部分。在单元内我们对单元的位移进行插值，也就是说单元内的位移可e以表示为weNede

e

(7.4.1)其中N称为单元形函数矩阵，它需要满足一些特别的性质(见有限元的相关著作)。当然，这样得到的位移在单元内是连续的。因此整个区域内的位移w也可以表示为总体节点位移d的一个插值，也就是说wNd其中N为形函数矩阵。对应的应变为TET()wE()Nd(7.4.2)根据本构关系可以得到应力应变关系为D其中D为弹性矩阵，是对称正定。整个区域的总势能为(7.4.3)T这里，12DdVuTfdVuTpdSuTGi(xxi)dVif,p,Gi分别为体积力，表面力和作用在位置xi的集中力，(x-xi)是脉冲函数。把TT

1dKddF2位移表达式代入总势能中可以得到(7.4.4)T

其中K

12TTE()NDE()NdV

称为总体刚度矩阵，它是对称。而FNTfdVNTpdSNTGi(x-xi)dVi称为等效节点力。根据最小势能原理，弹性力学的解应使总势能取极小值，现在总势能是关于节点位移d的函数，因此有(d)TKd(d)TF0也就是说KdF(7.4.)这就是有限元的平衡方程。在代入边界条件（也就是边界上节点位移）后可以得到一个非奇异的线性代数方程组，求解该方程组后得到所有节点的位移，进而求得每个单元内的应变和应力。当然在实际求解的时候，往往是先在单元内进行插值wNd，再求得单元的刚度矩阵K和等效节点力向量F，然后再把单元刚度矩阵组装成总体刚度矩阵K，把单元等效节点力向量F组装成总体等效节点力向量F。因此说，有限元法可以看成是里兹法的一个推广。里兹法在整个区域内采用同一个插值e

e

e

e

e

e

函数，而有限元方法则把整个区域分成一个个子区域，然后在每个区域内对位移分别进行插值。由于有限元插值的区域较小（每个单元），因此插值函数的形式相对可以比较简单，不需要很高的精度。以上我们讨论的是根据最小势能原理得到的位移有限元方法，同样可以得到基于最小余能原理的应力有限元法。在实际应用中还可以得到基于广义变分原理的混合有限元法，他们的特点是包含两种或者两种以上的场变量。同样可以放松有限元在边界上的连续条件，通过引进交界面上的场变量来建立修正的广义变分原理，得到相应的杂交有限元(见7.6节)。关于有限元方法的详细讨论可以参考相关的著作(如辛可维奇的有限元著作)。7.5有限元方法的收敛性

先看一个简单的例子。细绳在拉力N作用下的（小）挠度u可以用二阶微分方程来描述d2u

N2p(x),u0atx=0,ldx

这里p(x)为横向分布载荷。上述问题也可以化为下列泛函的极小问题来描述l1du

(u){N()2p(x)u}dx

02dx

(7.5.1)(7.5.2)假定极小问题的解为u，r(ur)min(u)1uH0(7.5.3)1

H0

是下列函数的集合∶在[0,l]上存在平方可积的一阶导数,且u(0)u(l)0。在广的意义下，(7.5.1)和(7.5.3)是等价的。现在用有限元方法求(7.5.3)的近似解。首先，我们指出，要使(7.5.2)有意义的一个充分条件是u在[0,l]上连续、且存在分段连续的一阶导数。符合这一条件的有限元，称为协调元。将[0,l]分为n个线段[xi1,xi],

这里i1,2,,n

x00,x(n)l。每个线段上用线性插值函数u(x)

xixxxi1

ui1ui,i1,2,,n

xixi1xixi1

(7.5.4)这样得到的函数满足上述的协调要求。取上式中的ui为（7.5.3）中极小问题的真解在节点上的值xixrxxi1rui1ui,i1,2,,nxixi1xixi1u()

(xxi1)(xxi)(xi-xi1)2,xi1xxi2!(7.5.5)uiur(xi)uir,(x)u有(x)ur(x)u

代入(7.5.2))(ur)(d2),(udmaxxixi1i(7.5.6)将(7.5.4)代入(7.5.3),则把一个泛函求极值问题化为函数求极值问题，记此时的解为uiuie,ue(x)xixexxi1eui1ui,i1,2,,nxixi1xixi1e(u)。这样由势能的极小原理可得其相应的势能为)(ur)(d2)(ur)(ue)(u

即(7.5.7)lim(ue)(ur)

d0

1

H0意味着有限元的解（在上）收敛于真实的解。这里我们用到了(7.5.8)(1)泛函(7.5.2)是正定的（即解是极小问题的解）；(2)插值函数估计式(7.5.5)成立(完备性)；(3)泛函(7.5.2)的积分可以化为各单元的积分之和（协调性）。一般来说，以上三条是保证有限元收敛性的充分条件（注意∶不是必要条件），不过具体形式可以随问题的不同而不同。下面看看薄板弯曲中的有限元收敛性的条件。由(7.6.2)可知22222222DwwwwwU222v22(1v)dxdy22xyxyxy显然其总势能是正定泛函，满足上述第1条正定性的要求。3

(d)，则满足完备性的要w如果我们在单元内的插值函数具有二次精度、即余项为求。这个要求也可以等价于∶每个单元内的插值函数包含一个完备的二次多项式。而前面例子提到的拉格朗日插值，只包含一个完备的一次多项式，所以线性插值函数不能用于薄板弯曲问题。为了满足协调性要求，插值函数仅仅在整个区域上连续是不够的，还需要一阶导数在整个区域上连续。换言之，插值函数不仅在单元边界上需连续，而且插值函数的法向导数在单元边界上需连续。这个要求给构造板的有限元带来了很大困难，从而吸引了很多从事这方面的工作。工作分两方面进行∶一是寻找满足协调性的单元，这往往带来很繁琐的公式；另一方面放松协调性的要求，构造不满足协调要求、但仍能保证收敛的非协调元，因为协调性仅仅是收敛的充分条件、而非必要条件。7.6应力杂交元

为了解决上述满足协调性要求的困难，卞学璜提出了应力杂交元。其基本思想是∶单元内部以满足平衡方程的应力作为基本变量，而单元边界上以假定位移作为未知量。譬如在薄板弯曲问题中，设单元为四边形或三角形单元，每个节点上取定在单元边界上位移w是三次多项式（可由边两端的六个参数是线性函数（可由边两端的四个参数w,wx,wy

为待定参数。假w,wx,wy

表示），法向导数wn

wx,wy表示），这样保一证了相邻单元的位移和导数是连续的。此外在每个单元内部，内力满足2Mxy2My2Mx

2q022

xxyy

MxyMyMxMxy

Qx,Qyxyxy

板的两类变量的广义余能为(7.6.1)2Mxy2My2Mxw2[V(22q)w]dxdyMdsnxxyy2nC1C1C2(MnswMnsQn)wds(MnMn)ds(QnPz)wdssnsC2C3C3wMnswVdxdy[M(Q)w]dsMdsnnnnsneeC1C2C3eC1MnswMns(Qn)wds(MnMn)ds(QnPz)wdssnsC1C2C2C3C3式中(7.6.2)C1(固支边):ww,ww

nn；MnMn；MnsQnPzs。C2(简支边):ww,C3(自由边):我们将MnMn,2分成两部分∶一部分是内力变分引起的，另一部分是由单元边界上的位移和20

位移法向导数的变分引起的。由于单元内的内力变量与周围单元是的，所以由得到的内力变分部分可由每个单元上得到e{VdxdyeeC1C2C3[MnwMns(Qn)w]ds}0,foreveryens也就是将每个单元视为已知位移边界的板弯曲问题，然后用最小余能原理，将内力用单元边界上的位移和位移法向导数表示。将每个单元的内力代入(7.6.2)，可以得到用单元边界上的位移和位移的法向导数表示(可进一步用每个节点的位移参数法。w,wx,wy

表示)的总余能（实质上已变成势能），最后利用变分原理可求得每个节点的位移参数。这就是应力杂交元