一种复合材料层合结构三维振动分析方法

摘要

本发明涉及的是一种应用于工程力学和振动工程领域的复合材料层合结构三维振动分析方法。本发明包括：提取结构的几何、材料和边界条件参数，对结构每一层沿厚度方向设置J个非均匀分布的计算平面；配置结构在计算平面上的位移；对每个计算平面上的结构位移进行面内分量展开并通过边界特征函数施加边界条件，得到计算平面上结构位移；计算第l层第j计算平面上结构的面内方向应变；求得结构在第l层任意位置的应变和应力；建立结构能量方程；求得结构特征方程。本发明通过空间分层取面把结构分解成多个计算平面，一方面降低结构维度，从而提高计算速度，节约计算成本，另一方面把结构化整为零，便于并行计算，从而提高计算效率。

说明书

技术领域

本发明涉及的是一种应用于工程力学和振动工程领域的复合材料层合结构三维振动分析 方法。

背景技术

复合材料层合结构具有质量轻、比刚度高、比强度大，隔热，隔音和优良的减振降噪性 能而被广泛应用于航空航天、军事装备、和科技建筑等领域。复合材料层合结构动力学分析 一直是很多学者关注和探讨的重点。与常规结构相比，复合材料层合结构的组成材料复杂、 铺层方式多样，因此其动力学行为更为复杂。目前，国内外的绝大部分研究都还是把三维的 复合材料层合结构通过ESL方法简化成一维或者二维的经典各向异性结构来处理。这种处理 降低了研究难度，对比较薄的层合结构来说其计算结果的精度是可以接受。但是这种方法忽 略了结构内部铺层之间在厚度方向上的正应力和剪切应力的不连续性，因而当结构的厚度比 比较高或者不同铺层之间材料属性差异较大时，其计算结果相差甚远。因此研究和建立一种 能够适用任意厚度复合材料层合结构的振动分析方法具有十分重要的意义。

本发明提供了一种基于Chebyshev多项式与空间分层取面的复合材料层合结构三维振动 分析方法。这种方法具有适用任意边界条件和任意厚度、精度高、收敛快、计算成本低、计 算方法简单等特点。

发明内容

本发明的目的在于提供一种用以求解任意厚度复合层合结构在任意边界条件下的三维振 动问题的基于Chebyshev多项式与空间分层取面的一种复合材料层合结构三维振动分析方 法。

本发明的目的是这样实现的：

(1)提取结构的几何、材料和边界条件参数，对结构每一层沿厚度方向设置J个非均匀 分布的计算平面；其中第1个和第J个计算平面分别选取为该层的下表面和上表面；第2到 第J-1个计算平面在厚度方向上的位置选取为移位勒让德多项式P_J-2(z)在区间[0,h_l]的实零点； P_J-2(z)表达式如下：

$\left(\begin{array}{cc}{P}_{J-2}\left(z\right)=\frac{1}{2^{J-2}(J-2)!}\frac{d^{J-2}}{{\mathrm{dz}}^{J-2}}{\left[{\left(\frac{2z}{h_l}-1\right)}^2-1\right]}^{J-2}=0,& J>2\end{array}\right)$

其中，l指的是层合结构的第l层，h_l为该层厚度；

(2)配置结构在计算平面上的位移，即结构在第l层第j个计算平面上的结构位移设置 为其中α，β为结构空间坐标系面内坐标，i＝1,2,3分别为结构位移在α，β和z 方向上的分量；利用拉格朗日插值将结构任意位置的结构位移设定为

$u_{il}(\alpha ,\beta ,z)=\underset{j=1}{\overset{J}{\Sigma }}\underset{k\ne j}{\overset{J}{\Pi }}\frac{z-z_l^k}{z_l^j-z_l^k}{u}_{il}^j(\alpha ,\beta )U_{il}^j$

其中，为步骤(1)中得到的多项式P_J-2(z)在区间[0,h_l]的实零点；

