随着材料科学的发展, 具有多物理场耦合特性的电磁弹性复合材料已经广泛用于智能器件的开发与制造中^[1-2]。该类材料往往由两相或多相压电和压磁材料按照一定的方式复合而成。由于不同材料相间的失配性, 电磁弹性复合材料在制造和使用过程中会不可避免地出现界面裂纹。这些裂纹会在外荷载作用下产生应力集中, 进而引发结构失效破坏, 造成安全事故。因此, 研究电磁弹性复合材料的界面断裂问题具有重要的理论和实际意义。

在现有研究中, 解析方法往往仅限于求解无限大或半无限大电磁弹性材料中的裂纹问题, 如积分方程方法^[3]、Stroh变换方法^[4]。对于有限几何尺寸的含裂纹电磁弹性材料还主要依赖于数值方法, 如有限元方法^[5]、边界元法^[6]和无网格方法^[7]。然而, 大部分数值方法仅适用于单相电磁弹性材料, 无法直接应用于双相电磁弹性复合材料的界面断裂分析。为解决上述问题, 本文提出一种适用于电磁弹性复合材料反平面界面断裂分析的辛离散有限元方法。该方法将一类基于哈密顿体系的解析方法^[8]与传统有限元方法相结合, 能够简单、高效地获得相关断裂参数, 并同时获得裂纹尖端附近奇异物理场的显式表达式。目前, 该方法已经成功应用于弹性材料与压电材料的断裂分析中^[9]。

本文提出的辛离散有限元方法可以分为2步：第一, 根据裂纹尖端的位置将整体结构分为2类区域, 即包含奇异性的近场和无奇异性的远场；第二, 通过在裂纹尖端附近引入解析的辛本征解函数, 实现该区域内的未知量变换, 将数量庞大的节点未知量转化为少量的辛本征解待定系数。该方法可以避免传统有限元方法在计算断裂参数时的网格敏感性和路径依赖性, 直接提高计算效率和计算精度。

1 电磁弹性材料反平面有限元列式

考虑如图 1所示的含界面裂纹的电磁弹性复合材料。上下层材料分别记为材料1(M₁)和材料2(M₂), z轴选取为电磁弹性介质的极化方向, 坐标原点位于裂纹尖端。结构承受反平面剪应力τ₀, 平面内电位移D₀和磁感应强度B₀。曲线Γ₀将整体结构划分为2类区域, 即近场区域和远场区域。

在直角坐标系下, 电磁弹性复合材料在反平面荷载作用下的基本方程可以表示为:

本构方程

(1)

几何方程

(2)

平衡方程

(3)

式中，w⁽ⁱ⁾, ϕ⁽ⁱ⁾和φ⁽ⁱ⁾分别为反平面位移、电势和磁势; σ_k⁽ⁱ⁾, D_k⁽ⁱ⁾和B_k⁽ⁱ⁾是应力、电位移和磁感应强度分量; ε_kz⁽ⁱ⁾, E_k⁽ⁱ⁾和H_k⁽ⁱ⁾是应变、电场强度和磁场强度分量; C₄₄⁽ⁱ⁾, e₁₅⁽ⁱ⁾, h₁₅⁽ⁱ⁾, g₁₁⁽ⁱ⁾为弹性、压电、压磁和电磁耦合常数; κ₁₁⁽ⁱ⁾, χ₁₁⁽ⁱ⁾为介电常数和磁通率; F⁽ⁱ⁾, Q⁽ⁱ⁾和I⁽ⁱ⁾是体力、点电荷和电流密度; k=x, y, 上标i表示材料区域编号, M⁽ⁱ⁾为材料常数矩阵, 见附录(A-1)。

本文使用的电磁弹性材料单元为八节点四边形单元。单元内任意点未知量由节点未知量表示

(4)

式中，N_j是形函数。单元势能可以表示为

(5)

式中，L⁽ⁱ⁾为应变能密度函数

(6)

式中，σ⁽ⁱ⁾={σ_xz⁽ⁱ⁾, σ_yz⁽ⁱ⁾}^T, D⁽ⁱ⁾={D_x⁽ⁱ⁾, D_y⁽ⁱ⁾}^T, B⁽ⁱ⁾={B_x⁽ⁱ⁾, B_y⁽ⁱ⁾}^T, ε⁽ⁱ⁾={ε_xz⁽ⁱ⁾, ε_yz⁽ⁱ⁾}^T, E⁽ⁱ⁾={-E_x⁽ⁱ⁾, -E_y⁽ⁱ⁾}^T, H⁽ⁱ⁾={-H_x⁽ⁱ⁾, -H_y⁽ⁱ⁾}^T。

将(1)式、(2)式和(4)式带入(5)式可得

(7)

式中，P为等效节点载荷; u_e⁽ⁱ⁾为单元节点位移向量u_e⁽ⁱ⁾={w₁⁽ⁱ⁾, …, w₈⁽ⁱ⁾, ϕ₁⁽ⁱ⁾, …, ϕ₈⁽ⁱ⁾, φ₁⁽ⁱ⁾, …, φ₈⁽ⁱ⁾}^T。

根据最小势能原理δΠ⁽ⁱ⁾=0, 可得

(8)

式中，K_e是单元刚度矩阵。根据节点编号, 组集整体刚度矩阵、节点位移向量和整体载荷向量, 可得电磁反平面有限元列式为

(9)

式中，K为刚度矩阵, f为荷载向量, 下标N和F分别代表近场和远场区域。

2 电磁弹性材料Ⅲ型界面断裂辛方法

在近场区域, 定义全状态向量为

(10)

则电磁弹性复合材料Ⅲ型界面断裂问题的哈密顿控制方程可以表示为

(11)

