土木与环境工程学报  2020, Vol. 42 Issue (4): 175-184   doi: 10.11835/j.issn.2096-6717.2020.023   PDF    
基于变分渐近法的双周期梯形波纹板等效模型
邓兵 , 钟轶峰 , 席森彪 , 禹辽     
重庆大学 土木工程学院; 山地城镇建设与新技术教育部重点实验室, 重庆 400045
收稿日期：2019-11-24
基金项目：国家自然科学基金(51778088)；重庆市自然科学基金(cstc2017jcyjBX00036)
作者简介：邓兵(1993-), 男, 主要从事轻质高强结构与材料研究, E-mail:20171602025t@cqu.edu.cn.
通信作者：钟轶峰(通信作者), 男, 教授, 博士生导师, E-mail:zhongjy58@sina.com.
摘要：双周期梯形波纹板是沿两个平面方向均呈周期性梯形波纹变化的新型轻质结构，现阶段对其力学性能与结构参数之间关系的研究较少，限制了该结构的推广应用。针对该问题，以凸起为中心的典型代表单元单胞为研究对象，利用变分渐近法建立单胞等效刚度的数值计算模型。再通过均匀化技术将双周期梯形波纹板转换为具有相同刚度特性的正交异性板进行分析。基于该方法分别计算出不同结构参数下的等效刚度、宏观变形和屈曲模态，通过与三维有限元模拟结果对比验证等效模型的精确性。计算结果表明：由于形貌变化，双周期梯形波纹板的等效刚度较平板的拉伸刚度降低，弯曲刚度增加；随着凸起间距的增大和板高度的减小，波纹板的拉伸刚度逐渐增大，弯曲刚度逐渐降低；而随着板厚度和腰部长度的增加，弯曲刚度呈上升的趋势。构建模型适用于波纹周期远小于结构尺寸的情况，由于等效板模型的近似能量与原三维波纹板能量尽可能接近(通过对能量泛函变分主导项的渐近扩展分析加以保证)，可用于计算不同材料波纹板等效弯矩和等效最大拉压应力等，能有效减少计算量和节省计算资源。
关键词：变分渐近法    梯形波纹板    等效刚度    等效模型    降维分析    
Equivalent model of bi-periodic trapezoidal corrugated plate based on variational asymptotic method
Deng Bing , Zhong Yifeng , Xi Senbiao , Yu Liao     
School of Civil Engineering; Key Laboratory of New Technology for Construction of Cities in Mountain Area, Ministry of Education, Chongqing University, Chongqing 400045, P. R. China
Received: 2019-11-24
Foundation item: National Natural Science Foundation of China (No. 51778088); Chongqing Natural Science Foundation (No. cstc2017jcyjBX00036)
Author brief: Deng Bing(1993-), main research interests: lightweight and high strength structure and materials, E-mail: 20171602025t@cqu.edu.cn.
Corresponding author: Zhong Yifeng(corresponding author), professor, doctorial supervisor, E-mail: zhongjy58@sina.com.
Abstract: The bi-periodic trapezoidal corrugated plate is a new type of lightweight structure with periodic trapezoidal corrugated changes along two plane directions. At present, there is little research on the relationship between its mechanical properties and structural parameters, limiting the popularization and application of this structure. In this article, a numerical model of the represent element (or called unit cell) is established by using the variational asymptotic method to calculate the equivalent stiffness. Then, the bi-periodic trapezoidal corrugated plate is converted to orthotropic plate with the same stiffness characteristics by homogenization technology. Based on this method, the equivalent stiffness with different structural parameters, global deformation and buckling modes are calculated, and the accuracy of the equivalent model is verified by comparing with the three-dimensional finite element simulation results. The results show that the equivalent stiffness of the bi-periodic trapezoidal corrugated plate is lower than that of the plate due to the change of the shape, and the bending stiffness is increased. The tensile stiffness increase and the bending stiffness decrease with increasing interval length and decreasing plate height. While with the increase of plate thickness and web length, the bending stiffness shows an upward trend. The constructed model is suitable for the case that the period of corrugation is far less than the structure size. The approximate energy of the equivalent plate model is as close as that of the original three-dimensional corrugation plate, which is guaranteed by the asymptotic expansion analysis of the leading terms of energy functional. The equivalent plate model can be used to calculate the equivalent bending moment and the equivalent maximum tensile/compressive stress of corrugated plates with different material properties, which can effectively reduce the calculation amount and save the calculation resources.
Keywords: variational asymptotic method    trapezoidal corrugated panel    equivalent stiffness    equivalent model    dimensional reduction analysis    

双周期梯形波纹板是沿两个平面方向均呈周期性梯形波纹变化的新型轻质结构，具有较强吸能质量比，较好抗弯曲性能。与蜂窝等结构相似，双周期梯形波纹结构也可以作为夹层结构的芯层使用，由于纵横波周期分布，双周期梯形波纹板具有很大的惯性矩和抗弯刚度，有效提高了整个结构的强度和承载能力，并大大降低自重而不影响其他使用性能。该结构既满足了强度和刚度的力学性能要求，又能保证材料对吸能和隔热等功能要求[1-5]。

由于波纹板复杂的形貌特征，宏观刚度特性难以通过弹塑性理论方法得到[6]。虽然有限元法可用于评估结构的刚度，但梯形波纹结构形状复杂，数值模拟需要精细网格划分，导致较高的计算成本。由于这类板壳结构的构造和材料分布在面内具有周期性，且周期的尺寸相比整个板的宏观尺寸小很多，所以在宏观尺寸上可以看作是均匀分布的，可使用等效模型以简化复杂波纹板结构的建模[7-9]。

学者们为获得复杂波纹板结构在各种工况载荷下力学行为的等效模型进行了大量的研究。Briassoulis[10]研究了矩形波纹板的等效抗弯刚度。Samanta等[11]在考虑拉伸和弯曲刚度的基础上，对梯形波纹板进行了静态和动态分析。Yokozeki等[12]通过实验和分析研究了由碳环氧复合材料制成的波纹层压板的性能。Peng等[13]通过无网格迦辽金法研究了正弦波形和梯形波纹板的等效弹性性质。Liew等[14]使用该方法进行波纹板的几何非线性分析。张勇等[15]提出了一种将波纹板等效为正交异性板的方法，根据等刚度公式可以将板件的材料参数进行等效。吴存利等[16]将复合材料层合板理论拓展到复合材料波纹板，得到拉伸、耦合和弯曲刚度的解析表达式。高轩能等[17]对槽型波纹板局部屈曲进行了研究，得出整体结构的极限承载力取决于其在半跨荷载作用下槽型波纹截面板的局部屈曲承载力的结论。

