一种电磁材料层状壳体电磁弹耦合仿真模拟方法

摘要

本发明提供了一种电磁材料层状壳体电磁弹耦合仿真模拟方法，利用Hamilton扩展原则建立含电-磁-弹耦合效应的三维能量方程；基于变分渐近法将三维能量渐近扩展为系列二维递归能量，并利用壳体固有的小参数渐近修正二维递归能量中主导变分项，从而得到与原三维能量尽可能接近的修正模型，并转换为工程常用的Reissner-Mindlin模型形式；基于得到的二维全局响应和各阶翘曲函数推导了三维场变量重构关系。该模型不需先验性假设，可准确预测多场作用下结构的电磁弹耦合性能，计算量小，计算效率高于高阶层合理论和三维有限元解，占用计算机资源少。

说明书

技术领域

本发明涉及材料力学性能分析领域，具体涉及一种基于变分渐近法的电磁材料层状壳 体电磁弹耦合性能的仿真模型，能有效预测多场作用下结构的电-磁-弹耦合性能。

背景技术

由压电相和压磁相组成的智能复合材料会产生单相压电及压磁材料所没有的磁电耦合 效应。利用此特性，可改进压电与压磁单相材料的小频宽、小致动的缺点，因此广泛应用 于航空、汽车、电子、土木、医学等领域。

Van Suchtelen于1972年提出了压电、压磁的结合会产生新的材料特性即磁电耦合效 应。20世纪90年代中期，Nan、Huang和Kuo从理论上提出了压电、压磁材料的细观力学 模型来估算其耦合效应。Benveniste利用均匀场的概念得到了纤维状压电、压磁弹性体耦 合效应中的不同部分之间的精确关系。20世纪90年代末期，Li进一步发展了细观力学来 研究对于非均匀的耦合介质在静力场、电场、磁场之间同样存在着耦合，同时得到了无限 压电、压磁弹性体的双圆柱夹杂和各向异性体问题的解。Wu和Huang研究出压电、压磁复 合材料磁电耦合效应的封闭解。在国内，对于压电、压磁复合材料研究的起步较晚，2001 年，刘金喜等研究了压电、压磁材料二维问题Green函数的重要特性。2003年，王建国等 研究了横观各向同性无限压电、压磁的状态变量解。2006年，姚伟岸等研究了电磁弹性固 体的边界元方法；同年，周振功等研究了电磁复合材料的裂纹对弹性波的散射问题。

由上述分析可知，对于电磁材料所组成的智能结构分析中，需要有数学、力学、电学 和电磁学等多学科的知识，并且由于电、磁和弹性介质间耦合的复杂性，对问题求解带来 了很大的困难，虽然已有了一些有益的工作，但仍有许多方面需进一步完善。如现有文献 大多数在简化模型中假设线性或者高阶的位移和电、磁势分布，忽略或部分忽略了局部变 形的非均匀性，因而不能反映电磁器件与结构结合部的局部应力和局部电、磁场。

发明内容

针对现有技术中存在的上述不足，本发明提供一种计算量小，占用计算机资源少，且 效率高的基于渐近变分法电磁材料层状壳体电磁弹耦合性能仿真模拟方法。

为解决上述技术问题，实现发明目的，本发明采用的技术方案如下：

一种电磁材料层状壳体电磁弹耦合仿真模拟方法，包括以下步骤：

1）基于Hamilton原则和旋转张量分解概念建立电磁材料层状壳体中含电-磁-弹耦合 效应的三维能量方程；

2）基于渐近变分法将三维能量方程降维分析得到二维能量泛函，并利用电磁材料层状 壳体中固有的小参数渐近修正二维能量泛函中含翘曲项的主导变分项，得到修正模型，将 修正模型转换为Reissner-Mindlin模型形式；

3）基于推导的三维场重构关系，利用Reissner-Mindlin模型分析得到的二维壳面的 全局响应和降维分析得到的翘曲函数重构场变量沿厚度方向分布，由本构关系重构三维应 力、电位移、磁感应势；对电磁材料层状壳体电磁弹耦合进行仿真模拟。

进一步，所述步骤1具体为：基于Hamilton原则和旋转张量分解概念建立电磁材料层 状壳体中含电-磁-弹耦合效应的三维能量方程：

${\int }_{t_1}^{t_2}[\delta (K_{2D}+K^{*}-U)+\delta {\bar{W}}_{2D}+\delta {\bar{W}}^{*}]\mathrm{dt}=0;$

其中，t₁，t₂为任意两固定时刻点；K^*为动载产生的广义余能，K_2D为动载产生的二维结 构动能；为载荷、电/磁场所做的二维虚功；为载荷、电/磁场所做的剩余虚功，上 划线用于表明该虚功不需要对泛函精确变分，δ为变分符号，U为载荷、电/磁场在壳体 内产生的内能，其具体形式为：

