有限大八次准晶薄板刚性夹杂问题的应力分析

Stress Analysis of Plane Problems in Two-Dimensional Octagonal Quasicrystals Containing an Elliptic Rigid Inclusion

摘要: 准晶因其性能良好而广泛应用于发动机等设备的表面涂层中。由于准晶材料非常脆，因而对准晶部件应力集中处的应力场分析已引起广泛关注。本文利用复边界元法研究八次对称二维准晶中的椭圆形刚性夹杂问题。首先，通过推广的Stroh公式推导出在集中点力作用下，含椭圆刚性夹杂八次对称二维准晶平面弹性问题的Green函数。其次，基于不计体力的平衡方程和椭圆形刚性夹杂问题所对应的边界条件，构建边界积分方程，最后，通过Guass数值积分公式离散该边界积分方程并求解，分别获得了声子场和相位子场的孔边应力值。数值实例讨论了椭圆形刚性夹杂所引起的孔边应力变化，与含孔洞问题相比，内部夹杂的存在使基体孔边应力集中处的应力值减弱。

Abstract: Quasicrystals are widely used in surface coatings for devices such as engines due to their excellent properties. Because quasicrystal materials are very brittle, the analysis of the stress field at stress concentration areas of quasicrystal components has attracted widespread attention. The plane elastic problem of octagonal symmetric two-dimensional quasicrystals containing an elliptical rigid inclusion is considered by using the complex boundary element method. First, Green’s functions are obtained utilizing the extended Stroh formalism under concentrated force. Second, based on the equilibrium equation satisfying body force free and the boundary conditions corresponding to the elliptical rigid inclusion problem, a boundary integral equation is constructed, which is discretized and solved using the Guass’s formula of numerical integration. The stresses of the phonon field and the phason field on the boundary of an elliptic hole are obtained. The effect of the elliptical rigid inclusion on the stress is also discussed by comparing with the problem of containing an elliptic hole, and the presence of the internal inclusion weakens the stress values at the stress concentration points.

1. 引言

随着材料科学技术的飞速进步，高强度合金、先进的陶瓷材料、纤维增强复合材料以及颗粒增强复合材料，因其高硬度、耐磨、耐久和耐高温等卓越的性能，在军事、航空航天、汽车制造、电子、能源和生物医学等领域得到了广泛应用。在这些复合材料中，刚性夹杂作为一种常见的微观结构特征，对材料的整体力学性能产生显著影响，夹杂问题已引起广泛关注和研究[1]-[11]。

文献[1]利用保角变换技术，分析了含椭圆夹杂在远场受声子场和相子场的反平面机械载荷和平面内电载荷的无限一维六方压电准晶的线性压电弹性问题，获得了复势的封闭解以及基体和夹杂中的所有场量。文献[2]研究含双周期分布的圆形刚性夹杂在无穷远受纵向剪切的弹性平面问题，采用复变函数方法，导出了圆形刚性夹杂双周期分布的界面应力解析表达式。文献[3]研究了含椭圆形刚性夹杂的颗粒增强复合材料的反平面问题，构造复应力函数，采用坐标变换和复变函数的保角映射方法，推导出在无穷远双向均匀剪切，椭圆刚性夹杂的界面应力解析表达式。文献[4]首先利用Riemann-Schwartz解析延拓技术和复势函数的奇性主部分析法，研究了十次对称二维准晶在无穷远拉伸作用和集中力作用下含刚性圆弧夹杂的平面问题，得到在无穷远拉伸作用和集中力作用下含圆形夹杂界面刚性线问题的一般解答，分析了夹杂半径对热应力集中现象的影响。文献[5]应用复变函数的Faber级数展开方法，研究了在无穷远拉伸作用下含椭圆形刚性夹杂的压电材料平面问题，给出了问题的封闭解。文献[6]利用Muskhelishvili复变量法和保角映射技术，研究了含椭圆夹杂热电材料的二维问题，给出了电势、温度和应力分量的封闭解。文献[7]利用Stroh公式、保角变换和叠加原理，研究了在远场任意均匀力、电、热流载荷作用下，带椭圆夹杂的无限大各向异性压电体问题，得到了基体和夹杂体中温度场和应力函数的解析解。