(3)利用Chebyshev多项式对每个计算平面上的结构位移进行面内分量展开并通过边界 特征函数施加边界条件，得到计算平面上结构位移：

$u_{il}^j\left(\alpha ,\beta \right)U_{il}^j=F_i\left(\alpha ,\beta \right)\underset{m=0}{\overset{M}{\Sigma }}\underset{n=0}{\overset{N}{\Sigma }}{U}_{ilmn}^j\cos \left[m\arccos \left(\alpha \right)\right]\cos \left[n\arccos \left(\beta \right)\right];$

其中F_i(α,β)＝(1-α)^a(1+α)^b(1-β)^c(1+β)^d为边界特征函数；a,b,c,d为边界特征函数系数，这些 系数由边界条件确定；M，N为截断级数；

(4)由步骤(2)和(3)计算第l层第j计算平面上结构的面内方向应变和横向应变和横向剪切应变和分别为：

$\left(\begin{array}{cc}{ϵ}_{l\alpha }^j=\frac{\partial u_{1lmn}^j\left(\alpha ,\beta \right)}{\partial \alpha }{U}_{1lmn}^j,& {ϵ}_{l\beta }^j=\frac{\partial u_{2lmn}^j\left(\alpha ,\beta \right)}{\partial \beta }{U}_{2lmn}^j\end{array}\right)$

${\gamma }_{l\alpha \beta }^j=\frac{\partial u_{1lmn}^j(\alpha ,\beta )}{\partial \beta }{U}_{1lmn}^j+\frac{\partial u_{2lmn}^j(\alpha ,\beta )}{\partial \alpha }{U}_{2lmn}^j$

${ϵ}_{lz}^j=u_{3lmn}^j(\alpha ,\beta )\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{3lmn}^j$

${\gamma }_{l\alpha z}^j=u_{1lmn}^j\left(\alpha ,\beta \right)\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{1lmn}^j+\frac{\partial u_{3lmn}^j\left(\alpha ,\beta \right)}{\partial \alpha }{U}_{3lmn}^j$

${\gamma }_{l\beta z}^j=u_{2lmn}^j\left(\alpha ,\beta \right)\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{2lmn}^j+\frac{\partial u_{3lmn}^j\left(\alpha ,\beta \right)}{\partial \beta }{U}_{3lmn}^j$

(5)由步骤(2)和(4)求得结构在第l层任意位置的应变和应力表达式为

$\left[{ϵ}_{l\alpha },{ϵ}_{l\beta },{ϵ}_{lz},{ϵ}_{l\alpha \beta },{ϵ}_{l\alpha z},{ϵ}_{l\beta z}\right]=\underset{j=1}{\overset{J}{\Sigma }}\underset{k\ne j}{\overset{J}{\Pi }}\frac{z-z_l^k}{z_l^j-z_l^k}\left[{ϵ}_{l\alpha }^j,{ϵ}_{l\beta }^j,{ϵ}_{lz}^j,{ϵ}_{l\alpha \beta }^j,{ϵ}_{l\alpha z}^j,{ϵ}_{l\beta z}^j\right]$

σ_l＝Cε_l；σ_l＝[σ_lα,σ_lβ,σ_lz,τ_lαβ,τ_lαz,τ_lβz]^T；ε_l＝[ε_lα,ε_lβ,ε_lz,γ_lαβ,γ_lαz,γ_lβz]^T

其中C为结构材料系数矩阵；

(6)根据步骤(5)建立结构能量方程(U，T)：

$\left(\begin{array}{cc}U=\frac{1}{2}\underset{l}{\Sigma }\int \int \int {ϵ}_l^T{C}_l{ϵ}_{l}dzd\beta d\alpha ,& T=\frac{1}{2}\underset{l}{\Sigma }\int \int \int {\rho }_l\left\{\frac{{\partial }^2{u}_{1l}}{\partial t^2}+\frac{{\partial }^2{u}_{2l}}{\partial t^2}+\frac{{\partial }^2{u}_{3l}}{\partial t^2}\right\}dzd\beta d\alpha \end{array}\right)$

