![](https://lxxb.cstam.org.cn/fileLXXB/richhtml/0459-1879/rich_nlm3_0/images/REemail.gif)
A QUADRILATERAL QUADRATURE PLATE ELEMENT BASED ON REDDY'S HIGHER-ORDER SHEAR DEFORMATION THEORY AND ITS APPLICATION1)
ShenZhiqiang![](https://lxxb.cstam.org.cn/fileLXXB/richhtml/0459-1879/rich_nlm3_0/images/REemail.gif)
中图分类号:O343.1
文献标识码:A
通讯作者:
收稿日期:2018-07-10
网络出版日期:2018-09-18
版权声明:2018力学学报期刊社力学学报期刊社 所有
基金资助:
展开
摘要
关键词:
Abstract
Keywords:
-->0
PDF (753KB)元数据多维度评价相关文章收藏文章
本文引用格式导出EndNoteRisBibtex收藏本文-->
引 言
板结构是工程实践中常见的一类结构形式,特别是近年来随着各类高性能复合材料和功能梯度材料的出现,由这些新型材料所构成的板结构在航空航天、土木工程以及军事防护领域得到了 广泛的应用[1-2].该类板结构的显著特点是其材料性能沿板厚变化,需合理考虑其横向剪切变形,而经典板理论假设横向剪切应变为零、一阶剪切变形理论则假设横向剪切应变沿板厚方向为常量,往往不适用于该类新型材料板结构的力学分析.在这种情况下,高阶剪切变形理论[2-3] (higher-order shear deformation theory)越来越受到学术界关注.在一般的高阶剪切变形理论中,板的位移场被表示为板厚坐标$z$多项式或其他函数形式,可利用板上下表面横向剪应变为零的条件确定各项系数间的关系,横向剪应力的计算不需要像一阶剪切变形理论那样人为地引入剪切修正因子.在各类高阶剪切变形理论中,由Reddy等[4]提出的位移描述形式(式(1))不仅具有一般高阶剪切变形理论的优点,同时用于描述位移场的独立自由度物理意义明确,且计算成本较低,被广泛地应用于各类复合材料、功能梯度材料中厚板的分析.由于传统的解析方法仅适用于特定载荷、几何形式及边界条件问题的求解,国内外****基于高阶剪切变形理论,利用各类数值方法对众多新型材料结构如碳纤维增强复合材料层合板、陶瓷$\!$-$\!$-$\!$金属功能梯度板、功能梯度碳纳米管增强复合材料板的相关力学问题进行了大量研究,这些数值方法中最为常见的是有限元法. 观察式(1)不难发现,$w_{,x} $项和$w_{,y}$项的存在要求单元$C^{1}$连续. 针对这一要求,一方面可通过构造较为复杂的协调板元或者非协调板元予以解决.如Singh等[5]和Liu等[6]利用Hermite插值,分别构造了矩形四节点协调元;Phan等[7]和Ren等[8]则构造了矩形四节点非协调元. 另一方面,不少****通过变量替换,将原问题转化为$C^{0}$连续以简化单元的构造.如Kant等[9]和Pervez等[10]令$\theta _x^\ast = - \left( {4 / {3h^2}} \right)\left( {\phi + w_{,x} }\right)$, Shankara等[11]和Ansari 等[12]则令$\theta _x^\ast = w_{,x}$,Nayakc等[13-14]和Sheikh 等[15-16]则令$\theta _x^\ast = \left( {\phi +w_{,x} } \right)$, 分别构造了$C^{0}$型矩形单元或三角形单元,这样做虽可避开$C^{1}$连续要求,但也牺牲了单元的计算精度.同时,也有****通过构造位移$\!$-$\!$-$\!$力混合单元[17]或阶谱单元[18]的方式来解决这一问题,这在一定程度上增加了问题的复杂度.此外,微分求积法[19]、无网格方法[20-21]等新型数值方法及商业有限元软件[22]也被用于高阶板的力学分析.目前这一领域的研究仍十分活跃,满足$C^{1}$连续的单元相对较少,对于简单高效数值计算方法的需求依然迫切.
弱形式求积元法(weak form quadrature element method) 简称 求积元法,是一种基于待求解问题的弱形式,利用统一的单元节点计算问题控制泛函中积分和微分的新型高阶数值方法[23],具有准确高效、步骤直接、前后处理过程简便等优点,在处理复杂几何形状、非均匀材料、强几何非线性问题上相比传统低阶有限元具有独特的优势[24-28].利用该方法,Zhong和Yue[26]基于经典板理论对等厚度三角形板、方板及变厚度环板进行了弯曲和自由振动分析;Zhong和Wang[27]对均质材料高阶梁的弯曲、自由振动及特征值屈曲进行了分析;作者则基于Euler梁理论,对考虑界面滑移效应的楔形变截面组合梁弯曲和自由振动问题进行了分析[29].
本文基于Reddy高阶剪切变形理论,建立$C^{1}$连续的四边形求积元板单元.对典型均质材料、复合材料、功能梯度材料构成的等厚度矩形板、变厚度矩形板及等厚度斜板的弯曲和自由振动问题进行计算分析,通过与现有文献中相应计算结果进行对比,考察求积元板单元的计算效率.
1 四边形求积元板单元的建立
1.1 高阶剪切变形理论
考虑空间中某任意四边形板,板的各边边长分别为$a$,$b$,$c$和$d$,厚度为$h\left( {x,y}\right)$,同时建立参考坐标系$O - XYZ$,$XOY$平面与板几何中性面重合,如图1所示.![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/thumbnail/img_1.png)
图1空间中的四边形板
-->Fig.1A quadrilateral plate in physical domain
-->
依据Reddy提出的高阶剪切变形理论,板的位移场可表示为
$${\pmb d} = \left[\!\! \begin{array}{c} {u + z\phi - \left( {4/{3h^2}} \right)z^3\left( {\phi + w_{,x} } \right)} {v + z\varphi - \left( {4 /{3h^2}} \right)z^3\left( {\varphi + w_{,y} } \right)} \hfill w \end{array} \!\! \right] (1)$$
式中,$u$, $v$和$w$分别为板几何中性面上物质点沿$X$, $Y$和$Z$向的位移,$\phi $, $ - \varphi$分别为该点处中性面法线绕$Y$, $X$轴转角. $w_{,x} $, $w_{,y}$表示$w$对$x$和$y$的偏导数,下同. 则板的应变场即可依据线性小应变理论表示为
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_2.png)
并与相应的应力分量通过本构关系加以联系,即
$$\left[ \!\!\begin{array}{c} {\sigma _{xx} } {\sigma _{yy} } {\sigma _{xy} } \end{array} \!\! \right] = {\pmb E}_{\rm p} \left[\!\!\begin{array}{c} {\varepsilon _{xx} } {\varepsilon _{yy} } {\gamma _{xy} }\end{array} \!\! \right]\,, \quad\left[\!\!\begin{array}{c} {\sigma _{yz} } {\sigma _{xz} }\end{array} \!\!\right] = {\pmb E}_{\rm t} \left[\!\!\begin{array}{c} {\varepsilon _{yz} } {\gamma _{xz} }\end{array} \!\!\right] (3)$$
式中,${\pmb E}_{\rm p} $和${\pmb E}_{\rm t}$是由弹性模量等常数构成的矩阵,其具体形式与板的构成材料有关. 如对于均质各向同性材料
$$\left. \!\!{\pmb E}_{\rm p} = \dfrac{E}{1 - v^2}\left[\!\!\begin{array}{ccc} 1 & v & 0 v & 1 & 0 0 & 0 & (1 - v) / 2 \end{array}\!\!\right] {\pmb E}_{\rm t} = \dfrac{E}{1 - v^2}\left[\!\!\begin{array}{cc} {(1 - v)} / 2 & 0 0 & {{(1 - v)} /2}\end{array}\!\!\right] \right\} (4)$$
式中,$E$, $\nu $分别为材料的杨氏模量与泊松比. 而对于其他如正交各向异性材料、复合材料及功能梯度材料,${\pmb E}_{\rm p} $和${\pmb E}_{\rm t}$略微复杂,也可沿板厚变化,具体的表达形式可参阅相关教材或文献[1].
1.2 求积元板单元
物理域中四边形求积元板单元如图2(a)所示,4个角点按逆时针方向依次编号(图中编为1$\sim$4),与之相应标准域板单元如图2(b)所示.![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/thumbnail/img_3.png)
图2四边形求积元板单元几何变换及自由度设置
-->Fig. 2The geometric mapping and degrees of freedom for quadrilateral quadrature plate element
-->
该单元$e$的应变能$U^{(e)}$可表示为
$$U^{(e)} = \int_\varOmega \dfrac{1}{2}{\pmb\varepsilon}^{\rm T}{\pmb E}{\pmb\varepsilon} d xd y (5)$$
式中, $\varOmega $为物理域板单元中性面区域. ${\pmb\varepsilon } $和${E}$为广义应变向量与弹性模量矩阵,分别为
$$ {\pmb \varepsilon} = [ {u_{,x} } \ \ {v_{,y} } \ \ {u_{,y} + v_{,x} } \ \ {\phi _{,x} } \ \ {\varphi _{,y} } \ \ {\phi _{,y} + \varphi _{,x} } \qquad {w_{,xx} } \ \ {w_{,yy} } \ \ {2w_{,xy} } \ \ {\phi _y + w_{,y} } \ \ {\phi _x + w_{,x} } ]^{\rm T} (6)$$
及
$${\pmb E} = \int_{ - h /2}^{h /2} \left[ \!\!\begin{array}{cccc} {{\pmb E}_{\rm p} } & {\alpha _1 {\pmb E}_{\rm p} } & {\alpha _2 {\pmb E}_{\rm p} } & {\bf 0} {\alpha _1 {\pmb E}_{\rm p} } & {\alpha _1^2 {\pmb E}_{\rm p} } & {\alpha _1 \alpha _2 {\pmb E}_{\rm p} } & {\bf 0} {\alpha _2 {\pmb E}_{\rm p} } & {\alpha _1 \alpha _2 {\pmb E}_{\rm p} } & {\alpha _2^2 {\pmb E}_{\rm p} } & {\bf 0} {\bf 0} & {\bf 0} & {\bf 0} & {\alpha _3^2 {\pmb E}_{\rm t} } \end{array} \!\!\right] d z (7)$$
式中
$$\alpha _1 = z - \dfrac{4}{3h^2}z^3\,, \ \ \alpha _2 = - \dfrac{4}{3h^2}z^3\,, \ \ \alpha _3 = 1 -\dfrac{4}{h^2}z^2 (8)$$
对于本文涉及的直边四边形板,可利用简单的线性变换,将问题的求解从物理域转换到标准域,即
$$x = \sum_{k = 1}^4 {S_k \left( {\xi ,\eta } \right)x_k } \,, \ y = \sum_{k = 1}^4 {S_k \left( {\xi ,\eta } \right)y_k } (9)$$
式中,($x_k , y_k $)为物理域四边形角点坐标,如图2(a)所示,$k = 1, 2, 3, 4$,且
$$\left. \begin{array}{l} S_1 \left( {\xi ,\eta } \right) = {\left( {1 - \xi } \right)\left( {1 - \eta } \right)} / 4 S_2 \left( {\xi ,\eta } \right) = {\left( {1 + \xi } \right)\left( {1 - \eta } \right)} /4 S_3 \left( {\xi ,\eta } \right) = {\left( {1 + \xi } \right)\left( {1 + \eta } \right)} / 4 {S_4 \left( {\xi ,\eta } \right) = {\left( {1 - \xi } \right)\left( {1 + \eta } \right)} / 4}\end{array} \,, \ { - 1 \leqslant \xi ,\eta \leqslant 1} \right \} (10)$$
并将${\pmb\varepsilon} $中变量对相应物理域坐标的导数转化到标准域,即
$${\pmb\varepsilon} ={\pmb D}\bar {\varepsilon } (11)$$
其中,标准域中的应变向量定义为
$$ \bar {\pmb\varepsilon } = [ {u_{,\xi } } \ \ {u_{,\eta } } \ \ {v_{,\xi } } \ \ {v_{,\eta } } \ \ \phi \ \ {\phi _{,\xi } } \ \ {\phi _{,\eta } } \ \ \varphi \ \ {\varphi _{,\xi } } \ \ {\varphi _{,\eta } } \qquad {w_{,\xi } } \ \ {w_{,\eta } } \ \ {w_{,\xi \xi } } \ \ {w_{,\eta \eta } } \ \ {w_{,\xi \eta} } ]^{\rm T} (12)$$
矩阵${\pmb D}\left( {11\times 15} \right)$的非零元素为
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_4.png)
式中所涉及的对物理域的各阶导数可利用式(9)并结合链式法则显式计算.
利用Gauss-Lobatto 数值积分[30]计算单元应变能式(5)
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_5.png)
其中,$\bar {\varOmega }$为标准域内板单元中性面区域. ${Nl}$为板单 $\xi $ 向积分点数,$Nm$为 $\eta $ 向积分点数,变量右下标$ij$表示该变量在该积分点处的值,$H_{i}$, $H_{j}$为相应积分点对应的积分权系数,${\pmb J}$为线性变换(9)的 Jacobi 矩阵.
利用微分求积法则和广义微分求积法则[31],将式(14)中$\bar{\varepsilon }_{ij}$所含积分点处的导数表示为积分点处基本自由度的线性加权代数和. 轴向位移$u$, $v$和转动$\phi $, $\varphi $的导数采用微分求积法则离散,如
$$\left( {u_{,\xi } } \right)_{ij} = \sum_{k = 1}^{Nl} {C_{ik}^{\left( 1 \right)} u_{kj} } \,, \ \ \left( {u_{,\eta } } \right)_{ij} = \sum_{k = 1}^{Nm} {C_{jk}^{\left( 1 \right)} u_{ik} } (15)$$
考虑到问题$C^{1}$连续的特点,横向位移对$\xi $向和$\eta $向的导数采用广义微分求积法则离散,对$\xi $向和$\eta $向二阶混合导数采用广义微分求积法结合传统微分求积法计算
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_6.png)
式(15)和式(16)中,$C_{ij}^{(m)} $, $G_{ij}^{(m)}$分别为$m$阶微分求积及广义微分求积系数.
上述过程可以简洁的表示为
$$ \bar {\pmb \varepsilon }_{ij} = {\pmb B}_{ij} \bar{\pmb d}^{(e)} (17) $$
式中,${\pmb B}_{ij} $是由微分求积及广义微分求积系数构成的矩阵. $\bar{\pmb d}^{(e)}$为标准域中单元自由度向量,每个积分点均设置$u,v,\phi ,\varphi,w$自由度,标准域角点增设$w_{,\xi } ,w_{,\eta } $自由度,除角点外的边点增设$w_{,\xi } (\xi = \pm 1) $或$w_{,\eta } (\eta = \pm 1)$自由度,如图2所示. 所有自由度经合理排列后构成标准域自由度向量$\bar {\pmb d}^{(e)}$.
为便于单元拼接和施加边界条件,物理域单元自由度向量${\pmb d}^{(e)}$按如下方式设置:每个积分点均设置$u,v,\phi ,\varphi ,w$自由度,物理域角点增设$w_{,x} ,w_{,y}$自由度,除角点外的边点增设$w_{,n} $自由度,${\pmb n}$为预先设定的单元边界法线方向. 角点处自由度$w_{,\xi } ,w_{,\eta } $可利用线性变换(9)的Jacobi矩阵与$w_{,x} ,w_{,y} $相联系,而边点自由度$w_{,\xi } (\eta = \pm 1) $或$w_{,\eta } (\xi = \pm 1)$可利用方向导数的计算方式与$w_{,n} $相联系,即
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_7.png)
式 中,$n_{x}$, $n_{y}$分别为单元边界法线关于$X$轴、$Y$轴的方向余弦,等式右端$w_{,\xi } $, $w_{,\eta } $的计算也需运用微分求积法则. 最终$\bar {\pmb d}^{(e)}$可表示为
$$\bar{\pmb d}^{(e)} = {\pmb T}{\pmb d}^{(e)} (19)$$
利用式(17),则式(14)可进一步表示为
$ U^{(e)} = \dfrac{1}{2}\bar {\pmb d}^{(e) \rm T}\left( {\sum_{i = 1}^{Nl} \sum_{j = 1}^{Nm} H_i H_j {\pmb B}_{ij}^{\rm T} {\pmb D}_{ij}^{\rm T} {\pmb E}_{ij} {\pmb D}_{ij} {\pmb B}_{ij} \left| {\pmb J} \right|_{ij} } \right)\bar{\pmb d}^{(e)} =$
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_8.png)
类似地,单元的外力势能$V^{(e)} $、动能$T^{(e)} $均可采用相同的步骤进行离散.在单元水平运用Hamilton原理
$$\delta \int_{t_1 }^{t_2 } {\left( {T^{(e)} - U^{(e)} - V^{(e)}} \right)} d t = 0 (21)$$
并将$U^{(e)} $, $V^{(e)} $, $T^{(e)} $代入,可得
$$\bar{\pmb M}^{\,(e)}{\mathop{\bar{\pmb d}^{\,(e)}}\limits^{\cdot\cdot \ \ \ }}+\bar{\pmb K}^{\,(e)}\bar{d}^{\,(e)} =\bar{\pmb F}^{\,(e)} (22)$$
相应地,物理域中的单元运动方程
$${\pmb M}^{\,(e)}\ddot {\pmb d}^{\,(e)} + {\pmb K}^{\,(e)} {\pmb d}^{\,(e)} ={\pmb F}^{\,(e)} (23)$$
可利用式(19)得到. 式(23)中,$\ddot{\pmb d}^{(e)}$表示单元自由度向量对时间的二阶导数,${\pmb M}^{(e)}$, ${\pmb K}^{(e)} $和${\pmb F}^{(e)}$分别为单元质量矩阵、刚度矩阵以及载荷向量.
如有必要,可进行单元拼接,得到总体刚度矩阵. 对于弯曲问题,式(23)转化为线性代数方程组. 对于自由振动问题,式(23)转化为广义代数特征值问题.
1.3 边界条件
求积元板单元的建立过程基于问题的弱形式描述,因而只需考虑位移边界条件. 后文中提及的边界条件定义如下. 简支边界$$w = 0\,, \ \ u_s = 0\,, \ \ \phi _s = 0 (24)$$
固支边界
$$w = w_n = u_s = u_n = \phi _s = \phi _n = 0 (25)$$
自由边界则不约束任何自由度. 式(24)和式(25)中,$u_s $和$u_n$分别表示沿单元边线切向和法向的位移,$\phi _s $和$\phi _n$分别表示中性面法线绕单元边线法向和切向的转角. 单元边线的切向${\pmb s}$和法向${\pmb n}$应预先合理定义. 具体计算中,针对不同的边界条件利用下述关系
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_9.png)
将${\pmb d}^{(e)}$中相应自由度进行转换,方可准确施加边界条件. 式(26)中,$s_{x}$, $s_{y}$分别为单元边界切线关于$X$轴、$Y$轴的方向余弦.
2 数值算例
通过对典型均质材料、复合材料、功能梯度材料构成的等/变厚度矩形板及等厚度斜板进行受荷弯曲和自由振动分析,并与文献中相关结果进行比较,考察求积元板单元的计算效率.2.1 等厚度矩形板的弯曲和自由振动分析
(1)均质材料板考察一等厚各向同性均质材料矩形板(边长分别为$a$,$b$;泊松比$\nu =0.3$)受均布载荷$q_{0}$的弯曲和自由振动问题,将整块板划分为1个求积元单元进行计算.
表1和表2分别给出了四边简支条件下,不同高跨比($h/a$)时,方板($a=b$)板中的无量纲挠度与弯矩.为便于比较,表中同时给出 了基于相同力学模型的有限元计算结果[15]与Reddy给出的Navier解[32]. 表中 "QEM 7$\times $7''表示,采用1个求积元板单元,积分点数${Nl}= {Nm}=7$,下同.
Table 1
表1
表1方板板中无量纲挠度$100{wEh^3} /{q_0 a^4(1 - \nu ^2)}$
Table 1The dimensionless deflection at the center of the plate
A/a | 0.001 | 0.01 | 0.10 | 0.20 | 0.25 |
---|---|---|---|---|---|
QEM 7x7 | 0.406 4 | 0.406 6 | 0.4274 | 0.490 3 | 0.537 3 |
QEM 9x9 | 0.406 2 | 0.406 5 | 0.4273 | 0.490 3 | 0.537 5 |
QEM 1x11 | 0.406 2 | 0.406 4 | 0.4273 | 0.490 3 | 0.537 3 |
QEM3x13 | 0.406 2 | 0.406 4 | 0.4273 | 0.490 3 | 0.537 4 |
Sheikh[15] | 0.406 4 | 0.406 6 | 0.4274 | 0.490 5 | 0.537 5 |
Reddy[32] | 0.406 2 | 0.406 4 | 0.4273 | 0.490 3 | 0.537 4 |
新窗口打开
Table 2
表2
表2方板板中无量纲弯矩$10{M_x } / {q_0 a^2}$
Table 2The dimensionless moment at the center of the plate
h/a | 0.001 | 0.01 | 0.10 | 0.20 | 0.25 |
---|---|---|---|---|---|
QEM 7x7 | 0.4785 | 0.478 5 | 0.478 4 | 0.478 2 | 0.478 2 |
QEM 9x9 | 0.4789 | 0.478 9 | 0.478 8 | 0.478 9 | 0.478 8 |
QEM 11x11 | 0.4789 | 0.478 9 | 0.478 9 | 0.478 8 | 0.478 8 |
QEM 13x13 | 0.4789 | 0.478 9 | 0.478 9 | 0.478 9 | 0.478 9 |
Sheikh[15] | 0.478 8 | 0.478 8 | 0.478 8 | 0.478 8 | 0.478 8 |
Reddy[32] | 0.4789 | 0.478 9 | 0.478 9 | 0.478 9 | 0.478 9 |
新窗口打开
观察表1和表2可知,随着积分点数的增加,位移和弯矩的求积元解均能较快地收敛于解析解,随着板厚的增大,收敛速度略有降低,采用9个积分点即能得到满意的计算结果. Sheikh等[15]利用512个的6节点(每节点7自由度)三角形有限元单元未能得到与解析解一致的计算结果.
表3给出了两对边简支($x=0$,$x=a$),$y=0$边固支,$y=b$边自由的边界条件下,不同高跨比($h/a$)和长宽比($b$/$a$)时板中的无量纲挠度$10{wEh^3}/ {q_0 a^4(1 - \nu ^2)}$、自由边板中弯矩${M_x } / {q_0 a^2}$、固支边板中弯矩${M_y }/{q_0 a^2}$.观察表3可知,对于挠度,求积元计算结果与有限元较为一致;而对于弯矩,求积元的计算结果略高于有限元. 原因除单元理论不同外,两种方法对于应力(内力)的恢复计算(后处理)方式也存在差异.求积元利用高阶插值的方式得到单元域内任意点应变,进而通过本构关系计算应力;而传统有限元在利用形函数计算单元应变及应力后,也会采取平滑化处理技术改善应力分布.
Table 3
表3
表3矩形板无量纲挠度与弯矩
Table 3The dimensionless deflection and moment of the rectangular plate
b/a | h/a | Method | $\bar{w}$ | $\bar{ Mx }$ | $\bar{ My }$ |
---|---|---|---|---|---|
1.5 | 0.01 | QEM 9x9 | 0.142 6 | 0.123 3 | 0.1236 |
1.5 | 0.01 | QEM 11x11 | 0.141 6 | 0.123 3 | 0.1237 |
1.5 | 0.01 | QEM 13x13 | 0.141 6 | 0.123 2 | 0.1237 |
1.5 | 0.01 | Sheikh[15] | 0.141 7 | 0.123 1 | 0.1177 |
1.5 | 0.1 | QEM 9x9 | 0.1478 | 0.1202 | 0.1188 |
1.5 | 0.1 | QEM 11x11 | 0.148 1 | 0.1198 | 0.1187 |
1.5 | 0.1 | QEM 13x13 | 0.148 1 | 0.1198 | 0.1187 |
1.5 | 0.1 | Sheikh[15] | 0.148 1 | 0.1170 | 0.115 8 |
2.0 | 0.01 | QEM 9x9 | 0.1496 | 0.1305 | 0.1246 |
2.0 | 0.01 | QEM 11x11 | 0.1496 | 0.1305 | 0.1246 |
2.0 | 0.01 | QEM 13x13 | 0.1496 | 0.1304 | 0.1246 |
2.0 | 0.01 | Sheikh[15] | 0.1500 | 0.131 4 | 0.115 2 |
2.0 | 0.1 | QEM 9x9 | 0.155 2 | 0.1279 | 0.1201 |
2.0 | 0.1 | QEM 11x11 | 0.155 6 | 0.1273 | 0.1198 |
2.0 | 0.1 | QEM 13x13 | 0.155 8 | 0.1271 | 0.1196 |
2.0 | 0.1 | Sheikh[15] | 0.155 8 | 0.1262 | 0.1136 |
新窗口打开
表4给出了四边简支方板,高跨比$h/a=0.1$时的横向无量纲自振频率.为便于比较,表中同时列出了基于相同力学模型的有限元(64个9节点矩形单元,每节点7个自由度)计算结果[13]、Reddy给出的Navier解[34]以及Shufrin给出的扩展Kantorovich法解[33].模态栏"1,1''表示$X$向,$Y$向的半波数,下同.
Table 4
表4
表4方板无量纲自振频率$\omega h\sqrt {{2\rho (1 + \nu )} / E} $
Table 4The dimensionless frequencies of the square plate
Mode | QEM 9x9 | QEM 10 x 10 | QEM 11 x 11 | Reddy [34] | Shufrin[33] Nayak [13] | |
---|---|---|---|---|---|---|
1,1 | 0.093 0 | 0.0930 | 0.093 0 | 0.093 1 | 0.0930 | 0.0930 |
1,2 | 0.222 0 | 0.2220 | 0.2220 | 0.2222 | 0.2220 | 0.2222 |
2,2 | 0.3406 | 0.3406 | 0.3406 | 0.341 1 | 0.3406 | 0.341 2 |
1,3 | 0.4090 | 0.4151 | 0.4151 | 0.4158 | 0.4151 | 0.4169 |
2,3 | 0.511 2 | 0.5208 | 0.5208 | 0.5221 | 0.5208 | 0.5231 |
1,4 | 0.569 4 | 0.6455 | 0.6526 | 0.6545 | — | 0.6599 |
新窗口打开
观察表4可知,求积元单元的计算结果与现有文献吻合较好,与文献[33]计算结果一致.随着积分点数的增加,各阶自振频率均能 较快地收敛,对于基频的计算,采用9个积分点即能得到满意的计算结果.
(2)复合材料层合板
考虑一四边简支的正方形复合材料层合板,其边长为$a$. 铺层层数为4,铺层方式为(0$^\circ$/90$^\circ$)$_{\rm s}$,各层材料均为相同的正交各向异性材料且层厚相等,材料常数关系为
$$ \left. {E_1 }/{E_2 } = 25\,, \ \ {G_{12} } / {E_2 } = 0.5\,, \ \ {G_{23} }/ {E_2 } = 0.2 \nu _{12} = 0.25\,, \ \ G_{12} = G_{13}\,, \ \ \rho = 1 \right\} (27) $$
该板承受正弦分布载荷,$q = q_0 \sin \left( {{\pi x} / a} \right)\sin\left( {{\pi y}/ a} \right)$.
考虑到问题的对称性以及与现有文献进行直接对比,分别用1/4几何模型和全模型进行计算,两种模型均划分为1个求积元单元.表5给出了不同高跨比时,板中挠度及相应位置处应力分量的无量纲计算值. 无量纲量计算方式为
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/img_10.png)
为便于比较,表中同时列出了基于相同力学模型的有限元协调元[6]和非协调元[7]计算结果以及Reddy给出的Navier解[32],其中协调元为16个4节点矩形单元,每节点8个自由度;非协调元为25个4节点矩形单元,每节点7个自由度.
Table 5
表5
表5复合材料层合板无量纲挠度与应力
Table 5The dimensionless deflection and stress of the laminated plate
h/a | Method | $\bar{w}$ | $\bar{\sigma_1}$ | $\bar{\sigma_2}$ | $\bar{\sigma_3}$ |
---|---|---|---|---|---|
0.1 | 1/4 model QEM 7 x 7 | 0.7156 | 0.5470 | 0.3893 | 0.0268 |
0.1 | 1/4 model QEM 9 x 9 | 0.7147 | 0.545 6 | 0.388 8 | 0.0268 |
0.1 | 1/4 model QEM 11 x 11 | 0.7147 | 0.545 6 | 0.388 8 | 0.0268 |
0.1 | full model QEM 7 x 7 | 0.7147 | 0.545 5 | 0.3888 | 0.0268 |
0.1 | full model QEM 9 x 9 | 0.7147 | 0.545 6 | 0.3888 | 0.0268 |
0.1 | full model QEM 11 x 11 | 0.7147 | 0.545 6 | 0.3888 | 0.0268 |
0.1 | Reddy[32] | 0.7147 | 0.545 6 | 0.3888 | 0.0268 |
0.1 | 1/4 model Phan[7] | 0.7161 | 0.5427 | 0.385 5 | 0.0266 |
0.1 | 1/4 model Liu[6] | 0.7176 | 0.541 3 | 0.3873 | 0.0266 |
0.01 | 1/4 model QEM 7 x 7 | 0.4343 | 0.5387 | 0.2708 | 0.021 3 |
0.01 | full model QEM 7 x 7 | 0.4349 | 0.5401 | 0.271 1 | 0.021 5 |
0.01 | full model QEM 9 x 9 | 0.4343 | 0.5387 | 0.2708 | 0.021 3 |
0.01 | full model QEM 11 x 11 | 0.4343 | 0.5387 | 0.2708 | 0.021 3 |
0.01 | Reddy[32] | 0.4343 | 0.5387 | 0.2708 | 0.021 3 |
0.01 | 1/4 model Phan[7] | 0.4320 | 0.5301 | 0.2664 | 0.021 1 |
0.01 | 1/4 model Liu[6] | 0.4408 | 0.5349 | 0.2670 | 0.021 1 |
新窗口打开
观察表5可知,随着积分点的增加,挠度和应力的求积元的计算结果均能较快地收敛于解析解[32],采用9个积分点即能得到满意的计算结果. 14几何模型和全模型没有显著差异,下文中均采用全模型进行计算. 同时,两类有限元的计算结果未能收敛至解析解.
将该板的铺层方式变为$(45^\circ/-45^\circ )_{4}$,层数为8,各层材料层厚相等. 表6给出了不同高跨比时,板中挠度及相应位置处应力分量的无量纲计算值,具体计算方式为
$$\left. \bar {w} = 100w\left( {a /2,a / 2} \right){E_2 h^3}/{q_0 }a^4 \bar {\sigma }_{xz} = \sigma _{xz} \left( {0,a /2,z} \right)h/{q_0 }a \right\} (29)$$
表中同时列出基于相同力学模型的有限元协调元[5]和$C^{0}$单元[9]计算结果,其中协调元为4个4节点矩形单元,每节点14个自由度;$C^{0}$单元为16个9节点矩形单元,每节点7个自由度.观察表6可知,对于挠度而言,求积元的计算结果与协调元吻合度较好;而对于应力,求积元的计算结果略低于协调元.
Table 6
表6
表6复合材料层合板无量纲挠度与应力
Table 6The dimensionless deflection and stress of the laminated plate
![]() |
新窗口打开
考察一四边简支矩形复合材料软芯三明治层合板的自由振动问题,其层数为7,边长分别为$a$和$b$. 铺层方式为0$^\circ$90$^\circ$0$^\circ$/芯层0$^\circ$90$^\circ$0$^\circ$,除芯层外各层材料均为相同的正交各向异性聚酯树脂材料,材料常数为:$E_{1}=24.51$,GPa,$E_{2}=7.77$,GPa,$G_{12}=G_{13}=3.34$,GPa,$G_{23}=1.34$,GPa, $\nu_{12}=0.078$, $\rho =1,800$,kg/m$^{3}$. 芯层材料为聚氯乙烯泡沫,$E_{\rm c}=103.63$,MPa,$G_{\rm c}=50$,MPa, $\nu _{\rm c}=0.32$, $\rho_{\rm c}=130$,kg/m$^{3}$. 芯层厚度$h_{\rm c}=0.88h$,其余各层厚度相同.
表7给出了不同高跨比、长宽比时,前三阶无量纲自振频率的计算值. 表中同时列出了基于相同力学模型的有限元$C^{0}$单元[13]和阶谱单元的[18]计算结果,其中$C^{0}$单元为36个9节点矩形单元,每节点7个自由度;阶谱单元为1个矩形单元,阶谱项数为14. 观察表7可知,求积元方法的计算结果与阶谱单元近乎一致,而$C^{0}$单元的计算结果与前二者差异较大,最大相对差异为5.6%.
Table 7
表7
表7矩形板无量纲自振频率$\omega \sqrt {{\rho _{\rm c} }/ {E_{\rm c} }} a{^2} / h$
Table 7The dimensionless frequencies of the rectangular laminated plate
Mode | h/a = 0.2 | h/a = 0.05 | |||||
---|---|---|---|---|---|---|---|
a/b = 0.5 | a/b = 1.0 | a/b = 1.5 | a/b = 0.5 | a/b = 1.0 | a/b = 1.5 | ||
QEM | 1 | 7.50 | 9.71 | 12.83 | 13.94 | 19.36 | 29.22 |
2 | 9.71 | 16.31 | 18.80 | 19.36 | 42.11 | 51.86 | |
13 X 13 | 3 | 10.21 | 16.81 | 20.41 | 29.22 | 45.52 | 72.21 |
1 | 7.50 | 9.71 | 12.83 | 13.94 | 19.36 | 29.22 | |
Ref.[18] | 2 | 9.72 | 16.31 | 18.80 | 19.36 | 42.11 | 51.86 |
3 | 10.21 | 16.81 | 20.41 | 29.22 | 45.52 | 72.21 | |
1 | 7.29 | 9.42 | 12.36 | 13.85 | 19.23 | 28.97 | |
Ref.[13] | 2 | 9.43 | 15.56 | 17.74 | 19.24 | 41.70 | 51.12 |
3 | 10.20 | 15.96 | 20.41 | 29.16 | 44.88 | 71.17 |
新窗口打开
2.2 变厚度矩形板的自由振动
考察一各向同性均质材料正方形板(边长为$a$, 泊松比 $\nu =0.3$)的自由振动问题,该板厚度沿$X$方向以四类方式(式(30))单向变化,如图3所示.$$\left. \begin{array}{l} {\rm I}: h\left( x \right) = h_0 \left( {1 - x /{2a}} \right) {\rm II}:h\left( x \right) = h_0 \left( {1 - {x^2}/ {2a^2}} \right) {\rm III}: h\left( x \right) = h_0 \left({1 + {x^2}/ {2a^2} - x / a} \right) {\rm IV}: h\left( x \right) = h_0 \left( {1 + 2{x^2} /{a^2} - 2x /a} \right) \end{array}\!\!\right\} (30)$$
表8给出了四边简支方形板($h_{0}$/$a$=0.2)在I类厚度变化条件下的前6阶横向无量纲自振频率. 表8同时列出了Shufrin等基于一阶剪切变形理论(FSDT)[35]以及Bacciocchi等基于一般高阶剪切变形理论(HSDT)[36]得出的计算结果. 观察表8可知,求积元单元仅需9个积分点即可使前四阶自振频率收敛,其计算结果高于一阶剪切变形理论,低于一般高阶剪切变形理论.
![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/thumbnail/img_12.png)
图3板厚沿$X$方向变化的4种类型
-->Fig.34 different thickness profiles of the plate along $X$ direction
-->
Table 8
表8
表8变厚度方板无量纲自振频率$\omega {a^2}/{\pi ^2}\sqrt {{12\rho \left({1 - \nu ^2} \right)}/{Eh_0^2 }} $
Table 8The dimensionless frequencies of the square plate with variable thickness (Type I)
Mode | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
QEM | 1.374 9 | 3.1150 | 3.1331 | 4.6741 | 5.451 5 | 5.5566 |
9x9 | ||||||
QEM | 1.374 9 | 3.1150 | 3.1331 | 4.6741 | 5.5050 | 5.5829 |
10 x 10 | ||||||
QEM | 1.374 9 | 3.1150 | 3.1331 | 4.6741 | 5.5047 | 5.5839 |
11 x 11 | ||||||
Ref.[35] | 1.373 8 | 3.1096 | 3.1276 | 4.661 3 | 5.4883 | 5.5657 |
Ref.[36] | 1.378 7 | 3.1324 | 3.1500 | 4.6923 | 5.5509 | 5.631 0 |
新窗口打开
表9给出了四边固支方形板($h_{0}/a =0.1$)在II、III、IV类厚度变化条件下的前6阶横向无量纲自振频率. 同时也给出了基于相同力学模型的扩展Kantorovich法解[35]以及基于一阶剪切变形理论的广义微分求积法解[36]. 观察表9可知,求积元法与文献[35]计算结果吻合较好,最大相对误差0.2%,二者计算结果高于一阶剪切变形理论[36].
Table 9
表9
表9变厚度方板无量纲自振频率$\omega {a^2}/{\pi ^2}\sqrt {{12\rho \left({1 - \nu ^2} \right)}/{Eh_0^2 }} $
Table 9The dimensionless frequencies of the square plate with variable thickness (Type II, III, IV)
Mode | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
Type II | ||||||
QEM 13x13 | 2.7450 | 5.3082 | 5.421 8 | 7.6561 | 9.0295 | 9.1302 |
Ref.[35] | 2.7454 | 5.3031 | 5.4186 | 7.6483 | 9.0192 | 9.1099 |
Ref.[36] | 2.7392 | 5.2879 | 5.4007 | 7.6143 | 8.9747 | 9.0768 |
Type III | ||||||
QEM 13x13 | 2.2846 | 4.3393 | 4.5317 | 6.4331 | 7.3198 | 7.7467 |
Ref.[35] | 2.2850 | 4.3375 | 4.5294 | 6.4291 | 7.311 8 | 7.7358 |
Ref.[36] | 2.281 8 | 4.3298 | 4.5207 | 6.4107 | 7.2945 | 7.7157 |
Type IV | ||||||
QEM 13x13 | 2.4382 | 4.2238 | 4.800 9 | 6.5089 | 6.9839 | 8.0295 |
Ref.[35] | 2.4384 | 4.2222 | 4.7987 | 6.5044 | 6.9782 | 8.021 2 |
Ref.[36] | 2.4356 | 4.2163 | 4.789 4 | 6.4882 | 6.9645 | 7.9971 |
新窗口打开
2.3 等厚度斜板的弯曲和自由振动分析
平行四边形斜板如图4所示,相邻两边边长为$a$,板厚为$h$,相邻边一边与$X$轴平行,另一边与$Y$轴夹角为$\alpha $.![](https://lxxb.cstam.org.cn/fileLXXB/journal/article/html/0459-1879-50-5-1093/thumbnail/img_13.png)
图4斜板的几何尺寸与坐标系
-->Fig. 4The geometry and configuration of a skew plate
-->
考虑一承受均布载荷$q_{0}$的复合材料层合板,其层数为3,铺层方式为0$^\circ$/90$^\circ$/0$^\circ$,各层材料均为相同的正交各向异性材料,材料常数关系见式(27).仍将整块板划分为一个求积元单元,随着夹角 $\alpha $ 的增大,所需积分点数逐步增加. 当 $\alpha=60^\circ$时,为得到较为理想的计算结果,需使用15个积分点.表10给出了不同高跨比及夹角条件下四边简支板板中挠度及应力的无量纲计算值,具体计算方式为
$$\left. \bar {w} = 100w [a\left( {1 + \sin \alpha } \right) /2,{a\cos \alpha } /2] {E_2 h^3} / (q_0 a^4) \bar {\sigma }_1 = \sigma _{xx} [a\left( {1 + \sin \alpha } \right) /2,{a\cos \alpha } /2,h/2] {h^2} /(q_0 a^2) \bar {\sigma }_2 = \sigma _{xx} [a\left( {1 + \sin \alpha } \right) / 2,{a\cos \alpha } /2,h /6] {h^2}/ (q_0 a^2) \right\} (31)$$
同时,表中给出了Sheikh等[15]利用512个的6节点(每节点7自由度)三角形有限元单元的计算结果. 观察表10可知,对于挠度而言,求积元的计算结果与有限元吻合度较好;而对于应力,求积元的计算结果略高于有限元.
Table 10
表10
表10斜板板中无量纲挠度与应力
Table 10The dimensionless deflection and stress at the center of the skew laminated plate
h/a | α | Method | $\bar{w}$ | $\bar{sigma_{1}}$ | $\bar{sigma_{2}}$ |
---|---|---|---|---|---|
0.01 | 30° | QEM 15x15 | 0.546 5 | 0.6592 | 0.2802 |
0.01 | 30° | Ref.[15] | 0.5452 | 0.6444 | 0.2629 |
0.01 | 45° | QEM 15x15 | 0.3633 | 0.4488 | 0.3197 |
0.01 | 45° | Ref.[15] | 0.3631 | 0.4421 | 0.3007 |
0.01 | 60° | QEM 15x15 | 0.1464 | 0.2025 | 0.2723 |
0.01 | 60° | Ref.[15] | 0.1455 | 0.201 1 | 0.2572 |
0.10 | 30° | QEM 15x15 | 0.861 6 | 0.6749 | 0.3896 |
0.10 | 30° | Ref.[15] | 0.8621 | 0.661 7 | 0.3709 |
0.10 | 45° | QEM 15x15 | 0.5708 | 0.4602 | 0.3980 |
0.10 | 45° | Ref.[15] | 0.5707 | 0.4543 | 0.3786 |
0.10 | 60° | QEM 15x15 | 0.2503 | 0.2059 | 0.3152 |
0.10 | 60° | Ref.[15] | 0.2505 | 0.2058 | 0.3023 |
新窗口打开
考虑一四边简支的功能梯度碳纳米管增强复合材料板的自由振动问题. ${\pmb E}_{\rm p}$和${\pmb E}_{\rm t}$按照文献[37]介绍的扩展混合律模型进行计算,基体材料参数为:$E_{\rm m}=2.1$\,GPa,$\rho_{\rm m}=1\,150$\,kg/m$^{3}$,$\nu_{\rm m}=0.34$. 碳纳米管材料参数为:$E_{11}=5.646\,6$\,TPa,$E_{22}=7.080\,0$\,TPa,$G_{12}=1.944\,5$\,TPa,$\nu_{12}=0.175$,$\rho =1\,400$\,kg/m$^{3}$. 表11和表12分别给出了高跨比$h/a=0.1$,碳纳米管平均体积含量为0.11时,碳纳米管沿厚度均匀分布(UD)与O型分布(FG-O)时的前6阶无量纲自振频率. 表中同时给出了Ansari等[12]利用256个9节点(每节点7自由度)矩形有限元单元的计算结果.
Table 11
表11
表11斜板无量纲自振频率(UD) $\omega {a^2} / h\sqrt {{\rho _m }/{E_m }} $
Table 11The dimensionless frequencies of the FG-CNTRC skew plate (UD)
Mode | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
α | =15° | |||||
QEM 15 x 15 | 13.746 | 18.278 | 18.327 | 22.576 | 28.358 | 33.426 |
Ref.[12] | 14.088 | 18.320 | 18.895 | 22.618 | 29.333 | 33.641 |
α | = 30° | |||||
QEM 15 x 15 | 14.519 | 18.668 | 20.719 | 28.286 | 32.165 | 34.834 |
Ref.[12] | 14.898 | 18.830 | 21.372 | 28.432 | 33.108 | 35.118 |
α | = 45° | |||||
QEM 15 x 15 | 16.707 | 20.197 | 26.207 | 36.689 | 38.864 | 39.160 |
Ref.[12] | 17.186 | 20.512 | 26.990 | 37.014 | 39.632 | 39.657 |
新窗口打开
Table 12
表12
表12斜板无量纲自振频率(FG-O)$\omega {a^2} / h\sqrt {{\rho _m }/ {E_m }} $
Table 12The dimensionless frequencies of the FG-CNTRC skew plate (FG-O)
Mode | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
α | = 15° | |||||
QEM 15 x 15 | 11.548 | 16.756 | 18.326 | 22.635 | 27.174 | 29.610 |
Ref.[12] | 11.225 | 16.857 | 18.273 | 22.561 | 27.702 | 28.690 |
α | = 30° | |||||
QEM 15 x 15 | 12.379 | 18.718 | 19.190 | 28.361 | 30.548 | 31.116 |
Ref.[12] | 12.113 | 18.775 | 19.364 | 28.356 | 30.314 | 30.780 |
α | = 45° | |||||
QEM 15 x 15 | 14.661 | 20.253 | 24.547 | 35.709 | 36.336 | 36.789 |
Ref.[12] | 14.520 | 20.439 | 24.733 | 35.270 | 36.212 | 36.901 |
新窗口打开
观察表11和表12可知,求积元的计算结果与有限元吻合较好,最大相对差异小于3.5%.由于文献[12]中提供的信息有限,无 法保证两类方法计算过程中除单元外的其他步骤完全一致,如对式(7)中沿厚度积分的具体计算方式,这些因素也将导致计算结果的差异.
3 结 论
本文基于Reddy高阶剪切变形理论,建立了$C^{1}$连续的四边形求积元板单元. 利用该单元对典型均质材料、复合材料、功能梯度材料构成的等厚度矩形板、变厚度矩形板及等厚度斜板的线弹性弯曲和自由振动问题进行了计算分析,并与现有文献中相应的计算结果进行了对比. 主要结论如下:(1)基于高阶剪切变形理论的四边形求积元板单元具有较高的计算效率和良好的适应性,对于本文中涉及的各类材料等变厚度矩形板,等厚度斜板均只需一个求积元单元即可得到理想的计算结果.
(2)对于等/变厚度矩形板,使用1个求积元单元9$\times$9个积分点即可得到满意的计算结果. 对于等厚度斜板,随着夹角$\alpha $的增大,所需积分点的数目逐渐增多至15$\times$15.
(3)本文建立的四边形求积元板单元可进一步应用于新型复合材料板的非线性分析.
The authors have declared that no competing interests exist.
参考文献 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子
[1] | |
[2] | . . |
[3] | . . |
[4] | . |
[5] | . |
[6] | . |
[7] | . |
[8] | . |
[9] | . |
[10] | . |
[11] | . |
[12] | . |
[13] | . |
[14] | . |
[15] | . |
[16] | . |
[17] | . |
[18] | . |
[19] | . |
[20] | . . |
[21] | . |
[22] | . . |
[23] | . |
[24] | . |
[25] | |
[26] | |
[27] | . |
[28] | . |
[29] | . . |
[30] | . |
[31] | |
[32] | . |
[33] | . |
[34] | . |
[35] | . |
[36] | . |
[37] | . |