求解弹性地基上自由矩形中厚板弯曲问题的辛-叠加方法*

李 锐, 田 宇, 郑新然, 王 博

(大连理工大学 工程力学系; 国际计算力学研究中心; 工业装备结构分析国家重点实验室(大连理工大学), 辽宁 大连 116024)

(第五届“钱令希计算力学奖(青年奖)”特邀论文)

引 言

弹性地基上的矩形中厚板是土木工程、机械工程、航空航天及海洋工程中的一种常见的结构形式,主要用于刚性路面、桥面板、舰船甲板等承重结构由于Kirchhoff薄板理论忽略了板在弯曲时的横向剪切变形,其应用于实际工程中较厚的板时往往会产生较大误差于是学者们相继建立了若干中厚板理论,其中Mindlin中厚板理论是应用最为广泛的理论之一[1],自其建立以来,不断有学者致力于求解相应的弯曲问题然而,由于中厚板的偏微分控制方程阶次较高,其边值问题的解析求解成为一类难题:除了半逆法容易处理的对边简支情况以外,其余边界条件下矩形中厚板的弯曲问题都不易获得解析解,其中尤以含自由边的板难以求解上述状况引得学者们采用多种数值方法求解相关问题,例如Rayleigh-Ritz法[2]、有限差分法[3]、有限元法[4]、样条元法[5]、边界元法[6]、微分求积法[7]、微分容积法[8]、离散奇异卷积法[9]、小波配点法[10]、状态空间法[11]

各类有效的数值方法无疑可以获得工程上可接受的中厚板问题的解,但这并未动摇解析解的地位首先,由于能够精确反映各参量之间的关联,解析解在力学研究中具有重要价值例如,解析解可以作为检验各类近似/数值方法的基准,能为快速参数分析和优化提供有力的工具,也是高效指导复杂实验设计的重要理论基础此外,有限元等数值方法本身就是一种近似解,其在求解板壳问题时也存在着一些问题,例如边界条件不能总是精确满足、应力等位移的高阶导数结果不够精确等可以说,解析解的获得不仅对板壳力学求解理论的发展和完善具有重要价值,而且对工程结构分析与设计具有重要的指导意义

近年来,基于钟万勰院士、姚伟岸教授等开创的辛弹性力学求解体系[12],李锐等提出了一种板壳力学求解新方法:辛-叠加方法(symplectic superposition method),获得了若干矩形薄板、中厚板弯曲、振动和稳定性问题的新解析解[13-21]辛-叠加方法的基本思路是:将待求问题导入Hamilton体系,形成Hamilton对偶方程边值问题;将原问题拆分为若干子问题,对这些子问题,运用辛数学方法构建相应的本征值问题,解析地获得本征解,进而利用辛本征展开等手段,精确表征子问题中板的挠度(弯曲问题)、振型(振动问题)或屈曲模态(稳定性问题)等控制变量;通过叠加子问题的解,获得原问题的解通过该方法,对于静力问题,将直接得到结构的位移场解析解;对于振动问题,将解析获得结构的固有频率和相应的振型;对于稳定性问题,将解析得到结构的临界屈曲载荷和相应的屈曲模态辛-叠加方法兼备了辛数学方法无需事先假定解形式的优点以及叠加法程式化的优点,具有统一的数学结构,不受问题具体类型的限制该方法一方面克服了其他传统方法求解过程中遇到的分离变量、半逆法失效等难题,同时规避了运用传统辛数学方法所出现的本征值方程无法解析求解等问题,在解析求解板壳力学问题中展现出独特的优势

本文首次应用辛-叠加方法解析求解了基于Mindlin理论的弹性地基上四边自由矩形中厚板的弯曲问题由于边界约束最弱,因此该类问题实际上是各类边界条件下矩形板问题中最难求解的情况之一最终的数值算例表明,无论是所关心的位移还是内力解,本文结果均与精细有限元分析结果吻合得很好,从而证明了该方法对于求解中厚板问题的适用性

1 控制方程的Hamilton体系表达

xOy直角坐标系下,弹性地基上中厚板的平衡方程为

(1)

(2)

(3)

其中,q为横向分布外载荷;K为地基刚度;板内单位长度的内力包括弯矩MxMy,扭矩Mxy,以及剪力QxQy上述各内力可以由广义位移(挠度W、转角ψxψy)表达为

(4)

(5)

(6)

(7)

(8)

其中,板的抗弯刚度D和剪切刚度C与弹性模量E、厚度h以及Poisson比ν存在以下关系:

D=Eh3/[12(1-ν2)],C=5Eh/[12(1+ν)],

剪切修正系数取为5/6

由式(2)、(3)、(7)和(8),可得

(9)

(10)

(11)

(12)

可得

(13)

(14)

将式(13)和(14)代入式(1),可得

(15)

将式(7)代入式(13),可得

(16)

将式(8)代入式(14),可得

(17)

对式(16)两边关于x求导,得

(18)

对式(17)两边关于y求导,得

(19)

利用式(11)和(15),对式(18)和(19)求和,可给出

C(2W-M)=KW-q

(20)

对式(16)两边关于y求导,得

(21)

对式(17)两边关于x求导,得

(22)

利用式(12),将式(21)和(22)相减,可给出

(23)

综上,求解弹性地基上中厚板的弯曲问题可归结为求解关于控制函数MWΨ的方程组(15)、(20)和(23):所有关键力学量均可经式(16)、(17)连同式(4)~(8),由上述3个函数表达

为将所求问题导入Hamilton体系,令

(24)

(25)

(26)