(7)在步骤(6)基础上建立结构拉格朗日能量泛函L＝U-T，然后利用Ritz法求得结构 特征方程：

(K-ω²M)＝0；

其中ω为圆频率；

(8)最后应用Arnoldi算法建立MATLAB求解器输出结构的振动特征数据并判定计算精 度，若满足精度要求则输出振动特征数据的固有频率、模态，不满足则继续优化空间分层数 量和增加面内位移展开式级数截取量。

本发明的优势在于：通过空间分层取面把结构分解成多个计算平面，一方面降低结构维 度，从而提高计算速度，节约计算成本，另一方面把结构化整为零，便于并行计算，从而提 高计算效率。

附图说明

图1是本发明的流程图；

图2是复合材料层合矩形板结构及其空间分层取面示意图。

具体实施方式

下面结合附图对本发明做进一步描述。

本发明包括：

(1)在结构的每一层沿厚度方向均设置了J个非均匀分布的计算平面，并且第1个和第 J个计算平面分别选取为该层的下表面和上表面。第2到第J-1个计算平面在厚度方向上的位 置为移位勒让德多项式P_J-2(z)在区间[0,h_l]的实零点，即：

$\left(\begin{array}{cc}{P}_{J-2}\left(z\right)=\frac{1}{2^{J-2}(J-2)!}\frac{d^{J-2}}{{\mathrm{dz}}^{J-2}}{\left[{\left(\frac{2z}{h_l}-1\right)}^2-1\right]}^{J-2}=0,& J>2\end{array}\right)$

在区间[0,h_l]的解。其中，l指的是层合结构的第l层，h_l为该层厚度。

(2)每个计算平面上的结构位移在面内方向展开成Chebyshev多项式和边界特征函数乘 积的形式。

(3)每个计算平面上结构的面内方向应变ε_lα，ε_lβ和γ_lαβ由所在计算平面上设置的位移 直接求偏导得到，而横向应变和横向剪切应变ε_lz，γ_lαz和γ_lβz则由所在层所有计算平面上的位 移的偏导数加权叠加而成。其中α，β和z为结构空间坐标系面内坐标，

(4)结构相邻层之间共用一个计算平面来满足层间的连续性条件。

(5)边界条件由边界特征函数F_i(α,β)＝(1-α)^a(1+α)^b(1-β)^c(1+β)^d表述。其中a,b,c,d为边 界特征函数系数，这些系数由边界条件确定。i＝1,2,3分别为边界特征函数在α，β和z方向 位移上的分量。

下面结合图2，以计算下述层合板边界条件为四边自由(F-F-F-F)时的固有频率(Hz) 为实例，进行方法说明。

(1)提取结构的几何、材料和边界条件参数，对结构每一层沿厚度方向设置J个非均匀 分布的计算平面。其中第1个和第J个计算平面分别选取为该层的下表面和上表面。第2到 第J-1个计算平面在厚度方向上的位置选取为移位勒让德多项式P_J-2(z)在区间[0,h_l]的实零点。 P_J-2(z)表达式如下：

$\left(\begin{array}{cc}{P}_{J-2}\left(z\right)=\frac{1}{2^{J-2}(J-2)!}\frac{d^{J-2}}{{\mathrm{dz}}^{J-2}}{\left[{\left(\frac{2z}{h_l}-1\right)}^2-1\right]}^{J-2}=0,& J>2\end{array}\right)$

其中，l指的是层合结构的第l层，h_l为该层厚度。

(2)配置结构在计算平面上的位移，即结构在第l层第j个计算平面上的结构位移设置 为其中α，β为结构空间坐标系面内坐标，i＝1,2,3分别为结构位移在α，β和z 方向上的分量。同时，利用拉格朗日插值将结构任意位置的结构位移设定为如下形式：