现阶段对波纹板的研究主要集中在单周期结构(仅一个方向有周期波纹变化)上，双周期梯形波纹板的力学性能分析缺少必要的研究。笔者利用波纹形状的周期远小于整个板结构尺寸的特点，以凸起为中心的典型代表单元——单胞为研究对象，利用变分渐近法[18-20]建立单胞等效刚度的数值计算模型。再通过均匀化技术将双周期梯形波纹板等效为具有相同刚度特性的正交异性板进行分析。基于该方法分别计算研究不同结构参数下的等效刚度和变形响应，并与三维有限元模拟结果进行对比，以验证等效刚度的精确性。

本文利用结构最小构建单元——单胞作为黑箱，通过渐近均匀化过程得到等效板模型的等效刚度，大大降低了三维有限元求解的繁琐程度，便于工程应用，尤其在求解大规模双周期梯形波纹板刚度中，具有计算简便的优点。

1 双周期梯形波纹板降维分析

对于构造上正交各向异性板的等效建模一般采用理论推导的方法，利用弹性力学理论得到该简化板在各方向的刚度，进而得到等效材料参数[13-14]。本文从能量概念角度建立双周期梯形波纹板的等效模型。由于双周期梯形波纹板沿平面两个坐标轴方向均具有异质性，其典型单胞为三维结构(图 1(b))。板的参考面可视为二维连续介质，该连续介质的每个材料点都有一个相关的单胞作为其微观结构。

图 1 双周期梯形波纹板降维分析及坐标系 Fig. 1 Dimensional reduction analysis and coordinate system of bi-periodic trapezoidal corrugated plate

将典型三维单胞分割为l1、l2、l3、l4、l55部分(其中l4=l2，l5=l1)，波纹板高度为h，板厚度为t，如图 1(b)所示。引入两组坐标系：全局坐标系xi和局部坐标系yi。对于降维后的板结构，场变量是定义在参考面上的(x1, x2)的函数，x3消失(如图 1所示)。全局和局部坐标的关系可表示为yi=xi/η(η是小参数，下标i, j=1, 2, 3；α=1, 2)。在降维模型中，原三维结构的场变量可用剩余宏观坐标xα和局部坐标yj表示，其偏导数为

$ \begin{array}{l} \frac{{\partial f({x_\alpha };{y_j})}}{{\partial {x_\alpha }}} = {\left. {\frac{{\partial f({x_\alpha };{y_j})}}{{\partial {x_\alpha }}}} \right|_{{y_j} = {\rm{const}}}} + \frac{1}{\eta }{\left. {\frac{{\partial f({x_\alpha };{y_j})}}{{\partial {y_j}}}} \right|_{{x_\alpha } = {\rm{const}}}} \equiv \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} f{,_\alpha } + \frac{1}{\eta }{f_{;j}} \end{array} $ (1)

为使用变分渐近法进行双周期梯形波纹板分析，首先需用二维应变量表示原结构的三维位移场

$ \begin{array}{*{20}{c}} {{u_1}({x_\alpha };{y_i}) = v\underline {_1({x_\alpha }) - \eta {y_3}{v_{3,1}}({x_\alpha })} + \eta {w_1}({x_\alpha };{y_i}),}\\ {{u_2}({x_\alpha };{y_i}) = \underline {{v_2}({x_\alpha }) - \eta {y_3}{v_{3,2}}({x_\alpha })} + \eta {w_2}({x_\alpha };{y_i}),}\\ {{u_3}({x_\alpha };{y_i}) = \underline {{v_3}({x_\alpha })} + \eta {w_3}({x_\alpha };{y_i})} \end{array} $ (2)

式中：ui、vi分别表示原三维板和二维板的位移；wi是待求波动函数，下划线项为参考面(中面)产生的位移，需满足条件

$ h{v_3}({x_\alpha }) = [{u_3}],h{v_\alpha }({x_\alpha }) = [{u_\alpha }] + [\eta {y_3}]{v_{3,2}} $ (3)

式中：[·]表示对单胞体积域积分。

式(3)意味着对波动函数存在约束

$ \langle {w_i}\rangle = 0 $ (4)

三维线弹性理论的应变场可表示为

$ {\varGamma _{ij}} = \frac{1}{2}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) $ (5)

将式(5)代入式(2)，并根据变分渐近法忽略高阶小项，得到三维应变场为

$ \begin{array}{*{20}{c}} {{\varGamma _{11}} = {\varepsilon _{11}} + \eta {y_3}{\kappa _{11}} + {w_{1,1}}}\\ {2{\varGamma _{12}} = 2{\varepsilon _{12}} + 2\eta {y_3}{\kappa _{12}} + {w_{1,2}} + {w_{2,1}}}\\ {{\varGamma _{22}} = {\varepsilon _{22}} + \eta {y_3}{\kappa _{22}} + {w_{2,2}}}\\ {2{\varGamma _{13}} = {w_{1,3}} + {w_{3,1}}}\\ {2{\varGamma _{23}} = {w_{2,3}} + {w_{3,2}}}\\ {{\varGamma _{33}} = {w_{3,3}}} \end{array} $ (6)

式中：二维面内拉伸应变εαβ与弯曲曲率καβ可分别表示为

$ {\varepsilon _{\alpha \beta }}({x_1},{x_2}) = ({v_{\alpha ,\beta }} + {v_{\beta ,\alpha }})/2,{\kappa _{\alpha \beta }}({x_1},{x_2}) = - {v_{3,\alpha \beta }} $ (7)

三维应变场可表示为矩阵形式

$ \begin{array}{*{20}{c}} {{\varGamma _{\rm{e}}} = {{[{\varGamma _{11}}{\varGamma _{22}}2{\varGamma _{12}}]}^{\rm{T}}} = \varepsilon + {x_3}\kappa + {\partial _{\rm{e}}}{w_{\left| {\ \ } \right|}},}\\ {2{\varGamma _{\rm{s}}} = {{[2{\varGamma _{13}}2{\varGamma _{23}}]}^{\rm{T}}} = {w_{\left| {\ \ } \right|,3}} + {\partial _{\rm{t}}}{w_3},}\\ {{\varGamma _{\rm{t}}} = {\varGamma _{33}} = {w_{3,3}}} \end{array} $ (8)
$ {\rm{式中}}:{\left( {\ \ } \right)_{\left| {\ \ } \right|}} = {\left[ {{{\left( {\ \ } \right)}_1}{{\left( {\ \ } \right)}_2}} \right]^{\rm{T}}}, $
$ \varepsilon = {[{\varepsilon _{11}}2{\varepsilon _{12}}{\varepsilon _{22}}]^{\rm{T}}},\kappa = {[{\kappa _{11}}{\kappa _{12}} + {\kappa _{21}}{\kappa _{22}}]^{\rm{T}}}, $
$ {\partial _{\rm{e}}} = \left[ {\begin{array}{*{20}{c}} {{{\left( {{\kern 1pt} {\kern 1pt} } \right)}_{,1}}}&0\\ {{{\left( {{\kern 1pt} {\kern 1pt} } \right)}_{,2}}}&{{{\left( {{\kern 1pt} {\kern 1pt} } \right)}_{,1}}}\\ 0&{{{\left( {{\kern 1pt} {\kern 1pt} } \right)}_{,2}}} \end{array}} \right],{\partial _{\rm{t}}} = \left\{ {\begin{array}{*{20}{c}} {{{\left( {{\kern 1pt} {\kern 1pt} } \right)}_{,1}}}\\ {{{\left( {{\kern 1pt} {\kern 1pt} } \right)}_{,2}}} \end{array}} \right\} $

