1. 引言

随着高强材料与薄壁结构在土木工程中的应用，板壳在面内压力作用下的屈曲破坏作为一个重要的力学问题被广泛研究。王忠民等 [1] 用Levy法和有限差分法研究了两对边简支、另两对边为三种典型支撑任意组合的六种弹性矩形薄板在切向均布随从力作用下的动力稳定性问题。杨端生等 [2] 求得了正交异性矩形薄板屈曲位移函数微分方程的一般解，可用以求解任意边界矩形板的稳定问题。甘立飞等 [3] 对矩形薄板受面内非均匀分布载荷的稳定性问题进行了有限元分析，提出即使在传统意义的薄板范围内，当相对厚度较大时，厚度对临界荷载的影响有时已不容忽略。Eisenberger和Alexandrov [4] 给出了受到双轴压缩且厚度在平行于两侧的方向上变化的矩形薄板分叉屈曲载荷的精确解，采用扩展的Kantorovich方法分析了具有可变抗弯刚度的受压构件的稳定性。Shufrin和Eisenberger [5] 使用多项扩展的Kantorovich方法计算了四边受剪切荷载的弹性矩形薄板的屈曲荷载。Shufrin等 [6] 提出了矩形板弹性非线性稳定性分析的半解析方法，考虑了任意边界条件和一般平面外和平面内荷载，利用薄板理论、非线性von Kármán应变和变分多项扩展Kantorovich方法，导出了弹性矩形板的几何非线性公式。Ruocco等 [7] 提出了一种在Kirchhoff-Love假设下研究任意边界条件下矩形薄板弹性稳定性的有限元分析方法，该方法可以得到不依赖于单元数目的屈曲载荷和模态位移。Ge等 [8] 基于弹性理论和Galerkin方法，建立了板内参数化窄带噪声激励下矩形薄板的动力学模型。

能量法是求解矩形薄板屈曲问题的一种经典方法，其使用三角级数展开式逼近薄板挠度，求解结果精度随着三角函数取值项数的增加而增加。龙志勤等 [9] 采用能量法和有限元法对简支矩形薄板的弹性失稳问题进行了求解并进行了非线性分析研究，结果表明，特征值屈曲的临界荷载是薄板实际压屈的上限，非线性分析的临界荷载低于特征值屈曲的结果。张永存等 [10] 分析了负泊松比材料对矩形平板在单轴压缩屈曲行为中的增益机理，考虑了通常忽略的非加载边弹性面内平移约束和扭转约束，利用能量法给出了单轴压缩矩形平板屈曲的显示解析表达式。Larionov等 [11] 采用能量法估算了纵向力作用下钢筋混凝土构件的稳定性。Popa等 [12] 分析了径向肋圆板在中心集中力F = 77,000 N作用下产生的变形，确定了解析计算得到的挠度，并与有限元法得到的位移进行了比较。Huang和Peng [13] 引入无惩罚项的改进板深能量法，研究了不规则Kirchhoff板的弯曲、自由振动和屈曲问题。

本文采用能量法计算了纵向集中荷载下四边简支矩形薄板的屈曲临界荷载，并利用有限元软件ABAQUS对简支薄板的线性及非线性压屈问题进行了计算，将计算结果与能量法的求解结果对比，验证了分析结果的正确性。

2. 矩形薄板屈曲临界荷载理论解

当薄板在特定纵向荷载的分布下处于平面平衡状态时，我们可以通过以下方式来判断该状态是否稳定：如果薄板在受到横向干扰力后进入附近的某一弯曲状态，在干扰力消失后，它是否能恢复到原来的平面状态。为了判断这一点，我们需要观察薄板从平面状态进入弯曲状态时势能的变化。具体而言，我们需要考虑以下情况：当薄板从平面状态进入弯曲状态时，势能是增加还是减少。如果势能增加，表示平面状态的势能是极小的，对应于稳定平衡；如果势能减少，表示平面状态的势能是极大的，对应于不稳定平衡；如果势能保持不变，表示平面状态的平衡是稳定平衡的极限，而相应于这一极限状态的纵向荷载就是临界荷载 [14] 。势能保持不变的原因在于荷载势能的减少恰好等于形变势能的增加，而荷载所做的功也等于荷载势能的减少。因此，从能量的角度来看，可以通过以下条件求解临界荷载：当薄板从平面状态进入附近的弯曲状态时，纵向荷载所做的功W等于形变势能的增加。