$u_{il}(\alpha ,\beta ,z)=\underset{j=1}{\overset{J}{\Sigma }}\underset{k\ne j}{\overset{J}{\Pi }}\frac{z-z_l^k}{z_l^j-z_l^k}{u}_{il}^j(\alpha ,\beta )U_{il}^j$

其中，为步骤(1)中得到的多项式P_J-2(z)在区间[0,h_l]的实零点。

(3)利用Chebyshev多项式对每个计算平面上的结构位移进行面内分量展开并通过边界 特征函数施加边界条件，得到计算平面上结构位移的详细表达式如下所示：

$u_{il}^j\left(\alpha ,\beta \right)U_{il}^j=F_i\left(\alpha ,\beta \right)\underset{m=0}{\overset{M}{\Sigma }}\underset{n=0}{\overset{N}{\Sigma }}{U}_{ilmn}^j\cos \left[m\arccos \left(\alpha \right)\right]\cos \left[n\arccos \left(\beta \right)\right]$

其中F_i(α,β)＝(1-α)^a(1+α)^b(1-β)^c(1+β)^d为边界特征函数。a,b,c,d为边界特征函数系数，这些 系数由边界条件确定。M，N为截断级数。

(4)由步骤(2)和(3)计算第l层第j计算平面上结构的面内方向应变和横向应变和横向剪切应变和分别为：

$\left(\begin{array}{cc}{ϵ}_{l\alpha }^j=\frac{\partial u_{1lmn}^j\left(\alpha ,\beta \right)}{\partial \alpha }{U}_{1lmn}^j,& {ϵ}_{l\beta }^j=\frac{\partial u_{2lmn}^j\left(\alpha ,\beta \right)}{\partial \beta }{U}_{2lmn}^j\end{array}\right)$

${\gamma }_{l\alpha \beta }^j=\frac{\partial u_{1lmn}^j\left(\alpha ,\beta \right)}{\partial \beta }{U}_{1lmn}^j+\frac{\partial u_{2lmn}^j\left(\alpha ,\beta \right)}{\partial \beta }{U}_{2lmn}^j$

${ϵ}_{lz}^j=u_{3lmn}^j(\alpha ,\beta )\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{3lmn}^j$

${\gamma }_{l\alpha z}^j=u_{1lmn}^j\left(\alpha ,\beta \right)\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{1lmn}^j+\frac{\partial u_{3lmn}^j\left(\alpha ,\beta \right)}{\partial \alpha }{U}_{3lmn}^j$

${\gamma }_{l\beta z}^j=u_{2lmn}^j\left(\alpha ,\beta \right)\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{2lmn}^j+\frac{\partial u_{3lmn}^j\left(\alpha ,\beta \right)}{\partial \beta }{U}_{3lmn}^j$

(5)由步骤(2)和(4)求得结构在第l层任意位置的应变和应力表达式为

$\left[{ϵ}_{l\alpha },{ϵ}_{l\beta },{ϵ}_{lz},{ϵ}_{l\alpha \beta },{ϵ}_{l\alpha z},{ϵ}_{l\beta z}\right]=\underset{j=1}{\overset{J}{\Sigma }}\underset{k\ne j}{\overset{J}{\Pi }}\frac{z-z_l^k}{z_l^j-z_l^k}\left[{ϵ}_{l\alpha }^j,{ϵ}_{l\beta }^j,{ϵ}_{lz}^j,{ϵ}_{l\alpha \beta }^j,{ϵ}_{l\alpha z}^j,{ϵ}_{l\beta z}^j\right]$

σ_l＝Cε_l；σ_l＝[σ_lα,σ_lβ,σ_lz,τ_lαβ,τ_lαz,τ_lβz]^T；ε_l＝[ε_lα,ε_lβ,ε_lz,γ_lαβ,γ_lαz,γ_lβz]^T

其中C为结构材料系数矩阵。

(6)根据步骤(5)建立结构能量方程(U，T)：

