IMPROVED LOCAL RADIAL BASIS FUNCTION COLLOCATION METHOD FOR BAND STRUCTURE ANALYSIS OF CRACKED PHONONIC CRYSTALS
-
摘要: 基于直接法局部径向基函数配点法(local radial basis function collocation method, LRBFCM)提出了含裂纹声子晶体反平面弹性波的能带结构计算方法, 分析了其能带结构特性, 并通过有限元法计算对比验证该数值结果的准确性和有效性. 在边界上采用直接法沿方向取局部节点, 解决了边界求解不稳定问题, 讨论了数值方法的处理技巧对结果的影响. 并通过考虑不同声阻抗比、散射体形状(方形、圆形)和裂纹情况, 验证了该方法计算含裂纹声子晶体能带结构的适用性. 研究了形状参数、配点数对计算结果的影响. 最后深入分析了不同长度和位置的裂纹对声子晶体能带结构的影响特性, 并进行了对比分析. 文章的创新性在于用直接法局部径向基函数配点法解决了含裂纹声子晶体能带结构计算的问题, 极大增加了声子晶体的应用价值. 研究结果显示: 随着裂纹的扩展, 能带结构带隙逐渐变窄; 当裂纹扩展到一定程度时, 金散射体的声子晶体带隙数量会增加, 且新增带隙会随着裂纹扩展而变宽, 但铝散射体的声子晶体带隙数量会减少; 直接法局部径向基函数配点法可以大大提高含裂纹声子晶体能带结构的计算效率和计算精度.Abstract: Based on the direct local radial basis function collocation method (LRBFCM), the band structure calculation algorithm for the anti-plane elastic wave of the cracked phononic crystals is proposed. The band structure characteristics are analyzed, and the accuracy and validity of the numerical results are verified with the comparison of the finite element method analysis. For the boundary collocation nodes, the direct method is adopted to select local nodes along the normal direction to solve the problems of stability, and the influence of the processing techniques of the numerical methods on the results is discussed. The applicability of the method for calculating the band structure of cracked phononic crystals is verified by considering different acoustic impedance ratios, scatterer shapes (square, circular) and crack conditions. The effects of shape parameters and the number of point numbers on the calculation results are investigated. Finally, the effects of crack of different length and its location on the band structure of phononic crystals are analyzed comprehensively, and a comparative analysis is carried out. The innovation of this paper is to solve the problem of calculating the band structure of cracked phononic crystals by the direct method of local radial basis function collocation method, which greatly increases the application value of phononic crystals. The results in this work show that due to the crack propagation, the bandgap of the phononic crystals structure gradually narrows; When the crack expands to a certain extent, the number of bandgaps in the phononic crystals of aurum scatterers increases, and the new bandgaps widen with the expansion of the crack. However, the number of bandgaps in the phononic crystals of aluminum scatterers decreases; The direct local radial basis function collocation method can improve the calculation efficiency and the computational accuracy significantly for the band structure of cracked phononic crystals.
-
引 言
声子晶体是由不同材料按结构周期性组合排列形成的人工复合材料. 弹性波在声子晶体中传播时, 受到其内部周期结构的影响, 形成特殊的能带结构. 两条能带结构之间弹性波/声波无法穿透的频率, 简称带隙, 具有独特的性质. 理论上带隙频率范围的弹性波/声波传播被抑制, 只有其他频率范围的弹性波/声波可以通过声子晶体传播[1]. 因此这种特性被应用于波的调控、能量捕获等多个方面[2-3]. 通过声子晶体的周期结构的缺陷、材料特性和几何参数等设计, 可以调整能带结构中的禁带带宽的大小[4], 有效达到控制噪声和振动的目的[5]. 因此, 研究声子晶体中裂纹和能带结构之间的关系非常重要. 而实验研究成本较大, 材料裂纹加工也存在困难, 建立含裂纹的声子晶体的能带计算模型, 并通过数值方法研究弹性波在含裂纹的声子晶体中传播特性很有必要.
近期许多计算声子晶体带隙的数值方法被开发用于研究声子晶体的能带结构, 如平面波展开法(PWE)[6-7]、多重散射理论(MST)[8]、时域有限差分法(FDTD)[9]、Dirichlet-to-Neumann映射法(DtN-map)[10]、边界元法(BEM)[11]和有限元法(FEM)[12-13]等. 然而, 每一种数值方法都有其优缺点. 如PWE和FDTD方法在处理散射体和基体之间存在较大声阻抗比时不够准确, 尤其是在处理流固界面条件时. MST和DtN-map方法理论推导较为复杂, 算法适用性不强, 且由于太多的矩阵, 造成很大的计算量. BEM求解复杂问题时采用的基本解限制了其应用, 特别是对三维一般问题. FEM在流固相互作用界面上, 又必须使用大量的单元以满足精度要求, 会导致计算效率大大降低. 因此, 还需要寻找一些更合适且更有效的数值方法研究声子晶体.
由于径向基函数配点法(radial basis function collocation method, RBFCM)采用欧几里得距离作为基函数, 径向基函数配点法具有无网格求解高维问题等优点[14-17]. 该方法在20世纪90年代被Kansa提出求解偏微分方程, 也被称为Kansa方法. 2000年左右受到学者们的广泛关注, 其中Chen[18]和Zhang等[19]都分别对该方法做出了改进, Chen等[20]提出了一种加权径向基函数配点法, 该方法解决了边界和内部节点矩阵元素数值差异较大的问题, 增加了该方法的适用性. 但该方法仍然对形状参数选取较为敏感. 最近Chen等[21]提出了一种源点和配点分离的技术, 该方法可以让径向基函数配点法收敛性大大提高, 并且对形状参数也不敏感, 但还是无法解决径向基函数配点法满秩矩阵的问题.
对于径向基函数配点法的满秩矩阵问题, 国内外学者提出了很多方法解决这一难题. 其中方法之一是更换径向基核函数, 最著名的是吴宗敏[22]所提出的紧支径向基函数(compact support RBF), 和Wendland[23]提出的Wendland径向基函数. 这类径向基核函数对较远距离的节点影响较小, 被应用于配点法[24-25]. 第2种方法是局部径向基函数配点法(local radial basis function collocation method, LRBFCM), 也叫差分格式的径向基函数配点法 (radial basis function finite difference (RBF-FD) method)[26]. 其关键特征是配点发生在重叠的局部域上, 可以在不降低精度的前提下减小配置矩阵的大小, 求解效率得到大大提高. 但是局部径向基函数配点法在求解第二类边界时, 存在着不稳定性. 有限线法[27]是一种类似直接法的半解析方法, 首先对每个配置点建立相交的线系, 线系形成后, 在其中的每一条线上定义若干个节点进行函数逼近. 有限线法利用了线单元的概念, 把高维问题变成若干一维问题, 简化了积分处理, 提高了稳定性, 但对线系的依赖性较强. Zheng等[28-29]提出了3种不同的数值技术来计算边界条件导数, 并成功应用于频域空间求解声子晶体的能带结构, 显示出了很高的计算精度及效率, 其中沿方向取局部节点的直接法, 使用方向上的一些相邻节点的位置信息进行计算, 这种方法的优势在于可以将二维或者三维问题情况下的边界导数求解计算简化成一维情况, 使得稳定性和精确度大大提高. 直接法的思路与有限线法类似, 但直接法是基于局部径向基函数配点法提出的, 往往仅应用于边界导数的处理. Yan等[30]通过Houbolt方法离散时间研究了局部径向基函数配点法应用于二维声子晶体在平面内的传播问题. 对于含裂纹问题, Wen等[31]提出了含有奇异性的径向基函数, 较为精确地确定裂纹尖端的应力强度因子, 有效地扩充了RBF的应用. 但含裂纹声子晶体的径向基函数配点法目前研究较少见.
由于裂纹问题中应力含有奇异性, 因此传统的数值方法需要在裂纹尖端周围采用较密的网格, 计算效果并不理想. 本文利用沿方向取点的直接法来处理裂纹面含导数边界条件问题结合局部径向基函数配点法, 创新性地提出了一种快速有效解决单胞含裂纹声子晶体的能带结构计算的方法. 并在考虑不同材料及不同散射体形状情况下, 充分讨论和研究了裂纹对反平面声子晶体能带结构的影响.
本文分别选取金(Au)和铝(Al)作为散射体, 形状分别考虑为方形和圆形的含裂纹声子晶体. 利用直接法局部径向基函数配点法计算该声子晶体在裂纹长度和位置改变情况下的能带结构, 并在此基础上深入分析了裂纹对声子晶体能带结构的影响特性.
1. 局部径向基函数配点法
径向基函数(RBF)是一类以欧氏距离(点$x$到点${x_i}$的距离)为自变量的函数. 在多变量逼近理论中, 用RBF进行插值已经成为了一种非常有效的数值计算工具, 其广泛应用于求解偏微分方程(PDE)边值问题. 常用的基函数类型有MQ型($ \sqrt {{r^2} + {c_s}^2} $)、三次样条型($ {r^3} $)、薄板样条型($ {r^2}\ln r $)、四次B样条型等.
径向基函数的选取和问题有关, 当求解问题的解是线性光滑时, 样条型径向基函数理论上得到较好的结果. 如果问题比较复杂, 尤其是有应力集中的情况下, 那么高阶径向基函数在同样节点的情况下可以得到较好的结果. 本文考虑的是含有裂纹的声子晶体, 选取具有较高收敛率的MQ径向基函数[32]
$$ \psi = \sqrt {{r^2} + {c_s}^2} $$ (1) 其中, $ {c_s} $是形状参数, $r = \left\| {{{\boldsymbol{x}}_i} - {\boldsymbol{x}}} \right\|$是${{\boldsymbol{x}}_i}$和${\boldsymbol{x}}$之间的欧氏距离, ${\boldsymbol{x}}$表示${{\boldsymbol{x}}_i}$附近的$ {N_s} $个节点的坐标, 即图1中的虚线圈内的节点, 另外蓝色点表示内部点, 黑色点表示第二类边界(Neumann boundary)点, 橙色点表示第一类边界(Dirichlet boundary)点.
考虑以下边值问题
$$\qquad\qquad Lu({\boldsymbol{x}}) = f({\boldsymbol{x}}), \quad {\boldsymbol{x}} \in \varOmega $$ (2) $$\qquad\qquad Bu({\boldsymbol{x}}) = g({\boldsymbol{x}}),\quad {\boldsymbol{x}} \in {\varGamma _1} $$ (3) $$\qquad\qquad u({\boldsymbol{x}}) = h({\boldsymbol{x}}), \quad {\boldsymbol{x}} \in {\varGamma _2} $$ (4) 其中, $ u $为计算域的解集, $ L $和$ B $是微分算子, $ \varOmega $, $ {\varGamma _1} $和$ {\varGamma _2} $分别对应于图1中的计算域和计算边界, 式(2)为控制方程, 式(3)为Neumann边界条件, 式(4)为Dirichlet边界条件.
如图1所示, 局部径向基函数配点法以插值基点${x_i}$为中心, 选取附近的$ {N_s} $个节点构成局部计算域, 中心点的位移向量$ u $的通解可近似为
$$ u({{\boldsymbol{x}}_i}) = \sum\limits_{n = 1}^{{N_s}} {\psi (\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_n}} \right\|)} {\alpha _n} $$ (5) 其中, $ {N_s} $为局部域内的节点总数, $ n $为局部域内节点的序号, 将式(5)代入式(2) ~ 式(4)可写为
$$ Lu({\boldsymbol{x}}) = \sum\limits_{n = 1}^{{N_s}} {L\psi (\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_n}} \right\|)} {\alpha _n},\quad {\boldsymbol{x}} \in \varOmega $$ (6) $$ Bu({\boldsymbol{x}}) = \sum\limits_{n = 1}^{{N_s}} {B\psi (\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_n}} \right\|)} {\alpha _n},\quad {\boldsymbol{x}} \in {\varGamma _1} $$ (7) $$ u({\boldsymbol{x}}) = \sum\limits_{n = 1}^{{N_s}} {\psi (\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_n}} \right\|)} {\alpha _n},\quad {\boldsymbol{x}} \in {\varGamma _2} $$ (8) 由式(5)可得
$$ {{\boldsymbol{\alpha }}_{{N_s} \times 1}} = {{\boldsymbol{\varPsi}} ^{ - 1}}_{{N_s} \times {N_s}}{{\boldsymbol{u}}_{{N_s} \times 1}} $$ (9) 其中, $ {{\boldsymbol{u}}_{{N_s} \times 1}} = {\left[ {{u_1}, {u_2}, \cdots , {u_{{N_s}}}} \right]^{\mathrm{T}}} $是局部域内节点的位移向量, ${{\boldsymbol{\alpha}} _{{N_s} \times 1}} = {\left[ {{\alpha _1}, {\alpha _2}, \cdots , {\alpha _{{N_s}}}} \right]^{\mathrm{T}}}$是未知系数, $ {{\boldsymbol{\varPsi}} _{{N_s} \times {N_s}}} = {[\psi \left( {\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_n}} \right\|} \right)]_{1 \leqslant i, n \leqslant {N_s}}} $是大小为$ {N_s} \times {N_s} $的矩阵. 将式(9)代入式(5)可得
$$ u({{\boldsymbol{x}}_i}) = \sum\limits_{n = 1}^{{N_s}} {\psi (\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_n}} \right\|)} {\alpha _n} = {{\boldsymbol{\varTheta}} _{1 \times {N_s}}}{\boldsymbol{\varPsi}} _{{N_s} \times {N_s}}^{ - 1}{{\boldsymbol{u}}_{{N_s} \times 1}} $$ (10) 其中, $ {{\boldsymbol{\varTheta}} _{1 \times {N_s}}} = \left[ {\psi \left( {\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_1}} \right\|} \right),\psi \left( {\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_2}} \right\|} \right), \cdots , \psi \left( {\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_n}} \right\|} \right)} \right] $. 为方便起见, 定义$ {\bar {\boldsymbol{\varPsi}} _{1 \times {N_s}}} = {{\boldsymbol{\varTheta}} _{1 \times {N_s}}}{\boldsymbol{\varPsi}} _{{N_s} \times {N_s}}^{ - 1}{\text{ = }}{\delta _{in}} $, 式中$ {\delta _{in}} $为狄拉克函数, 则式(6) ~ 式(8)可写成
$$\qquad\quad \begin{split} & Lu({\boldsymbol{x}}) = L{{\boldsymbol{\varTheta}} _{1 \times {N_s}}}{\boldsymbol{\varPsi}} _{{N_s} \times {N_s}}^{ - 1}{{\boldsymbol{u}}_{{N_s} \times 1}} = \\ &\qquad L{{\bar {\boldsymbol{\varPsi}} }_{1 \times {N_s}}}{{\boldsymbol{u}}_{{N_s} \times 1}} = f({\boldsymbol{x}}), \quad {\boldsymbol{x}} \in \varOmega \end{split} $$ (11) $$\qquad\quad \begin{split} & Bu({\boldsymbol{x}}) = B{{\boldsymbol{\varTheta}} _{1 \times {N_s}}}{\boldsymbol{\varPsi}} _{{N_s} \times {N_s}}^{ - 1}{{\boldsymbol{u}}_{{N_s} \times 1}} = \\ &\qquad B{{\bar {\boldsymbol{\varPsi}} }_{1 \times {N_s}}}{{\boldsymbol{u}}_{{N_s} \times 1}} = g({\boldsymbol{x}}),\quad {\boldsymbol{x}} \in {\varGamma _1} \end{split} $$ (12) $$\qquad\quad \begin{split} & u({\boldsymbol{x}}) = {{\boldsymbol{\varTheta}} _{1 \times {N_s}}}{\boldsymbol{\varPsi }}_{{N_s} \times {N_s}}^{ - 1}{{\boldsymbol{u}}_{{N_s} \times 1}} = \\ &\qquad {{\bar {\boldsymbol{\varPsi}} }_{1 \times {N_s}}}{{\boldsymbol{u}}_{{N_s} \times 1}} = h({\boldsymbol{x}}),\quad {\boldsymbol{x}} \in {\varGamma _2}\end{split} $$ (13) 令除去局部域内计算点外其余值为0, 则可得到计算域内全部点的数值解稀疏矩阵, 用矩阵表示式(11) ~ 式(13), 可得
$$ Lu({\boldsymbol{x}}) = L{{\bar {\boldsymbol{\varPsi}} }_{1 \times {N_s}}}{{\boldsymbol{u}}_{{N_s} \times 1}} = L{{\tilde {\boldsymbol{\varPsi}} }_{1 \times N}}{{\tilde {\boldsymbol{u}}}_{N \times 1}} = f({\boldsymbol{x}}),\;\; {\boldsymbol{x}} \in \varOmega $$ (14) $$ Bu({\boldsymbol{x}}) = B{{\bar {\boldsymbol{\varPsi}} }_{1 \times {N_s}}}{{\boldsymbol{u}}_{{N_s} \times 1}} = B{{\tilde {\boldsymbol{\varPsi}} }_{1 \times N}}{{\tilde {\boldsymbol{u}}}_{N \times 1}} = g({\boldsymbol{x}}), \quad {\boldsymbol{x}} \in {\varGamma _1} $$ (15) $$ u({\boldsymbol{x}}) = {{\bar {\boldsymbol{\varPsi}} }_{1 \times {N_s}}}{{\boldsymbol{u}}_{{N_s} \times 1}} = {{\tilde {\boldsymbol{\varPsi}} }_{1 \times N}}{{\tilde {\boldsymbol{u}}}_{N \times 1}} = h({\boldsymbol{x}}),\quad {\boldsymbol{x}} \in {\varGamma _2} $$ (16) 将上述微分方程写为矩阵形式, 则有
$$ {\boldsymbol{A}}\tilde {\boldsymbol{u}} = {\boldsymbol{b }}$$ (17) 其中$ {\boldsymbol{A}} = {\left[ \begin{gathered} {(L\tilde {\boldsymbol{\varPsi}} )_{{N_i} \times N}} \\ {(B\tilde {\boldsymbol{\varPsi}} )_{{N_{b1}} \times N}} \\ {{\tilde {\boldsymbol{\varPsi}} }_{{N_{b2}} \times N}} \\ \end{gathered} \right]_{N \times N}} $, $ {\boldsymbol{b}} = {\left[ \begin{gathered} f{({\boldsymbol{x}})_{{N_i} \times 1}} \\ g{({\boldsymbol{x}})_{{N_{b1}} \times 1}} \\ h{({\boldsymbol{x}})_{{N_{b2}} \times 1}} \\ \end{gathered} \right]_{N \times 1}}. $
最后可以直接由矩阵运算得到数值解
$$ \tilde {\boldsymbol{u}}({\boldsymbol{x}}) = {{\boldsymbol{A}}^{ - 1}}{\boldsymbol{b}} $$ (18) 为了解决局部径向基函数配点法在边界导数求解不稳定的问题, 在边界上沿着方向取局部节点的“直接法”往往被广泛采用, 如图2(b)所示直接法使用n方向上的一些相邻节点的位置信息进行计算, 左图所示为选取临近节点信息, 具体内容参见文献[29-31].
![]()
图 2 (a)传统法和(b)直接法选点示意图Figure 2. (a) The traditional method and (b) the direct method of node selection为验证局部径向基函数配点法求解反平面裂纹问题的有效性, 考虑一个含边裂纹的反平面问题(如图3所示), 弹性归一化参数$E = 1$, $\nu = 0. 3$, 长度为$b = 2$, 宽度为$h = 2$, 剪切载荷为$ {\tau _0} $, 方向垂直纸面, 形状参数${c_s} = 1$. 根据对称性, 可取模型的一半进行计算分析, 由于各向同性弹性体反平面裂纹问题$xOy$平面上的位移$ {u_x} = {u_y} = 0 $, 且$z$方向上的位移$u = u(x, y) \ne 0$, 仅仅是$x$与$y$的函数, 因此, 反平面裂纹问题的控制方程选用拉普拉斯方程如下
$$ \frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}} = 0 $$ (19) 边界条件为
$$ \left. \begin{split} & y = 0,\;\; x > a,\;\; u = 0 \\ & y = 0,\;\; 0 \leqslant x \leqslant a,\;\; {\sigma _{yz}} = 0 \\ & 0 \leqslant y \leqslant \frac{h}{2}, \;\; x = 0\& x = b,\;\; {\sigma _{xz}} = 0 \\ & y = \frac{h}{2},\;\; 0 \leqslant x \leqslant b, \;\;{\sigma _{zy}} = \mu \frac{{\partial u}}{{\partial n}} = {\tau _0} \end{split} \right\} $$ (20) 其中$ \mu = \dfrac{E}{{2(1 + \nu )}} $为剪切模量.
图4(a)给出了传统方法选点和直接法选点计算的位移结果与有限元(COMSOL)裂纹面上的位移结果对比, 可以看出用直接法计算的位移结果与有限元结果更为吻合. 采用COD方法, 对裂纹处的应力强度因子进行计算, 将位移代入式(21)算得的应力强度因子与有限元结果进行对比, 用线性拟合分析, 得出裂纹尖端处的应力强度因子分布如图4(b)所示. 有限元计算归一化应力强度因子结果为1.1812, 而应力强度因子手册[33]给出该问题理论解为1.1808, 误差精度达3.39 × 10−4. 由于有限元法计算结果具有很高的精度, 因此本文算例均采用有限元结果作为参考解. 用传统方法计算出的归一化裂纹强度因子为1.2475, 与有限元结果的相对误差为5.61 × 10−2, 而直接法计算出的结果为1.1578, 相对误差为1.98 × 10−2. 对于裂纹问题2%左右的误差在允许的范围内
$$ {K_{{\mathrm{III}}}} = \frac{E}{{8(1 + \nu )}}\sqrt {\frac{{2\text{π} }}{l}} \Delta u $$ (21) 其中$l$为计算位置到裂纹尖端处的距离, $ \Delta u $为裂纹面的位移.
接下来考虑不同裂纹长度, 分别再计算占模型边长1/4, 3/4裂纹长度的归一化应力强度因子如表1所示. 最后作出两种取局部点方法与有限元结果的相对误差如图5所示, 从相对误差可以得出直接法将边界节点法向导数计算简化为一维情况, 提高了计算精度. 因此选用直接法能够提升局部径向基函数配点法在计算反平面裂纹问题时的精度.
表 1 归一化应力强度因子结果Table 1. Normalized stress intensity factor resultsPercentage of crack lengths 1/4 1/2 3/4 RBF (traditional) 0.5619 1.2475 1.3943 RBF (direct) 0.5143 1.1578 1.2872 FEM 0.5287 1.1812 1.3178 2. 含裂纹声子晶体以及数值离散
2.1 含周期裂纹声子晶体
二维声子晶体是由晶格常数为a的正方形或三角形周期排列的无限圆柱体组成, 如图6所示. 内部散射体的横截面可以是任意的. 如果弹性波在垂直于柱体轴线(z轴)的平面(xOy面)内传播, 则出现反平面横波. 描述反平面弹性波运动的控制方程[28]可以表示为
![]()
图 6 二维声子晶体结构: (a)正方形晶格类型, (b)正方形单元格和(c)第一布里渊区Figure 6. 2D phononic crystal structure: (a) the square lattices, (b) corresponding square unit-cell and (c) the first Brillouin zone$$ {\mu _j}\left[\frac{{{\partial ^2}{u_j}({\boldsymbol{x}})}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_j}({\boldsymbol{x}})}}{{\partial {y^2}}}\right] = - {\rho _j}{\omega ^2}{u_j}({\boldsymbol{x}}),\;\; {\text{ }}j = 0, 1 $$ (22) 其中, $ {\mu _j} $和$ {\rho _j} $分别为剪切模量和密度; $ {u_j} $是沿$z$方向的位移, 下标为“0”的量表示散射体, 下标为“1”的量表示基体, $\omega $为特征频率.
根据Bloch定理, 弹性波在声子晶体结构中可以表示为周期函数. 在本文中, 由于声子晶体的周期特性, 单元格间的边界是施加周期性条件. 基体与散射体交界面的连续边界条件可以表示为以下方程
$$ {u_0}({{\boldsymbol{x}}_{{\varGamma _0}}}) = {u_1}({{\boldsymbol{x}}_{{\varGamma _0}}}),\quad {{\boldsymbol{x}}_{{\varGamma _0}}} \in {\varGamma _0} $$ (23) $$ {T_0}({{\boldsymbol{x}}_{{\varGamma _0}}}) = {T_1}({{\boldsymbol{x}}_{{\varGamma _0}}}), \quad{{\boldsymbol{x}}_{{\varGamma _0}}} \in {\varGamma _0} $$ (24) 其中, $ {T_0}({\boldsymbol{x}}) $和$ {T_1}({\boldsymbol{x}}) $是面力向量, 而在反平面弹性问题中的$ T({\boldsymbol{x}}) $可以表示成
$$ T({\boldsymbol{x}}) = \mu \frac{{\partial u({\boldsymbol{x}})}}{{\partial {\boldsymbol{n}}}} $$ (25) 其中, ${\boldsymbol{n}} = {({n_x}, {n_y})^{\mathrm{T}}}$是垂直于界面的单位法向量.
利用Bloch定理, 反平面弹性波周期边界条件可以表示成
$$ u({\boldsymbol{x}} + {\boldsymbol{a}}) = {{\mathrm{e}}^{{\mathrm{i}}{\boldsymbol{ka}}}}u({\boldsymbol{x}}) $$ (26) $$ T({\boldsymbol{x}} + {\boldsymbol{a}}) = {{\mathrm{e}}^{{\mathrm{i}}{\boldsymbol{ka}}}}T({\boldsymbol{x}}) $$ (27) 其中, ${\boldsymbol{k}} = ({k_x}, {k_y})$是Bloch波矢向量. $ u({\boldsymbol{x}}) $和$ T({\boldsymbol{x}}) $是满足Bloch周期条件的位移和面力. 且${\mathrm{i}} = \sqrt { - 1} $, ${\boldsymbol{a}} = {m_1}{{\boldsymbol{a}}_1} + {m_2}{{\boldsymbol{a}}_2}, {\boldsymbol{m}} = ({m_1}, {m_2}) \in {Z^2}$, ${{\boldsymbol{a}}_1}$和${{\boldsymbol{a}}_2}$是晶格的基本平移向量.
考虑第一布里渊区边界上的波矢向量${\boldsymbol{k}}$. 在单元格的边界上, 应用式(26)和式(27)的Bloch周期条件, 可得
$$ u({{\boldsymbol{x}}_{{\varGamma _1}}}) = u({{\boldsymbol{x}}_{{\varGamma _3}}}){{\mathrm{e}}^{{\mathrm{i}}{\boldsymbol{ka}}}},\quad u({{\boldsymbol{x}}_{{\varGamma _2}}}) = u({{\boldsymbol{x}}_{{\varGamma _4}}}){{\mathrm{e}}^{{\mathrm{i}}{\boldsymbol{ka}}}} $$ (28) $$ T({{\boldsymbol{x}}_{{\varGamma _1}}}) = T({{\boldsymbol{x}}_{{\varGamma _3}}}){{\mathrm{e}}^{{\mathrm{i}}{\boldsymbol{ka}}}},\quad T({{\boldsymbol{x}}_{{\varGamma _2}}}) = T({{\boldsymbol{x}}_{{\varGamma _4}}}){{\mathrm{e}}^{{\mathrm{i}}{\boldsymbol{ka}}}} $$ (29) 其中, $ {{\boldsymbol{x}}_{{\varGamma _i}}} \in {\varGamma _i}\;(i = 1, 2, 3, 4) $.
对于图5(b)所示的单胞晶格模型, 裂纹面在图示的红色位置. 裂纹面为自由边界, 考虑反平面自由边界条件可得
$$ {T_0}({{\boldsymbol{x}}_{{\varGamma _b}}}) = 0,\quad {T_0}({{\boldsymbol{x}}_{{\varGamma _c}}}) = 0 $$ (30) 其中, $ {T_0}({{\boldsymbol{x}}_{{\varGamma _b}}}) $和$ {T_0}({{\boldsymbol{x}}_{{\varGamma _c}}}) $分别是左右裂纹面的面力向量.
2.2 数值离散
对裂纹面自由边界条件, 通过式(14) ~ 式(16)对其进行数值离散, 可将式(30)表示为如下形式
$$ {\mu _0} \frac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{\varGamma _b}}})}}{{\partial {\boldsymbol{n}}}}\tilde {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _b}}}) = 0 $$ (31) $$ {\mu _0} \frac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{\varGamma _c}}})}}{{\partial {\boldsymbol{n}}}}\tilde {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _c}}}) = 0 $$ (32) 同理再将控制方程式(22)、基体与散射体交界面上的面力连续条件式(24)以及周期边界条件式(29)离散, 形成如下矩阵形式
$$\begin{split} &\left[ {\begin{array}{*{20}{c}} {\dfrac{{{\partial ^2}\tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{D_0}}})}}{{\partial {x^2}}} + \dfrac{{{\partial ^2}\tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{D_0}}})}}{{\partial {y^2}}}}&0 \\ 0&{\dfrac{{{\partial ^2}\tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{D{}_1}})}}{{\partial {x^2}}} + \dfrac{{{\partial ^2}\tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{D{}_1}})}}{{\partial {y^2}}}} \\ 0&{\dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{\varGamma _1}}})}}{{\partial {\boldsymbol{n}}}} - \dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{\varGamma 3}})}}{{\partial {\boldsymbol{n}}}}{{\mathrm{e}}^{ - {\mathrm{i}}{k_y}a}}} \\ 0&{\dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{\varGamma _2}}})}}{{\partial {\boldsymbol{n}}}} - \dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{\varGamma _4}}})}}{{\partial {\boldsymbol{n}}}}{{\mathrm{e}}^{ - {\mathrm{i}}{k_x}a}}} \\ {{\mu _0}\dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({\boldsymbol{x}}_{{\varGamma _0}}^0)}}{{\partial {\boldsymbol{n}}}}}&{ - {\mu _1}\dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({\boldsymbol{x}}_{{\varGamma _0}}^1)}}{{\partial {\boldsymbol{n}}}}} \\ {\dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{\varGamma _b}}})}}{{\partial {\boldsymbol{n}}}}}&0 \\ {\dfrac{{\partial \tilde {\boldsymbol{\varPsi}} ({{\boldsymbol{x}}_{{\varGamma _c}}})}}{{\partial {\boldsymbol{n}}}}}&0 \end{array}} \right]\cdot \\ &\qquad \left[ \begin{gathered} \tilde {\boldsymbol{u}}({{\boldsymbol{x}}_0}) \\ \tilde {\boldsymbol{u}}({{\boldsymbol{x}}_1}) \\ \end{gathered} \right] = - {\omega ^2}\left[ {\begin{array}{*{20}{c}} {\dfrac{{{\rho _0}}}{{{\mu _0}}}\delta _{{\mathrm{in}}}^0}&0 \\ 0&{\dfrac{{{\rho _1}}}{{{\mu _1}}}\delta _{{\mathrm{in}}}^1} \\ 0&0 \\ 0&0 \\ 0&0 \\ 0&0 \\ 0&0 \end{array}} \right]\left[ \begin{gathered} \tilde {\boldsymbol{u}}({{\boldsymbol{x}}_0}) \\ \tilde {\boldsymbol{u}}({{\boldsymbol{x}}_1}) \\ \end{gathered} \right] \end{split}$$ (33) 可以将式(33)简写为以下矩阵形式
$$ {\boldsymbol{KU}} = - {\omega ^2}{\boldsymbol{PU}} $$ (34) 其中, ${\boldsymbol{K}}$是和RBF相关的离散矩阵. ${\boldsymbol{P}}$是与给定初始条件以及质量密度$\rho $和剪切模量$\mu $相关的向量. 由于离散矩阵${\boldsymbol{K}}$不是方阵, 所以式(34)中的线性方程组不能直接求解, 因此我们将基体与散射体的交界面的位移连续条件和位移周期边界条件及裂纹位置的位移边界条件通过矩阵${\boldsymbol{K}}$的加减解析得到满足, 具体如下
$$ {\boldsymbol{K}} = [{{{\boldsymbol{K}}}_{{{D}_{0}}}},{{{\boldsymbol{K}}}_{{{D}_{1}}}},{{{\boldsymbol{K}}}_{{{\varGamma }_{1}}}},{{{\boldsymbol{K}}}_{{{\varGamma }_{2}}}},{{{\boldsymbol{K}}}_{{{\varGamma }_{3}}}}, {{{\boldsymbol{K}}}_{{{\varGamma }_{4}}}},{\boldsymbol{K}}_{{{\varGamma }_{0}}}^{0},{\boldsymbol{K}}_{{{\varGamma }_{0}}}^{1},{\boldsymbol{K}}_{{{\varGamma }_{b}}}^{0},{\boldsymbol{K}}_{{{\varGamma }_{c}}}^{0}] $$ (35) $$ {\boldsymbol{P}} = [{{\boldsymbol{P}}_{{D_0}}}, {{\boldsymbol{P}}_{{D_1}}}, {{\boldsymbol{P}}_{{\varGamma _1}}}, {{\boldsymbol{P}}_{{\varGamma _2}}}, {{\boldsymbol{P}}_{{\varGamma _3}}}, {{\boldsymbol{P}}_{{\varGamma _4}}}, {\boldsymbol{P}}_{{\varGamma _0}}^0, {\boldsymbol{P}}_{{\varGamma _0}}^1, {\boldsymbol{P}}_{{\varGamma _b}}^0, {\boldsymbol{P}}_{{\varGamma _c}}^0] $$ (36) $$ \begin{split} & {\boldsymbol{U}} = [{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_0}}}), {\boldsymbol{u}}({{\boldsymbol{x}}_{{D_1}}}), {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _1}}}), {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _2}}}), {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _3}}}), \\ &\qquad {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _4}}}), {\boldsymbol{u}}({\boldsymbol{x}}_{_{{\varGamma _0}}}^0), {\boldsymbol{u}}({\boldsymbol{x}}_{_{{\varGamma _0}}}^1), {\boldsymbol{u}}({\boldsymbol{x}}_{_{{\varGamma _b}}}^0), {\boldsymbol{u}}({\boldsymbol{x}}_{_{{\varGamma _c}}}^0){]^{\mathrm{T}}} \end{split} $$ (37) 其中, $ {{\boldsymbol{K}}_{{D_0}}} $, $ {{\boldsymbol{K}}_{{D_1}}} $, $ {{\boldsymbol{K}}_{{\varGamma _1}}} $, $ {{\boldsymbol{K}}_{{\varGamma _2}}} $, $ {{\boldsymbol{K}}_{{\varGamma _3}}} $, $ {{\boldsymbol{K}}_{{\varGamma _4}}} $, $ {\boldsymbol{K}}_{{\varGamma _0}}^0 $, $ {\boldsymbol{K}}_{{\varGamma _0}}^1 $, $ {\boldsymbol{K}}_{{\varGamma _b}}^0 $和$ {\boldsymbol{K}}_{{\varGamma _c}}^0 $分别是与$ {\boldsymbol{u}}({{\boldsymbol{x}}_{{D_0}}}) $, $ {\boldsymbol{u}}({{\boldsymbol{x}}_{{D_1}}}) $, $ {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _1}}}) $, $ {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _2}}}) $, $ {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _3}}}) $, $ {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _4}}}) $, $ {\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^0) $, $ {\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^1) $, $ {\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _b}}^0) $和$ {\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _c}}^0) $相关的矩阵列, ${\boldsymbol{P}}$矩阵同理可以拆分开成与${\boldsymbol{K}}$相似的矩阵. 因此式(34)的左边及右边可以写成
$$ \begin{split} & {\boldsymbol{KU}} = {{\boldsymbol{K}}_{{D_0}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_0}}}) + {{\boldsymbol{K}}_{{D_1}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_1}}}) + {{\boldsymbol{K}}_{{\varGamma _1}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _1}}}) +\\ &\qquad {{\boldsymbol{K}}_{{\varGamma _2}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _2}}}) + {{\boldsymbol{K}}_{{\varGamma _3}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _3}}}) + {{\boldsymbol{K}}_{{\varGamma _4}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _4}}}) + {\boldsymbol{K}}_{{\varGamma _0}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^0) +\\ &\qquad {\boldsymbol{K}}_{{\varGamma _0}}^1{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^1) + {\boldsymbol{K}}_{{\varGamma _b}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _b}}^0) + {\boldsymbol{K}}_{{\varGamma _c}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _c}}^0) \end{split} $$ (38) $$ \begin{split} & {\boldsymbol{PU}} = {{\boldsymbol{P}}_{{D_0}}}{\boldsymbol{u}}({x_{{D_0}}}) + {{\boldsymbol{P}}_{{D_1}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_1}}}) + {{\boldsymbol{P}}_{{\varGamma _1}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _1}}}) +\\ &\qquad {{\boldsymbol{P}}_{{\varGamma _2}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _2}}}) + {{\boldsymbol{P}}_{{\varGamma _3}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _3}}}) + {{\boldsymbol{P}}_{{\varGamma _4}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _4}}}) + {\boldsymbol{P}}_{{\varGamma _0}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^0) + \\ &\qquad {\boldsymbol{P}}_{{\varGamma _0}}^1{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^1) + {\boldsymbol{P}}_{{\varGamma _b}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _b}}^0) + {\boldsymbol{P}}_{{\varGamma _c}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _c}}^0) \end{split} $$ (39) 式(23)的基体-散射体界面上的位移连续性条件和式(28)中的位移周期性条件通过加减列来考虑, 其可以重写为以下形式
$$ \begin{split} & \tilde{{\boldsymbol{K}}}\tilde{{\boldsymbol{U}}} = {{\boldsymbol{K}}_{{D_0}}}{\boldsymbol{u}}({x_{{D_0}}}) + {{\boldsymbol{K}}_{{D_1}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_1}}}) + ({{\boldsymbol{K}}_{{\varGamma _1}}} +\\ &\qquad {{\boldsymbol{K}}_{{\varGamma _3}}}{{\mathrm{e}}^{{\mathrm{i}}{{\boldsymbol{k}}_y}{\boldsymbol{a}}}}){\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _1}}}) + ({{\boldsymbol{K}}_{{\varGamma _2}}} + {{\boldsymbol{K}}_{{\varGamma _4}}}{{\mathrm{e}}^{{\mathrm{i}}{{\boldsymbol{k}}_x}{\boldsymbol{a}}}}){\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _2}}}) + \\ &\qquad ({\boldsymbol{K}}_{{\varGamma _0}}^0 + {\boldsymbol{K}}_{{\varGamma _0}}^1){\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^0) + {\boldsymbol{K}}_{{\varGamma _b}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _b}}^0) + {\boldsymbol{K}}_{{\varGamma _c}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _c}}^0)\end{split} $$ (40) $$ \begin{split} & \tilde{{\boldsymbol{P}}}\tilde{{\boldsymbol{U}}} = {{\boldsymbol{P}}_{{D_0}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_0}}}) + {{\boldsymbol{P}}_{{D_1}}}{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_1}}}) + ({{\boldsymbol{P}}_{{\varGamma _1}}} + \\ &\qquad {{\boldsymbol{P}}_{{\varGamma _3}}}{{\mathrm{e}}^{{\mathrm{i}}{k_y}{\boldsymbol{a}}}}){\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _1}}}) + ({{\boldsymbol{P}}_{{\varGamma _2}}} + {{\boldsymbol{P}}_{{\varGamma _4}}}{{\mathrm{e}}^{{\mathrm{i}}{k_x}{\boldsymbol{a}}}}){\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _2}}}) + \\ &\qquad(K_{{\varGamma _0}}^0 + {\boldsymbol{P}}_{{\varGamma _0}}^1){\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^0) + {\boldsymbol{P}}_{{\varGamma _b}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _b}}^0) + {\boldsymbol{P}}_{{\varGamma _c}}^0{\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _c}}^0) \end{split} $$ (41) 然后, 可以将式(34)改写为如下方程
$$ \tilde {\boldsymbol{K}}\tilde {\boldsymbol{U}} = - {\omega ^2}\tilde {\boldsymbol{P}}\tilde {\boldsymbol{U}} $$ (42) 其中
$$ \begin{split} & \tilde {\boldsymbol{K}} = \left[ {{\boldsymbol{K}}_{{D_0}}}, {{\boldsymbol{K}}_{{D_1}}}, {{\boldsymbol{K}}_{{\varGamma _1}}} + {{\boldsymbol{K}}_{{\varGamma _3}}}{{\mathrm{e}}^{{\mathrm{i}}{k_y}{\boldsymbol{a}}}}, {{\boldsymbol{K}}_{{\varGamma _2}}} + {{\boldsymbol{K}}_{{\varGamma _4}}}{{\mathrm{e}}^{{\mathrm{i}}{k_x}{\boldsymbol{a}}}},\right. \\ &\qquad \left.{\boldsymbol{K}}_{{\varGamma _0}}^0 + {\boldsymbol{K}}_{{\varGamma _0}}^1, {\boldsymbol{K}}_{{\varGamma _b}}^0, {\boldsymbol{K}}_{{\varGamma _c}}^0 \right] \\ & \tilde {\boldsymbol{P}} = \left[ {{\boldsymbol{P}}_{{D_0}}}, {{\boldsymbol{P}}_{{D_1}}}, {{\boldsymbol{P}}_{{\varGamma _1}}} + {{\boldsymbol{P}}_{{\varGamma _3}}}{{\mathrm{e}}^{{\mathrm{i}}{k_y}{\boldsymbol{a}}}}, {{\boldsymbol{P}}_{{\varGamma _2}}} + {{\boldsymbol{P}}_{{\varGamma _4}}}{{\mathrm{e}}^{{\mathrm{i}}{k_x}{\boldsymbol{a}}}},\right.\\ &\qquad \left. {\boldsymbol{P}}_{{\varGamma _0}}^0 + {\boldsymbol{P}}_{{\varGamma _0}}^1, {\boldsymbol{P}}_{{\varGamma _b}}^0, {\boldsymbol{P}}_{{\varGamma _c}}^0 \right] \\ & \tilde {\boldsymbol{U}} = {\left[ {{\boldsymbol{u}}({{\boldsymbol{x}}_{{D_0}}}), {\boldsymbol{u}}({{\boldsymbol{x}}_{{D_1}}}), {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _1}}}), {\boldsymbol{u}}({{\boldsymbol{x}}_{{\varGamma _2}}}), {\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _0}}^0), {\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _b}}^0), {\boldsymbol{u}}({\boldsymbol{x}}_{{\varGamma _c}}^0)} \right]^{\mathrm{T}}} \end{split} $$ 通过这种方式可以得到$\tilde {\boldsymbol{K}}$并且解析满足位移的连续条件, 显著提高了直接法局部径向基函数配点法的稳定性. 在求解特征值问题中, 我们直接调用Matlab自带的“eigs”程序, 实现了对特征频率的求解.
3. 数值结果与讨论
形状参数${c_s}$是影响局部径向基函数配点法结果的一个非常重要的因素, 它与布点的距离和问题的光滑程度有关系, 可以通过样本函数去确定形状参数的有效性[34], 本文中形状参数${c_s} = 1$, 数值计算均采用均匀节点分布.
本节将环氧树脂基体和金(Au)散射体所组成的声子晶体及由环氧树脂基体和铝(Al)散射体所组成的声子晶体作为研究对象. 对裂纹长度和裂纹位置以及不同散射体形状对声子晶体带隙结构的影响进行研究. 将计算结果与有限元软件(COMSOL)声学模块的计算结果对比, 验证了所使用局部径向基函数配点法的准确性和高效性.
首先我们研究嵌在环氧树脂基体中的金(Au)散射体声子晶体. 组分材料的密度和剪切模量如表2所示, 所考虑的声子晶体的声阻抗比为$Z = {\rho _1}{c_1}/ ({\rho _0}{c_0}) = 17. 64$. 其中${c_j} = \sqrt {{\mu _j}/{\rho _j}} $表示的是材料的波速.
表 2 材料参数Table 2. Material parametersMedium $\rho /({\mathrm{kg}}\cdot{{\mathrm{m}}^{-3}})$ $\mu /{\mathrm{Pa}}$ Au 19500 2.99 × 1010 Al 2799 2.68 × 1010 epoxy 1180 1.59 × 109 基于声子晶体模型, 我们利用局部径向基函数配点法计算完整声子晶体无裂纹情况下的能带结构, 其相应的布点情况如图7(a), 其中蓝色点表示的是基体、红色点表示的是散射体、绿色点表示的是基体与散射体的交界面. 角点因为采用了直接法, 让它作为其中一条边界的节点即可, 在数值离散的时候局部点选取和定义一致. 比较了数值结果与有限元软件(COMSOL)计算的结果, 并得出能带结构对比图(如图7(b)所示), 可以得出, 该数值结果与有限元结果吻合很好, 因此用该方法计算声子晶体是可靠的.
![]()
图 7 方形晶格中方形散射体(金/环氧树脂): (a)节点分布和(b)能带结构Figure 7. (a) Node distribution and (b) band structure of square scatterer (Au/epoxy) in the square lattice3.1 方形散射体
3.1.1 金/环氧树脂
本节研究方形散射体, 材料为金(Au). 方形散射体的面积为0.4 × 0.4, 单胞面积为1 × 1, 填充率为0.16. 布点方式如图8所示, 黑色点为裂纹面位置, 晶格中的节点数为51 × 51. 内部点用传统方法选取邻近的9个相邻节点, 裂纹面使用直接法沿法向选取5个相邻节点, 基体与散射体交界面的连续边界条件和基体的周期性边界条件同样采用直接法并沿法向选取5个相邻节点. 当裂纹从散射体的中间位置扩展且长度占散射体边长的1/4时, 相应的有限元模型具有2410个三角单元、4525个自由度.
![]()
图 8 不同裂纹长度和位置方形散射体声子晶体能带结构(金/环氧树脂)Figure 8. Band structure of phononic crystals of square scatterers with different crack lengths and positions (Au/epoxy)裂纹从散射体的中间向两边扩展, 裂纹的长度分别占散射体边长的1/4, 1/2, 3/4和整个散射体边长. 从图8能带图可以得出, 局部径向基函数配点法与有限元所计算出的能带结构吻合良好, 无频率缺失现象.
研究在同种材料及散射体形状下, 考虑裂纹位置从基体与散射体交界面开始扩展. 裂纹的长度分别占散射体边长的1/4, 1/2, 3/4和整个散射体边长. 从图8对应能带结构的结果可以看出数值结果也与有限元结果吻合良好.
随后分析上述裂纹对完整声子晶体能带结构的影响情况. 结合图8和图9可以看出, 随着裂纹从中间位置向外逐渐扩展, 第1能带带隙开始变窄, 第2能带带隙也逐渐变窄. 当裂纹扩展至3/4的散射体边长时第2能带消失, 直至裂纹扩展贯穿整个散射体时能带结构中出现第3条能带带隙. 而当裂纹从基体与散射体交界面开始扩展时, 第1和第2能带带隙同样会逐渐变窄. 当裂纹扩展至3/4散射体边长时出现第3带隙. 可以看出原胞裂纹长度的变化对低频能带带隙的影响较小, 但对中低频能带带隙的影响较大. 另外新增带隙随着裂纹扩展而变宽.
![]()
图 9 不同裂纹长度对比完整声子晶体能带结构(金/环氧树脂)Figure 9. Comparison of complete phononic crystal band structure with different crack lengths (Au/epoxy)3.1.2 铝/环氧树脂
本节再用铝(Al)散射继续体验证该方法的有效性, 所考虑的声子晶体的声阻抗比为$Z = {\rho _1}{c_1}/({\rho _0}{c_0}) = 6. 32$, 其他条件和(1)中算例相同. 从图10可以看出, 不论是裂纹从中间往两边扩展, 还是从基体与散射体交界面扩展, 该方法计算出的数值结果与有限元结果对比吻合度都非常好, 并且未出现明显的频率缺失现象. 与上一节相比, 可以看出散射体材料的变化对数值结果影响不大.
![]()
图 10 不同裂纹长度和位置方形散射体声子晶体能带结构(铝/环氧树脂)Figure 10. Band structure of phononic crystals of square scatterers with different crack lengths and positions (Al/epoxy)接着分析裂纹对铝/环氧树脂声子晶体能带结构的影响. 与金散射体的情况进行对比, 由于声子晶体声阻抗比的大幅减小, 对应声子晶体的能带带隙也明显减少且急剧变窄, 但不论裂纹从中间还是交界面扩展, 第1能带带隙会逐渐变窄, 而当裂纹扩展至贯穿散射体时第1能带完全消失.
3.2 圆形散射体
3.2.1 金/环氧树脂
本节考虑圆形散射体、材料为金(Au). 圆形散射体半径为0.3, 单胞面积为1 × 1, 填充率为0.283. 采用均匀节点分布(41 × 41). 裂纹位置使用直接法沿法向选取5个相邻节点. 当裂纹从中间扩展且长度占散射体边长的1/4时, 相应有限元模型具有2224个三角单元和4889个自由度.
从图11可以看出, 在散射体形状从方形变成圆形时, 局部径向基函数配点法的结果与有限元的结果仍然有非常高的吻合度. 与3.1.1节的算例相比, 可以看出散射体形状对计算结果的符合程度几乎没有影响. 随着裂纹扩展, 声子晶体的能带带隙范围逐渐减小, 第1、第2和第3带隙逐渐变窄; 当裂纹扩展至占1/2散射体边长时出现新增带隙, 并且新增带隙会随着裂纹的扩展逐渐变宽.
![]()
图 11 不同裂纹长度和位置圆形散射体声子晶体能带结构(金/环氧树脂)Figure 11. Band structure of phononic crystals of circular scatterers with different crack lengths and positions (Au/epoxy)3.2.2 铝/环氧树脂
现在分析圆形散射体, 材料为铝(Al), 研究不同裂纹长度及扩展位置的影响并与有限元结果进行比对. 从图12可以得出, 裂纹不管是从中间向外扩展还是从基体与散射体交界面扩展, 局部径向基函数配点法结果和有限元结果吻合良好, 这也更进一步说明不同材料对该数值方法的适用性没有影响. 并且随着裂纹的扩展, 第1带隙和第2带隙都在逐渐变窄, 直到裂纹扩展至3/4的散射体边长时, 第2带隙开始消失.
将3.1.1节和3.2.1节、3.1.2节和3.2.2节的算例结果进行对比可知, 随着裂纹的扩展, 金散射体的声子晶体带隙数量会增加, 且新增带隙会随着裂纹的扩展变宽, 而铝散射体的声子晶体能带带隙数量会相应减少.
![]()
图 12 不同裂纹长度和位置圆形散射体声子晶体能带结构(铝/环氧树脂)Figure 12. Band structure of phononic crystals of circular scatterers with different crack lengths and positions (Al/epoxy)3.3 计算效率分析
本节主要分析和讨论局部径向基函数配点法的计算精度及效率, 表3给出了不同散射体形状以及裂纹位置的数值结果与有限元结果的相对误差, 相对误差定义为$ \displaystyle\sum {\left| {{E_r} - {E_f}} \right|} \Bigr/\displaystyle\sum {{E_f}} $, 其中, $ {E_r} $表示局部径向基函数配点法结果, $ {E_f} $表示有限元结果. 从表中数据可以得出, 数值结果和有限元结果的相对误差在10−4 ~ 10−3之间, 说明该数值方法具有较为良好的精度.
表 3 计算效率和精度比较Table 3. CPU time and accuracy comparisonScatterer form
Crack locationSquare Circular 
RBF number of nodes 2601 2601 1771 1771 time cost/s 5.27 5.18 3.39 3.37 FEM degrees of freedom 4525 4261 4889 4986 time cost/s 96 85 54 67 Comparison errors/10−4 9.80 1.34 21.2 25.2 time saving/% 94.5 93.9 93.7 95.0 表3还给出了不同方法表现的比较, 包括计算时间等, 无论是方形散射体还是圆形散射体, 无网格局部径向基函数配点法都可以大大提高计算效率, 在保持足够的精度下还能节省94%左右的时间.
从图13(a)不同形状参数的相对误差中可以看出, 当形状参数取在0.3 ~ 5之间时对计算精度影响可以忽略不计, 而从图13(b)不同节点数的相对误差可以看出, 当节点数超过1000时可以获得较好的计算精度.
4. 结论
本文基于直接法局部径向基函数配点法计算含裂纹声子晶体的能带结构, 并给出了该方法的控制方程、周期性边界条件、交界面连续性边界条件和裂纹面自由边界条件的数值计算全过程. 通过与有限元法结果比对, 验证了直接法局部径向基函数配点法在计算含裂纹声子晶体能带结构的有效性, 并分析了裂纹对声子晶体能带结构的影响. 得出主要结论如下.
(1)利用直接法处理边界条件法向导数的计算可以大大提高局部径向基函数配点法计算裂纹声子晶体能带结构的稳定性.
(2)通过对不同材料声阻抗比和散射体形状(方形、圆形)的结果分析, 含裂纹声子晶体能带结构的改进局部径向基函数配点法有较强的适用性.
(3)改进局部径向基函数配点法对不同的形状参数以及节点数有较好的收敛性, 且相比有限元法可以大大提升含裂纹声子晶体能带结构的计算效率.
(4)随着裂纹的扩展, 能带结构带隙逐渐变窄, 当扩展到一定程度时, 金散射体的声子晶体带隙数量会增加, 新增带隙会随着裂纹的扩展而变宽; 但铝散射体的声子晶体能带带隙数量会减少. 对设计可调制的声子晶体带隙声学器件具有一定的理论参考价值和实际的指导意义.
局部径向基函数配点法在复杂边界条件下或者复杂几何形状裂纹中, 理论上有较好的精度, 且该方法可以有效模拟裂纹扩展的问题[31]. 在非均匀布点的情况下, 每个局部区域可能有非常近的点, 导致计算结果不稳定. 需要通过调整形状参数, 或者舍去一些临近点等其他技术来实现其稳定性, 这个也是一个研究热点问题. 未来进一步的工作还将继续研究局部径向基函数配点法在声子晶体复杂边界条件下或者复杂裂纹情况下计算结果的精度和效率, 并考虑非均匀布点情况下计算的结果, 充分发挥无网格法没有网格畸变和网格依赖特性的优势[35].
-
![]()
图 2 (a)传统法和(b)直接法选点示意图
Figure 2. (a) The traditional method and (b) the direct method of node selection
![]()
图 6 二维声子晶体结构: (a)正方形晶格类型, (b)正方形单元格和(c)第一布里渊区
Figure 6. 2D phononic crystal structure: (a) the square lattices, (b) corresponding square unit-cell and (c) the first Brillouin zone
![]()
图 7 方形晶格中方形散射体(金/环氧树脂): (a)节点分布和(b)能带结构
Figure 7. (a) Node distribution and (b) band structure of square scatterer (Au/epoxy) in the square lattice
![]()
图 8 不同裂纹长度和位置方形散射体声子晶体能带结构(金/环氧树脂)
Figure 8. Band structure of phononic crystals of square scatterers with different crack lengths and positions (Au/epoxy)
![]()
图 9 不同裂纹长度对比完整声子晶体能带结构(金/环氧树脂)
Figure 9. Comparison of complete phononic crystal band structure with different crack lengths (Au/epoxy)
![]()
图 10 不同裂纹长度和位置方形散射体声子晶体能带结构(铝/环氧树脂)
Figure 10. Band structure of phononic crystals of square scatterers with different crack lengths and positions (Al/epoxy)
![]()
图 11 不同裂纹长度和位置圆形散射体声子晶体能带结构(金/环氧树脂)
Figure 11. Band structure of phononic crystals of circular scatterers with different crack lengths and positions (Au/epoxy)
![]()
图 12 不同裂纹长度和位置圆形散射体声子晶体能带结构(铝/环氧树脂)
Figure 12. Band structure of phononic crystals of circular scatterers with different crack lengths and positions (Al/epoxy)
表 1 归一化应力强度因子结果
Table 1 Normalized stress intensity factor results
Percentage of crack lengths 1/4 1/2 3/4 RBF (traditional) 0.5619 1.2475 1.3943 RBF (direct) 0.5143 1.1578 1.2872 FEM 0.5287 1.1812 1.3178 
下载: 导出CSV
表 2 材料参数
Table 2 Material parameters
Medium $\rho /({\mathrm{kg}}\cdot{{\mathrm{m}}^{-3}})$ $\mu /{\mathrm{Pa}}$ Au 19500 2.99 × 1010 Al 2799 2.68 × 1010 epoxy 1180 1.59 × 109
下载: 导出CSV
表 3 计算效率和精度比较
Table 3 CPU time and accuracy comparison
Scatterer form
Crack locationSquare Circular RBF number of nodes 2601 2601 1771 1771 time cost/s 5.27 5.18 3.39 3.37 FEM degrees of freedom 4525 4261 4889 4986 time cost/s 96 85 54 67 Comparison errors/10−4 9.80 1.34 21.2 25.2 time saving/% 94.5 93.9 93.7 95.0
下载: 导出CSV
-
[1] Kushwaha MS, Halevi P, Dobrzynski L, et al. Acoustic band structure of periodic elastic composites. Physical Review Letters, 1993, 71(13): 2022-2025 doi: 10.1103/PhysRevLett.71.2022
[2] 卢一铭, 曹东兴, 申永军等. 局域共振型声子晶体板缺陷态带隙及其俘能特性研究. 力学学报, 2021, 53(4): 1114-1123 (Lu Yiming, Cao Dongxing, Shen Yongjun, et al. Study on the bandgaps of defect states and application of energy harvesting of local resonant phononic crystal plate. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(4): 1114-1123 (in Chinese) doi: 10.6052/0459-1879-20-436 Lu Yiming, Cao Dongxing, Shen Yongjun, et al. Study on the bandgaps of defect states and application of energy harvesting of local resonant phononic crystal plate. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(4): 1114-1123 (in Chinese) doi: 10.6052/0459-1879-20-436
[3] 邱克鹏, 陈智谋, 张建刚等. 基于形状记忆合金声子晶体的带隙优化设计. 力学学报, 2023, 55(6): 1278-1287 (Qiu Kepeng, Chen Zhimou, Zhang Jiangang, et al. Bandgap optimization design of phononic crystals based on shape memory alloy. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(6): 1278-1287 (in Chinese) doi: 10.6052/0459-1879-23-024 Qiu Kepeng, Chen Zhimou, Zhang Jiangang, et al. Bandgap optimization design of phononic crystals based on shape memory alloy. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(6): 1278-1287 (in Chinese) doi: 10.6052/0459-1879-23-024
[4] Chen LY, Guo YJ, Yi H. Optimization study of bandgaps properties for two-dimensional chiral phononic crystals base on lightweight design. Physics Letters A, 2021, 388: 127054 doi: 10.1016/j.physleta.2020.127054
[5] 韩东海, 张广军, 赵静波等. 新型 Helmholtz 型声子晶体的低频带隙及隔声特性. 物理学报, 2022, 71(11): 114301 (Han Donghai, Zhang Guangjun, Zhao Jingbo, et al. Low-frequency bandgaps and sound isolation characteristics of a novel Helmholtz-type phononic crystal. Acta Physica Sinica, 2022, 71(11): 114301 (in Chinese) doi: 10.7498/aps.71.20211932 Han Donghai, Zhang Guangjun, Zhao Jingbo, et al. Low-frequency bandgaps and sound isolation characteristics of a novel Helmholtz-type phononic crystal. Acta Physica Sinica, 2022, 71(11): 114301 (in Chinese) doi: 10.7498/aps.71.20211932
[6] 丁兰, 丁彪, 吴巧云等. 含双自由度周期振子的平行并联梁弯曲振动带隙特性. 工程力学, 2023, 40(10): 1-10, 57 (Ding Lan, Ding Biao, Wu Qiaoyun, et al. Flexural vibration bandgap characteristics of double beams in parallel with oscillators with two degrees of freedom. Engineering Mechanics, 2023, 40(10): 1-10, 57 (in Chinese) doi: 10.6052/j.issn.1000-4750.2022.01.0086 Ding Lan, Ding Biao, Wu Qiaoyun, et al. Flexural vibration bandgap characteristics of double beams in parallel with oscillators with two degrees of freedom. Engineering Mechanics, 2023, 40(10): 1-10, 57 (in Chinese) doi: 10.6052/j.issn.1000-4750.2022.01.0086
[7] 金亚斌, 温治辉. 二维共振结构板波动调控研究进展. 力学季刊, 2022, 43(1): 1 (Jin Yabin, Wen Zhihui. Research progress in wave modulation of two-dimensional resonant structured plates. Chinese Quarterly of Mechanics, 2022, 43(1): 1 (in Chinese) Jin Yabin, Wen Zhihui. Research progress in wave modulation of two-dimensional resonant structured plates. Chinese Quarterly of Mechanics, 2022, 43(1): 1 (in Chinese)
[8] Kafesaki M, Economou EN. Multiple-scattering theory for three-dimensional periodic acoustic composites. Physical Review B, 1999, 60(17): 11993-12001 doi: 10.1103/PhysRevB.60.11993
[9] Cebrecos A, Krattiger D, Sánchez-Morcillo VJ, et al. The finite-element time-domain method for elastic band-structure calculations. Computer Physics Communications, 2019, 238: 77-87 doi: 10.1016/j.cpc.2018.12.016
[10] Li FL, Zhang CZ, Wang YS. Analysis of the effects of viscosity on the SH-wave band-gaps of 2D viscoelastic phononic crystals by Dirichlet-to-Neumann map method. International Journal of Mechanical Sciences, 2021, 195: 106225 doi: 10.1016/j.ijmecsci.2020.106225
[11] 李凤莲, 吕梅. 基于边界元法的非完好界面三角晶格声子晶体带隙计算. 人工晶体学报, 2020, 49(1): 27-32 (Li Fenglian, Lyu Mei. Band-gap calculations of phononic crystals in a triangular lattice with imperfect interfaces based on the boundary element method. Journal of Synthetic Crystals, 2020, 49(1): 27-32 (in Chinese) doi: 10.3969/j.issn.1000-985X.2020.01.004 Li Fenglian, Lyu Mei. Band-gap calculations of phononic crystals in a triangular lattice with imperfect interfaces based on the boundary element method. Journal of Synthetic Crystals, 2020, 49(1): 27-32 (in Chinese) doi: 10.3969/j.issn.1000-985X.2020.01.004
[12] 王翌伟, 徐晓美, 林萍. 基于COMSOL局域共振声子晶体薄板振动带隙研究. 噪声与振动控制, 2022, 42(3): 73-79 (Wang Yiwei, Xu Xiaomei, Lin Ping. Study on vibration band gaps of local resonant phononic crystal thin plates based on COMSOL. Niose and Vibration Contral, 2022, 42(3): 73-79 (in Chinese) doi: 10.3969/j.issn.1006-1355.2022.03.013 Wang Yiwei, Xu Xiaomei, Lin Ping. Study on vibration band gaps of local resonant phononic crystal thin plates based on COMSOL. Niose and Vibration Contral, 2022, 42(3): 73-79 (in Chinese) doi: 10.3969/j.issn.1006-1355.2022.03.013
[13] 曹蕾蕾, 朱旺, 武建华等. 基于人工神经网络的声子晶体逆向设计. 力学学报, 2021, 53(7): 1992-1998 (Cao Leilei, Zhu Wang, Wu Jianhua, et al. Inverse design of phononic crystals by artificial neural networks. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1992-1998 (in Chinese) doi: 10.6052/0459-1879-21-142 Cao Leilei, Zhu Wang, Wu Jianhua, et al. Inverse design of phononic crystals by artificial neural networks. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1992-1998 (in Chinese) doi: 10.6052/0459-1879-21-142
[14] An XY, Fan HL, Zhang CZ. Elastic wave and vibration bandgaps in planar square metamaterial-based lattice structures. Journal of Sound and Vibration, 2020, 475(6): 1-15
[15] Chen W, Tanaka M. A meshless, integration-free, and boundary-only RBF technique. Computers & Mathematics with Applications, 2002, 43(3-5): 379-391
[16] Golberg MA, Chen CS, Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM. Engineering Analysis with Boundary Elements, 1999, 23(4): 285-296 doi: 10.1016/S0955-7997(98)00087-3
[17] Fu ZJ, Chen W, Yang HT. Boundary particle method for Laplace transformed time fractional diffusion equations. Journal of Computational Physics, 2013, 235: 52-66 doi: 10.1016/j.jcp.2012.10.018
[18] Chen W. New RBF collocation methods and kernel RBFs with applications//Griebel M, Schweitzer MA, eds. Meshfree Methods for Partial Differential Equations. Berlin, Heidelberg: Springer Berlin Heidelberg, 2003: 75-86
[19] Zhang X, Song KZ, Lu MW, et al. Meshless methods based on collocation with radial basis functions. Computational Mechanics, 2000, 26: 333-343 doi: 10.1007/s004660000181
[20] Chen JS, Wang LH, Hu HY, et al. Subdomain radial basis collocation method for heterogeneous media. International Journal for Numerical Methods in Engineering, 2009, 80(2): 163-190 doi: 10.1002/nme.2624
[21] Chen CS, Karageorghis A, Zheng H. Improved RBF collocation methods for fourth order boundary value problems. Communications in Computational Physics, 2020, 27(5): 1530-1549 doi: 10.4208/cicp.OA-2019-0163
[22] 吴宗敏. 径向基函数、散乱数据拟合与无网格偏微分方程数值解. 工程数学学报, 2002, 19(2): 1-12 (Wu Zongmin. Radial basis function scattered data interpolation and the meshless method of numerical solution of PDEs. Journal of Engineering Matehematics, 2002, 19(2): 1-12 (in Chinese) Wu Zongmin. Radial basis function scattered data interpolation and the meshless method of numerical solution of PDEs. Journal of Engineering Matehematics, 2002, 19(2): 1-12 (in Chinese)
[23] Wendland H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Advances in Computational Mathematics, 1995, 4: 389-396 doi: 10.1007/BF02123482
[24] 韩向科, 钱若军, 苏波等. 基于紧支径向基函数的流固交互作用数据传递. 同济大学学报: 自然科学版, 2011, 39(1): 48-52 (Han Xiangke, Qian Ruojun, Su Bo, et al. Data exchange method for fluid-structure interaction based on interpolation algorithm adopting compactly supported radial based function. Journal of Tongji Unversity: Nature Science, 2011, 39(1): 48-52 (in Chinese) Han Xiangke, Qian Ruojun, Su Bo, et al. Data exchange method for fluid-structure interaction based on interpolation algorithm adopting compactly supported radial based function. Journal of Tongji Unversity: Nature Science, 2011, 39(1): 48-52 (in Chinese)
[25] Gelas A, Bernard O, Friboulet D, et al. Compactly supported radial basis functions based collocation method for level-set evolution in image segmentation. IEEE Transactions on Image Processing, 2007, 16(7): 1873-1887 doi: 10.1109/TIP.2007.898969
[26] Šarler B, Vertnik R. Meshfree explicit local radial basis function collocation method for diffusion problems. Computers and Mathematics with Applications, 2006, 51(8): 1269-1282 doi: 10.1016/j.camwa.2006.04.013
[27] 刘华雩, 高效伟, 范伟龙. 分区有限线法及其在复合结构热应力分析中的应用. 力学学报, 2023, 55(6): 1394-1406 (Liu Huayu, Gao Xiaowei, Fan Weilong. Zonal finite line method and its applications in analyzing thermal stress of compositestructures. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(6): 1394-1406 (in Chinese) doi: 10.6052/0459-1879-23-003 Liu Huayu, Gao Xiaowei, Fan Weilong. Zonal finite line method and its applications in analyzing thermal stress of compositestructures. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(6): 1394-1406 (in Chinese) doi: 10.6052/0459-1879-23-003
[28] Zheng H, Zhang CZ, Wang YS, et al. A meshfree local RBF collocation method for anti-plane transverse elastic wave propagation analysis in 2D phononic crystals. Journal of Computational Physics, 2016, 305(1): 997-1014
[29] Zheng H, Zhou CB, Yan DJ, et al. A meshless collocation method for band structure simulation of nanoscale phononic crystals based on nonlocal elasticity theory. Journal of Computational Physics, 2020, 408(5): 1-17
[30] Yan XB, Zheng H, Yan DJ. Analysis of the band structure of transient in-plane elastic waves based on the localized radial basis function collocation method. Applied Mathematical Modelling, 2024, 125: 468-484 doi: 10.1016/j.apm.2023.09.002
[31] Wen PH, Aliabadi MH. Meshless method with enriched radial basis functions for fracture mechanics. Structural Durability & Health Monitoring, 2007, 3(2): 107
[32] Hardy RL. Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Research, 1971, 76(8): 1905-1915 doi: 10.1029/JB076i008p01905
[33] Sih GC. Handbook of Stress-Intensity Factors. Bethlehem, Pennsylvania: Lehigh University, 1973
[34] Zheng H, Yao G, Kuo LH, et al. On the selection of a good shape parameter of the localized method of approximated particular solutions. Advances in Applied Mathematics and Mechanics, 2018, 10(4): 896-911 doi: 10.4208/aamm.OA-2017-0167
[35] 王莉华, 阮剑武. 配点型无网格法理论和研究进展. 力学季刊, 2021, 42(4): 613-632 (Wang Lihua, Ruan Jianwu. Theory and research progress of the collocation-type meshfree methods. Chinese Quarterly of Mechanics, 2021, 42(4): 613-632 (in Chinese) Wang Lihua, Ruan Jianwu. Theory and research progress of the collocation-type meshfree methods. Chinese Quarterly of Mechanics, 2021, 42(4): 613-632 (in Chinese)
计量
- 文章访问数: 133
- HTML全文浏览量: 23
- PDF下载量: 42
