样条有限元条法

1、第6章 样条有限元条法 6.1绪论20世纪60年代在R.W.CIough命名后迅速发展起来的有限单元法，被公认为结构分析 最强有力的工具，理论上可适用于所有的结构。但对于规则区域的结构，有限元法的解题效 率不如差分法；在计算机技术取得突破性进展之前，解决大型空间结构工程设计问题，往往需要大型计算机和昂贵的费用，使有限元法受到种种限制；有些结构问题用有限元法也难以解决，需要创造一些新的数值方法，于是又产生了有限条法、边界元法、样条元法以及新的 加权残数(余量)法等等。适用于几何形状规则和边界条件简单的土建结构的一种行之有效的特殊有限元法一一 半解析有限条法，是由YKCheung (张佑启)和 E

2、.L.Wilson在1968年提出的17。与有限元法相比，主要不同点在于所取得位移函数，一般是以多项式和正交级数乘积的形式给出， 使得弹性力学问题降维，从而使总刚度矩阵大大降阶，既省机时，精度也高。但对于非简支 端边界条件，级数耦合，计算繁杂，尤其是对于集中荷载作用和内部支承情况，所取项数多，收敛慢。对于处理沿跨向材料突变或是变截面问题，结果也不能令人满意。1979年石钟慈提出样条有限元法18，用三次B样条变分方法解规则区域上板梁组合弹 性结构的平衡问题，导出适用于各种边界条件的统一计算格式，比通常的有限元法计算量少，精度高，便于在小型计算机上实现。稍后，秦荣提出以样条函数、梁振动函数及能量变

3、分为 基础的样条有限点法和后来的样条子域法19。1982年Y.K.Cheung等人又提出结构分析的样条有限条法，以克服经典有限条法的缺点。作者在分析空间悬挂结构吊桥时，采用元条相结合的方法20，将主索刚度叠加于桥面板条单元的结线上，分析索结构的动力特性取得成果。 6.2 B样条与样条函数样条函数是现代函数逼近的一个十分活跃的分支，是计算方法的一个重要基础，已得 到广泛应用。样条函数来源于实际生产中的样条曲线，由于样条函数是一个分段多项式，利用它去逼近任意函数，具有更大的灵活性和适应性。样条函数可以利用基本样条(如B样条) 函数的线性组合来构造，在计算力学中，常用B样条函数来构造位移函数及应力函

4、数。 B样条及其性质B样条函数是基本样条函数(简称 B样条)，它与广义函数一一3函数有密切的内在联 系，并可以逼近3函数， 所以也称3样条函数。B样条函数系可作成样条函数空间的一组基。 n次B样条函数可写成统一的表达式n 1k / n 1n #n(x)( 1) (k )(x Xk) /n!(6.2.1)k 0式中，(n 1)是二项式系数的组合表示，样条结点xk k 2, k 0,1,2, n 1 。2吊桥动力分析主要应用一次和三次B样条，其分段表达式和曲线形状如下：X 1X1,0(X)1 XX0, 1(6.2.2)0X13(X)(X 2)3(x1(2(22)3x)34(x 1)334(1 x)