$\left(\begin{array}{cc}U=\frac{1}{2}\underset{l}{\Sigma }\int \int \int {ϵ}_l^T{C}_l{ϵ}_{l}dzd\beta d\alpha ,& T=\frac{1}{2}\underset{l}{\Sigma }\int \int \int {\rho }_l\left\{\frac{{\partial }^2{u}_{1l}}{\partial t^2}+\frac{{\partial }^2{u}_{2l}}{\partial t^2}+\frac{{\partial }^2{u}_{3l}}{\partial t^2}\right\}dzd\beta d\alpha \end{array}\right)$

(7)在步骤(6)基础上建立结构拉格朗日能量泛函L＝U-T，然后利用Ritz法求得结构 特征方程：

(K-ω²M)＝0

其中ω为圆频率。

(8)最后应用Arnoldi算法建立MATLAB求解器输出结构的振动特征数据并判定计算精 度，若满足精度要求则输出振动特征数据(固有频率，模态等)，不满足则继续优化空间分层 数量和增加面内位移展开式级数截取量。

层合板长L＝0.11179m，宽b＝0.0127m，总厚度H＝0.0038m,铺层形式为[45°/45°]且各层厚 度和材料均相等。材料参数如下：杨氏模量E₁＝37.41GPa，E₂＝E₃＝13.67GPa，剪切模量为 G₁₂＝5.478GPa,G₁₃＝6.03GPa，G₂₃＝6.666GPa，泊松比为μ₁₂＝μ₁₃＝μ₂₃＝0.3，密度为ρ＝1968.9kg/m³。 具体步骤如下：

(1)提取层合板的几何和材料参数并根据板结构特征选择以下坐标参数：α＝x，β＝y。 同时，设置结构的面内和横向位移为u₁(x,y,z)，u₂(x,y,z),和u₃(x,y,z)。

(2)对结构每一层沿厚度方向设置J个非均匀分布的计算平面。其中第1个和第J个计 算平面分别选取为该层的下表面和上表面。第2到第J-1个计算平面在厚度方向上的位置选 取为移位勒让德多项式P_J-2(z)在区间[0,h_l]的实零点，即：

第一层:$z_1^1=0,z_1^J=\frac{H}{2},z_1^j$为$P_{J-2}\left(z\right)=\frac{1}{2^{J-2}(J-2)!}\frac{d^{J-2}}{{\mathrm{dz}}^{J-2}}{\left[{\left(\frac{4z}{H}-1\right)}^2-1\right]}^{J-2}$

在区间的实零点

第二层:$z_2^1=\frac{H}{2},z_2^J=H,z_2^j-\frac{H}{2}$为$P_{J-2}\left(z\right)=\frac{1}{2^{J-2}(J-2)!}\frac{d^{J-2}}{{\mathrm{dz}}^{J-2}}{\left[{\left(\frac{4z}{H}-1\right)}^2-1\right]}^{J-2}$

在区间的实零点

同时，将第l层第j个计算平面上的结构位移设定为并将结构第l层任意位置的结 构位移设定为如下形式：

$u_{il}(x,y,z)=\underset{j=1}{\overset{J}{\Sigma }}\underset{k\ne j}{\overset{J}{\Pi }}\frac{z-z_l^k}{z_l^j-z_l^k}{u}_{il}^j(x,y)U_{il}^j$

(3)利用Chebyshev多项式对每个计算平面上的结构位移进行面内分量展开并通过边界 特征函数施加边界条件，得到计算平面上结构位移的详细表达式如下所示：

$u_{il}^j(x,y)U_{il}^j={(1+x)}^a{(1-x)}^b{(1+y)}^c{(1-y)}^d\underset{m=0}{\overset{M}{\Sigma }}\underset{n=0}{\overset{N}{\Sigma }}{U}_{ilmn}^{j}cos\left[m\arccos \left(x\right)\right]cos\left[n\arccos \left(y\right)\right]$

由自由边界条件求得边界特征函数系数a＝b＝c＝d＝0，即

