基于精确几何的曲梁分析新方法

, 　,

(长江科学院 材料与结构研究所，武汉　430010)

1　研究背景

梁、板、壳是工程结构的常见形式，其中，曲梁和曲壳由于其优美的流线造型和良好的受力特性，在工程中有着广泛的应用。

目前，以有限元法为代表的数值计算方法是梁板壳力学分析的主要方法[1]。梁板壳分析一般基于2个基本假设：①对于原先垂直于中面的横截面，在变形过程中始终保持为平面；②沿厚度方向产生的直接应变可以忽略。对于细长梁和薄板壳，可进一步忽略剪切变形，在基本假设①的基础上考虑“变形后的截面仍垂直于中面”的“直法线”假设。具体分析中，大都以中面位移为未知量，考虑位移沿厚度方向的线性变化，由此推导出不同于实体分析的控制方程，一般要求的C1连续性给构造近似函数带来一些困难；另一种方式是位移和转角独立插值，或直接从实体单元退化得到梁板壳单元，仅要求C0连续，并能考虑厚度较大时的剪切变形，但在求解细长梁或薄板壳时易出现剪切自锁。

对于曲梁和曲壳而言，由于存在初始曲率，其受力分析比直梁和平板要复杂得多。参考文献[1]—文献[4]，曲梁和曲壳的有限元分析主要有2种途径：

(1)一种简单方式就是用直梁或平板单元近似地模拟曲梁或曲壳。随着单元数量的增加，在几何上和变形上都会逼近真实的曲梁或曲壳。这种方式固然简单，但为保证几何相似性会使用较多单元，而且难以很好地反映曲梁或曲壳实际存在的一系列耦合问题，造成求解过程中收敛慢和计算误差大。

(2)另一途径是直接构造曲梁和曲壳单元。但由于曲率的存在，需要对任意形状(常曲率、变曲率)下的应变-位移关系进行精确描述，还要考虑轴力、弯矩、剪力、扭矩等作用的耦合，甚至要涉及三维空间的非欧几何，几何方程、物理方程和平衡微分方程都比较复杂，还存在不同假设或简化条件下的壳体理论，因此，曲梁和曲壳单元的研究难度远大于直梁和平板单元。

近年来兴起的等几何分析方法[5]，几何模型和分析模型采用统一的NURBS(Non-Uniform Rational B-Splines，非均匀有理B样条)描述方式，实现了几何模型的保形性，计算精度高，适合于对几何精确性要求高的曲梁和曲壳分析。但等几何分析方法也存在基函数相对复杂、本质边界条件不易施加、提高阶次也不能完全避免薄壳计算中的剪切自锁等问题[6]。

考虑到曲梁和曲壳有着相似的变形特征，曲梁可看作曲壳在空间上的退化形式，因此，本文仅讨论二维曲梁，但本质上是着眼于曲梁和曲壳分析的一般性方法。

笔者在文献[7]中基于“独立覆盖流形法”，提出梁的独立覆盖分析方法，采用实体计算退化方式，使某些多项式基函数不参与计算，就能模拟梁的基本假设，而且仅要求近似函数的C0连续并从根本上解决了剪切自锁等问题。在此基础上，本文将这种方法推广到二维曲梁，直接在精确的中面曲线参数方程基础上计算曲梁的变形，为三维曲梁和下一步的曲壳分析探索新的路径。为简便起见，本文以单位宽度的矩形截面梁为例，且暂时只考虑中面用一条曲线描述的情况。

2　独立覆盖流形法简介

图1　任意形状的独立覆盖及其条形连接Fig.1　Independent coverswith arbitrary shapeconnected through strips

图2　曲梁的独立覆盖划分及其条形连接Fig.2A curved beam divided by independent coversconnected through strips

2012年，作者在石根华博士的建议下，在流形法中首次引入“独立覆盖”[10-11]，即单位分解函数φi=1、近似函数V就是给定级数Vi的覆盖区域，独立覆盖之间为窄条形的覆盖重叠区域，其单位分解函数取为有限元线性形函数以实现覆盖函数的线性过渡。这样，每个覆盖就包含1个独立覆盖及其条形重叠区域。首先研究了矩形独立覆盖，2013年提出任意形状的独立覆盖[12]，如图1所示，3个任意形状的独立覆盖用条形(包括条形交叉处的过渡单元)连接，这样，覆盖的划分不仅能适应求解域的物理边界，还能严格施加本质边界条件。比如，将图2所示的细长曲梁划分为5个独立覆盖，其间用条形相连，条形上的单位分解函数就是一维有限单元的形函数，同时，在梁的端部通过条形与定义了边界条件的结点相连。