上述这些工作为进一步分析含刚性夹杂的无限平面问题提供了可行的方法。但在实际工程应用中，构件的尺寸是有限的，无限平面问题的解析解不能满足有限平面问题的边界条件，且有限尺寸刚性夹杂问题尚未获得解析解。有限差分法、有限元法和边界元法是解决工程实际问题的常用数值方法，有限元法和边界元法更适合求解带复杂边界的边值问题，如果能获得线性微分方程的基本解，则边界元法是首选[11]，因为边界元法引入边界积分方程，使求解问题的维度降低一维，其将解析的基本解和数值离散相结合，数值结果具有较高的精度，而只需将区域的表面离散化，大大减少了未知数的个数，简化了前处理工作，且数值解的误差仅来源于在边界进行的数据插值和数值积分。文献[12]和[13]利用实边界元法分别研究了一维六方热弹准晶复合材料的界面断裂力学问题和一维六方热准晶复合材料的界面断裂力学问题，讨论了各种材料常数对裂纹尖端应力强度因子的影响。文献[14]简述了求解含椭圆刚性夹杂各向异性弹性板的断裂力学问题的复边界元法。文献[15]基于扩展的Stroh公式，利用复边界元方法研究了含椭圆孔有限大二十面体准晶板平面弹性问题，分析了在板端垂向拉伸作用下，板的大小、孔口尺寸及倾斜角度对孔边应力的影响。

准晶是一种复杂的金属合金，因其良好的性能而广泛应用于发动机等设备的表面涂层中，发动机的热障涂层在节能减排方面发挥着重要的作用。八次对称准晶是二维准晶的一种，根据某些具体构型，即使其弹性问题可以化成一个平面弹性问题与一个反平面声子场弹性问题的叠加，平面问题的最终控制方程仍很复杂，求解难度很大。文献[16]引进应力势函数，利用发展了的Lekhnitskii方法解决了八次对称二维准晶的椭圆缺口问题和位错问题。文献[17]应用Fourier变换与对偶积分方程理论得到了无限大八次对称二维准晶材料Ⅱ型Griffith裂纹的解析解。文献[18]用有限元方法求解了有限大平面八次对称准晶受内压圆筒问题以及I型Griffith裂纹问题。文献[19]用复边界元法对含椭圆孔的有限大八次对称二维准晶平面弹性问题进行了求解，并将椭圆孔退化到裂纹，计算了裂纹尖端的应力强度因子。

本文在不考虑夹杂相对基体旋转的情况下，利用复边界元法，研究有限大八次对称二维准晶椭圆刚性夹杂问题，即，考虑z轴方向极薄的等厚度板，刚性夹杂的构型如图1所示，板端施加平行于板面且不沿厚度变化的拉伸力。

2. 基本方程和Stroh公式

2.1. 基本方程

不妨假定八次对称二维准晶的周期方向是z轴，准周期平面是xoy平面，夹杂不随周期方向变化，即$\frac{\partial }{\partial z}=0$ 。由文献[16]，点群8 mm八次对称二维准晶的广义Hooke定律为

$\{\begin{array}{l}{\sigma }_{xx}=L\left({\epsilon }_{xx}+{\epsilon }_{yy}\right)+2M{\epsilon }_{xx}+R\left({\omega }_{xx}+{\omega }_{yy}\right),\\ {\sigma }_{yy}=L\left({\epsilon }_{xx}+{\epsilon }_{yy}\right)+2M{\epsilon }_{yy}-R\left({\omega }_{xx}+{\omega }_{yy}\right),\\ {\sigma }_{xy}={\sigma }_{yx}=2M{\epsilon }_{xy}+R\left({\omega }_{yx}-{\omega }_{xy}\right),\\ {\tau }_{xx}=K_1{\omega }_{xx}+K_2{\omega }_{yy}+R\left({\epsilon }_{xx}-{\epsilon }_{yy}\right),\\ {\tau }_{yy}=K_1{\omega }_{yy}+K_2{\omega }_{xx}+R\left({\epsilon }_{xx}-{\epsilon }_{yy}\right),\\ {\tau }_{xy}=\left(K_1+K_2+K_3\right){\omega }_{xy}+K_3{\omega }_{yx}-2R{\epsilon }_{xy},\\ {\tau }_{yx}=\left(K_1+K_2+K_3\right){\omega }_{yx}+K_3{\omega }_{xy}+2R{\epsilon }_{xy},\end{array}$ (1)

变形几何方程为

${\epsilon }_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right),{\omega }_{ij}=v_{i,j},$ (2)