$U=\frac{1}{2}{\int }_v[{\Gamma }^T{C}^{E,H}\Gamma -E^T{d}^{\Gamma ,H}E-H^T{\mu }^{\Gamma ,E}H-2E{e}^H\Gamma -2H{q}^E\Gamma -2E{a}^{\mathrm{\Gamma H}}]\mathrm{dv}$

式中：Γ,E,H分别为载荷、电/磁场下壳体产生的应变强度、电场强度和磁场强度；e^H,q^E,a^Γ分别为磁场强度恒定时的压电常数矩阵、电场强度恒定时的压磁常数矩阵和应变强度恒定 时的磁电常数矩阵；C^E,H,d^Γ,H,μ^Γ,E分别为电场强度和磁场强度恒定时的弹性常数矩阵、应 变强度和磁场强度恒定时的介电常数矩阵以及应变强度和电场强度恒定时的磁导率常数 矩阵；v为壳体空间体积。

进一步，所述步骤2具体为：基于渐近变分法将三维能量降维分析得到二维能量泛函： 利用电磁材料层状壳体中固有的小参数渐近修正二维能量泛函中含翘曲项的主导变分项， 得到零阶和一阶渐近修正模型：

2Π₀=∈^TA∈；

$2{\Pi }_1={\in }^T{A}_R\in +{\in }_{;\alpha }^T{B}_{\mathrm{\alpha \beta }}{\in }_{;\beta }+2{\in }^{T}F;$

其中，Π₀,Π₁分别为零阶，一阶近似得到的能量泛函；∈为广义二维应变量；A为二维刚 度矩阵；A_R,B_αβ为考虑壳体初始曲率修正后的刚度矩阵；α=β=1,2， A_α为壳体中二维壳面基向量的拉梅参数，x_α为在电磁材料层状壳体上建立的平面坐标系， x₁，x₂轴分别为沿壳体参考面的曲面方向、长度方向；F为荷载相关项，上标T表示转置 矩阵；

将一阶渐近修正模型转换为Reissner-Mindlin模型形式：

其中，为Reissner-Mindlin模型的应变量；$\gamma ={\left(\begin{array}{cc}2{\gamma }_{13}& 2{\gamma }_{23}\end{array}\right)}^T,$γ₁₃和 γ₂₃为横向剪切量；G为剪切刚度矩阵；分别为转换为Reissner-Mindlin模型的刚度 矩阵和荷载相关项。

进一步，所述步骤3具体为：基于推导的三维场重构关系，利用Reissner-Mindlin模 型分析得到的二维壳面的全局响应和降维分析得到的翘曲函数重构场变量沿厚度方向分 布：

$U_i=u_i+x_3{\left(\begin{array}{ccc}{C}_{31}& C_{32}& C_{33}\end{array}\right)}^T+S(V_0+{V_1}_{});$

其中，U_i，u_i分别为三维壳体变形和二维壳体变形列阵；C_ij为全局旋转张量；x₃为沿壳 体参考面的厚度方向上的坐标，V₀,V₁为零阶和一阶渐近修正翘曲节点值，S为形函数；

由本构关系重构三维应力、电位移、磁感应势：

$\left(\begin{array}{c}\sigma \\ D\\ B\end{array}\right)=\begin{array}{c}\left(\begin{array}{ccc}C& -e& -q\\ {-e}^T& -d& -a\\ {-q}^T& {-a}^T& -\mu \end{array}\right)\Gamma \end{array};$

其中，σ,D,B分别为载荷、电/磁场下壳体的三维应力、电位移、磁感应势；C,e,q,d,a,μ 分别为含弹性、压电、压磁、介电、磁电和导磁率常数的矩阵；Γ为广义应变量矩阵。

相比于现有技术，本发明具有如下优点：

1、本发明提供的电磁材料层状壳体电磁弹耦合仿真模拟方法，通过与三维有限元精 确解对比分析表明：采用变分渐近法来求解由非均质、各向异性双相电磁材料构成的壳体 结构的三维场分布精确度高。

2、本发明建立的电磁材料层状壳体电磁弹耦合模型属于单层壳模型，计算量小，计 算效率高于高阶层合理论和三维有限元解。

3、本发明提供的电磁材料层状壳体电磁弹耦合仿真模拟方法，只需建立单层壳模型， 划分单元和节点较三维有限元大为减少，占用计算机资源少。

4、本发明提供的电磁材料层状壳体电磁弹耦合仿真模拟方法，不需任何先验性假设， 具有数学上的严密性。

5、本发明提供的电磁材料层状壳体电磁弹耦合仿真模拟方法，可模拟完全各向异性 电磁材料层壳体的电-磁-弹耦合性能，突破了有限元只能处理宏观正交各向异性材料的局 限。

附图说明

图1为本发明提供的电磁材料层状壳体电磁弹耦合仿真模拟方法的流程图。

