混合能量原理求解矩形板在静水压力作用下的弯曲

陈英杰 阙春发 郭敦

(燕山大学建筑工程与力学学院,河北秦皇岛066004)

矩形板作为结构基本构件之一，我们可以把矩形板理论分析方法归纳分为精确解法和近似解法，且这两类方法均以能量原理为基础。能量法又称为变分法，可以追溯到17世纪末，有着悠久的历史。变分法可以用于处理函数变量的问题。截止到目前，国内外学者普遍非常重视变分原理[1-2]的研究，而且已经有了丰硕的成果。例如：王根会等[3]基于能量变分法推导了新型组合箱梁的总势能式与微分控制方程，为新型组合箱梁的应用推广提供了理论和技术支持；曾祥勇等[4]以能量变分法的最小势能原理为基础分析了矩形筏板的受弯问题。混合变量的极值变分原理[5-6]包含混合变量的最小势能原理[7-9]与混合变量的最小余能原理[10]。混合变量的极值变分原理，其容许位移与容许内力分别被弱容许位移与弱容许内力所代替，与传统理论[11-12]相比增强了等价方程，而且需要预先满足的条件减弱了。

本文将应用变分原理中的混合变量的最小势能原理讨论矩形板在静水压力作用下的弯曲问题，并最终得到矩形板受静水压力的数值精确解。应用Matlab软件计算出矩形板挠度的精确值，应用ANSYS模拟软件得到矩形板挠度模拟值，并对两项数值进行对比分析，说明本文方法具有更好的准确性。

1 三边简支一边固定矩形板

1.1 挠曲线方程

现考虑受静水压力的矩形板，边界条件如图1(a)所示。图中q0为三角形分布的均布载荷，b为矩形板沿y轴方向的长度，a为矩形板沿x轴方向的长度，w为挠度。以分布弯矩My0代替固定边的弯曲约束，创建如图1(b)所示的变形后的等效图。

图1 受静水压力的三边简支一边固定的矩形板

分布载荷

假设固定端弯矩

式中Cm为待定系数，αm=mπ/a，同时假设板的弱容许挠度为

式中，βn=nπ/b。

图1(b)对应的矩形板混合变量总势能

式中D=Eh3/[12(1−ν2)]，为板的抗弯刚度，E为弹性模量，ν为泊松比。

将式(1)和式(2)代入式(4)经过积分运算并且对Amn取变分极值得

将式(5)代入式(3)得

在推导本边界条件挠曲线方程的过程中，式(6)是用混合变量最小势能原理导出的以正弦双重三角级数表示的弱容许挠度。而对于非齐次挠度和弯矩边界条件，正弦双重三角级数将会在边界上出现第二类间断点。为避免出现这种情况，并且能够加快级数收敛速度，需要将式(6)转换成为在边界上连续可微的挠度，而该挠度即为本次边界条件下的挠度表达式。应用文献[5]中附录式(A92)和式(A47)分别对式(6)各项进行转化，则得

1.2 边界条件

固定端转角为0，故式(7)应符合式(8)的边界条件，其他边界条件已自动满足

将方程式(7)代入式(8)中，经过运算整理可得执行方程

1.3 数值计算

取矩形板的各项参数：a=b=1 m，泊松比ν=0.3,弹性模量E=2.0×1011Pa，板的厚度h=0.01 m，静水压力q0=1 MPa。其中含有未知项Cm。其次，应用Matlab软件编程计算，可求得具体挠度值。对计算程序循环50次，便可保障计算结果的收敛性。

同时，应用模拟软件对本次边界条件下的矩形板进行建模分析，作为参考解与本文解进行对比。此种边界条件选用模型shell63单元，该单元具有弯曲能力和薄膜效应，并且忽略剪切变形，可以很好地应用于矩形板变形过程中的数值模拟。弯曲矩形板在x/a=0.1，x/a=0.3，x/a=0.5和x/a=0.8与y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8线上的模拟数值与本文计算数值两项数据绘制成图和表，如图2与表1所示。

图2 x/a=0.1,x/a=0.3,x/a=0.5和x/a=0.8处挠度分布曲线图

1.4 结果分析

图2为矩形板沿x轴不同位置处的本文计算值和模拟值的挠度对比曲线图，可以直观看到,由于结构受三角形载荷作用，矩形板在静水压力作用下的挠度并非对称分布变化，且符合实际受力变化规律，计算结果更为贴近实际，这也间接说明混合变量的最小势能原理解决矩形板受静水压力问题的适用性。

表1给出了矩形板在y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8位置处沿x轴方向变化的挠度计算值和模拟参考值。通过对两项数值的对比分析，可以得到两者挠度的最大相对差值分别为：4.0%，3.8%，3.7%，本文研究方法计算的挠度值略大于有限元模拟的参考解，表明本方法可以更好地保证结构安全，反映结构在载荷作用下的变形规律，对矩形板弯曲解拥有更好的有效性和适用性，可以更精准地计算矩形板的挠度问题。

表1 y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8处挠度沿x分布值(10−7 m)

2 两邻边固定两邻边简支矩形板

2.1 挠曲线方程

现考虑受静水压力的矩形板，边界条件如图3(a)所示。以分布弯矩Mx0，My0代替固定边的弯曲约束，创建如图3(b)所示的变形后的等效图。

图3 受静水压力的两邻边固定两邻边简支矩形板

分布载荷

假设固定端弯矩

式中，An为待定系数。同时假设板的弱容许挠度为