不计体力的平衡方程为

${\sigma }_{ij,j}=0,{\tau }_{ij,j}=0,$ (3)

其中$i=1,2$ ，$j=1,2$ ；${\sigma }_{ij},{\epsilon }_{ij},u_i$ 表示声子场的应力、应变和位移；${\tau }_{ij},{\omega }_{ij},v_i$ 是相位子场的应力、应变和位移；L和M表示声子场的弹性常数；$K_1$ 、$K_2$ 及$K_3$

表示相位子场的弹性常数；R表示声子场和相位子场耦合的弹性常数。下标中的“，”表示偏微分，即

$u_{i,j}=\frac{\partial u_i}{\partial x_j},v_{i,j}=\frac{\partial v_i}{\partial x_j};$

指标符号$x_1$ 和$x_2$ 分别表示空间直角坐标x和y。用${\sigma }_{12}$ 简记${\sigma }_{xy}$ ，${\sigma }_{xy}$ 表示作用在y面(其外法线沿着y方向)上沿x方向的应力，其他符号的定义是类似的。同一个下标符号在一个单项中出现两次表示对该下标符号的所有可能取值求和，如

${\sigma }_{ij,j}={\sigma }_{i1,1}+{\sigma }_{i2,2},{\tau }_{ij,j}={\tau }_{i1,1}+{\tau }_{i2,2},$

上式中i是自由指标，j是哑标。

2.2. Stroh公式[16]

假设薄板各点沿z轴方向所有的场变量都为零，我们引入广义位移向量

$U=\xi g\left(z\right),$ (4)

其中，$U={\left(u_1,u_2,v_1,v_2\right)}^{\text{T}}$ ，$g\left(z\right)$ 为任意的解析函数，$z=x_1+\lambda x_2$ ，上标“T”表示转置，$\lambda$ 是复数，$\xi$ 是列向量。$U$ 可以由边界条件确定。把式(1)代入式(3)并结合式(4)和式(2)得

$\Lambda \xi =0,$ (5)

其中

$\begin{array}{l}\Lambda =D+\lambda \left(N+N^{\text{T}}\right)+{\lambda }^{2}S,D=\left[\begin{array}{cc}{Q}_1& R_1\\ R_1& T_1\end{array}\right],N=\left[\begin{array}{cc}{Q}_2& R_2\\ R_2^T& T_2\end{array}\right],S=\left[\begin{array}{cc}{Q}_3& -R_1\\ -R_1& T_3\end{array}\right],\\ Q_1=\left[\begin{array}{cc}L+2M& 0\\ 0& M\end{array}\right],Q_2=\left[\begin{array}{cc}0& L\\ M& 0\end{array}\right],Q_3=\left[\begin{array}{cc}M& 0\\ 0& L+2M\end{array}\right],R_1=\left[\begin{array}{cc}R& 0\\ 0& R\end{array}\right],\\ R_2=\left[\begin{array}{cc}0& R\\ -R& 0\end{array}\right],T_1=\left[\begin{array}{cc}{K}_1& 0\\ 0& K_1+K_2+K_3\end{array}\right],T_2=\left[\begin{array}{cc}0& K_2\\ K_3& 0\end{array}\right],T_3=\left[\begin{array}{cc}{K}_1+K_2+K_3& 0\\ 0& K_1\end{array}\right].\end{array}$ (6)

若式(5)存在非平凡解$\xi$ ，须有

$\text{det}\left(\Lambda \right)=0.$ (7)