当薄板从平面状态进入弯曲状态时，和它受横向荷载作用而弯曲时一样，挠度w是从零开始的，所以形变势能的增加也就是薄板全部的弯曲形变势能V_ε。于是，我们有功能方程W = V_ε，即：

$V_{\epsilon }-W=0$ (1)

而其中的弯曲形变势能V_ε如公式(2)或(3)所示。

在等厚度的薄板中，薄板的抗弯刚度D是常量。以矩形薄板一角点为原点O，薄板平面上角点相邻两边作x轴、y轴，垂直xOy平面作z轴，建立空间直角坐标系，有：

$V_{\epsilon }=\frac{D}{2}{\displaystyle \iint \left\{{\left({\nabla }^{2}w\right)}^2-2\left(1-\mu \right)\left[\frac{{\partial }^{2}w}{\partial x^2}\frac{{\partial }^{2}w}{\partial y^2}-{\left(\frac{{\partial }^{2}w}{\partial x\partial y}\right)}^2\right]\right\}\text{d}x\text{d}y}$ (2)

其中， $D=\frac{E{\delta }^3}{12\left(1-{\mu }^2\right)}$ ，E为板的弹性模量， $\delta$ 为板厚， $\mu$ 为泊松比。

如果薄板的全部边界都是夹支边，则不论边界的形状如何，在边界上都有 $\frac{\partial w}{\partial x}=0$ 。于是：

$V_{\epsilon }=\frac{D}{2}{\displaystyle \iint {\left({\nabla }^{2}w\right)}^2\text{d}x\text{d}y}$ (3)

如果一个矩形薄板没有自由边，而只有夹支边和简支边，则在x为常量的边界上有 $\text{d}x=0$ 及 $\frac{{\partial }^{2}w}{\partial y^2}=0$ ，在y为常量的边界上有 $\text{d}y=0$ 及 $\frac{\partial w}{\partial x}=0$ 。式(3)仍然成立。

纵向荷载所做的功W就等于纵向荷载乘以作用边上力的作用点距离的缩短，具体计算临界荷载时，可先设定薄板压曲以后的、满足位移边界条件的挠度表达式，然后按照这些表达式，用公式(2)或(3)求出 $V_{\epsilon }$ ，并求出纵向荷载所做的功W。最后命 $V_{\epsilon }=W$ ，即得出用纵向荷载表示的压曲条件，从而求得薄板的临界荷载。

为了使得设定的挠度较好地符合临界荷载下的挠度，从而求得较精确的临界荷载，可以设定挠度的表达式为：

$w={\displaystyle {\sum }_m{C}_m{w}_m}$ (4)

其中， $w_m$ 是满足位移边界条件的函数， $C_m$ 是互不依赖的待定系数。选择 $C_m$ 时，可以应用最小势能原理。以薄板在平面状态下的形变势能及荷载势能均为零，则薄板在压曲状态下的形变势能为 $V_{\epsilon }$ ，荷载势能为 $V=-W$ ，而总势能为 $V_{\epsilon }+V$ ，也就是 $V_{\epsilon }-W$ 。于是，由最小势能原理得到：

$\frac{\partial }{\partial C_m}\left(V_{\epsilon }-W\right)=0$ (5)

这将给出 $C_m$ 的m个齐次线性方程。为了w具有非零解，就必须 $C_m$ 具有非零解，因而这个齐次线性方程组的系数行列式必须等于零。这样就得出求解临界荷载的方程。

3. 实例分析

3.1. 能量法求临界荷载