2015年笔者在文献[13]中正式命名“基于独立覆盖的数值流形方法”，简称“独立覆盖流形法”，其中，独立覆盖和条形都是基本的计算单元(或称为覆盖网格)，但强调独立覆盖的主体作用。

3　曲梁的独立覆盖函数

如图3所示的整体直角坐标系x-y下，在第i个独立覆盖中，中面坐标(x0i,y0i)由曲线参数方程描述，以(x0i,y0i)为原点建立沿中面曲线变化的局部正交坐标系xi-yi，其中xi和yi分别沿梁的轴线方向(即中面曲线的切向)和厚度方向(即中面曲线的法向)，xi轴与整体坐标x轴的夹角为αi。xi考虑为曲线坐标，可由弧长s或曲线参数方程中的参数来描述，yi表示梁上的点到中面的距离。若中面由多条曲线连接而成，相同曲线的独立覆盖可以合并使用同一条曲线定义的曲线坐标。

图3　整体直角坐标系下的曲梁及局部坐标Fig.3　A curved beam and its local coordinates underthe global rectangular coordinate system

与文献[7]类似，将梁作为二维实体考虑，设局部坐标xi和yi方向的位移分别为ui和vi，参考文献[1]，对应于文献[7]中的Timoshenko直梁，设曲梁各点在局部坐标下的位移为

(1)

梁中面的轴向和厚度方向位移分别设为f1(xi)和f3(xi)，假设xi=s，表明梁中面位移仅仅随弧长s变化，与yi没有关系。根据基本假设②——沿厚度yi方向的应变可忽略，因而对于固定了xi坐标的梁上任一点，有vi(xi)=f3(xi)，再根据基本假设①，用yif2(xi)描述轴向纤维沿厚度方向呈线性变化的伸长量(相对于中面位移f1(xi))，其中f2(xi)代表独立定义的随xi=s变化的截面转角θ。

考虑关于局部坐标的完全多项式，则有

(2)

(a) 轴线方向

(b) 厚度方向

图4多项式排列
Fig.4Polynomialorder

图5　曲梁的独立覆盖及其连接Fig.5　Independent covers and their connections oncurved beams

4　引入局部坐标的曲梁分析公式

如图5所示，以相邻的2个独立覆盖为例。在第i(i=1，2)个独立覆盖内，引入转换矩阵Li(各项为局部坐标轴的方向余弦)，转换到整体坐标系下的位移为

(3)

在独立覆盖1和独立覆盖2的条形重叠区域，比如图5的A点和B点之间，有

(4)

设

Tink=wiLiNink。

(5)

(6)

(7)

上式：

(8)

其中：

(9)

与文献[7]相比，考虑到αi随坐标的变化，因此式(7)多了局部坐标系的方向余弦关于整体坐标的导数。

如图3所示，由简单的几何关系，得到各点的整体直角坐标与局部坐标的转换公式为

(10)

式(10)再结合中面(x0i,y0i)的曲线表达式，就可以将式(9)、式(8)乃至式(7)计算出来。

按最小势能原理推导单元刚度矩阵的子矩阵为

(11)

式中：Ω为计算单元的区域；i=1，2；j=1，2；q，t依次对第i个和第j个多项式覆盖函数中参与计算的所有项进行循环；D为平面应力的弹性矩阵，其表达式为

(12)

式中:E为弹性模量;μ为泊松比。

集中荷载的子向量为

(13)

式中Fx和Fy分别为水平向和竖向的集中荷载。分布荷载也有类似的采用积分形式的公式。

将上述单元刚度矩阵和单元荷载向量分别组集成整体刚度矩阵和整体荷载向量，然后加上边界条件，就可以求解线性方程组得到多项式的系数，过程与文献[7]类似。

5　积分方式及具体计算过程

∬ΩF(x(φ,r),y(φ,r))dxdy=∬ΩF(φ,r)dsdr=

∬ΩF(φ,r)s′(φ)dφdr。

(14)

式中，ds=s′(φ)dφ。

令

(15)

式中：r1(φ)和r2(φ)分别表示沿厚度方向积分的下限和上限(即-h/2到h/2，h为厚度，可以考虑h随φ的变化)。采用一维Gauss数值积分[2]计算,即

(16)

积分点的个数nl与覆盖函数关于xi=φ的阶次有关，Hl为积分权值。对于固定的积分点φl，式(15)的F1(φl)一般也采用一维Gauss积分，由于式(1)中的覆盖函数只出现关于yi=r的常数项和一次式，因此，只用2点Gauss积分就可以满足积分精度要求。更进一步有可能在式(15)中用解析积分的方式消掉r，推导出简化的F1(φ)表达式，这样就只需一次一维Gauss积分。