图2为具体实施例中三层电磁材料壳体几何构型。

图3为具体实施例中在三层电磁材料壳体上建立的坐标系。

图4为具体实施例中工况一中壳体沿壳体厚度方向的电位势分布图。

图5为具体实施例中工况一中壳体沿壳体厚度方向的磁位势分布图。

图6为具体实施例中工况一中壳体沿壳体厚度方向的磁感应势分布图。

图7为具体实施例中工况一中壳体沿壳体厚度方向的电位移分布图。

图8为具体实施例中工况一中壳体沿壳体厚度方向的横法向应力分布图。

图9为具体实施例中工况一中壳体沿壳体厚度方向的弹性位移分布图。

图10为具体实施例中工况二中壳体沿壳体厚度方向的电位势分布图。

图11为具体实施例中工况二中壳体沿壳体厚度方向的磁位势分布图。

图12为具体实施例中工况二中壳体沿壳体厚度方向的磁感应势分布图。

图13为具体实施例中工况二中壳体沿壳体厚度方向的电位移分布图。

图14为具体实施例中工况二中壳体沿壳体厚度方向的横法向应力分布图。

图15为具体实施例中工况二中壳体沿壳体厚度方向的弹性位移分布图。

具体实施方式

如图2、3所示的三层电磁材料壳体（图2中三角形标志处设置有固定支座，圆形标 志处设置有活动支座），取壳体中面为参考面，x₁，x₂，x₃轴分别为沿壳体参考面的曲面 方向、长度方向和厚度方向，沿x₁轴的初始曲率为k₁₁，h，R，l分别为壳体的厚度、曲率 半径和长度值。电磁材料层状壳体电磁弹耦合仿真模拟方法的流程如图1所示。

1、三维能量方程。

因涉及运算符号较为复杂，先做如下表述：下标i=1,2,3，j=1,2,3，α=1,2，β=1,2， A_α为与二维壳面基向量有关的拉梅参数，<>表示沿壳 体厚度方向的定积分，表示作用在壳体顶、底面的物理量。

基于Hamilton原理（即哈密顿原理），电磁材料壳体的弹性动力性能可表示为：

${\int }_{t_1}^{t_2}[\delta (K-U)+\delta \bar{W}]\mathrm{dt}=0---\left(1\right)$

式中：t₁，t₂为任意两固定时刻；K,U分别为动载产生的结构动能和载荷、电/磁场在壳体 内产生的内能，为载荷、电/磁场所做虚功，上划线表示不需对该项精确变分，δ为变 分符号。

对处在多物理场中的壳体结构，其内能泛函为压电、压磁弹性介质材料的广义势能以 及压电能、压磁能、电磁能组成的能量泛函：

$U=\frac{1}{2}{\int }_v[{\Gamma }^T{C}^{E,H}\Gamma -E^T{d}^{\Gamma ,H}E-H^T{\mu }^{\Gamma ,E}H-2E{e}^H\Gamma -2H{q}^E\Gamma -2E{a}^{\mathrm{\Gamma H}}]\mathrm{dv}---\left(2\right)$

式中：Γ,E,H分别为载荷、电/磁场下壳体产生的应变、电场和磁场强度；e^H,q^E,a^Γ分别 为磁场强度恒定时的压电常数矩阵、电场强度恒定时的压磁常数矩阵和应变恒定时的磁电 常数矩阵；C^E,H,d^Γ,H,μ^Γ,E分别为电场和磁场强度恒定时的弹性常数矩阵、应变和磁场强度 恒定时的介电常数矩阵以及应变和电场强度恒定时的磁导率常数矩阵；v为壳体空间体积。

基于旋转张量分解概念，式(2)中应变矢量Γi_j可定义为：

Γ_ij=(F_ij+F_ji)/2-△_ij   (3) 式中：△i_j为Kronecker符号（即克罗内克符号）；F_ij、F_ji均为变形梯度矢量的混合基分 量。

式(2)中电场和磁场强度可由电位势φ_ij和磁位势Ψ_ij确定。若不计体电荷和体电流，则 电位势和磁位势可用如下形式的变量变换表示：

$\left(\begin{array}{c}\phi \left(x_i\right)=\phi (x_1,x_2)+w_{\phi }(x_1,x_2,x_3),\\ \psi \left(x_i\right)=\psi (x_1,x_2+)w_{\psi }(x_1,x_2,x_3)\end{array}\right)---\left(4\right)$式中：φ为二维电势，取值为三维电势φ沿厚度方向的平均值；Ψ为二维磁势，取值为三 维电势Ψ沿厚度方向的平均值。这意味着电翘曲函数w_φ和磁翘曲函数w_Ψ需满足：

<w_φ(x₁,x₂,x₃)>=0,<w_ψ(x₁,x₂,x₃)>=0   (5)