设有四边简支的矩形薄板，长边长度a = 60 mm，短边长度b = 40 mm，厚度为t = 1 mm，采用Q345钢材( $E=2.06\times {10}^5\text{MPa}$ ， $\mu =0.3$ )，在薄板两长边三分点处施加F = 1000 N的集中荷载，如图1所示。

1) 挠度应满足位移边界条件，取压曲以后的挠度表达式为：

$w=A\sin \frac{\pi x}{a}\sin \frac{\pi y}{b}$ (6)

将挠度w代入公式(3)求得：

$V_{\epsilon }=\frac{{\pi }^{4}abD}{8}{A}^2{\left(\frac{1}{a^2}+\frac{1}{b^2}\right)}^2$ (7)

直接计算纵向荷载在薄板压曲时所做的功，这个功就等于力F乘以两对边上力的作用点距离的缩短，即：

$\begin{array}{c}W=F{\displaystyle {\int }_0^b\frac{1}{2}}{\left(\frac{\partial w}{\partial y}\right)}_{x=\frac{a}{3}}^2\text{d}y+F{\displaystyle {\int }_0^b\frac{1}{2}}{\left(\frac{\partial w}{\partial y}\right)}_{x=\frac{2a}{3}}^2\text{d}y\\ =\frac{F}{2}{\displaystyle {\int }_0^b\frac{{\pi }^2}{b^2}}{\left(A\sin \frac{\pi x}{a}\cos \frac{\pi y}{b}\right)}_{x=\frac{a}{3}}^2\text{d}y+\frac{F}{2}{\displaystyle {\int }_0^b\frac{{\pi }^2}{b^2}}{\left(A\sin \frac{\pi x}{a}\cos \frac{\pi y}{b}\right)}_{x=\frac{2a}{3}}^2\text{d}y\\ =\frac{3{\pi }^{2}F}{4b^2}{A}^2{\displaystyle {\int }_0^b{\cos }^2\frac{\pi y}{b}\text{d}y}=\frac{3{\pi }^{2}F}{8b}{A}^2\end{array}$ (8)

于是，由 $\frac{\partial }{\partial A}\left(V_{\epsilon }-W\right)=0$ 得到：

$\frac{{\pi }^{4}abD}{8}2A{\left(\frac{1}{a^2}+\frac{1}{b^2}\right)}^2-\frac{3{\pi }^{2}F}{8b}2A=0$ (9)

命这一方程中A的系数等于零，即得该板的临界荷载：

$F=F_C=\frac{{\pi }^{2}D}{3a}{\left(\frac{b}{a}+\frac{a}{b}\right)}^2$ (10)

可得：

$F_C=\frac{{3.14}^2}{3\times 0.06}\times {\left(\frac{0.04}{0.06}+\frac{0.06}{0.04}\right)}^2\times \frac{2.06\times {10}^{11}\times {0.001}^3}{12\times \left(1-{0.3}^2\right)}\text{N}=4850.8\text{N}$ 。

2) 若取压曲以后的挠度表达式为：

$w=A\sin \frac{\pi x}{a}\sin \frac{\pi y}{b}+B\sin \frac{5\pi x}{a}\sin \frac{\pi y}{b}$ (11)

将挠度w代入公式(3)求得：

$V_{\epsilon }=\frac{{\pi }^{4}abD}{8}\left[A^2{\left(\frac{1}{a^2}+\frac{1}{b^2}\right)}^2+B^2{\left(\frac{25}{a^2}+\frac{1}{b^2}\right)}^2\right]$ (12)

同样直接计算纵向荷载在薄板压曲时所做的功，这个功就等于力F乘以两对边上力的作用点距离的缩短，即：