可得图3(b)对应的矩形板混合变量总势能式

将式(10)∼式(12)代入式(14)经过积分运算并且对Amn取变分极值得

将式(15)代入式(13)得

同理，需要将式(16)转换成双曲函数和三角级数混合表示的基本解。应用文献[5]中附录式(A47)和式(A92)分别对式(16)各项进行转化，得

2.2 边界条件

固定端应满足转角为零，故式(17)应符合式(18)和式(19)的边界条件，其他边界条件已自动满足

将方程式(17)代入式(18)和式(19)中，运算整理得执行方程

2.3 数值计算

通过之前的公式推导，已经得到本边界条件的执行方程。且矩形板的计算参数与1.3相同。本节中所需要求解的未知数是An和Cm，对计算程序循环50次，运行结果不再变化，得到本边界条件下的挠度值。应用模拟软件对本次边界条件下的矩形板进行建模分析，单元类型采用Shell63，将弯曲矩形板在x/a=0.1，x/a=0.3，x/a=0.5和x/a=0.8与y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8线上的模拟数值与本文计算数值分别绘制成图和表，如图4与表2所示。

图4 x/a=0.1，x/a=0.3，x/a=0.5和x/a=0.8处挠度分布曲线图

2.4 结果分析

图4给出了矩形板沿x轴不同位置处的本文计算值和模拟值的挠度对比曲线图，可以直观看出，本文研究方法计算的挠度值略大于有限元模拟的参考解，表明本方法可以更好地保证结构安全，反映结构在载荷作用下的变形规律，对矩形板弯曲解具有更好的有效性和适用性，可以更精准地计算矩形板的挠度问题。由于受到载荷分布位置与边界条件不同的影响，矩形板的挠度曲线均沿坐标轴先增大后减小，且在x/a=0.6的位置达到了最大值。表明计算出的矩形板在静水压力作用下的挠度变化是符合受力变化规律的。

表2给出了矩形板在y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8位置处沿x轴方向变化的挠度公式计算值和模拟值。通过对两项数值的对比分析，可以得到两者挠度的最大相对差值分别为：3.5%、3.4%、3.0%，本文研究方法计算的挠度值略大于有限元模拟的参考解，表明本方法可以更好地保证结构安全，反映结构在载荷作用下的变形规律，使矩形板挠度问题的计算更为精确。

表2 矩形板在y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8处挠度沿x分布值(10−7 m)

3 三边固定一边简支矩形板

3.1 挠曲线方程

现考虑受静水压力的矩形板，边界条件如图5(a)所示。以分布弯矩Mx0，My0，Mxa代替固定边的弯曲约束，创建如图5(b)所示的变形后的等效图。

分布载荷

假设固定端弯矩

假设板的弱容许挠度为

可得图5(b)对应的矩形板混合变量总势能式

图5 受静水压力的三边固定一边简支弯曲矩形板

将式(22)∼式(24)代入式(26)经过积分运算并且对Amn取变分极值得

将式(27)代入式(25)得

同理，需要将式(28)转换成双曲函数和三角级数混合表示的基本解。应用文献[5]中附录式(A47)和式(A92)分别对式(28)各项进行转化，可得

3.2 边界条件

固定端应满足转角为零，故式(29)应符合式(30)和式(31)的边界条件，其他边界条件已自动满足

把方程式(29)代入式(30)和式(31)中，经过运算整理得执行方程。

3.3 数值计算

通过之前的公式推导，已经得到本边界条件的执行方程。且矩形板的计算参数与1.3相同。本节中所需要求解的未知数是An和Cm，同理为了保障计算数值的收敛性，对计算程序循环50次，得到本边界条件下的挠度值，并通过模拟软件模拟分析，将模拟结果作为参考值。将弯曲矩形板在x/a=0.1，x/a=0.3，x/a=0.5和x/a=0.8与y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8线上的模拟数值与本文计算数值两项数据分别绘制成图和表，如图6与表3所示。

3.4 结果分析

图6给出了矩形板沿y轴不同位置处的本文计算值和模拟值的挠度对比曲线图，可以直观看出，由于受到载荷分布位置与边界条件不同的影响，矩形板的挠度曲线都沿坐标轴先增大后减小，且在x/a=0.6的位置达到了最大值。表明计算出的矩形板在静水压力作用下的挠度变化符合受力变化规律。

图6 x/a=0.1，x/a=0.3，x/a=0.5和x/a=0.8处挠度分布曲线图

表3给出了矩形板在y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8位置处沿x轴方向变化的挠度公式计算值和模拟值。通过对两项数值的对比分析，可以得到两者挠度的最大相对差值分别为：4.2%、3.9%、3.2%，本文研究方法计算的挠度值略大于有限元模拟的参考解，表明本方法可以更好地保证结构安全，并且对求解矩形板弯曲解拥有更好的有效性和适用性。

表3 y/b=0.1，y/b=0.3，y/b=0.5和y/b=0.8处挠度沿x分布值(10−7 m)

4 结论

本文依据混合变量的最小势能原理推导出了三边简支一边固定、两邻边固定两邻边简支、三边固定一边简支三种不同边界条件的弯曲矩形板在静水压力作用下的边界应力函数的表达式和封闭解析解，数值计算软件求得的数值解与模拟软件的模拟值进行归纳分析，表明本文方法的正确性，说明了本文研究方法可以使矩形板的弯曲问题的求解简单化和精确化。