方程(7)有4对共轭特征根${\lambda }_k$ 和${\overline{\lambda }}_k$ ，由应变能的正定性可知，特征根$\lambda$ 是非实的，不失一般性，我们取$\lambda =\left\{{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4\right\}$ 是虚部为正的特征根，相应的特征向量$\xi =\left\{{\xi }_k\right\}$ 由方程(5)给出。方程(3)的基本解可以表示为

$U=2\mathrm{Re}\left\{Ag\left(z\right)\right\},$ (8)

其中，$A={\left({\xi }_1,{\xi }_2,{\xi }_3,{\xi }_4\right)}^{\text{T}}$ ，$g\left(z\right)={\left(g_1\left(z_1\right),g\left(z_2\right),g_3\left(z_3\right),g_4\left(z_4\right)\right)}^{\text{T}}$ ，$z_k=x_1+{\lambda }_k{x}_2$ 。

取$\varsigma =\left(N^{\text{T}}+\lambda S\right)\xi$ ，声子场和相位子场应力可以表示为

${\sigma }_{j1}=-\lambda {\varsigma }_j{g}^{\prime }\left(z\right),{\sigma }_{j2}={\varsigma }_j{g}^{\prime }\left(z\right),{\tau }_{j1}=-\lambda {\varsigma }_{j+2}{g}^{\prime }\left(z\right),{\tau }_{j2}={\varsigma }_{j+2}{g}^{\prime }\left(z\right),j=1,2.$ (9)

引入应力函数P，声子场和相位子场应力可以表示为

${\sigma }_{j1}=-P_{j,2},{\sigma }_{j2}=P_{j,1},{\tau }_{j1}=-P_{j+2,2},{\tau }_{j2}=P_{j+2,1},j=1,2.$ (10)

则可求得应力函数

$P=2\text{Re}\left\{Bg\left(z\right)\right\},$ (11)

其中，$B=\left[{\varsigma }_1,{\varsigma }_2,{\varsigma }_3,{\varsigma }_4\right]$ 。

若无特殊说明，本文取$k=1,2,3,4$ 。

2.3. 基本解

下面考虑含刚性夹杂的无限大准晶平面弹性问题。假设薄板在$z^0=\left(x_1^0,x_2^0\right)$ 处受一集中点力$q=\left(q_1,q_2,q_3,q_4\right)$ 作用，在无限远处不受力，如图1。

Figure 1. Configuration and concentration point

图1. 构型及集中点

如果刚性夹杂相对基体不发生旋转，则相应的边界条件表示为[14]

$U=0$ ，沿着椭圆边界；

${\displaystyle {\oint }_{C_1}\text{d}P}=q,{\displaystyle {\oint }_{C_2}\text{d}U}=\text{0}$ ，$C_1$ 表示包含点的封闭曲线，$C_2$ 表示任意封闭曲线；(12)

${\sigma }_{ij}\to 0,{\tau }_{ij}\to 0$ ，在无穷远处。

式(12)表明，如果点力q非零，应力函数$g_k$ 应该是多值函数，要得到该多值函数，需要选择一个与点$z^0$ 有关的多值函数$g_k\left(z_k,z_k^0\right)$ 。参考文献[15]，我们选择

$g_k\left(z_k,z_k^0\right)=\ln \left(z_k-z_k^0\right){\eta }_k,$ (13)

其中${\eta }_k$ 是未知复常数，由边界条件确定。结合式(13)，由式(8)和式(11)得到的基本解可以写为

$U=2\text{Re}\left\{A\langle \ln \left(z_k-z_k^0\right)\rangle \eta \right\},P=2\text{Re}\left\{B\langle \ln \left(z_k-z_k^0\right)\rangle \eta \right\},$ (14)

其中，$\eta ={\left({\eta }_1,{\eta }_2,{\eta }_3,{\eta }_4\right)}^{\text{T}}$ ，$\langle \text{ln}\left(z_k-z_k^0\right)\rangle$ 表示4阶对角矩阵。对式(13)进行积分，则有

${\displaystyle {\oint }_C\text{d}}g\left(z\right)=2\text{πi}.$

把式(14)代入边界条件式(12)的第二组条件，可得

$2\text{πi}\left\{A\eta -\overline{A\eta }\right\}=0\text{, }2\text{πi}\left\{B\eta -\overline{B\eta }\right\}=q.$(15)

利用正交性

$\left(\begin{array}{cc}{B}^{\text{T}}& A^{\text{T}}\\ {\overline{B}}^{\text{T}}& {\overline{A}}^{\text{T}}\end{array}\right)\left(\begin{array}{cc}A& \overline{A}\\ B& \overline{B}\end{array}\right)=\left(\begin{array}{cc}I& 0\\ 0& I\end{array}\right),$

由式(15)可得$\eta =A^{\text{T}}q/2\text{πi}$ 。

2.4. Green函数

边界元法的基础是Green函数。为了方便问题求解，利用保角变换

${\chi }_k=\frac{z_k+\sqrt{z_k^2-a^2-{\lambda }_k^2{b}^2}}{a-\text{i}{\lambda }_{k}b}$ (16)

将不同倾斜角椭圆孔的外部映射到单位圆孔的外部。根据解析延拓方法[14]，设椭圆刚性夹杂问题的Green函数为

$U=2\text{Re}\left\{A\left[g^{\text{I}}\left(\chi \right)+g^{\text{R}}\left(\chi \right)\right]\right\},P=2\text{Re}\left\{B\left[g^{\text{I}}\left(\chi \right)+g^{\text{R}}\left(\chi \right)\right]\right\},$ (17)

其中$g^{\text{I}}\left(\chi \right)=\langle \text{ln}\left({\chi }_k-{\chi }_k^0\right)\rangle {\eta }_{0}q$ ，${\eta }_0=A^{\text{T}}/2\text{πi}$ ，$g^{\text{R}}\left(\chi \right)$ 是对应于刚性夹杂问题摄动场的解析函数，由边界条件确定。$U$ 是无限大八次准晶平面弹性问题的基本解，如式(14)。把式(17)代入式(12)的第一个条件，经过计算得

$g^{\text{R}}\left(\chi \right)=A^{-1}\overline{A}\langle \text{ln}\left({\chi }^{-1}-{\overline{{\chi }^0}}_k\right)\rangle {\overline{A}}^{\text{T}}q/2\text{πi}={\displaystyle \sum }_{j=1}^4\langle \text{ln}\left({\chi }_k^{-1}-{\overline{{\chi }^0}}_k\right)\rangle A^{-1}\overline{A}{E}_j{\overline{A}}^{\text{T}}q/2\text{πi}\text{,}$

其中$E_j={\left(e_{mn}\right)}_{4\times 4}$ ，$e_{mn}=\{\begin{array}{l}1,m=n=j\\ 0,其它\end{array}$ 。

3. 边界元法

3.1. 边界条件

含刚性夹杂有限尺寸准晶薄板如图2所示，椭圆刚性夹杂的中心位于坐标原点，椭圆孔边界用参数方程

$x_1=a\cos \theta ,x_2=b\sin \theta$ (18)

表示，参数$\theta \in \left[0,2\text{π}\right]$ ，a和b分别为椭圆的长半轴和短半轴。W为薄板宽，H为薄板高，板的边界用$\Gamma$ 表示，椭圆中心与薄板中心重合。基体的椭圆孔边没有加载，仅在试样两端施加拉伸力。含椭圆孔有限尺寸准晶薄板如图3所示[19]。

Figure 2. Rigid elliptical inclusion quasicrystalline thin plate

图2. 含椭圆刚性夹杂准晶薄板

Figure 3. Quasicrystalline thin plate with elliptical orifice

图3. 含椭圆孔口准晶薄板

在体力忽略不计且基体孔边不受力的情况下，边界条件为

$f_k={\sigma }_{kj}{n}_j={\tilde{f}}_k,\left(x_1,x_2\right)\in {\Gamma }_1;u_k={\tilde{u}}_k,\left(x_1,x_2\right)\in {\Gamma }_2,$ (19)

其中，$k=1,2,3,4$ ；$j=1,2$ ；${\text{Γ}}_1\cup {\text{Γ}}_2=\text{Γ}$ 。

3.2. 边界积分方程及离散化

引入记号

${\sigma }_{3j}\triangleq {\tau }_{1j},{\sigma }_{4j}\triangleq {\tau }_{2j},{\epsilon }_{3j}=\frac{1}{2}{u}_{3,j}\triangleq \frac{1}{2}{\omega }_{1j},{\epsilon }_{4j}=\frac{1}{2}{u}_{4,j}\triangleq \frac{1}{2}{\omega }_{2j},$

则式(2)和式(3)简记为

${\epsilon }_{kj}=\left(u_{k,j}+u_{j,k}\right)/2$

和

${\sigma }_{kj,j}=0.$ (20)

当k > 2时，$u_{j,k}=0$ 。式(19)和式(20)组成偏微分方程边值问题，利用加权余量法，其定解问题的弱形式表示为

${\displaystyle {\iint }_{\Omega }{\sigma }_{kj,j}{\widehat{u}}_k\text{d}\Omega }-{{\displaystyle \int }}_{{\Gamma }_1}\left(f_k-{\tilde{f}}_k\right){\widehat{u}}_k\text{d}\Gamma +{{\displaystyle \int }}_{{\Gamma }_2}\left(u_k-{\tilde{u}}_k\right){\widehat{p}}_k\text{d}\Gamma =0.$(21)

其中${\widehat{u}}_k$ 和${\widehat{p}}_k$ 为权函数。

对式(21)的第一项应用高斯公式，结合式(19)并利用功的互等定理，则有

$\begin{array}{c}{\displaystyle {\iint }_{\Omega }{\sigma }_{kj,j}{\widehat{u}}_k\text{d}\Omega }={{\displaystyle \int }}_{\Omega }{\left({\sigma }_{kj}{\widehat{u}}_k\right)}_{,j}\text{d}\Omega -{{\displaystyle \int }}_{\Omega }{\sigma }_{kj}{\widehat{u}}_{k,j}\text{d}\Omega \\ ={{\displaystyle \int }}_{\Gamma }{\sigma }_{kj}{\widehat{u}}_k{n}_j\text{d}\Gamma -{{\displaystyle \int }}_{\Omega }{\sigma }_{kj}{\widehat{u}}_{k,j}\text{d}\Omega \\ ={{\displaystyle \int }}_{\Gamma }{f}_k{\widehat{u}}_k\text{d}\Gamma -{{\displaystyle \int }}_{\Omega }{\sigma }_{kj}{\widehat{u}}_{k,j}\text{d}\Omega \\ ={{\displaystyle \int }}_{\Gamma }{f}_k{\widehat{u}}_k\text{d}\Gamma -{{\displaystyle \int }}_{\Omega }{\widehat{\sigma }}_{kj}{u}_{k,j}\text{d}\Omega \\ ={{\displaystyle \int }}_{\Gamma }{f}_k{\widehat{u}}_k\text{d}\Gamma -{{\displaystyle \int }}_{\Omega }{\left({\widehat{\sigma }}_{kj}{u}_k\right)}_{,j}\text{d}\Omega +{{\displaystyle \int }}_{\Omega }{\widehat{\sigma }}_{kj,j}{u}_k\text{d}\Omega \\ ={{\displaystyle \int }}_{{\Gamma }_1}{f}_k{\widehat{u}}_k\text{d}\Gamma +{{\displaystyle \int }}_{{\Gamma }_2}{f}_k{\widehat{u}}_k\text{d}\Gamma -{{\displaystyle \int }}_{{\Gamma }_1}{\widehat{p}}_k{u}_k\text{d}\Gamma -{{\displaystyle \int }}_{{\Gamma }_2}{\widehat{p}}_k{u}_k\text{d}\Gamma +{{\displaystyle \int }}_{\Omega }{\widehat{\sigma }}_{kj,j}{u}_k\text{d}\Omega ,\end{array}$(22)

将式(22)代入式(21)得

${{\displaystyle \int }}_{\Omega }{\widehat{\sigma }}_{kj,j}{u}_k\text{d}\Omega +{{\displaystyle \int }}_{{\Gamma }_1}{\tilde{f}}_k{\widehat{u}}_k\text{d}\Gamma +{{\displaystyle \int }}_{{\Gamma }_2}{f}_k{\widehat{u}}_k\text{d}\Gamma -{{\displaystyle \int }}_{{\Gamma }_1}{u}_k{\widehat{p}}_k\text{d}\Gamma -{{\displaystyle \int }}_{{\Gamma }_2}{\tilde{u}}_k{\widehat{p}}_k\text{d}\Gamma =0.$(23)

若权值${\widehat{u}}_k$ 和${\widehat{f}}_k$ 为方程(20)的位移和应力基本解，则与${\widehat{u}}_k$ 对应的${\widehat{\sigma }}_{kj}$ 应满足

${\widehat{\sigma }}_{kj,j}+{\delta }_{kL}\Delta \left(z^0,z\right)=0,$ (24)

其中“δ”是Kronecker符号，“Δ”是Dirac符号，$L=1,2,3,4$ 。边界上统一记号，${\tilde{f}}_k$ 和$f_k$ 均由${\stackrel{⌣}{f}}_k$ 表示，${\tilde{u}}_k$ 和$u_k$ 均由${\stackrel{⌣}{u}}_k$ 表示，为了避免混淆，哑标k改写为$\alpha$ ，$\alpha =1,2,3,4$ 。由式(23)和式(24)可得区域Ω内的积分方程

${\stackrel{⌣}{u}}_L\left(z^0\right)={{\displaystyle \int }}_{\Gamma }{\stackrel{⌣}{f}}_{\alpha }{\widehat{u}}_{\alpha L}d\Gamma -{{\displaystyle \int }}_{\Gamma }{\stackrel{⌣}{u}}_{\alpha }{\widehat{f}}_{\alpha L}d\Gamma ,$ (25)

其中${\widehat{u}}_{\alpha L}$ 和${\widehat{p}}_{\alpha L}$ 为已知函数。记长度为4的向量$q_L={\left(q_k\right)}_4$ ，不失一般性，取

$q_k=\{\begin{array}{l}1,k=L,\\ 0,k\ne L.\end{array}$

${\widehat{u}}_{\alpha L}$ 由式(17)中的$U$ 直接得到，${\widehat{f}}_{\alpha L}$ 由式(17)中的$\frac{\partial P}{\partial m}$ 得到。式(25)表明，若已知边界上的位移$u_L$ 和面力$f_L$ ，即可求得区域内任意一点的声子场位移及相位子场位移。将区域内的点$z^0$ 边界逼近直至位于边界上，利用区域延拓法[16]可得边界积分方程

$D_{\alpha L}\left(z^0\right)u_{\alpha }\left(z^0\right)=\underset{\Gamma }{\overset{ }{{\displaystyle \int }}}{p}_{\alpha }\left(z\right)u_{\alpha L}\left(z^0,z\right)\text{d}\Gamma \left(z\right)-\underset{\Gamma }{\overset{ }{{\displaystyle \int }}}{u}_{\alpha }\left(z\right)f_{\alpha L}\left(z^0,z\right)\text{d}\Gamma \left(z\right),$ (26)

其中

$D_{\alpha L}\left(z^0\right)=-\underset{\Gamma }{\overset{ }{{\displaystyle \int }}}{f}_{\alpha L}\left(z^0,z\right)\text{d}\Gamma \left(z\right).$ (27)

式(26)中的$u_{\alpha L}\left(z^0,z\right)$ 和$f_{\alpha L}\left(z^0,z\right)$ 分别表示点$z^0$ 上$x_{\alpha }$ 方向的单位点力诱导出点z处$x_L$ 方向的位移和牵引力，$x_1$ 和$x_2$ 分别对应声子场的x方向和y方向，$x_3$ 和$x_4$ 分别对应相位子场的x方向和y方向，$D_{\alpha L}\left(z^0\right)$ 是依赖于$z^0$ 的系数。利用式(17)把$z^0$ 和z分别变换为${\chi }^0$ 和$\chi$ ，则边界积分方程式(26)的矩阵形式可写为

$D\left({\chi }^0\right)u\left({\chi }^0\right)=\underset{\Gamma }{\overset{ }{{\displaystyle \int }}}f\left(\chi \right)U\left({\chi }^0,\chi \right)\text{d}\Gamma \left(\chi \right)-\underset{\Gamma }{\overset{ }{{\displaystyle \int }}}u\left(\chi \right)F\left({\chi }^0,\chi \right)\text{d}\Gamma \left(\chi \right),$ (28)

待求解的物理量是$f$ 或$u$ 。由式(17)的格林函数可得$U\left({\chi }^0,\chi \right)$ 和$F\left({\chi }^0,\chi \right)$ 如下

$U\left({\chi }^0,\chi \right)=\frac{1}{\text{π}}\text{Im}\left\{A\left[\ln \left({\chi }_{\alpha }-{\chi }_{\alpha }^0\right)A^{\text{T}}+{{\displaystyle \sum }}_{i=1}^4\ln \left({\chi }_{\alpha }^{-1}-{\overline{\chi }}_i^0\right)A^{-1}\overline{A}{E}_i{\overline{A}}^{\text{T}}\right]\right\},$

$\begin{array}{c}F\left({\chi }^0,\chi \right)=\frac{1}{\text{π}}\mathrm{Im}\left\{B\frac{2{\chi }_{\alpha }^2\left(\partial x_1/\partial m+p_k\partial x_2/\partial m\right)}{\left({\chi }_{\alpha }-{\chi }_{\alpha }^0\right)\left[\left(a-\text{i}{\lambda }_{\alpha }b\right){\chi }_{\alpha }^2-\left(a+\text{i}{\lambda }_{\alpha }b\right)\right]}{A}^{\text{T}}\right\}\\ +\frac{1}{\text{π}}{{\displaystyle \sum }}_{i=1}^4\mathrm{Im}\left\{B\frac{-2{\chi }_{\alpha }\left(\partial x_1/\partial m+p_{\alpha }\partial x_2/\partial m\right)}{\left(1-{\chi }_{\alpha }{\overline{\chi }}_i^0\right)\left[\left(a-\text{i}{\lambda }_{\alpha }b\right){\chi }_k^2-\left(a+\text{i}{\lambda }_{\alpha }b\right)\right]}{A}^{-1}\overline{A}{E}_i{\overline{A}}^{\text{T}}\right\}.\end{array}$

式(27)中$f$ 和$u$ 均是长度为4的向量$D,U$ 和$F$ 均是四阶方阵。由式(27)可知，$D$ 为对角阵，平滑边界上$D$ 的主对角线元素都是1/2，内点处$D$ 是单位阵。

利用Guass数值积分公式对边界积分方程式(27)的积分进行离散，有

$D\left({\chi }^0\right)u\left({\chi }^0\right)+{{\displaystyle \sum }}_{n=1}^N{\tilde{Y}}_n\left({\chi }^0\right)u_n={{\displaystyle \sum }}_{n=1}^N{G}_n\left({\chi }^0\right)f_n$,

整理为

${{\displaystyle \sum }}_{n=1}^N{Y}_{in}{u}_n={{\displaystyle \sum }}_{n=1}^N{G}_{in}{f}_n,i=1,2,\cdots ,N,$ (29)

其中

$Y_{in}={\tilde{Y}}_{in},i\ne n;Y_{in}=Y_{in}+D_i,i=n.$

由式(29)可以计算出各边界节点处的位移或牵引力、应力及应变值。如果${\chi }^0$ 是内部点利用同样的数值积分公式，对区域内积分方程式(25)进行离散并整理，有

$u_i\left({\chi }^0\right)={{\displaystyle \sum }}_{n=1}^N{G}_{in}\left({\chi }^0\right)f_n-{{\displaystyle \sum }}_{n=1}^N{\widehat{Y}}_{in}\left({\chi }^0\right)u_n,i=1,2,\cdots ,N.$(30)

利用边界的物理量$f$ 和$u$ ，由式(30)可以计算出区域内任意点相应的物理量。

4. 数值实例

取八次对称二维准晶材料的弹性常数L = 234.3 GPa，M = 57.41GPa，K₁ = 122 GPa，K₂ = 24 GPa，K₃ = 12 GPa。有刚性夹杂和没有夹杂(如图3)的椭圆孔口大小相同，基体孔边声子场和相位子应力随$\theta$ 的变化(逆时针)如图4所示。

观察图4可以看到，由于内部刚性夹杂的存在，基体孔边应力值减小。当板端受到上下拉伸力时，声子场孔边应力${\sigma }_{22}$ 和${\sigma }_{11}$ 及相位子场孔边应力$H_{22}$ 和$H_{11}$ 关于$\theta =180˚$ 轴对称；A点和D点的应力值最大；B点和C点的应力值最小。

(1) 声子场孔边应力${\sigma }_{12}$ (2) 声子场孔边应力${\sigma }_{21}$ (3) 声子场孔边应力${\sigma }_{11}$

(4) 声子场孔边应力${\sigma }_{22}$ (5) 相位子场孔边应力$H_{22}$

(6) 相位子场孔边应力$H_{12}$ (7) 相位子场孔边应力$H_{21}$ (8) 相位子场孔边应力$H_{11}$

Figure 4. A comparison of the variation of phonon field and phase site stress with $\theta$ in elliptic rigid inclusion and aperture problem

图4. 椭圆刚性夹杂与孔口问题声子场和相位子应力随$\theta$ 的变化对比图

5. 结论

当材料中某点的应变超过最大弹性应变时，可以在损伤区域插入椭圆夹杂，由于内部夹杂的存在，基体孔边应力值减弱，损伤区应变值减小。椭圆夹杂属于常用的细观力学模型范畴，本文基于工程实际，利用复边界元法研究有限大二维准晶试样的夹杂问题，这些研究对提高材料的性能和可靠性，促进新型材料的发展有一定的指导意义。

基金项目

内蒙古自治区高校科研项目(No. NJZZ21003)。