控制梯形波纹板线弹性行为的应变能可表示为

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {U = \frac{1}{2}[\int_{ - {l_3}/2}^{{l_3}/2} {\int_{ - {l_3}/2}^{{l_3}/2} {\int_0^h {\varGamma _3^{\rm{T}}} } } {\mathit{\boldsymbol{D}}_3}{\varGamma _3}{\rm{d}}{y_3}{\rm{d}}{y_2}{\rm{d}}{y_1} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2\int_{{l_3}/2}^{{l_3}/2 + {l_2}} {\int_{{l_3}/2}^{{l_3}/2 + {l_2}} {\int_0^h {\varGamma _2^{\rm{T}}} } } {\mathit{\boldsymbol{D}}_2}{\varGamma _2}{\rm{d}}{y_3}{\rm{d}}{y_2}{\rm{d}}{y_1} + } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2\int_{{l_3}/2 + {l_2}}^{{l_3}/2 + {l_2} + {l_1}} {\int_{{l_3}/2 + {l_2}}^{{l_3}/2 + {l_2} + {l_1}} {\int_0^t {\varGamma _1^{\rm{T}}} } } {\mathit{\boldsymbol{D}}_1}{\varGamma _1}{\rm{d}}{y_3}{\rm{d}}{y_2}{\rm{d}}{y_1}] \equiv \\ \frac{1}{2}\langle {\varGamma ^{\rm{T}}}\mathit{\boldsymbol{D}}\varGamma \rangle = \frac{1}{2}\left\langle {{{\left\{ {\begin{array}{*{20}{l}} {{\varGamma _{\rm{e}}}}\\ {2{\varGamma _{\rm{s}}}}\\ {{\varGamma _{\rm{t}}}} \end{array}} \right\}}^{\rm{T}}}\left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{D}}_{\rm{e}}}}&{{\mathit{\boldsymbol{D}}_{{\rm{es}}}}}&{{\mathit{\boldsymbol{D}}_{{\rm{et}}}}}\\ {\mathit{\boldsymbol{D}}_{{\rm{es}}}^{\rm{T}}}&{{\mathit{\boldsymbol{D}}_{\rm{s}}}}&{{\mathit{\boldsymbol{D}}_{{\rm{st}}}}}\\ {\mathit{\boldsymbol{D}}_{{\rm{et}}}^{\rm{T}}}&{\mathit{\boldsymbol{D}}_{{\rm{st}}}^{\rm{T}}}&{{\mathit{\boldsymbol{D}}_{\rm{t}}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} {{\varGamma _{\rm{e}}}}\\ {2{\varGamma _{\rm{s}}}}\\ {{\varGamma _{\rm{t}}}} \end{array}} \right\}} \right\rangle \end{array} $ (9)

式中：De、Des、Det、Ds、Dst、Dt是三维6×6阶材料矩阵的对应分块矩阵。

外载所做的虚功可表示为

$ \delta {\bar W_{3\mathit{\boldsymbol{D}}}} = \delta {{\rm{\bar W}}_{2\mathit{\boldsymbol{D}}}} + \delta {{\rm{\bar W}}^*} $ (10)

式中：$\delta {\overline {\rm{W}} _{2\mathit{\boldsymbol{D}}}}, \delta {{\bar W}^*} $分别为与波动函数wi无关和相关的虚功，其表达式分别为

$ {\delta {{\bar W}_{2\mathit{\boldsymbol{D}}}} = \langle {p_i}\delta {v_i} + {q_{\rm{ \mathsf{ α} }}}\delta {v_{3,{\rm{ \mathsf{ α} }}}}\rangle } $ (11)
$ {\delta {{\bar W}^*} = \langle \langle {f_i}\delta {w_i}\rangle + {\tau _i}\delta w_i^ + + {\beta _i}\delta w_i^ - \rangle } $ (12)

其中：$ {{(\bullet )}^{+}}={{\left. (\bullet ) \right|}_{{{x}_{3}}=h/2}}, {{(\bullet )}^{-}}={{\left. (\bullet ) \right|}_{{{x}_{3}}=-h/2}}$表示作用在板顶、底面上的量；Ω表示二维板面域；τi、βi、fi分别表示作用板顶、底面的面力集度和体力；$ p_{i}=\left\langle f_{i}\right\rangle+\alpha_{i}+\beta_{i}, q_{\alpha}=\frac{h}{2}\left(\beta_{\alpha}-\alpha_{\alpha}\right)-\left\langle x_{3} f_{\alpha}\right\rangle$。

式(10)中的虚功$ \delta {\overline W _{2\mathit{\boldsymbol{D}}}}$在求解以波动函数为未知量的问题时可忽略，总势能密度可表示为

$ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \delta \Pi = \delta U - \delta {W^*} = }\\ {\frac{1}{2}\delta \langle {\varGamma ^{\rm{T}}}\mathit{\boldsymbol{D}}\varGamma \rangle - \langle \langle {f_i}\delta {w_i}\rangle + {\tau _i}\delta w_i^ + + {\beta _i}\delta w_i^ - \rangle = 0} \end{array} $ (13)

式中：仅wi是变化的。

2 梯形波纹板的变分渐近降维分析

为使用变分渐近法求解式(13)中的未知翘曲函数wi，首先需评估各项的阶数。由于梯形波纹板的厚度与边长之比很小($\eta=h / L \ll 1 $)，可利用该小参数进行降维分析。式(13)中各项的阶数评估为