由式(23),注意到2C/[D(1-ν)]=10/h2,可得

(27)

由式(15)和(20),可得

(28)

(29)

式(24)~(29)可以写成如下矩阵形式:

(30)

其中

矩阵H满足HT=JHJ,其中是单位辛矩阵,I3是三阶单位矩阵,因此H是一个Hamilton算子矩阵[12],而式(30)即为弹性地基上中厚板弯曲问题的Hamilton对偶方程容易看出,该对偶方程形式非常简洁,因此对后续求解非常有利

需要指出,式(30)仅是所求问题控制方程的一种Hamilton体系表达,采用其他推导方法(例如基于Hellinger-Reissner变分原理)也可以导出其他形式的Hamilton对偶方程当然,此时的状态向量Z所包含的元素也将有所不同

2 弹性地基上对边滑支矩形中厚板的辛解析解

本节利用辛数学方法求解弹性地基上对边滑支矩形中厚板的弯曲问题,从而为子问题的求解奠定基础

方程(30)的齐次方程为

(31)

在辛数学求解框架下,分离变量法有效[12],因此令

Z=X(x)Y(y),

(32)

其中

X(x)=[M(x),W(x),Ψ(x),β(x),α(x),θ(x)]T

将式(32)代入式(31),得

HX(x)=μX(x),

(33)

(34)

其中μ为本征值,X(x)为本征向量

与式(33)对应的特征方程为

(35)

其中k=-K/Dδ=D/C展开式(35),得

(36)

方程(36)的根为

(37)

其中为i为虚数单位,

(38)

因此可以写出方程通解:

(39)

将其代入式(33),得到各常数之间的关系:

(40)

其中

(41)

对于x=0和x=a边滑支的矩形中厚板,其x方向的边界条件为

(42)

将式(39)和(40)代入式(42),令所得方程组的系数矩阵行列式为0,即得到弹性地基上对边滑支矩形中厚板弯曲问题的本征值超越方程:

sin(α1a)sinh(α2a)sinh(α3a)=0

(43)

求解得到两组本征值,即

(44)

(45)

其中αm=mπ/a(m=1,2,3,…)

本征值μ1μ4对应的本征向量分别为

(46)

本征值μm1μm6对应的本征向量分别为

(47)

本征向量组(46)及(47)满足共轭辛正交关系[12]

非齐次方程(30)的解可写成

Z=X(x)Y(y),

(48)

其中

X(x)=[X1(x),X2(x),X3(x),X4(x),…,Xm1(x),Xm2(x),

Xm3(x),Xm4(x),Xm5(x),Xm6(x),…],

Y(y)=[Y1(y),Y2(y),Y3(y),Y4(y),…,Ym1(y),Ym2(y),

Ym3(y),Ym4(y),Ym5(y),Ym6(y),…]T

将式(48)代入式(30),得到

X(x)dY(y)/dy=HX(x)Y(y)+f

(49)

注意到

HX(x)=X(x)M

(50)

其中M=diag(P0,…,Qm,…),P0=diag(μ1,μ2,μ3,μ4),Qm=diag(μm1,μm2,μm3,μm4,μm5,μm6)(m=1,2,3,…)f可以按辛本征向量展开为

f=X(x)G,

(51)

其中

G=[g1,g2,g3,g4,…,gm1,gm2,gm3,gm4,gm5,gm6,…]T

为展开系数列阵,各元素可通过对式(51)两边同时左乘X(x)TJdx并关于x从0到a积分,由共轭辛正交关系求出于是式(49)给出

dY(y)/dy-MY(y)=G

(52)

对于在(x0,y0)处作用一集中载荷p的板,可得

(53)

其中H(y-y0)为单位阶跃函数,c1c4以及cm1cm6为待求常数,由板在y方向的边界条件决定

3 弹性地基上四边自由矩形中厚板的辛-叠加解析解

如图1(a)所示,弹性地基上四边自由矩形中厚板在(x0,y0)处作用一集中载荷p该问题可由图1(b)、 (c)、 (d)这3类子问题叠加求解: 子问题1为集中载荷作用下的四边滑支板(图1(b));子问题2为四边滑支板在y=0和y=b边分别承受级数表示的转角作用(图1(c)); 子问题3为四边滑支板在x=0和x=a边分别承受级数表示的转角作用(图1(d)),其中βm=mπ/b

对于子问题1,y方向的边界条件为

ψy|y=0,b=0,Qy|y=0,b=0,Mxy|y=0,b=0

(54)

将式(46)、(47)及(53)代入式(48),求出ψy,Qy,Mxy的表达式,然后代入式(54),求出常数c1c4以及cm1cm6,最终得到集中载荷作用下弹性地基上四边滑支矩形中厚板弯曲问题的辛解析解:

(55)

其中

对于子问题2,代入相应的边界条件,同理可得问题的解:

(56)

其中

对于子问题3,其解为

(57)

其中

通过以上推导给出了3个子问题的挠度解,其余力学量容易类似求出

为满足x=0边自由的边界条件,叠加以上3个子问题在x=0边的弯矩,并令其等于零,化简后得到第一组方程:

(58)

[Em-Fmcos(iπ)]+

(59)

其中i=1,2,3,…

为满足x=a边自由的边界条件,叠加子问题在x=a边的弯矩,并令其等于零,得到第二组方程:

(60)

ϑβi[E0-cos(iπ)F0]+

[Em-Fmcos(iπ)]+

(61)

其中i=1,2,3,…

为满足y=0边自由的边界条件,叠加子问题在y=0边的弯矩,并令其等于零,得到第三组方程:

(62)

[G