5、(a) 一次B样条x)3x 2, 1x 1,0x 0,1x 1, 2(6.2.3)(b)三次B样条图6.2.1 B样条曲线n次B样条函数具有紧致性(即$ n(x)及其导数以区间此外，处处为零)和分段光滑性(为n次分段多项式，在(n 1T ,-oo,m)n 1亠一1 为紧致支撑集，2区间，具有 n-1阶连续导数)以及对称(即$n(-X)= $ n(x)和面积不变(亦即n(x)dx 1 )等特性。尤其以及最佳逼近性质，是奇次样条函数具有极小模性质(可用来量度变形能极小性质)很大的数学与力学的使用价值。 样条函数样条函数曲线来源于生产实践中的数学放样。定义一次样条函数为二阶广义微分方程显示了的解，即

6、y (x)j (x Xj)jy(x)式中，Xj为结点，下标“ +01Xj(Xj为截断符号，即定义：Xj)(6.2.4(X Xj )n(X Xj)nXj 0(6.2.5Xj 0可见，一次样条函数是连续， 它表征了缆索受吊杆作用力时的实际情况。我们也可以定义三次样条函数为四阶广义微分方程(4).y (x)且分段线性的折线函数,其力学意义是集中力作用下的弹性弦,j (xjXj)的解，其一般表达式为y(x) 01X1 j(x Xj)33!(6.2.6)j=0时，便退化为普通的三次多项式)显然，这是一个分段三次多项式(当B，有二阶连续导 数，三阶阶跃导数。如果与梁的挠曲线方程y(x)=q(x)对照，可知

7、三次样条函数相当于弹性梁在结点X处受集中荷载B j作用下的挠曲线。因此，可用以反映吊桥加劲梁受吊杆力作 用的实际情况。同样，我们可以定义更高次的样条函数，例如，在研究连续应力场和应变场时，定义了五次样条函数21。二次样条函数也同样有其力学背景。这里不赘述。根据展开定理，n次样条函数Sn(x)可表为n次B样条函数$ n(x)及其平移$ n(x-k)的线 性组合，即NSn(x)kk n(x k) , PP 1(6.2.7当n为奇数时，若在0, I的区间上取N等分划, 于是上式可表为则h=l/N，在第k结点,有 xk=kh,对于一次和三次而 V Span组基函数，张成Sn(x)xn(-xk計)(6.

8、2.8)B样条，有Si(x)S3(X)1(h k)ko,i,k)k)n 为0,1上以 Xh(kN+1维空间。类似的，V Span(6.2.9(6.2.100,1, N)为结点的一次B样条的一3(h k)k 1,0,N1 为三次B样条的一N+3维空间。组基函数，张成构造出适合结构位移边界条件的基函6.2.10 )式改写后表示位移函数(6.2.11在线性空间等价的前提下，可改换基函数的形式, 数，使统一处理各种边界条件得到方便。例如，将(N 1W(x) k k(x)k 1由于基函数k(x)的局部紧凑，若要计算某样条结点x的位移值，上式中最多只有三项不为零，即W(Xj)i 1 i 1 (Xj )i

9、i( Xj )i 1 i 1 (Xj )(6.2.12)在求出广义参数后，经上式计算，才得到待求的位移值(与有限元不同)。由于k(x)在结点上的值k(Xi)是一些简单的现成的数，故需要的计算量很小。这些“现成数”已分几种边界情况编入子程序中。k(x)是一个与三次B样条3(- k)有关的基函数，考虑常见的边界条件，写成：1(x)3(h1)o(X)3(h)4 3(X 1) hX1X1（X）3（1)- 3(mh2 h2（X）3C2)N 2（X）N 1（X）3(X N 2)x3(h N 1)h3(X N 1)1 3(f N)h2 hN(X)3(f N) 7 N 1)N 1(X)3(- N 1)以上的基

10、函数对应于自由边界。若为简支边界，则删去第一行（列）或末一行（列） 若为固定边界，则删去首两行（列）或末两行（列）。若结构对称，可简化取半计算，将半结构分成N等分划，根据滑动固支边界（ W=0）,仅删去对称轴一端的三行三列，代以新的二行二列，即N1（x）3（7 N 1）3CX N 1）hhN（X）3（ N）h反对称时，取半结构的对称轴一端仍按铰支处理。利用参数法19，还可构造出与三次 B样条函数有关的新的基函数，得到具有物理意义的位移参数和广义参数混合的位移函数，使处理位移边界条件和连续条件都很方便。经过对B样条函数的求导和积分，计算出系数矩阵提供编程使用。 B样条函数的求导和积分B样条函数的

11、求导和积分公式为(j)k n 1 n 1 n i(6.2.13)n(X)ko(1)(k)(X T k) J/(nj)!式中，当J为正整数时，表示求J次导数；当为负整数时，表示积分（-J）次；当J=0时,以B样条基函数（或其导数）即为（6.2.1 ）式。为元素的矩阵乘积的积分是一组很有实用意义的系数矩阵, 其具体形式可事先算出，列表待查。吊桥的动力分析将用到以下的系数矩阵，记为:1y0口dy ,11 2y。门dy3y0门dy ,1 4y1t0 T dy(6.2.14)5y0门dy ,1 6y1t 1(0 dy)T(0 dy)以I 4y为例（暂不考虑修改边界），说明矩阵乘积的积分过程。0门dy 0

12、 3i） 3（h k）dy令thi ,ck i则I 4 yN i hi3(t)3(tc)dt反复进行分部积分。由于Ni331)(t)(331)(tc)dt0于是I 4 y3h (j 01)j3j)(t)(j31)(tc)N ii(6.2.15)同理，可进行其余矩阵的积分。特殊的系数矩阵（如修改边界时，可能是$3的某种线性组合）根据问题的需要区别处理，随时算出。由于三次B样条的紧凑性，所生成的系数矩阵一般是半带宽为 4的带状矩阵。 6.3样条有限元条的位移模式 吊桥简化假定以吊桥为例，在利用样条有限元条法进行吊桥动力特性分析时，作如下假定：小振幅振动，所有材料遵循虎克定律。缆索忽略弯曲刚度。初始

13、静载全部由主索承受，索呈二次抛物线的理想索形。忽略细长加劲桥面结构的剪切变形和转动惯量影响。主索的重量沿跨长均布，各跨加劲结构等截面。吊杆不可伸长，在竖直扭转振动时保持垂直。振动引起的索的附加水平张力和塔顶的附加压力可忽略不计。 忽略加劲结构的初应力和初始曲率。 结构离散与位移模式根据上述的线性化假定， 仅考虑桥面结构的弯扭振型时，悬挂的空间结构被离散为索单元，塔单元和桥面壳条元体系。1. 索作为一维样条弦杆单元考虑，采用一次B样条（由于实际吊杆密布，有时索也采用三次B样条，以便统一处理）。仅在竖直平面内有一个自由度（W）。由刚性吊杆的假定，其有W=W,这里W为加劲桥面的竖向位移。索与桥面系在

14、顺桥向都取相同的N个等分划,即有N+1个结点。位移模式为 Wc cWci，而c 1 0 1 i n 1, 矩阵中各元素i（y）是与B样条函数有关的基函数。 WCi=W-1 WWWWN+1 T ,其中 W 是待定的与时间有关的竖向结点位移参数。2. 桥面系视为一般的正交异性折板结构（桥宽不大时，也可视为一维梁），包括了板梁、箱梁和闭合桁架等典型桥面结构，并进一步离散为样条有限条的矩形壳条单元（如图6.3.1所示），坐标轴x,y为其主方向。顺桥向采用N个等分划的三次 B样条，横桥向采用二结点低阶条的常规形函数17。假定不计壳条弯曲和薄膜效应之间的相互作用，每条结线（内设结点）上有四个自由度，用位移

15、参数u,v,w, B表示。特殊情况作为板条单元，只计弯曲效应，每条结线只有两个自由度，位移参数为w, 0 （ B为埃米特插值所需）。位移模式表为：式中，A是形函数矩阵，而 S 是位移参数列矩阵。00N,0N200 0其中N0 N10 N20 000(6.3.2)0 00 0N3 N4N5N6矩阵兀素N(x)为低阶条的常规形函数，其中:Ni1 xN2xN31 3X2 2x3N4x(1:2x x2)N53X22x3N6x(X2x)xx/bb为条宽uivi0ujvjwi0iwjj子矩阵ui1 0 12iN 1NN 1ui(631)Fs U V WT N A 式中，加“”记号的基函数一般是根据边界条件

16、局部修改过的。Nd ui0N2uj000000Nvi0N2vj00000000弘wiN4【iN5wjN6 ju iv iu jv jwiiw jj T子矩阵u iu 1u0u1u2uiuN 1uNTu N 1 ij1012iN 1NTN 1 j结线i3. 桥塔被离散为四个或更多的在竖直面内挠曲的悬臂梁（柱）单元，采用三次B样条,取N个等分划，每个内结点仅由一个自由度（V）,位移模式为：Vt tVtj，这里，t和Vtj的定义类似于索单元。 6.4动力特性矩阵样条有限元条法基本上是变分问题的直接法（Ritz-Galerkin 法）和以位移逼近的有限元法的综合。 结构振动的泛函根据哈密顿原理，可知结

17、构振动的泛函为t2J (U V T)dt(6.4.1)式中：t为时间，t1 t Wt2； U为结构的应变能；V为结构的外力功；T为结构的动能。 设结构振动时的位移函数为 f t=u,v,wT, p为板壳单位面积质量， t为应变， x * 为曲率（扭率），qt为分布外动荷载，则1 a bttU 2 0 0 ( T R t T D t)dxdy(6.4.2)a b t丄 dxdyV o o f t q t dxdy(6.4.3)(6.4.7)(6.4.4)把（6.4.5 ）式写成Jr ,由变分原理J 0转化为函数的极值问题,函数，在某一周期 T=2n / 3内对时间t积分，可得结构自由振动的泛函为

18、a bJ( T2 0 0、RT D2 f T f )dxdy(6.4.5)代入吊桥结构的位移模式A00f0c0Wc(6.4.6)0 0tVt研究无阻尼自由振动时，动荷载 qt=0。由于微幅谐振动，振动的位移函数是一个周期TWcTVtT T得到吊桥不计阻尼自由振动的动力矩阵方程:2(KM)r0求解该广义特征问题，得出结构的固有频率及其相应的振型。 动力特性矩阵在给出壳条或单元在平衡位置的势能和动能的表达式后， 条和元的刚度矩阵s和一致质量矩阵m。索的刚度、由影响系数的定义，容易得到 质量与条的刚度、质量在结线上叠力口；塔的刚度、质量可作为独立的矩阵块。通过坐标变换和根据边界条件局部修改样条函 数

19、后，集成结构总刚度矩阵K和总一致质量矩阵对于矩形扁壳条，广义应变为M。xyxyx其中，应变矩阵：2w2xuiviN2uj0vjuiviujvjwiwjwi2N3wi2N42N5wj2N6广义应力为xyMxMyMxy对于正交异性材料，弹性矩阵:CxC10C1Cy0其中：xyDxD10D1Dy0D xyDxExt3ExtCyEytC1xCyyCx , CxyGt,12(1DyEyt312(1以上各式中,Ex、Ey、x y )y和G是弹性常数，D1xD yxyt是条的厚度。当各向同性时,Gt3Ex=E=E,壳条应变能：Us i2 V(i)壳条刚度矩阵BT D BdAB dAB dxdysiiSI2S

20、21上式各分块矩阵中，S i2=S2i= 0sii其中:kiiki2k14ki3k58kiik33ki2k22ki3k23k33ki4k24k34k44k55k56k66443 y3y bk34C 11 5y2(k23)T2C-1 Cx 1 4yCxy15yy 3y6CylI 3yCx 1 4yb-Cxy(l 5y)TbC l3xy 3y3-丨xy 4y丄C丨xy 4ybk77-2d i.3 D x 1 4ybCxy(I 5y)T旦DyI35 yy 1 iy5bDxy|3yb3 D |105 y iy iib2 D |D y 1 1 y210 y y 13b24204b _ .Di l2yi5

21、 y7 Di1 2y5 iDiI54DxI4y b2 DxI b2色D I2 Dx I 4yb3 DxI 4yDi I 2yb5b2b .b y i5 yk88k78k67124y2yDy Iiy8b,Dxy I 3yi52 D Ixy 3y52D丨Dxy I 3y59b , 7oDy|iyb3 D |D y 1 iy140 y y以上各式中的系数矩阵，仅l5y为反对称矩阵，(2)壳条的一致质量矩阵l b壳条动能：TstS 2 0 024,Dxy 1 3y5b2b,Dxy 1 3yi5其余均为对称矩阵。FS T FS dA -S S 2dxdymlt0 bTo A A dxdy01m21m12

22、m22(6.4.10)各分块矩阵为:m12=m 21= 0 ,201015622b5413dtb201tb4b213d3b2mn-6 对20 I4y，mb2420对15622b 4y称2称4b2桥面系沿横向可以有不同方位布置的条元，如图641所示，条元的局部坐标系 y轴和结构整体坐标系 y轴平行(相邻条交线重合)为x和x轴之间的夹角(顺时针为正)。壳条坐标转换矩阵为TT11T12(6.4.11)T21T22cos0sin010其中T1 2221T2cossin0100y(v)图6.4.1坐标转换z(w)(3)索的弹性刚度矩阵将索的位移模式代入索张力变化所储存应变能表达式，导出索的弹性刚度矩阵为

23、l式中，a0EcAc /WTwi dy , Le(HW)23secW为单位长度的索重加桥面结构重,EcA / W、2 a ()I Hewdy为索的等效长度,Hw为索的静水平张力。(6412)I6y为一系数矩阵,(4)索的重力刚度矩阵索的动伸长对索挠度的影响与索竖向振幅相比是高阶微量，可忽略不计。代入索的位移模式和相容方程，索的重力势能可表达为1 11 2Vcg0W Wcdy-Hw0( Wc/ y)2dy11 TTHw 0 W iwiwi W idy则索的重力刚度矩阵为Scg(5)索的一致质量矩阵1索的动能表达式：Tc-c 2则索的一致质量矩阵为mcHw owiwi dyH w l3y1T11T

24、wco cc Wc dy -cW ic 0 iiTc owiwi dyc1 4yTwi(6413)wi W i dy(6414)(6)塔柱的刚度矩阵塔柱取等效系统分析， 假设其弯曲和扭转由各柱在竖直面内的挠曲形成。考虑塔柱的应变能和外力势能，得到塔柱弹性刚度矩阵为htSte0EJt tj tjdz Etlthz(6415)式中，l1Z(及下式I3Z、|4Z)为一系数矩阵，参照|1Y等。塔柱变刚度时，若惯性矩It为Z的连续函数，则可参与积分，给出新的系数矩阵；若刚度呈阶梯形变化，则应分段，各等分划 处理。塔柱的几何刚度矩阵考虑塔顶的静、动压力hTStg0 Pwtjtj dzRw和P(t)以及塔顶

25、等效弹簧刚度k,写成kPwl3z k(6.4.16)(7)塔柱的一致质量矩阵考虑单塔柱的动能，得到单塔柱的一致质量矩阵表达式：h Tmt t 0 tj tj dz tl4z(6.4.17)对于三跨简支吊桥，如果桥面结构的材料、构造都相同，可按跨分段，各取等分划。根据简支边界条件局部修改 B样条后，三跨的基函数集成对角分块矩阵， 其阶数可随各跨分划 数的不同而不同，样条矩形壳条的刚度矩阵形式同单跨， 系数矩阵按跨分块，集成对角带状 矩阵；索的弹性刚度矩阵也由三跨分块集成，对称局部满阵。位移参数列阵也由各跨集成。对于三跨连续的吊桥桥面结构，可在中间支承处分段。为方便起见，系数矩阵需经修改，构造出具

26、有混合参数的位移函数，在此不赘述。 6.5算例例1用样条有限元(按一维工程梁)计算如图6.5.1所示单跨吊桥的竖向弯曲振动。特性数据如下：加劲梁：L=853.44m, W=4246.016 kg /m,巳=203932.16MN/m； ls=7.695953m422索：f=70.714m, A c=0.123548m , E c=179129.6MN/m , L E=1219.2m, H w=53578kN.图6.5.1单跨吊桥对于一维梁单元，设 W=w刚度矩阵：SEc Ac匸WHT)21 6y致质量矩阵:mWT4y几种方法计算结果示于表6.5.1。比较表明，样条有限元法具有较高的计算精度。表6.5.1单跨吊桥竖向振动频率频率序自振频率(rad/sec )对称振型反对称振型解析解有限元S-12S-16解析解有限元S-12S-1611.4001.39751.39741.39741.33181.33301.33041.330422.6962.70472.70222.70224.49014.