三维电场和磁场强度E_i,H_i可定义为：

$\left(\begin{array}{c}{E}_{\alpha }=E_{2D_{\alpha }}-w_{\phi ;\alpha },E_3=-w_{\phi ,3}\\ H_{\alpha }=H_{2D_{\alpha }}-w_{\Psi ;\alpha },H_3=-w_{\Psi ,3}\end{array}\right)---\left(6\right)$式中：分别为二维电场和二维磁场，

由式(3)和(5)得到含h/R,(h/l)²阶次项的三维广义应变场为：

$\hat{\Gamma }={\Gamma }_h\hat{w}+{\Gamma }_{\in }+{\Gamma }_{\mathrm{Rh}}\hat{w}{\Gamma }_{R\in }+{\Gamma }_{l_{\alpha }}{\hat{w}}_{;\alpha }---\left(7\right)$

式中：

$\hat{\Gamma }={\left(\begin{array}{cccccccccccc}{\Gamma }_{11}& 2{\Gamma }_{12}& {\Gamma }_{22}& 2{\Gamma }_{13}& 2{\Gamma }_{23}& {\Gamma }_{33}& E_1& E_2& E_3& H_1& H_2& H_3\end{array}\right)}^T$

其中，w₁、w₂、w₃分别为荷载沿x₁、x₂、x₃轴产生的翘曲变形 值；$\in ={\left(\begin{array}{cccccccccc}{ϵ}_{11}& 2{ϵ}_{12}& {ϵ}_{22}& K_{11}& K_{12}+K_{21}& K_{22}& E_{2D_1}& E_{2D_2}& H_{2D_1}& H_{2D_2}\end{array}\right)}^T$其中

ε_αβ,K_αβ,(α,β=1,2)分别为壳体的面内应变和面外应变；Γ_h,Γ_∈,Γ_Rh,Γ_R∈为相应的积分 算子矩阵。

可将内能泛函表示为：

$U=\frac{1}{2}{\int }_{\Omega }<{\hat{\Gamma }}^T\left(\begin{array}{ccc}{C}^{E,H}& -e^H& -q^E\\ -{\left(e^H\right)}^T& -d^{\Gamma ,H}& -a^{\Gamma }\\ -{\left(q^E\right)}^T& -{\left(a^{\Gamma }\right)}^T& -{\mu }^{\Gamma ,E}\end{array},{\hat{\Gamma }}_{\mu }\right){\hat{\Gamma }}^{\mu }>\mathrm{d\Omega }=\frac{1}{2}{\int }_{\Omega }<{\hat{\Gamma }}^{T}D\hat{\Gamma }{\mu }^{}>\mathrm{d\Omega }---\left(8\right)$

式中：Ω表示参考壳面面域；μ为壳体几何修正系数，μ=1+x₂(k₁₁)+O(h²/R²)，O(h²/R²) 表示高于h²/R²阶项；D为12×12阶含弹性、压电、压磁、介电、磁电和导磁率常数的矩 阵，若为单斜对称材料，且绕自身中性轴旋转，则矩阵中的一些元素始终为零，得到：

式中：C^E,H表示矩阵$\left(\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{12}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right),$C_kl为电场强度E和磁场强度H恒定 时的弹性常数；-e^H表示矩阵$\left(\begin{array}{ccc}0& 0& e_{31}\\ 0& 0& e_{32}\\ 0& 0& e_{33}\\ 0& e_{24}& 0\\ e_{15}& 0& 0\\ 0& 0& 0\end{array}\right),$e_kl为磁场强度H恒定时的压电常数的负值； -q^E表示矩阵$\left(\begin{array}{ccc}0& 0& q_{31}\\ 0& 0& q_{32}\\ 0& 0& q_{33}\\ 0& q_{24}& 0\\ q_{15}& 0& 0\\ 0& 0& 0\end{array}\right),$q_kl为电场强度E恒定时的压磁常数的负值；-(e^H)^T表示矩 阵$\left(\begin{array}{cccccc}0& 0& 0& 0& e_{15}& 0\\ 0& 0& 0& e_{24}& 0& 0\\ e_{31}& e_{32}& e_{33}& 0& 0& 0\end{array}\right);$-d^Γ,H表示矩阵$\left(\begin{array}{ccc}{d}_{11}& 0& 0\\ 0& d_{22}& 0\\ 0& 0& d_{33}\end{array}\right),$d_kl为应变Γ和磁场强度H 恒定时的介电常数的负值；-a^Γ表示矩阵$\left(\begin{array}{ccc}{a}_{11}& 0& 0\\ 0& a_{22}& 0\\ 0& 0& a_{33}\end{array}\right),$a_kl为应变Γ恒定时的磁电常数的 负值；-(q^E)^T表示矩阵$\left(\begin{array}{cccccc}0& 0& 0& 0& q_{15}& 0\\ 0& 0& 0& q_{24}& 0& 0\\ q_{31}& q_{32}& q_{33}& 0& 0& 0\end{array}\right);$-(a^Γ)^T表示矩阵$\left(\begin{array}{ccc}{a}_{11}& 0& 0\\ 0& a_{22}& 0\\ 0& 0& a_{33}\end{array}\right);$-μ^Γ,E表示矩阵$\left(\begin{array}{ccc}{\mu }_{11}& 0& 0\\ 0& {\mu }_{22}& 0\\ 0& 0& {\mu }_{33}\end{array}\right),$μ_kl应变Γ和电场强度E恒定时的磁导率常数的负值；下标 k=l=1,2,3,4,5,6。