$\begin{array}{c}W=F{\displaystyle {\int }_0^b\frac{1}{2}}{\left(\frac{\partial w}{\partial y}\right)}_{x=\frac{a}{3}}^2\text{d}y+F{\displaystyle {\int }_0^b\frac{1}{2}}{\left(\frac{\partial w}{\partial y}\right)}_{x=\frac{2a}{3}}^2\text{d}y\\ =\frac{F}{2}{\displaystyle {\int }_0^b\frac{{\pi }^2}{b^2}}{\left(A\sin \frac{\pi x}{a}\cos \frac{\pi y}{b}+B\sin \frac{5\pi x}{a}\cos \frac{\pi y}{b}\right)}_{x=\frac{a}{3}}^2\text{d}y\\ +\frac{F}{2}{\displaystyle {\int }_0^b\frac{{\pi }^2}{b^2}}{\left(A\sin \frac{\pi x}{a}\cos \frac{\pi y}{b}+B\sin \frac{5\pi x}{a}\cos \frac{\pi y}{b}\right)}_{x=\frac{2a}{3}}^2\text{d}y\\ =\frac{3{\pi }^{2}F}{8b^2}{\left(A-B\right)}^2{\displaystyle {\int }_0^b{\cos }^2\frac{\pi y}{b}\text{d}y}=\frac{3{\pi }^{2}F}{8b}{\left(A-B\right)}^2\end{array}$ (13)

于是，由 $\frac{\partial }{\partial A}\left(V_{\epsilon }-W\right)=0$ 及 $\frac{\partial }{\partial B}\left(V_{\epsilon }-W\right)=0$ 得到：