综上所述，只需在积分点(包括沿厚度方向和沿轴线方向)上计算被积函数值，再乘以相应的积分权值，就可以得到式(11)的刚度子矩阵。在具体计算过程中，关键是式(9)中的局部坐标和式(7)中的方向余弦关于整体坐标的导数计算，以下给出计算步骤：

(1)列出中面曲线的参数方程,即

(17)

(2)根据式(17)，求出积分点处的φ。

(5)求α的正弦和余弦关于x或y的导数,即:

(18)

图6　一段圆形曲梁及其变形(一端固定，另一端受水平力作用)Fig.6　A circular curved beam and its deformation(with one end constrained and another subjected tohorizontal load)

6　圆形曲梁算例

考虑中面曲线的曲率不变，如图6所示，计算第1象限的1/4圆形曲梁。按上节的5个步骤计算如下：

(1)中面曲线参数方程为

式中:(x0,y0)为圆心坐标;R为半径，得到弧长的微分ds=Rdφ。

设圆的半径为1 m，梁的截面厚度为0.01 m。只用1个独立覆盖，左端通过一个窄条与固定端相连，右端承受水平力F=1 kN。梁的弹性模量E=108kN/m2，泊松比μ=0。

轴向采用6点Gauss积分，计算得到荷载作用点处的位移为：覆盖函数取5阶多项式(总共17个自由度)时，u=-0.093 9 m，v=-0.059 8 m；取6阶和7阶多项式(分别为20个和23个自由度)时，u=-0.094 1 m，v=-0.059 8 m。

将文献[14]的理论解用关于φ和r的多项式展开(考虑R=1，s=Rφ=φ)，得

(19)

可见，局部坐标xi=φ，yi=r的位移uφ和ur在中面处随参数方程的参数φ而变化，uφ沿厚度方向呈线性变化(关于r2项的系数很小)，ur沿厚度方向的变化可忽略，这与梁的基本假设一致。在底部φ=0处，理论解的位移与本文计算结果非常接近，不仅如此，计算得出的独立覆盖上关于φ和r多项式的各项系数，也与式(19)比较接近。

如图7(a)所示，只用1个独立覆盖，左端和右端分别通过窄条与固定端相连，在梁的中点作用合力为F=1 kN的法向集中力。梁的弹性模量E=105kN/m2，泊松比μ=0。

轴向采用10点Gauss积分，覆盖函数取10阶和11阶多项式(分别为32个和35个自由度)，计算得到荷载作用点处的位移u=v=-0.107 m，与细密网格的有限元解u=v=-0.112 m的相对误差<5%。但进一步升阶时计算不收敛，有可能是积分精度不够。重新划分为5个独立覆盖，如图7(b)所示，各独立覆盖只需采用4阶多项式(总自由度为70个)，轴向采用6点Gauss积分，得到u=v=-0.110 m。因此，实际计算中不一定要用很高的多项式阶次，在控制好计算误差的前提下，可采用升阶和加密覆盖的2种方式来提高精度。

图7　圆形曲梁及其变形(两端固定，中点受集中力)Fig.7　A circular curved beam and its deformation(with two ends constrained and the middle pointsubjected to concentrated load)

值得关注的是，当采用1个独立覆盖时，独立覆盖的单元弧长与截面厚度之比超过150∶1，但未出现有限元法的Timoshenko梁单元在求解细长梁时的剪切自锁现象。

图8　椭圆形曲梁及其变形(一端固定，另一端受集中力作用)Fig.8　An ellipse curved beam and its deformation(with one end constrained and another subjected toconcentrated load)

7　椭圆形曲梁算例

考虑中面曲线的曲率变化，如图8所示，计算第1象限的1/4椭圆形曲梁。具体计算过程如下：

(5)按式(18)求α的正弦和余弦关于x或y的导数。

设椭圆的半长轴a=2 m，半短轴b=1 m，梁的截面厚度为0.01 m。只用1个独立覆盖，左端通过一个窄条与固定端相连，右端分别承受水平力和竖向力F=1 kN。梁的弹性模量E=108kN/m2，泊松比μ=0。

采用划分100个单元的有限元结果作为对比，本文用10阶多项式(自由度为32)的计算结果如下：水平力作用下，有限元为u=-0.168 5 m，v=-0.223 5 m，本文为u=-0.168 3 m，v=-0.223 0 m，如图7(a)所示；竖向力作用下，有限元为u=-0.223 5 m，v=-0.326 0 m，本文为u=-0.223 0 m，v=-0.325 2 m，如图7(b)所示。可以看出，各截面上的轴力、弯矩、剪力等各种内力荷载组合作用下的变形都能够反映出来。