式(5)电磁翘曲约束加上3个力场翘曲约束，最终得到多场耦合下层状壳体的约束为：

$<{\Gamma }_c\hat{w}>=0---\left(10\right)$

式中：Γ_c为5×5阶单位矩阵。

荷载、电/磁场所作的虚功可表示为：

$\delta \bar{W}={\int }_{\Omega }(<P·\delta \hat{R}>+\tau ·\delta {\hat{R}}^{+}+\beta ·\delta {\hat{R}}^{-}-D^{\pm }\delta {\phi }^{\pm }-B^{\pm }\delta {\Psi }^{\pm })\mathrm{d\Omega }---\left(11\right)$

式中：P,τ,β分别为壳体的体力，壳体顶面的表面力，壳体底面的表面力。

P=P_iB_i、τ=τ_iB_i、β=β_iB_i P_i为体力的荷载集度，τ_i为壳体顶面表面力的荷载集度，β_i表示壳体底面表面力的荷载集 度，B_i为变形后的三元基向量；D^±为壳体顶、底面的电位移(电通密度)，B^±分别为壳体 顶、底面的磁感应势(磁通密度)；φ^±，Ψ^±别为壳体顶、底面的电位势和磁位势，为位 移场(含虚拟位移和旋转)的拉格朗日变分，为壳体顶面的位移场的拉格朗日变分， 为壳体底面的位移场的拉格朗日变分。

由于翘曲变形较小，可安全忽略中翘曲与虚拟旋转的乘积项，将式(11)简化为：

$\delta \bar{W}=\delta {\bar{W}}_{2D}+\delta {\bar{W}}^{*}---\left(12\right)$

式中，

$\delta {\bar{W}}^{*}={\int }_{\Omega }[<P_i{\mathrm{\delta w}}_i>+{\tau }_i\delta w_i^{+}+{\beta }_i\delta w_i^{-}{D}^{\pm }\delta w_{\phi }^{\pm }-B^{\pm }\delta w_{\psi }^{\pm }]\mathrm{d\Omega }---\left(14\right)$

式中，分别为壳体产生的虚拟位移和旋转；D⁺为壳体顶面的电位移，D^-为壳体 底面的电位移，B⁺为壳体顶面的磁感应势，B^-为壳体底面的磁感应势，为荷载沿x_i轴 在壳体顶面产生的翘曲变形值，为荷载沿x_i轴在壳体底面产生的翘曲变形值，为壳 体顶、底面的电翘曲变形值，为壳体顶、底面的磁翘曲变形值，广义力f_i和弯矩m_α可 定义为：

f_i=<P_i>+τ_i+β_i；m_α=<x₃P_α>+h(τ_α-β_α)/2   (15)

对于壳体的动能，可表示为：

$K=\frac{1}{2}{\int }_v\rho {\hat{V}}^T\mathrm{dv}=K_{2D}+K^{*}---\left(16\right)$

式中，为壳面上质点的绝对速度；

$\left(\begin{array}{c}{K}_{2D}=\frac{1}{2}{\int }_{\Omega }(\bar{\mu }{\hat{V}}^T+\hat{V}+2{\omega }^T\mu \bar{\xi }\hat{V}+{\omega }^T\mathrm{j\omega })\mathrm{d\Omega }\\ K^{*}=\frac{1}{2}{\int }_v\rho [{(\tilde{\omega }\tilde{w}+\partial \tilde{w}/\partial t)}^T+2{(\hat{V}+\tilde{\omega }\xi )}^T(\tilde{\omega }\hat{w}+\partial \hat{w}/\partial t\left)\mathrm{dv}\right]\end{array}\right)---\left(17\right)$

式中，ρ为质量密度，为惯性角速度ω的反对称矩阵；j为动态 分析常用的惯性常数：$\bar{\mu }=<\rho >,\mu \bar{\xi }={\left(\begin{array}{ccc}0& 0& <x_3\rho >\end{array}\right)}^T,j\left(\begin{array}{ccc}<x_3^2\rho >& 0& 0\\ 0& <x_3^2\rho >& 0\\ 0& 0& 0\end{array}\right)$