$ {\varGamma _{ij}}\backsim {\varepsilon _{\alpha \beta }}\backsim h{\kappa _{\alpha \beta }}\backsim n,{w_i}\backsim hn,{w_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }}\backsim {w_{3,\alpha }}\backsim \frac{h}{L}n $ (14)
$ {{w_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}}\backsim {w_{3,3}}\backsim n,h{f_\alpha }\backsim {\alpha _\alpha }\backsim {\beta _\alpha }\backsim \mu \frac{h}{L}n,} $
$ {h{f_3}\backsim {\alpha _3}\backsim {\beta _3}\backsim \mu {{\left( {\frac{h}{L}} \right)}^2}n} $

式中：h与L分别表示板的厚度和边长；n为最小应变量的阶数；μ为材料属性的阶数。

2.1 零阶近似

总势能密度的显式表达式为

$ \begin{array}{l} 2\Pi = \langle \underline{\underline {{{(\varepsilon + {x_3}\kappa )}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{e}}}(\varepsilon + {x_3}\kappa )}} + \underline {2{{(\varepsilon + {x_3}\kappa )}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{e}}}{\partial _{\rm{e}}}{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }}} + \underline {2{{({\partial _e}{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }})}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{e}}}{\partial _{\rm{e}}}{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }}} + 2{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{D_{{\rm{es}}}}{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \underline {2{{(\varepsilon + {x_3}\kappa )}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{es}}}}{\partial _{\rm{t}}}{w_{3,\alpha }}} + \underline {2{{({\partial _{\rm{e}}}{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }})}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{es}}}}({\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + {\partial _{\rm{t}}}{w_{3,\alpha }})} + 2{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{et}}}}{w_{3,3}} + \underline {2{{({\partial _e}{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }})}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{et}}}}{w_{3,3}}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} w_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{s}}}{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + \underline {2w_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{s}}}{\partial _{\rm{t}}}{w_{3,\alpha }}} + \underline {2{{({\partial _{\rm{t}}}{\mathit{\boldsymbol{w}}_{3,\alpha }})}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{s}}}{\partial _{\rm{t}}}{w_{3,\alpha }}} + 2\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{{\rm{st}}}}{w_{3,3}} + \underline {2{{({\partial _{\rm{t}}}{w_{3,\alpha }})}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{st}}}}{w_{3,3}} + {\mathit{\boldsymbol{D}}_{\rm{t}}}w_{3,3}^2} \rangle - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \underline {2([f_i^{\rm{T}}{w_i}] + \tau _i^{\rm{T}}w_i^{\rm{T}} + \beta _i^{\rm{T}}w_i^{\rm{T}})} \end{array} $ (15)

式中：下划线项为h/L阶或更高阶项，在零阶近似时可忽略；双下划线项为常量，不影响求解未知波动函数，也可忽略。

引入拉格朗日乘子λi以考虑对波动函数的约束，

$ \delta (\Pi + {\lambda _i}\langle {\mathit{\boldsymbol{w}}_i}\rangle ) = 0 $ (16)

求解波动函数的零阶近似变分表达式为

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \langle [{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{es}}}} + \mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{s}}} + \mathit{\boldsymbol{w}}_{3,3}^{\rm{T}}\mathit{\boldsymbol{D}}_{{\rm{st}}}^{\rm{T}}]\delta {\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + {\lambda _i}\delta {\mathit{\boldsymbol{w}}_i} + \\ [{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{et}}}} + \mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{{\rm{st}}}} + \mathit{\boldsymbol{w}}_{3,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_t}]\delta {\mathit{\boldsymbol{w}}_{3,3}}\rangle = 0 \end{array} $ (17)

对式(17)部分积分，得到相应的Euler-Lagrange方程

$ \begin{array}{l} [{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{es}}}} + \mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{s}}} + \mathit{\boldsymbol{w}}_{3,3}^{\rm{T}}\mathit{\boldsymbol{D}}_{{\rm{st}}}^{\rm{T}}]{,_{,3}} = {\lambda _{\left| {{\kern 1pt} {\kern 1pt} } \right|}}\\ [{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{et}}}} + \mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{{\rm{st}}}} + \mathit{\boldsymbol{w}}_{3,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{t}}}]{,_{,3}} = {\lambda _3} \end{array} $ (18)

式中：$ \lambda_{\|}=\left[\lambda_{1} \lambda_{2}\right]^{\mathrm{T}}$。

根据自由表面条件，式(18)方括号内的表达式在板顶、底面应为零，可定义板顶、底面边界条件为

$ \begin{array}{l} {[{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{es}}}} + \mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{s}}} + \mathit{\boldsymbol{w}}_{3,3}^{\rm{T}}\mathit{\boldsymbol{D}}_{{\rm{st}}}^{\rm{T}}]^{ + / - }} = 0\\ {[{(\varepsilon + {x_3}\kappa )^{\rm{T}}}{\mathit{\boldsymbol{D}}_{{\rm{et}}}} + \mathit{\boldsymbol{D}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{{\rm{st}}}} + \mathit{\boldsymbol{D}}_{3,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{t}}}]^{ + / - }} = 0 \end{array} $ (19)

式中：上标“+/-”表示作用在板顶、底面上的量。

将上述条件代入式(18)，可求解w||, w3

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} = {{\langle - (\varepsilon + {x_3}\kappa ){{\mathit{\boldsymbol{\bar D}}}_{{\rm{es}}}}{{({\mathit{\boldsymbol{D}}_{\rm{s}}})}^{ - 1}}\rangle }^{\rm{T}}},}\\ {{w_3} = \langle - (\varepsilon + {x_3}\kappa ){{\mathit{\boldsymbol{\bar D}}}_{{\rm{et}}}}{{({\mathit{\boldsymbol{D}}_{\rm{t}}})}^{ - 1}}\rangle } \end{array} $ (20)

式中：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar D}}}_{{\rm{es}}}} = {\mathit{\boldsymbol{D}}_{{\rm{es}}}} - {{\mathit{\boldsymbol{\bar D}}}_{{\rm{et}}}}{{({\mathit{\boldsymbol{D}}_{{\rm{st}}}})}^{\rm{T}}}{{({{\mathit{\boldsymbol{\bar D}}}_{\rm{t}}})}^{ - 1}},}\\ {{{\mathit{\boldsymbol{\bar D}}}_{{\rm{et}}}} = {\mathit{\boldsymbol{D}}_{{\rm{et}}}} - {\mathit{\boldsymbol{D}}_{{\rm{es}}}}{{({\mathit{\boldsymbol{D}}_{\rm{s}}})}^{ - 1}}{\mathit{\boldsymbol{D}}_{{\rm{st}}}},}\\ {{{\mathit{\boldsymbol{\bar D}}}_{\rm{t}}} = {\mathit{\boldsymbol{D}}_{\rm{t}}} - {{({\mathit{\boldsymbol{D}}_{{\rm{st}}}})}^{\rm{T}}}{{({\mathit{\boldsymbol{D}}_{\rm{s}}})}^{ - 1}}{\mathit{\boldsymbol{D}}_{{\rm{st}}}}} \end{array} $ (21)

