样条有限元条法

1、第6章 样条有限元条法§6.1 绪论20世纪60年代在R.W.Clough命名后迅速发展起来的有限单元法，被公认为结构分析最强有力的工具，理论上可适用于所有的结构。但对于规则区域的结构，有限元法的解题效率不如差分法；在计算机技术取得突破性进展之前，解决大型空间结构工程设计问题，往往需要大型计算机和昂贵的费用，使有限元法受到种种限制；有些结构问题用有限元法也难以解决，需要创造一些新的数值方法，于是又产生了有限条法、边界元法、样条元法以及新的加权残数（余量）法等等。适用于几何形状规则和边界条件简单的土建结构的一种行之有效的特殊有限元法半解析有限条法，是由Y.K.Cheung（张佑启）和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样条），它与广义函数函数有密切的内在联系，并可以逼近函数，所以也称样条函数。B样条函数系可作成样条函数空间的一组基。n次B样条函数可写成统一的表达式 (）式中，是二项式系数的组合表示，样条结点 。吊桥动力分析主要应用一次和三次B样条，其分段表达式和曲线形状如下： (） () (a) 一次B样条 (b) 三次B样条图 B样条曲线n次B样条函数具有紧致性（即n(x)及其导数以区间为紧致支撑集，此外，处处为零）和分段光滑性（为n次分段多项式，在（-，）区间，具有n-1 阶连续导数）以及对称（即n(-x)=n(x)）和面积不变（亦即）等特性。尤其是奇次

5、样条函数具有极小模性质（可用来量度变形能极小性质）以及最佳逼近性质，显示了很大的数学与力学的使用价值。 样条函数样条函数曲线来源于生产实践中的数学放样。定义一次样条函数为二阶广义微分方程 的解，即 (）式中，xj为结点，下标“+”为截断符号，即定义： (）可见，一次样条函数是连续，且分段线性的折线函数，其力学意义是集中力作用下的弹性弦，它表征了缆索受吊杆作用力时的实际情况。我们也可以定义三次样条函数为四阶广义微分方程 的解，其一般表达式为 ()显然，这是一个分段三次多项式（当j=0时，便退化为普通的三次多项式），有二阶连续导数，三阶阶跃导数。如果与梁的挠曲线方程y(4)(x)=q(x)对照，可

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

7、空间。在线性空间等价的前提下，可改换基函数的形式，构造出适合结构位移边界条件的基函数，使统一处理各种边界条件得到方便。例如，将（）式改写后表示位移函数 (）由于基函数的局部紧凑，若要计算某样条结点xi的位移值，上式中最多只有三项不为零，即 (）在求出广义参数后，经上式计算，才得到待求的位移值（与有限元不同）。由于k(x)在结点上的值是一些简单的现成的数，故需要的计算量很小。这些“现成数”已分几种边界情况编入子程序中。是一个与三次B样条有关的基函数，考虑常见的边界条件，写成： 以上的基函数对应于自由边界。若为简支边界，则删去第一行（列）或末一行（列）；若为固定边界，则删去首两行（列）或末两行（列

8、）。若结构对称，可简化取半计算，将半结构分成N等分划，根据滑动固支边界（W=0），仅删去对称轴一端的三行三列，代以新的二行二列，即 反对称时，取半结构的对称轴一端仍按铰支处理。 利用参数法19，还可构造出与三次B样条函数有关的新的基函数，得到具有物理意义的位移参数和广义参数混合的位移函数，使处理位移边界条件和连续条件都很方便。经过对B样条函数的求导和积分，计算出系数矩阵提供编程使用。 B样条函数的求导和积分B样条函数的求导和积分公式为 ()式中，当j为正整数时，表示求j次导数；当j为负整数时，表示积分(-j)次；当j=0时，即为（）式。以B样条基函数（或其导数）为元素的矩阵乘积的积分是一组很有

9、实用意义的系数矩阵，其具体形式可事先算出，列表待查。吊桥的动力分析将用到以下的系数矩阵，记为： ()以I4y为例（暂不考虑修改边界），说明矩阵乘积的积分过程。 令 则 ，反复进行分部积分。由于 于是 ()同理，可进行其余矩阵的积分。特殊的系数矩阵（如修改边界时，可能是3的某种线性组合）根据问题的需要区别处理，随时算出。由于三次B样条的紧凑性，所生成的系数矩阵一般是半带宽为4的带状矩阵。§6.3 样条有限元条的位移模式 吊桥简化假定以吊桥为例，在利用样条有限元条法进行吊桥动力特性分析时，作如下假定： 小振幅振动，所有材料遵循虎克定律。 缆索忽略弯曲刚度。 初始静载全部由主索承受，索呈二

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

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

12、每条结线只有两个自由度，位移参数为w,(为埃米特插值所需)。位移模式表为： ()式中，A是形函数矩阵，而是位移参数列矩阵。其中 (）矩阵元素Ni(x)为低阶条的常规形函数，其中： 子矩阵 式中，加“”记号的基函数一般是根据边界条件局部修改过的。 子矩阵 图 样条矩形壳条单元 3. 桥塔被离散为四个或更多的在竖直面内挠曲的悬臂梁（柱）单元，采用三次B样条，取N个等分划，每个内结点仅由一个自由度（Vt）,位移模式为：，这里，和Vtj的定义类似于索单元。§6.4 动力特性矩阵样条有限元条法基本上是变分问题的直接法（Ritz-Galerkin法）和以位移逼近的有限元法的综合。 结构振动的泛函

13、 根据哈密顿原理，可知结构振动的泛函为 ()式中：t为时间，t1tt2；U为结构的应变能；V为结构的外力功；T为结构的动能。设结构振动时的位移函数为ft =u，v，wT,为板壳单位面积质量，t为应变，t为曲率（扭率），qt为分布外动荷载，则 () () ()研究无阻尼自由振动时，动荷载qt=0。由于微幅谐振动，振动的位移函数是一个周期函数，在某一周期T=2/内对时间t积分，可得结构自由振动的泛函为 ()代入吊桥结构的位移模式 ()令 ()把（）式写成，由变分原理转化为函数的极值问题，即 ，得到吊桥不计阻尼自由振动的动力矩阵方程： ()求解该广义特征问题，得出结构的固有频率及其相应的振型。 动力

14、特性矩阵 在给出壳条或单元在平衡位置的势能和动能的表达式后，由影响系数的定义，容易得到条和元的刚度矩阵s和一致质量矩阵m。索的刚度、质量与条的刚度、质量在结线上叠加；塔的刚度、质量可作为独立的矩阵块。通过坐标变换和根据边界条件局部修改样条函数后，集成结构总刚度矩阵K和总一致质量矩阵M。对于矩形扁壳条，广义应变为其中，应变矩阵： 广义应力为 对于正交异性材料，弹性矩阵： 其中: 以上各式中，Ex、Ey、x、y和G是弹性常数，t是条的厚度。当各向同性时，Ex=Ey=E，。壳条应变能： （1）壳条刚度矩阵 ()上式各分块矩阵中，s12=s21= 0 其中：以上各式中的系数矩阵，仅I5y为反对称矩阵，

15、其余均为对称矩阵。（2）壳条的一致质量矩阵壳条动能：式中，是单位体积的重量，t为板厚，A=N。故壳条的一致质量矩阵为 ()各分块矩阵为： m12=m21= 0， ， 桥面系沿横向可以有不同方位布置的条元，如图所示，条元的局部坐标系y轴和结构整体坐标系y轴平行（相邻条交线重合），为x和x轴之间的夹角（顺时针为正）。壳条坐标转换矩阵为 (6.4.11)其中 图 坐标转换 （3）索的弹性刚度矩阵将索的位移模式代入索张力变化所储存应变能表达式，导出索的弹性刚度矩阵为 ()式中，为索的等效长度，为倾角，为一系数矩阵，为单位长度的索重加桥面结构重，为索的静水平张力。 （4）索的重力刚度矩阵索的动伸长对索挠

16、度的影响与索竖向振幅相比是高阶微量，可忽略不计。代入索的位移模式和相容方程，索的重力势能可表达为 则索的重力刚度矩阵为 () （5）索的一致质量矩阵索的动能表达式： 则索的一致质量矩阵为 () （6）塔柱的刚度矩阵塔柱取等效系统分析，假设其弯曲和扭转由各柱在竖直面内的挠曲形成。考虑塔柱的应变能和外力势能，得到塔柱弹性刚度矩阵为 ()式中，I1Z（及下式I3Z、I4Z）为一系数矩阵，参照I1Y等。塔柱变刚度时，若惯性矩It为z的连续函数，则可参与积分，给出新的系数矩阵；若刚度呈阶梯形变化，则应分段，各等分划处理。塔柱的几何刚度矩阵考虑塔顶的静、动压力PW和P(t)以及塔顶等效弹簧刚度k,写成 (

17、) （7）塔柱的一致质量矩阵考虑单塔柱的动能，得到单塔柱的一致质量矩阵表达式： ()对于三跨简支吊桥，如果桥面结构的材料、构造都相同，可按跨分段,各取等分划。根据简支边界条件局部修改B样条后，三跨的基函数集成对角分块矩阵，其阶数可随各跨分划数的不同而不同，样条矩形壳条的刚度矩阵形式同单跨，系数矩阵按跨分块，集成对角带状矩阵；索的弹性刚度矩阵也由三跨分块集成，对称局部满阵。位移参数列阵也由各跨集成。对于三跨连续的吊桥桥面结构，可在中间支承处分段。为方便起见，系数矩阵需经修改，构造出具有混合参数的位移函数，在此不赘述。 §6.5 算例 例1 用样条有限元（按一维工程梁）计算如图所示单跨吊

18、桥的竖向弯曲振动。特性数据如下：加劲梁：L=853.44m, W=4246.016/m, ES=203932.16MN/m2, IS=7.695953m4 索：f=70.714m, AC=0.123548m2, EC=179129.6MN/m2, LE=1219.2m, HW=53578kN.图 单跨吊桥对于一维梁单元，设W=w刚度矩阵: 一致质量矩阵：几种方法计算结果示于表。比较表明，样条有限元法具有较高的计算精度。表 单跨吊桥竖向振动频率频率序自 振 频 率 （rad/sec）对 称 振 型反 对 称 振 型解析解有限元S-12S-16解析解有限元S-12S-1611.4001.39751

19、.39741.39741.33181.33301.33041.330422.6962.70472.70222.70224.49014.48704.48854.484936.84726.85656.84369.71399.71639.76289.7073413.118713.270413.136217.020517.046217.456017.0376521.515322.456321.683826.412826.563628.360126.5569注：表中S-12和S-16分别表示样条有限元条法取12分划和16分划。例2 三跨简支吊桥旧Tacoma桥（该桥于1940年在不到20m/s的风速下发生颤振破坏）主要设计数据如下：ES=2.038318