利用式(12)和(16)，可将式(1)写为：

${\int }_{t_1}^{t_2}[\delta (K_{2D}+K^{*}-U)+\delta {\bar{W}}_{2D}+\delta {\hat{W}}^{*}]\mathrm{dt}=0---\left(1\right)$

式中：为动载产生的二维结构动能和载荷、电/磁场所做的二维虚功；分别 为动载产生的广义余能和载荷、电/磁场所做的剩余虚功。

至此，建立了基于参考面位移的能量泛函方程，但这仅是原三维多场耦合问题的另一 表达形式，直接求解难度较大。主要的难度在于未知翘曲函数现有文献通常根据二维 变量预先假设的变化形式，直接将原三维连续模型降维为二维壳面模型。但对于由非均 质、各向异性双相电磁材料构成的壳体结构，引入这种假设会产生明显的误差。本发明借 助变分渐近法，利用壳体固有的小参数渐近计算翘曲函数，无需先验性假设，并可使计算 得到简化，可模拟完全各向异性电磁材料层合壳的电-磁-弹耦合性能，突破了有限元只能 处理宏观正交各向异性材料的局限。

2降维分析。

多场耦合下壳体振动以低频振动为主，且荷载产生的翘曲变形相对较小，在降维分析 中可安全忽略式(19)中的K^*和仅对含内能和虚功的能量泛函进行渐近分析。根据 变分原理，未知翘曲函数可通过对能量泛函Π取驻值得到：

δΠ=0   (20) 为了能有效分析层状结构，并与有限元中的二维板壳求解器相衔接，应用有限元法将三维 翘曲函数离散为一维有限元形式：

w(x_i)=S(x₃)V(x₁,x₂)   (21) 式中：S为形函数；V为沿横法线方向的翘曲场节点值矩阵。

可得能量泛函的离散形式为：

$\left(\begin{array}{c}2\Pi =V^T\mathrm{EV}+{2V}^T(D_{h\in }+D_{\mathrm{Rh}\in }\in +D_{\mathrm{hRh}}V+D_{\mathrm{hR}\in }\in +D_{{\mathrm{hl}}_{\alpha }}{V}_{;\alpha })\\ +{\in }^T(D_{\in \in }+2D_{\in R\in })\in +V_{;\alpha }^T{D}_{l_{\alpha }{l}_{\beta }}{V}_{;\beta }+V_{;\alpha }^T{D}_{l_{\alpha }\in }+{2V}^{T}L\end{array}\right)---\left(22\right)$

式中：L为荷载相关项，L=-S^+Tτ-S^-Tβ-<S^TP>-S^±Tφ-S^±TΨ， S⁺为壳体顶面的形函数，S^-为壳体底面的形函数；新引入的与几何形状和材料属性有关 的矩阵为：

$\left(\begin{array}{c}E=<{\left[{\Gamma }_{h}S\right]}^{T}D\left[{\Gamma }_{h}S\right]\mu >,D_{h\in }=<{\left[{\Gamma }_{h}S\right]}^T\mathrm{D\Gamma }{\in }_{}\mu >\\ D_{{\mathrm{hl}}_{\alpha }}=<{\left[{\Gamma }_{h}S\right]}^{T}D\left[{\Gamma }_{l_{\alpha }}\right]>,D_{l_{\alpha }{l}_{\beta }}=<{\left[{\Gamma }_{l_{\alpha }}S\right]}^{T}D\left[{\Gamma }_{l_{\beta }}S\right]>\\ D_{\in \in }=<{{\Gamma }_{\in }}^T{\mathrm{D\Gamma }}_{\in }\mu >,D_{l_{\alpha }\in }=<{\left[{\Gamma }_{l_{\alpha }}S\right]}^T{\mathrm{D\Gamma }}_{\in }>\\ D_{\mathrm{hRh}}=<{\left[{\Gamma }_{h}S\right]}^{T}D\left[{\Gamma }_{\mathrm{Rh}}S\right]>,D_{\mathrm{hR}\in }=<{\left[{\Gamma }_{h}S\right]}^T{\mathrm{D\Gamma }}_{R\in }>\\ D_{\in R\in }=<{{\Gamma }_{\in }}^T{\mathrm{D\Gamma }}_{\in }>,D_{\mathrm{Rh}\in }=<{\left[{\Gamma }_{\mathrm{Rh}}S\right]}^T{\mathrm{D\Gamma }}_{\in }>\end{array}\right)---\left(23\right)$

其中，为积分算子矩阵。

2.1零阶近似。

为验证变分渐近法应用于智能材料壳体建模的有效性，首先基于变分渐近法构建电磁 材料层状壳体古典层合模型。

式(10)翘曲函数的离散形式为：

V^THΨ=0   (24) 式中：H=S^TS；Ψ为E₀(零初始曲率)的正交化核心矩阵，Ψ^THΨ=I，I为单位矩阵。 这样问题转化为式(24)约束下式(22)最小化问题。