类似地，关于r的表达式为低阶有理多项式，可用解析积分方式，但公式推导相对于圆形曲梁更为复杂。

8　结　语

本文采用独立覆盖流形法，在前期研究的关于直梁的独立覆盖分析方法基础上，只需借助随中面参数方程变化的局部坐标系，并考虑该坐标系的局部坐标以及方向余弦关于整体坐标的导数，就能实现精确几何描述下的曲梁分析。文中给出常曲率的一段圆形曲梁和变曲率的一段椭圆形曲梁2个算例，验证了方法的可行性。

本文方法的特点是：完全采用实体分析模式，只需使多项式覆盖函数中的某些项不参与计算，就能模拟梁的基本假设，相比较而言，通常的曲梁单元构造比较复杂，需要推导曲梁控制方程和相应的离散公式，考虑曲率对几何方程、平衡方程的影响以及各种内力耦合；同时，本文方法还保持了梁的独立覆盖分析方法的一些基本优点，如仅要求近似函数的C0连续性、无剪切自锁、不存在由于厚度远小于其他尺度而导致的数值病态、与实体连接方便等；在实现几何保形性方面，这是除了等几何分析方法之外的另一途径，与之相比，本文方法具有多项式基函数相对简单、严格施加本质边界条件、完全消除剪切自锁等优点。

正如引言所述，本文不限于讨论二维曲梁，而是希望找到解决曲梁和曲壳类问题的新途径。后续研究将紧紧围绕这一主题，探讨新方法对各种曲梁和曲壳的适应性：参考文献[7]，考虑到截面转角θ=dv/dxi，模拟细长曲梁的“直法线假设”；考虑中面为不同曲线参数方程的连接，并给出多种厚度以及变截面算例；最近完成了母线为曲线描述的旋转壳分析，将另文介绍；三维曲梁和曲壳分析的思路与本文一致，因此将按空间曲线或曲面参数方程建立随中面变化的曲线坐标系，并计算该坐标系的曲线坐标和方向余弦关于整体坐标的导数，最终在三维空间内实现精确几何描述下的曲梁和曲壳分析。

致谢：感谢石根华博士的指导！

参考文献：

[1]ZIENKIEWICZ O C, TAYLOR R L.有限元方法[M].5版.庄 茁，岑 松,译.北京: 清华大学出版社, 2006.

[2]王勖成．有限单元法[M]．北京：清华大学出版社，2003.

[3]潘科琪，刘锦阳. 柔性曲梁多体系统的研究现状和展望[J]. 力学进展，2011，41(6)：711-721.

[4]赵跃宇，康厚军，冯锐，等. 曲线梁研究进展[J]. 力学进展，2006，36(2)：170-186.

[5]HUGHES T J R，COTTRELL J A，BAZILEVS Y. Isogeometric Analysis：CAD，Finite Elements，NURBS，Exact Geometry and Mesh Refinement[J]. Computer Methods in Applied Mechanics and Engineering，2005,194(39/40/41): 4135-4195.

[6]李新康. 层合结构等几何分析研究[D]. 杭州：浙江大学，2015.

[7]苏海东，颉志强. 梁的独立覆盖分析方法[J]. 长江科学院院报, 2018,35(4):143-150.

[8]SHI G H. Manifold Method of Material Analysis[C]∥U.S. Army Research Office. Transactions of the Ninth Army Conference on Applied Mathematics and Computing. Minneapolis, Minnesota, U.S.A, June 18-21,1991: 51-76.

[9]BABUSKA I, MELENK J M. The Partition of Unity Method[J]. International Journal for Numerical Methods in Engineering, 1997, 40:727-758.

[10] 祁勇峰, 苏海东, 崔建华．部分重叠覆盖的数值流形方法初步研究[J].长江科学院院报, 2013,30(1):65-70.

[11] SU Hai-dong, QI Yong-feng, GONG Ya-qi,etal. Preliminary Research of Numerical Manifold Method Based on Partly Overlapping Rectangular Covers[C]∥DDA Commission of International Society for Rock Mechanics. Proceedings of the 11th International Conference on Analysis of Discontinuous Deformation (ICADD11), Fukuoka, Japan, August 27-29, 2013, London: Taylor & Francis Group, 2013: 341-347.

[12] 苏海东, 祁勇峰, 龚亚琦,等. 任意形状覆盖的数值流形方法初步研究[J]. 长江科学院院报, 2013, 30(12): 91-96.

[13] 苏海东，颉志强，龚亚琦,等．基于独立覆盖的流形法收敛性及覆盖网格特性[J]. 长江科学院院报，2016，33(2)：131-136.

[14] 王敏中，王炜，武际可．弹性力学教程(修订版)[M]．北京：北京大学出版社，2011.