将求解的波动函数代入式(17)，得到

$ \begin{array}{*{20}{l}} {{U_0} = \frac{1}{2}\langle {{(\varepsilon + {x_3}\kappa )}^{\rm{T}}}{{\mathit{\boldsymbol{\bar D}}}_e}(\varepsilon + {x_3}\kappa )\rangle = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}{{\left\{ {\begin{array}{*{20}{l}} \varepsilon \\ \kappa \end{array}} \right\}}^{\rm{T}}}\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{A}}&\mathit{\boldsymbol{B}}\\ \mathit{\boldsymbol{B}}&\mathit{\boldsymbol{D}} \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} \varepsilon \\ \kappa \end{array}} \right\}} \end{array} $ (22)

式中：A为3×3阶拉伸刚度子矩阵; D为3×3阶弯曲刚度子矩阵; B为3×3阶拉弯耦合刚度子矩阵，其具体表达式为

$ \mathit{\boldsymbol{A}} = \langle {\mathit{\boldsymbol{\bar D}}_{\rm{e}}}\rangle ,\mathit{\boldsymbol{B}} = [{x_3}{\mathit{\boldsymbol{\bar D}}_{\rm{e}}}],\mathit{\boldsymbol{D}} = \langle x_3^2{\mathit{\boldsymbol{\bar D}}_{\rm{e}}}\rangle $ (23)

式中：$\overline{\boldsymbol{D}}_{\mathrm{e}}=\boldsymbol{D}_{\mathrm{e}}-\overline{\boldsymbol{D}}_{\mathrm{es}} \boldsymbol{D}_{\mathrm{s}}^{-1} \boldsymbol{D}_{\mathrm{es}}^{\mathrm{T}}-\overline{\boldsymbol{D}}_{\mathrm{et}} \boldsymbol{D}_{\mathrm{et}}^{\mathrm{T}} / \overline{\boldsymbol{D}}_{\mathrm{t}} $。

零阶近似三维应变场可重构为

$ \varGamma _{\rm{e}}^0 = \varepsilon + {x_3}\kappa ,2\varGamma _{\rm{s}}^0 = - {\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}},\varGamma _{\rm{t}}^0 = {w_{3,3}} $ (24)

应力场重构为

$ \begin{array}{*{20}{c}} {\sigma _{\rm{e}}^0 = {{[\sigma _{11}^0\sigma _{12}^0\sigma _{22}^0]}^{\rm{T}}} = {{\bar D}_{\rm{e}}}(\varepsilon + {x_3}\kappa ),}\\ {\sigma _{\rm{s}}^0 = {{[\sigma _{13}^0\sigma _{23}^0]}^{\rm{T}}} = 0,}\\ {\sigma _{\rm{t}}^0 = \sigma _{33}^0 = 0} \end{array} $ (25)
2.2 一阶近似

由式(25)可知，零阶近似只能得到平面应力，为得到分析破坏至关重要的面外应力，需进行下阶近似。为此，可将零阶翘曲函数摄动为

$ {\mathit{\boldsymbol{w}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} = {\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}},{\mathit{\boldsymbol{w}}_3} = {\mathit{\boldsymbol{\bar v}}_3} + {\mathit{\boldsymbol{D}}_ \bot }\chi $ (26)

将式(26)代入式(13)，得到一阶近似主导项为

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2{\Pi _1} = \langle \mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{s}}}{{\mathit{\boldsymbol{\bar v}}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + {\mathit{\boldsymbol{D}}_{\rm{t}}}\mathit{\boldsymbol{\bar v}}_{3,3}^2 + 2\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}^{\rm{T}}{\mathit{\boldsymbol{C}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}}\chi {,_\alpha } + \\ 2\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}^{\rm{T}}{\mathit{\boldsymbol{D}}_{\rm{s}}}{\partial _{\rm{t}}}{\mathit{\boldsymbol{D}}_ \bot }{\chi _{,\alpha }} - 2\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}^{\rm{T}}{p_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} - 2\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}^{ + {\rm{T}}}{\tau _{\left| {{\kern 1pt} {\kern 1pt} } \right|}} - 2\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}^{ - {\rm{T}}}{\beta _{\left| {{\kern 1pt} {\kern 1pt} } \right|}}\rangle \end{array} $ (27)

式(27)的Euler-Lagrange方程为