$\{\begin{array}{l}\frac{{\pi }^{4}abD}{4}{\left(\frac{1}{a^2}+\frac{1}{b^2}\right)}^{2}A-\frac{3{\pi }^{2}F}{4b}\left(A-B\right)=0\\ \frac{{\pi }^{4}abD}{4}{\left(\frac{25}{a^2}+\frac{1}{b^2}\right)}^{2}B+\frac{3{\pi }^{2}F}{4b}\left(A-B\right)=0\end{array}$ (14)

简化以后得：

$\{\begin{array}{l}\left[{\pi }^{2}a{b}^{2}D{\left(\frac{1}{a^2}+\frac{1}{b^2}\right)}^2-3F\right]A+3FB=0\\ 3FA+\left[{\pi }^{2}a{b}^{2}D{\left(\frac{25}{a^2}+\frac{1}{b^2}\right)}^2-3F\right]B=0\end{array}$ (15)

命上列方程组的系数行列式等于零，得到：

$\left|\begin{array}{cc}{\pi }^{2}a{b}^{2}D{\left(\frac{1}{a^2}+\frac{1}{b^2}\right)}^2-3F& 3F\\ 3F& {\pi }^{2}a{b}^{2}D{\left(\frac{25}{a^2}+\frac{1}{b^2}\right)}^2-3F\end{array}\right|=0$ (16)

从而得到该板的临界荷载：

$F=F_C=\frac{{\pi }^{2}D}{3a}\frac{{\left(\frac{b}{a}+\frac{a}{b}\right)}^2{\left(\frac{25b}{a}+\frac{a}{b}\right)}^2}{{\left(\frac{b}{a}+\frac{a}{b}\right)}^2+{\left(\frac{25b}{a}+\frac{a}{b}\right)}^2}$ (17)

可得：

$F_C=\frac{{3.14}^2}{3\times 0.06}\times \frac{{\left(\frac{0.04}{0.06}+\frac{0.06}{0.04}\right)}^2{\left(\frac{25\times 0.04}{0.06}+\frac{0.06}{0.04}\right)}^2}{{\left(\frac{0.04}{0.06}+\frac{0.06}{0.04}\right)}^2+{\left(\frac{25\times 0.04}{0.06}+\frac{0.06}{0.04}\right)}^2}\times \frac{2.06\times {10}^{11}\times {0.001}^3}{12\times \left(1-{0.3}^2\right)}\text{N}=4782.8\text{N}$ 。

比较式(10)和式(17)可知，由两项三角级数挠度表达式导出的四边简支矩形薄板的临界荷载公式比由一项三角级数挠度表达式导出的四边简支矩形薄板的临界荷载公式在形式上更复杂，从理论上讲也更精确。这是因为在两项三角级数挠度表达式中，第二项的函数周期小于第一项的函数周期，与第一项叠加后能够反映更小范围内薄板挠度的变化，所以含有更小周期三角函数的挠度表达式能够更好地逼近薄板屈曲时真实的挠度曲线函数表达式，通过同样的理论推导方式得出的临界荷载也就越接近薄板屈曲时的真实临界荷载。通过两项三角级数挠度表达式推导出的临界荷载为4782.8 N，小于通过一项三角函数挠度表达式推导出的临界荷载(4850.8 N)，两者差值小于较小临界荷载的1.5%。随着三角级数挠度表达式取项的增加，推测临界荷载将收敛于4782.8 N附近某值。

在进行两项三角级数挠度表达式中第二项三角函数的选取时，并不是只要满足周期小于第一项三角函数这一条件，其还需根据荷载作用的位置进行选取。如本例子中F的作用位置为长边，由于长边平行于x轴，那么就需要在包含x的三角函数内添加系数；两F作用在长边的三分点处，则包含x的三角函数内添加的系数不能为3，否则导出第二项三角函数未知系数等于0；包含x的三角函数内添加的系数不能为偶数，否则仍导出第二项三角函数未知系数等于0。

本文尝试选择第二项为 $B\sin \frac{2\pi x}{a}\sin \frac{\pi y}{b}$ 、 $B\sin \frac{3\pi x}{a}\sin \frac{\pi y}{b}$ 、 $B\sin \frac{4\pi x}{a}\sin \frac{\pi y}{b}$ 进行推导时均失败，直至选择 $B\sin \frac{5\pi x}{a}\sin \frac{\pi y}{b}$ 时才得到正确的解。

3.2. 有限元法求临界荷载

采用ABAQUS软件对上节例子进行分析，以mm为单位建立板结构的有限元实体模型，纵向长度为60，横向长度为40，厚度为1。材料属性设置为Q345钢材， $E=2.06\times {10}^5\text{MPa}$ ， $\mu =0.3$ 。以N为单位在板两长边三分点处对称施加1000的纵向集中力，网格宽度为0.5。边界条件为四边简支，处理方式为板的下底面四边设置一般支撑，dz = 0。边界条件及荷载施加情况及网格划分情况如图2所示。

采用线性屈曲理论，对该桥进行特征值稳定分析，对板结构施加纵向集中力1000 N，计算结构前5阶屈曲模态分析结果，如图3所示。

一阶屈曲模态下临界荷载系数(即特征值)为5.0718，二阶屈曲模态下临界荷载系数为11.662，三阶屈曲模态下临界荷载系数为20.091，四阶屈曲模态下临界荷载系数为20.635，五阶屈曲模态下临界荷载系数为22.037。

由数值计算得到的一阶临界荷载为 $F_C=1000N\times 5.0718=5071.8\text{N}$ ，而使用能量法基于一项三角级数挠度表达式和两项三角级数挠度表达式计算的临界荷载分别为4850.8 N、4782.8 N，有限元法的计算结果与能量法的计算结果相差小于数值计算结果的6%，可以看出用能量法计算得到的结果与用数值方法计算得到的结果是十分接近的，验证了计算结果的准确性。然而，有限元特征值分析得到的结果比能量法的理论解要更大，这意味着在实际工程中选用能量法计算是更偏安全的。

4. 结论

本文分别使用能量法与有限元方法计算了四边简支矩形薄板的压屈临界荷载，并将两种方法的计算结果进行了对比，得出以下结论：

1) 在能量法中，挠度表达式一般是由三角级数表示，其取值项数对计算结果的精度有影响，计算发现，取两项时的计算结果小于取一项时的计算结果。

2) 在进行两项三角级数挠度表达式中第二项三角函数的选取时，需根据荷载作用的位置进行选取，选取不当会导出第二项三角函数未知系数为0的结果。

3) 在有限元分析中，将特征值屈曲计算得到的一阶临界荷载与能量法计算结果进行了对比，结果证明能量法在工程应用中更偏于安全。

参考文献