式中，H⁽ⁱ⁾为哈密顿矩阵, 见附录(A-2)。界面连续条件为

(12)

裂纹面边界条件为

(13)

式中，。

利用分离变量法求解哈密顿方程(11), 并结合(12)式和(13)式得到问题的辛本征值和本征解, 见附录(A-3~A-7)。因此, 双材料电磁反平面的位移, 电势和磁势可以表示为

(14)

式中，M是所取的本征解项数, c_{j, k}为待定系数。

3 辛离散有限元方法

由(14)式, 近场区域内节点未知量可以表示为

(15)

式中，c={c_{0, 1}, c_{0, 2}, c_{0, 3}, c_{1, 1}, …, c_{1, 6}, c_{2, 1}, …, c_{M, 6}}^T为辛展开函数的待定系数向量

(m=1, 2, …, N_ΩN; j=1, 2, …, 3M+6), N_ΩN为近场区域内总节点数, (r_m, θ_m)是第m个节点的极坐标。将(15)式代入(9)式, 则电磁反平面断裂问题的辛离散有限元表达式为

(16)

4 强度因子和能量释放率

由(1)式和辛本征解可知, 在裂纹尖端(r=0)^[10], 广义应力强度因子可以表示为

(17)

式中, K₃, K_D和K_B分别为应力强度因子、电位移强度因子和磁感应强度因子, 是由反平面剪应力τ₀, 平面内电位移D₀和磁感应强度B₀耦合作用的结果。能量释放率为

(18)

式中，E=[(R⁽¹⁾)^-1+(R⁽²⁾)^-1]/2。

5 数值算例

为验证本文提出方法的正确性与有效性, 数值算例计算了3种典型裂纹对应的断裂参数。本节中所有计算数据均采用无量纲形式^[10], 涉及的材料参数如表 1所示。

5.1 算例1

考虑1个受到均布荷载作用的含边裂纹电磁弹性材料矩形板, 其几何尺寸和有限元网格如图 2所示。取τ₀=1, D₀=1, B₀=1, d=0.5a, M₂=M₁, 表 2给出了不同高度和宽度下各奇异物理场的强度因子和能量释放率。由于结构关于x轴对称, 故有K₃^*=K_D^*=K_B^*=K^*[11]。从表中数据可以看出, 各强度因子与能量释放率均随高度和宽度的增加而减小, 并且逐渐趋近于稳定值, 即H→∞和W→∞对应的值^[12-13]。当W/a=5时, 最大误差仅为1.15%。这表明本文方法具有良好的计算精度, 能够有效计算含裂纹电磁弹性材料的断裂参数。

5.2 算例2

图 3为1个含有平行裂纹的电磁弹性材料矩形板, 其裂纹长度分别为2a₁和2a₂。该板两端受到均布荷载τ₀=1, D₀=1和B₀=1作用。令W=1, H/W=0.5, H₁/W=0.5, 表 3给出了不同裂纹长度对应的裂纹尖端A处的强度因子与能量释放率。从表中数据可以看出, 对于给定的a₂, 各强度因子与能量释放率总是随着的a₁增加而增大; 而对于给定的a₁, 断裂参数表现出相反的变化趋势。表 4给出了裂纹长度相同时(a₁=a₂)断裂参数随H/W的变化规律。从表中可以看到, 各断裂参数随着H/W单调变化, 并逐渐趋近于某一固定值。该现象可以通过圣维南原理解释, 当H/W增加到一定值时, 边缘效应不再对裂纹尖端物理场分布产生影响。

5.3 算例3

考虑如图 4所示的1个含有折线裂纹的电磁弹性材料矩形板。裂纹与水平方向的偏离角度分别为θ_A和θ_B。令W=H, a=0.4H, τ₀=1, D₀=1和B₀=1, 表 5和表 6分别给出了的裂纹尖端A和B处的强度因子和能量释放率随θ_A和θ_B的变化规律。从表 5中可以看出, 所有断裂参数随θ_A的增大均表现出先增大后减小的趋势。当θ_A=0时, 应力、电位移强度因子和能量释放率达到最大值。表 6表现出与表 5相似的变化趋势。从上述现象可以看出, 合理控制外部电磁场可以有效改变电磁弹性材料裂纹尖端的能量释放率, 从而阻滞裂纹扩展。

6 结论

本文将辛离散有限元方法成功应用于电磁弹性复合材料在反平面荷载作用下的界面断裂分析中。该方法将传统有限元方法与哈密顿体系辛方法有机结合, 可以直接计算出高精度的应力、电场、磁场强度因子和能量释放率, 并同时获得裂纹尖端附近区域各物理场的显式表达式。相比其他数值方法, 该方法有三方面优势:①在裂纹尖端奇异区域内引入具有辛展开形式的解析解, 将该区域内传统有限元对应的大量节点未知量转化为少量辛本征解待定系数, 解决了传统有限元方法中断裂分析的网格敏感性问题, 并同时大幅提高了计算效率; ②引入的辛本征解函数能够直接表征裂纹尖端的奇异性, 避免了传统有限元方法中断裂参数计算的路径依赖性问题, 直接提高了应力强度因子计算精度; ③无需额外引进新的单元以及后处理程序, 可以直接获得裂纹尖端附近奇异物理场的显式表达式, 有利于该类结构的前期优化设计与后期安全评估。

附录A

(A-1)

(A-2)

式中

(A-3)

(A-4)

(A-5)

(A-6)

(A-7)

式中，ϕ_j=r^μ_jsin(μ_jθ), φ_j=r^μ_jcos(μ_jθ); μ_j=j/2是辛本征值; ρ=(R⁽²⁾)^-1R⁽¹⁾, n=1, 2, 3, …。