$ ({\mathit{\boldsymbol{D}}_{\rm{s}}}{\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + {\mathit{\boldsymbol{D}}_{\rm{s}}}{\partial _{\rm{t}}}{\mathit{\boldsymbol{D}}_ \bot }{\chi _{,\alpha }}){,_{,3}} = {\mathit{\boldsymbol{C}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}}{\partial _{\rm{t}}}{\mathit{\boldsymbol{D}}_ \bot }{\chi _{,\alpha }} + {\mathit{\boldsymbol{g}}_{,3}} + {\mathit{\boldsymbol{\lambda }}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} $ (28)

式中：$ {\mathit{\boldsymbol{C}}_{||}} = - \partial _{\bf{e}}^{\bf{T}}\left[ {{\mathit{\boldsymbol{D}}_{||}}{x_3}{\mathit{\boldsymbol{D}}_{||}}} \right], {\mathit{\boldsymbol{g}}_{, 3}} = - {p_{||}}$。

$\bar{v}_{3} $与$ {{\mathit{\boldsymbol{\bar v}}}_{||}}$完全独立，由式(6)可知v3仅有平凡解。求解式(28)和式(6)，得到翘曲函数${\mathit{\boldsymbol{v}}_{||}} $解为

$ {\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} = ({\mathit{\boldsymbol{\bar C}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} + {\mathit{\boldsymbol{L}}_\alpha }){\chi _{,\alpha }} + \mathit{\boldsymbol{\bar g}} $ (29)

式中：

$ {{\mathit{\boldsymbol{\bar C}}}_{\left| {{\kern 1pt} {\kern 1pt} } \right| \cdot 3}} = \mathit{\boldsymbol{D}}_{\rm{s}}^{ - 1}{\mathit{\boldsymbol{C}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}},\langle {{\mathit{\boldsymbol{\bar C}}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}}\rangle = 0,{{\mathit{\boldsymbol{\bar g}}}_{,3}} = \mathit{\boldsymbol{D}}_{\rm{s}}^{ - 1}\mathit{\boldsymbol{\bar g}},[\mathit{\boldsymbol{\bar g}}] = 0, $
$ {\mathit{\boldsymbol{L}}_\alpha }{\mathit{\boldsymbol{\chi }}_{,\alpha }} = {\mathit{\boldsymbol{c}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}}/h,{{\mathit{\boldsymbol{\bar C}}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} = {\mathit{\boldsymbol{C}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}} + \frac{{{x_3}}}{h}\mathit{\boldsymbol{D}}_\alpha ^ \mp - \frac{1}{2}\mathit{\boldsymbol{D}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}^ \pm - {\mathit{\boldsymbol{D}}_{\rm{s}}}{e_\alpha }{\mathit{\boldsymbol{D}}_ \bot }, $
$ \mathit{\boldsymbol{\bar g}} = \mathit{\boldsymbol{g}} + \frac{{{x_3}}}{h}{\mathit{\boldsymbol{g}}^ \mp } - \frac{1}{2}{\mathit{\boldsymbol{g}}^ \pm } $

其中：()±=()++()-, ()=()--()+。

至此，修正到一阶的近似能量为

$ 2{\Pi _1} = {\chi ^{\rm{T}}}{\mathit{\boldsymbol{\bar D}}_{\rm{e}}}\chi + \chi _{,\alpha }^{\rm{T}}{\mathit{\boldsymbol{B}}_{\alpha \beta }}{\chi _{,\beta }} - 2{\chi ^{\rm{T}}}F $ (30)

式中：

$ {\mathit{\boldsymbol{B}}_{\alpha \beta }} = \langle {\mathit{\boldsymbol{D}}_{{{\rm{s}}_{\alpha \beta }}}}\mathit{\boldsymbol{D}}_ \bot ^{\rm{T}}{\mathit{\boldsymbol{D}}_ \bot } - \mathit{\boldsymbol{\bar C}}_\alpha ^{\rm{T}}\mathit{\boldsymbol{D}}_{\rm{s}}^{ - 1}{\mathit{\boldsymbol{\bar C}}_\beta }\rangle + L_\alpha ^{\rm{T}}\langle {C_{\beta ,3}}\rangle ; $ (31)
$ \mathit{\boldsymbol{F}} = [\mathit{\boldsymbol{D}}_ \bot ^{\rm{T}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {p_3}] - [\mathit{\boldsymbol{\bar C}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}^{\rm{T}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{D}}_{\rm{s}}^{ - 1}{\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{\bar g}}_{,\alpha }}] - {\mathit{\boldsymbol{L}}_\alpha }(\langle \bar p\rangle + \langle {p_{\left| {{\kern 1pt} {\kern 1pt} } \right|}}\rangle ){,_{,\alpha }} $
2.3 三维场重构

降维模型的可靠性取决于其对原三维结构位移、应力和应变场重构的准确性。为此，有必要提供重构关系以完善降维模型，即用二维变量和x3表征三维位移、应力和应变场。

由式(4)，可重构三维位移场为

$ {U_i} = {u_i} + {x_3}({C_{3i}} - {\delta _{3i}}) + {C_{ji}}{\mathit{\boldsymbol{w}}_j} $ (32)

式中：Ui、ui分别为三维位移和二维位移。

由式(12)重构三维应变场为

$ {\varGamma _{\rm{e}}} = \varepsilon + {x_3}\kappa ,2{\varGamma _{\rm{s}}} = {\mathit{\boldsymbol{\bar v}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + {\partial _{\rm{t}}}{\mathit{\boldsymbol{D}}_ \bot }{\chi _{,\alpha }},{\varGamma _t} = {\mathit{\boldsymbol{D}}_{ \bot ,3}}\chi $ (33)

三维应力场可重构为

$ \begin{array}{*{20}{c}} {{\sigma _{\rm{e}}} = {{[{\sigma _{11}}{\sigma _{12}}{\sigma _{22}}]}^{\rm{T}}} = {\mathit{\boldsymbol{D}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|}}(\varepsilon + {x_3}\kappa ) + {\mathit{\boldsymbol{D}}_{\rm{e}}}{\partial _{\rm{e}}}{{\mathit{\boldsymbol{\bar v}}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }},}\\ {2{\sigma _{\rm{s}}} = {{[{\sigma _{13}}{\sigma _{23}}]}^{\rm{T}}} = {\mathit{\boldsymbol{D}}_s}({{\bar v}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,3}} + {\partial _t}{\mathit{\boldsymbol{D}}_ \bot }{\chi _{,\alpha }}),}\\ {{\sigma _{\rm{t}}} = {\sigma _{33}} = \mathit{\boldsymbol{D}}_{{\rm{et}}}^{\rm{T}}{\partial _{\rm{e}}}{{\mathit{\boldsymbol{\bar v}}}_{\left| {{\kern 1pt} {\kern 1pt} } \right|,\alpha }}} \end{array} $ (34)

由上述推导过程可知，通过提取出双周期波纹板结构的单胞，分析其单胞结构上的受力特点并结合周期性边界条件，即可得到等效刚度特性(见图 2)；再通过均匀化技术将双周期梯形波纹板转换为具有相同刚度特性的正交异性板进行分析；最后，将等效板模型分析得到的全局部位移、应变代入所求位置的三维单胞内进行重构，可得到局部应力、位移场分布。

图 2 双周期波纹板等效刚度特性求解过程 Fig. 2 The solution process of equivalent stiffness of bi-periodic trapezoidal corrugated plate

3 模型验证与讨论

为讨论和分析典型双周期梯形波纹板的刚度特性，同时，对构建的模型进行验证，选取结构参数为(见图 3)：高度h=10 mm，厚度t=0.5 mm，周期长度l=30 mm；结构总长度为20个周期600 mm×600 mm；单个周期l1=5 mm，l2=5 mm，l3=10 mm。材料属性为各向同性：弹性模量E=206 GPa，泊松比v=0.3。

图 3 结构的单胞及网格划分 Fig. 3 Unit cell and meshing of the unit cells

考虑周期结构在x1和x2方向上均具有20个周期，采用ABAQUS有限元软件建立三维模型和等效模型，并赋予相应的结构参数。

3.1 模型验证

由于波纹方向与平面坐标一致，A13、A23、D13、D23消失；对于正交各向异性板，x1和x2方向的扭转刚度是不同的，D22取平均值；双周期波纹形状沿x1和x2方向相对坐标原点对称，则拉弯耦合刚度B11=B13=B22=B33=0；因为两个方向周期间距相同，A11=A22，D11=D22。基于构建的等效模型计算得到的等效刚度列于表 1，并与普通平板的刚度进行比较。由表 1可看出：双周期梯形波纹板的拉伸刚度A较平板分别下降至37%、15%、37%、10%；弯曲刚度D分别增加了2.96、2.23、2.96、15.84倍。这主要是因为双周期梯形波纹板的平面内应变过程需先将梯形拉直，再产生类似平板的拉伸，使得其面内拉伸刚度有所降低；而等效模型为符合Kirchhoff假设，其相对厚度较平板提高很多，等效弯曲刚度较平板有较大增加。在保持总质量不变的前提下，弯曲刚度有较大增加，可实现结构的轻型化。

表 1 双周期梯形波纹板与普通平板的刚度比较 Table 1 Stiffness comparison between bi-periodic trapezoidal corrugated plates and flat Plates

分别建立结构的实体模型和以等效刚度为基础的等效板模型，分别使用411 804个实体单元(C3D20R)和1 575个壳单元(S8R5)。图 4对比了两种模型在不同载荷(跨中施加集中载荷50 N和均匀载荷500 Pa)与边界条件(四边固结CCCC和两边自由两边固结FFCC)下的位移分布。可看出等效模型与三维模型在相同载荷下的位移分布吻合较好，四边固结(CCCC)施加集中力和均布载荷的最大误差分别为1%和1.2%；两边自由两边固结(FFCC)施加集中力和均布载荷的最大误差分别为3.17%和1.6%。验算比较最大误差在5%以内，等效模型具有较好的吻合性。因此，基于变分渐近法得到的等效刚度是准确的，可用于评估宏观刚度特性。

图 4 双周期梯形波纹板在不同载荷和边界条件下的位移比较 Fig. 4 Displacement comparison bi-periodic trapezoidal corrugated plate under different loads and boundary condtions

基于得到的等效刚度，同时对梯形波纹板和等效平板建模，通过线性屈曲分析求解结构的临界特征值列于表 2。

表 2 三维实体结构和等效板全局屈曲模态和临界荷载(N)比较(前5阶) Table 2 Comparison of first five buckling loads (N) and mode shapes between three-dimensional solid structures and equivalent plates

由表 2可以看出，等效板模型的全局屈曲模态和三维梯形波纹板的有限元结果相吻合，屈曲临界荷载误差整体较小；随着模态阶数的增加，误差整体出现增大的趋势，原因在于，等效模型的屈曲模态受屈曲波长的影响，但最大误差(第5阶3.51%)仍低于5%。因此，等效模型具有较高的精确性和有效性。

3.2 局部场重构

基于2.3节推导公式，将四边固结CCCC和集中载荷50 N下等效板模型分析得到板中点处(图 1(c)所示)的全局位移、应变代入对应三维单胞内进行重构，得到局部应力、位移分布(图 5、图 6)。值得注意的是，目前大多数简化模型无法预测出准确的局部场分布。

图 5 板中点处单胞内局部应力分量分布 Fig. 5 Distribution of local stress components within the unit cell at the center of the plate

图 6 板中点处单胞内局部位移分量分布 Fig. 6 Distribution of local displacements within the unit cell at the center of the plate

由图 5应力分布云图得知，在发生应变时，梯形板的应力分布不均匀，凸起没有完全承受荷载，在x=5 mm和x=25 mm处(平板和凸起的相交处)出现最值，应力在平板和凸起的相交处变化比较明显。σ11和σ22在凸起位置有明显应力增大现象且均为负值。σ33沿平板基本保持不变，但在交点处出现应力突变。这就合理解释了在交点处应力先破坏的现象。

根据图 6位移云图，在梯形凸起的顶点处出现极值，单胞突起的对称点呈现相反的位移变化。经过凸起位置，位移会发生突变，随后会保持平板的位移变化特征趋势。重构的U3最大位移与整体模型相一致，最大误差仅为1%。

3.3 几何参数影响分析

双周期梯形波纹板结构参数主要包括l1、l3、h和t。当材料属性不变时，利用图 3(b)所示的三维单胞模型讨论双周期梯形波纹板结构参数对等效刚度特性的影响规律，保持单胞凸起部分长度不变，改变参数l1、l3、h和t，研究等效刚度的变化规律。

图 7为其他参数不变，l1在6~18 mm范围内变化时结构等效刚度变化情况。由图 7可知，随着梯形板的结构参数l1即凸起间距的增大，拉伸刚度逐渐增大，弯曲刚度逐渐降低，且都呈现出非线性变化。主要原因在于，随着凸起间距的增大，板中凸起结构越来越稀疏，结合应力云图可知单胞凸起的组合特性愈发不明显，结构更趋近于普通平板。

图 7 凸起间隔长度对结构刚度特性的影响 Fig. 7 The influence of interval length on the structural stiffness

图 8为其他参数不变，高度h在6~18 mm范围内变化时结构等效刚度变化情况。由图 8可知，随着h增大，拉伸刚度逐渐减小，弯曲刚度逐渐增大，且弯曲刚度较拉伸刚度变化更加明显。主要原因在于随着梯形凸起高度的增加，结构的薄板特性越不明显，构造异性特征更加突出。

图 8 板高度h对结构刚度特性的影响 Fig. 8 The influence of height h on the structural stiffness

图 9为其他参数不变，厚度t在1.0~3.0 mm范围内变化时结构的等效刚度变化情况。由图 9可知，结构随着几何参数板厚t的增加，拉伸刚度和弯曲刚度都逐渐增大。拉伸刚度A随着板厚度t呈现出线性变化，且A22增长较快；弯曲刚度D呈现出非线性增长，增长速率逐渐加快。分析其原因是随着板的厚度增加，梯形板的构造正交异性特性逐渐模糊，呈现实体的特性。

图 9 板厚t对结构刚度特性的影响 Fig. 9 The influence of thickness t on the structural stiffness

腰部长度受到l2和l3的影响，为方便分析，只需要改变结构参数l3即可。图 10为l1=12 mm、h=12 mm、t=2mm、l3在4~12 mm范围内变化时梯形波纹板的等效刚度变化规律。由图 10可知，随着梯形板的结构参数l3增大即顶部变宽时，拉伸刚度逐渐减小，弯曲刚度逐渐增大，分析原因是，在顶部变宽的过程中，凸起部分承受水平荷载能力逐渐减小。l3对结构等效刚度的影响比l1对等效刚度影响小。

图 10 结构参数l3对刚度特性的影响 Fig. 10 The influence of l3 on the structural stiffness

4 结论

1) 基于变分渐近法建立了双周期梯形波纹板的等效模型，得到各刚度系数解析式和局部场重构关系，适用于波纹的周期远小于结构尺寸的情况。由于等效板模型的近似能量与原三维波纹板能量尽可能接近(通过对能量泛函变分主导项的渐近扩展分析加以保证)，可用于不同材料波纹板的计算等效弯矩和等效最大拉、压应力等。

2) 将构建的刚度计算方法以及建模方法应用于双周期梯型波纹钢板全局性能分析中，计算得到的全局位移和屈曲模态与Abaqus有限元结果相吻合，但计算效率大为提高。同时，将等效板模型分析得到的全局部位移、应变代入所求位置的三维单胞内进行重构，得到局部应力、位移分布。

3) 利用三维单胞模型得到了双周期梯形波纹板结构参数(板高度h、板厚度t和凸起间隔长度l1和腰部长度)对等效刚度的影响规律。随着凸起间距的增大和板高度的减小，板凸起的组合特性愈发不明显，更趋近于普通平板，导致拉伸刚度逐渐增大，弯曲刚度逐渐降低；而随着板厚度和腰部长度的增加，梯形波纹板的构造正交异性特性逐渐降低，弯曲刚度呈上升的趋势。未来可结合隔热、减振等性能，对双周期梯形波纹板的实际应用提供理论依据，使其在土木建筑、航空航天等领域得到更广泛的应用。

参考文献
[1]
ABED A M, ALGHOUL M A, SOPIAN K, et al. Design characteristics of corrugated trapezoidal plate heat exchangers using nanofluids[J]. Chemical Engineering and Processing:Process Intensification, 2015, 87(6): 88-103.
[2]
PEREL D, LIBOVE C. Elastic buckling of infinitely long trapezoidally corrugated plates in shear[J]. Journal of Applied Mechanics, 1978, 45(3): 579-582.
[3]
EMAMI F, MOFID M. On the hysteretic behavior of trapezoidally corrugated steel shear walls[J]. The Structural Design of Tall and Special Buildings, 2014, 23(2): 94-104.
[4]
BAHREBAR M, KABIR M Z, ZIRAKIAN T, et al. Structural performance assessment of trapezoidally-corrugated and centrally-perforated steel plate shear walls[J]. Journal of Constructional Steel Research, 2016, 122(11): 584-594.
[5]
SHON S, YOO M, LEE S. An experimental study on the shear hysteresis and energy dissipation of the steel frame with a trapezoidal-corrugated steel plate[J]. Materials, 2017, 10(3): 261.
[6]
王远, 臧勇, 管奔, 等. 基于单胞有限元的波纹板等效刚度特性[J]. 北京航空航天大学学报, 2018, 44(6): 1230-1238.
WANG Y, ZANG Y, GUAN B, et al. Equivalent stiffness property of dimpled sheet based on unit cell finite element[J]. Journal of Beijing University of Aeronautics and Astronautics, 2018, 44(6): 1230-1238. (in Chinese)
[7]
XIA Y, FRISWELL M I, FLORES E I S. Equivalent models of corrugated panels[J]. International Journal of Solids and Structures, 2012, 49(13): 1453-1462.
[8]
LI W, LI H X, LI G Q, et al. Numerical and experimental analysis of composite fouling in corrugated plate heat exchangers[J]. International Journal of Heat and Mass Transfer, 2013, 63: 351-360.
[9]
FAIZAL M, AHMED M R. Experimental studies on a corrugated plate heat exchanger for small temperature difference applications[J]. Experimental Thermal and Fluid Science, 2012, 36(1): 242-248.
[10]
BRIASSOULIS D. Equivalent orthotropic properties of corrugated sheets[J]. Computers & Structures, 1986, 23(2): 129-138.
[11]
SAMANTA A, MUKHOPADHYAY M. Finite element static and dynamic analyses of folded plates[J]. Engineering Structures, 1999, 21(3): 277-287.
[12]
YOKOZEKI T, TAKEDA S I, OGASAWARA T, et al. Mechanical properties of corrugated composites for candidate materials of flexible wing structures[J]. Composites Part A:Applied Science and Manufacturing, 2006, 37(10): 1578-1586.
[13]
PENG L X, LIEW K M, KITIPORNCHAI S. Analysis of stiffened corrugated plates based on the FSDT via the mesh-free method[J]. International Journal of Mechanical Sciences, 2007, 49(3): 364-378.
[14]
LIEW K M, PENG L X, KITIPORNCHAI S. Nonlinear analysis of corrugated plates using a FSDT and a meshfree method[J]. Computer Methods in Applied Mechanics and Engineering, 2007, 196(21-24): 2358-2376.
[15]
张勇, 刘锡良, 张福海. 银河金属拱型波纹屋顶的静力稳定承载力试验研究[J]. 建筑结构学报, 1997, 18(6): 46-54, 40.
ZHANG Y, LIU X L, ZHANG F H. Experimental study on static stability bearing capacity of milky way arched corrugated metal roof[J]. Journal of Building Structures, 1997, 18(6): 46-54, 40. (in Chinese)
[16]
吴存利, 段世慧, 孙侠生. 复合材料波纹板刚度工程计算方法及其在结构分析中的应用[J]. 航空学报, 2008, 29(6): 1570-1575.
WU C L, DUAN S H, SUN X S. Technique for calculation of composite corrugated plate stiffness and its application in structure finite element analysis[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(6): 1570-1575. (in Chinese)
[17]
高轩能, 王健. 槽型波纹板组结构的局部屈曲承载能力研究[J]. 南昌大学学报(工科版), 2004, 26(3): 68-71, 93.
GAO X N, WANG J. Study on local buckling of corrugated channel sections[J]. Journal of Nanchang University(Engineering & Technology), 2004, 26(3): 68-71, 93. (in Chinese)
[18]
钟轶峰, 周小平, 张亮亮. 压电圆柱棒基于变分渐近法的广义Timoshenko模型[J]. 工程力学, 2014, 31(10): 14-20.
ZHONG Y F, ZHOU X P, ZHANG L L. A new generalized Timoshenko model for piezoelectric cylindrical rods by using the variational asymptotic method[J]. Engineering Mechanics, 2014, 31(10): 14-20. (in Chinese)
[19]
钟轶峰, 刘森林, 陈约瑟, 等. 中厚度复合材料夹芯层板变分渐近精细模型[J]. 应用数学和力学, 2013, 34(7): 687-696.
ZHONG Y F, LIU S L, CHEN Y S, et al. Refined model for composite sandwich laminates of moderate thickness based on the variational asymptotic method[J]. Applied Mathematics and Mechanics, 2013, 34(7): 687-696. (in Chinese)
[20]
钟轶峰, 余文斌. 压电复合材料层合板的热压电弹性简化模型研究[J]. 工程力学, 2012, 29(6): 314-319.
ZHONG Y F, YU W B. Study on simplified thermopiezoelastic model for piezoelectric composite laminates[J]. Engineering Mechanics, 2012, 29(6): 314-319. (in Chinese)
X