应用变分渐近法，首先需根据泛函的不同阶数找到主导项。由于只有翘曲是变化的， 只需找到含翘曲的主导项。式(22)零阶近似后泛函的主导项为

${2\Pi }_0^{*}=V^T{E}_{0}V+{2V}^T{D}_{h\in 0}+{\in }^T{D}_{\in \in 0}\in ---\left(25\right)$

式中：E₀为式（23）中公式E中μ=1(无几何修正)的对应矩阵，Dhò0为式（23）中公式D_hò中μ=1(无几何修正)的对应矩阵，D_òò0为式（23）中公式D_òò中μ=1(无几何修正)的对应 矩阵。

最小化式(25)，可得到零阶翘曲函数为

$V={\hat{V}}_0\in {=V}_0---\left(26\right)$

式中零阶广义翘曲节点值，V₀为沿横法线方向的零阶渐近修正翘曲节点值矩阵。

将式(26)代回式(25)，可得到渐近修正到零阶的二维能量泛函Π₀为

2Π₀=∈^TA∈   (27) 式中：为二维刚度矩阵。

式(27)与古典壳体理论形式相同，但未引入Kirchhoff假设(变形前后中面的法线为直 线，且与中面相垂直)。

2.2一阶近似。

零阶近似得到的古典层合壳体模型可较好地预测薄壳结构的全局和面内变量分布，但 对中厚壳体，还需利用壳体固有小参数h/R,h/l进行更高阶近似，以准确预测面外变量分 布。

为考虑壳体初始曲率的影响，首先对能量泛函Π_R主导项进行h/R阶修正：

2Π_R=∈^TA_R∈   (28) 式中：

$A_R=A+{\hat{V}}_0^T{E}^{*}{\hat{V}}_0^T+D_{\in \in }^{*}+{2\hat{V}}_0^T(D_{h\in }^{*}+F_{\mathrm{hR}\in }+D_{\mathrm{Rh}\in })+2D_{\in R\in }+2{\hat{V}}_0^T{D}_{\mathrm{hRh}}{\hat{V}}_0---\left(29\right)$

式中：E^*为式（23）中公式E中μ=μ-1的对应矩阵，为式（23）中公式D_h∈中μ=μ-1 的对应矩阵，为式（23）中公式D_∈∈中μ=μ-1的对应矩阵。

其次，为准确描述壳体的横向剪切变形，需推导(h/l)²阶修正能量泛函。为此，将零 阶翘曲简单摄动为：

V=V₀+V₁   (30) 其中，V₁为沿横法线方向的一阶渐近修正翘曲节点值矩阵。

可得到一阶近似能量泛函主导项为

$2{\Pi }_1^{*}={V_1}^T{\mathrm{EV}}_1+2{V_1}^T{D}_{\alpha }{\in }_{;\alpha }+{{2V}_1}^{T}L---\left(31\right)$

式中：

与零阶近似类似，一阶翘曲可由下式求解

V₁=V_1α∈_;α+V_1L(32) 式中：V_1α,V_1L分别为一阶修正翘曲节点值和一阶荷载节点值。

式(32)代入式(31)，得到渐近修正到(h/l)²、(h/R)阶总能量泛函Π₁为

${2\Pi }_1={\in }^T{A}_R\in +{\in }_{;\alpha }^T{B}_{\mathrm{\alpha \beta }}{\in }_{;\beta }+2{\in }^{T}F---\left(33\right)$

式中：

$\left(\begin{array}{c}{B}_{\mathrm{\alpha \beta }}={\hat{V}}_0^T{D}_{l_{\alpha }{l}_{\beta }}{\hat{V}}_0+B_{1\alpha }^T{D}_{\beta },\\ F={\stackrel{Y}{V}}_0^{T}L-(D_{\alpha }^T{V}_{1L;\alpha }+V_{1\alpha }^T{L}_{;\alpha })/2\end{array}\right)---\left(34\right)$

2.3Reissner模型形式转换。

尽管式(33)渐近修正到二阶，但因含二维广义应变的导数∈_;α，难以在实际工程中有效 应用。为得到实用的能量泛函，可将近似能量转换为在实际工程中常用的 Reissner-Mindlin模型形式。为此，需增加2个横向剪切量$\gamma ={\left(\begin{array}{cc}2{\gamma }_{13}& 2{\gamma }_{23}\end{array}\right)}^T$作为独立自由 度。Reissner-Mindlin应变量与∈的关系为：

式中：为含古典广义应变量 的矩阵，其中，${ϵ}_{11}^{*},{ϵ}_{12}^{*},{ϵ}_{22}^{*},K_{11}^{*},K_{12}^{*},K_{22}^{*},E_{2D_1}^{*},E_{2D_2}^{*},H_{2D_1}^{*},H_{2D_2}^{*}$均为古典广义应变 量。