$u_{il}^j(x,y)U_{il}^j=\underset{m=0}{\overset{M}{\Sigma }}\underset{n=0}{\overset{N}{\Sigma }}{U}_{ilmn}^{j}cos\left[m\arccos \left(x\right)\right]cos\left[n\arccos \left(y\right)\right]$

(4)由结构理论和步骤(2)到(3)计算第l层第j计算平面上结构的面内应变和横向 应变为：

$\left(\begin{array}{cc}{ϵ}_{lx}^j=\frac{\partial u_{1lmn}^j\left(x,y\right)}{\partial x}{U}_{1lmn}^j,& {ϵ}_{ly}^j=\frac{\partial u_{2lmn}^j\left(x,y\right)}{\partial y}{U}_{2lmn}^j\end{array}\right)$

${\gamma }_{lxy}^j=\frac{\partial u_{1lmn}^j\left(x,y\right)}{\partial y}{U}_{1lmn}^j+\frac{\partial u_{2lmn}^j\left(x,y\right)}{\partial x}{U}_{2lmn}^j$

${ϵ}_{lz}^j=u_{3lmn}^j(x,y)\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{3lmn}^j$

${\gamma }_{lxz}^j=u_{1lmn}^j\left(x,y\right)\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{1lmn}^j+\frac{\partial u_{3lmn}^j\left(x,y\right)}{\partial x}{U}_{3lmn}^j$

${\gamma }_{lyz}^j=u_{2lmn}^j(x,y)\underset{i=1}{\overset{J}{\Sigma }}{M}_{li}\left(z_l^j\right)U_{2lmn}^j+\frac{\partial u_{3lmn}^j(x,y)}{\partial y}{U}_{3lmn}^j$

(5)由步骤(2)和(4)求得结构第l层任意位置的应变和应力表达式为

$\left[{ϵ}_{lx},{ϵ}_{ly},{ϵ}_{lz},{ϵ}_{lxy},{ϵ}_{lxz},{ϵ}_{lyz}\right]=\underset{j=1}{\overset{J}{\Sigma }}\underset{k\ne j}{\overset{J}{\Pi }}\frac{z-z_l^k}{z_l^j-z_l^k}\left[{ϵ}_{lx}^j,{ϵ}_{ly}^j,{ϵ}_{lz}^j,{ϵ}_{lxy}^j,{ϵ}_{lxz}^j,{ϵ}_{lyz}^j\right]$

σ_l＝Cε_l；σ_l＝[σ_lx,σ_ly,σ_lz,τ_lxy,τ_lxz,τ_lyz]^T；ε_l＝[ε_lx,ε_ly,ε_lz,γ_lxy,γ_lxz,γ_lyz]^T

(6)根据步骤(5)建立结构振动能量方程(U，T)：

$\left(\begin{array}{cc}U=\frac{1}{2}\underset{l=1}{\overset{2}{\Sigma }}\int \int \int {ϵ}_l^T{C}_l{ϵ}_{l}dzdydx,& T=\frac{1}{2}\underset{l=1}{\overset{2}{\Sigma }}\int \int \int {\rho }_l\left\{\frac{{\partial }^2{u}_{1l}}{\partial t^2}+\frac{{\partial }^2{u}_{2l}}{\partial t^2}+\frac{{\partial }^2{u}_{3l}}{\partial t^2}\right\}dzdydx\end{array}\right)$

(7)在步骤(6)的基础上建立结构的拉格朗日能量泛函：L＝U+U_s-T，然后利用Ritz法， 即

$\frac{\partial (U+U_s-T)}{\partial U_{ilmn}^j}=0,\left(\begin{array}{c}j=1,2,...J\\ m=1,2,...M\\ n=1,2,...N\end{array}\right)$

求得结构特征方程：(K-ω²M)＝0。

(8)应用Arnoldi算法建立MATLAB求解器输出层合板结构的固有频率。

计算所得结果如下表所示。从表中我们可以看出本发明的方法具有很好的计算精度。