将式(35)代入式(33)，可得到Reissner应变表示的修正到二阶的总能量泛函为：

为消除式(37)中变量的偏导数可利用式(15)的两个弯矩平衡方程推导出γ与的关系式为：

式中：G为剪切刚度矩阵，F_γ为剪切力，m₁,m₂分别为沿x₁,x₂方向的弯矩。

利用式(38)，可将式(37)改写为

式中：

广义Reissner-Mindlin模型的形式为

若式(39)中的

则式(39)与式(41)等效。令实际工程中U*并不为零， 但可通过最小二乘法等优化技术使U^*尽可能趋近于零，得到精确的计算结果，在此不再赘 述。

2.4重构关系。

构建的电磁层状壳体广义Reissner-Mindlin模型可准确预测壳面的全局响应，对于 更为关注的场变量沿厚度方向的分布情况，则需重构关系进行预测。为此，利用二维壳面 分析得到的全局响应和降维分析得到的翘曲函数重构场变量沿厚度方向分布：

$U_i=u_i+x_3{\left(\begin{array}{ccc}{C}_{13}& C_{32}& C_{33}-1\end{array}\right)}^T+S(V_0+V_1)---\left(43\right)$式中：U_i,u_i分别为三维壳体变形和二维壳体变形列阵；C_ij为全局旋转张量。

由本构关系可重构三维应力、电位移、磁感应势为

$\left(\begin{array}{c}\sigma \\ D\\ B\end{array}\right)=\begin{array}{c}\left(\begin{array}{ccc}C& -e& -q\\ {-e}^T& -d& -a\\ {-q}^T& {-a}^T& -\mu \end{array}\right)\Gamma \end{array}---\left(44\right)$

式中：σ,D,B分别为载荷、电/磁场下壳体的三维应力、电位移、磁感应势；Γ为广义应 变量矩阵，C、e、q、d、a、μ分别为含弹性、压电、压磁、介电、磁电和导磁率常数 的矩阵。

为了验证本发明建立的方法的有效性和准确性，利用柱形弯曲问题进行分析。电磁材 料层状壳体采用非均质、各向异性双相电磁材料构成的壳体结构，上下两层采用压电材料 BaTiO3(简记为B)；内层为磁性材料CoFe2O4(简记为F)。各层厚度均为0.1m，壳体半径 R=10mm，夹角采用的坐标系为x₂∈[0,l],x₃∈[-h/2,h/2]。分二种 工况：

工况一：壳体顶面作用正弦型载荷：无电场、磁场 作用。

工况二：压电层与压磁层接触面接地，压电层外表面作用电势 无机械场、磁场作用。电磁材料常数列于表1。

工况一中壳体沿壳体厚度方向的场变量分布图，如图4～图9所示，工况二中壳体沿 壳体厚度方向的场变量分布图如图10～图15所示。可以看出，采用本申请提供的电磁弹 耦合仿真模拟方法对上述两种不同的工况中的壳体进行仿真模拟时，重构的壳体的场变量 分布与三维有限元精确解吻合性很好。其中：①工况一的弹性位移u₃小于工况二，但两者 的弹性位移趋势完全相反；②二种工况下的磁位势Ψ曲线在交界处连续，但曲线倾角斜率 并不连续；③由于压电、压磁耦合效应在不同区域产生作用，两种工况下的横法向应力σ₃₃性质、变化趋势截然不同；④两种工况下的磁感应势B₃具有相同的变化趋势，即压电层磁 感应接近于零(压电层的压磁系数为零)，磁质材料层的磁感应呈非线性变化(磁电耦合效 应)，此时若仍采用线性假设会引起较大误差；⑤两种工况下压磁层的电位移D₃均为零(压 电系数为零)，但工况一的压电层电位移呈非线性变化(压电耦合效应)；工况二无规律可 循；⑥工况一压电层电位势呈非线性分布(压电耦合效应)，工况二压电层电位势呈线性分 布，压磁层电位势为零。

表2给出三维有限元解和本发明解的计算规模和计算时间。从表2中可看出：采用本 发明计算结构的响应时，需要数值求解的未知量个数远小于三维有限元法需要求解的未知 量个数。本发明解的计算量小，计算规模大大减少，提高了计算效率。

表1电磁材料常数

由此可见，本发明提供的电磁材料层状壳体电磁弹耦合仿真模拟方法，通过与三维有 限元精确解对比分析表明：采用变分渐近法来求解由非均质、各向异性双相电磁材料构成 的壳体结构的三维场分布精确度高。计算量小，计算效率高于高阶层合理论和三维有限元 解，并且，只需建立单层壳模型，划分单元和节点较三维有限元大为减少，占用计算机资 源少。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例 对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进 行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权 利要求范围当中。