摘要

以两个直角楔体组成的双材料半平面为模型，研究黏结界面端点的应力奇异性。采用艾里应力函数和分离变量法，由自由边与界面条件建立特征方程。两侧泊松比均取 0.3 时，数值结果显示，模量比越偏离 1，奇异阶次越大。两侧材料相同时，搜索区间内无奇异实根。互换两种材料所得阶次相同。模量比为 100 时，端点应力的奇异阶次约为 0.278。当模量比为 1.5、材料 1 的泊松比为 0.30 时，材料 2 的泊松比取 0.43 时，幂律奇异根消失。COMSOL 拉伸算例得到的界面应力趋势与理论一致。

一、引言

在涂层与基体、电子封装与基板等结构中，连接两种材料的界面往往延伸到外表面。沿界面向外走，到了结构边缘，共同传递面力的黏结界面与无面力的外表面相交。即使外表面是一条直线，端点附近也可能出现很陡的应力变化。用线弹性模型描述时，这种变化可以表现为应力随距离缩短而按幂律发散。

Bogy $^{[1]}$ 分析了两种材料组成的黏结直角楔体，指出奇异阶次与两侧的剪切模量和泊松比有关，随后将问题推广到任意楔角$^{[2]}$。本文取两个楔角均为直角，把端点附近看作一个复合半平面。通过固定泊松比、改变弹性模量比，考察外表面平直时产生奇异应力的条件，以及两侧刚度差异怎样改变应力随距离缩短而增长的速度。本文采用局部解，只需给出端点附近的材料和边界条件，即可比较各种模量比下的发散速度。在此基础上，本文固定模量比，计算泊松比变化时的临界条件，并以有限元拉伸算例比较临界线两侧的界面应力。

二、力学模型与边界条件

如图1 所示，取黏结界面端点 $O$ 为原点，以沿界面向内的方向为 $\theta = 0$。距离端点的长度为 $r$，逆时针量取角度 $\theta$。记直角为 $\theta_{0}$，上下两种材料的角度范围分别为

$$ \theta_{0} = \frac{\pi}{2},\quad\quad 0 < \theta < \theta_{0}\quad(i = 1),\quad\quad - \theta_{0} < \theta < 0\quad(i = 2) \tag{1} $$

材料 $i$ 的弹性模量和泊松比记为 $E_{i},\nu_{i}$。两种材料均匀、各向同性，采用小变形线弹性和平面应变假设，忽略体力。

在 $\theta = \theta_{0}$ 和 $\theta = - \theta_{0}$ 的自由边上，外法向沿环向，故法向应力 $\sigma_{\theta}$ 与切应力 $\tau_{r\theta}$ 均为零。在界面 $\theta = 0$ 上，沿同一环向比较，两侧这两个应力分量连续。完全黏结也要求两侧径向位移 $u_{r}$ 和环向位移 $u_{\theta}$ 连续。两条自由边各有两个条件，界面有四个条件。

图1　双材料黏结自由边模型及极坐标系

三、端点应力场与特征方程

1、艾里应力函数与角向方程

引入材料 $i$ 的艾里应力函数 $U_{i}$。在两种材料内部，$U_{i}$ 满足双调和方程 $\nabla^{4}U_{i} = 0$ $^{[3]}$。

边界和界面都是从端点出发的射线。沿同一射线，角度 $\theta$ 固定，距离 $r$ 可以变化，因此用径向幂函数乘一个角函数来表示一个模态。用 $\lambda$ 表示待求的特征指数，令

$$ U_{i}(r,\theta) = r^{\lambda + 1}f_{i}(\theta),\quad\quad r > 0 \tag{2} $$

这里 $f_{i}(\theta)$ 描述模态在材料 $i$ 内随角度的变化，两侧用同一个 $\lambda$，使界面两侧能够逐项比较相同的径向幂次。式（2） 中取指数 $\lambda + 1$，是因为求应力时会降两次幂，最后得到 $r^{\lambda - 1}$。

为计算角向方程，记 $s$ 为径向幂指数，$f$ 为任意角函数。撇号表示对 $\theta$ 求导。极坐标下对 $r^{s}f$ 作用一次拉普拉斯算子，得到

$$ \nabla^{2}(r^{s}f) = s(s - 1)r^{s - 2}f + sr^{s - 2}f + r^{s - 2}f^{''} = r^{s - 2}(s^{2}f + f^{''}) \tag{3} $$

第一次拉普拉斯对式（2） 的应力函数应用拉普拉斯算子，在式（3） 中取 $s = \lambda + 1$，得到

$$ \nabla^{2}U_{i} = r^{\lambda - 1}\left\lbrack (\lambda + 1)^{2}f_{i} + f_{i}^{''} \right\rbrack \tag{4} $$

式（4） 中的括号只依赖于 $\theta$，记为 $g_{i}(\theta)$。第二次拉普拉斯将式（3） 用于 $r^{\lambda - 1}g_{i}$，取 $s = \lambda - 1$，得到

$$ \nabla^{4}U_{i} = r^{\lambda - 3}\left\lbrack (\lambda - 1)^{2}g_{i} + g_{i}^{''} \right\rbrack = 0 \tag{5} $$

双调和方程要求式（5） 为零。因为 $r > 0$，约去 $r^{\lambda - 3}$，再把 $g_{i}$ 展开为 $f_{i}^{''} + (\lambda + 1)^{2}f_{i}$，得到方程

$$ \left\lbrack \frac{\mathrm{d}^{2}}{\mathrm{d}\theta^{2}} + (\lambda - 1)^{2} \right\rbrack\left\lbrack \frac{\mathrm{d}^{2}}{\mathrm{d}\theta^{2}} + (\lambda + 1)^{2} \right\rbrack f_{i} = 0 \tag{6} $$

记 $k_{+} = \lambda + 1$、$k_{-} = \lambda - 1$。两个二阶因子各对应一组正弦和余弦解。对 $0 < \lambda < 1$，两组频率的绝对值不同且均非零，四个基函数相互独立。因此

$$ f_{i +} = a_{i}\cos(k_{+}\theta) + b_{i}\sin(k_{+}\theta),\quad\quad f_{i -} = c_{i}\cos(k_{-}\theta) + d_{i}\sin(k_{-}\theta) \tag{7} $$

$$ f_{i} = f_{i +} + f_{i -} \tag{8} $$

式中 $a_{i},b_{i},c_{i},d_{i}$ 是材料 $i$ 的四个待定系数，正负下标用来标记两组频率。两种材料共有八个系数，正好对应第二节的八条自由边和界面条件。

2、由应力函数求应力

在极坐标中，$\sigma_{r}^{(i)}$ 和 $\sigma_{\theta}^{(i)}$ 分别为径向和环向正应力，$\tau_{r\theta}^{(i)}$ 为面内切应力。极坐标基向量随 $\theta$ 转动，因此对 $U_{i}$ 求导时会出现 $r^{- 1}$ 和 $r^{- 2}$。应力函数与三个分量的关系为 $^{[3]}$

$$ \begin{pmatrix} \sigma_{r}^{(i)} \\ \sigma_{\theta}^{(i)} \\ \tau_{r\theta}^{(i)} \end{pmatrix} = \begin{pmatrix} \frac{1}{r}\frac{\partial U_{i}}{\partial r} + \frac{1}{r^{2}}\frac{\partial^{2}U_{i}}{\partial\theta^{2}} \\ \frac{\partial^{2}U_{i}}{\partial r^{2}} \\ - \frac{\partial}{\partial r}\left( \frac{1}{r}\frac{\partial U_{i}}{\partial\theta} \right) \end{pmatrix} \tag{9} $$

把式（2） 代入式（9）。径向正应力中的两项分别成为 $(\lambda + 1)r^{\lambda - 1}f_{i}$ 和 $r^{\lambda - 1}f_{i}^{''}$。环向正应力是径向二阶导数，系数为 $\lambda(\lambda + 1)$。求切应力时，先在括号内得到 $r^{\lambda}f_{i}'$，再对 $r$ 求导并取负号。合写为：

$$ \begin{pmatrix} \sigma_{r}^{(i)} \\ \sigma_{\theta}^{(i)} \\ \tau_{r\theta}^{(i)} \end{pmatrix} = r^{\lambda - 1}\begin{pmatrix} (\lambda + 1)f_{i} + f_{i}^{''} \\ \lambda(\lambda + 1)f_{i} \\ - \lambda f_{i}' \end{pmatrix} \tag{10} $$

三个分量具有相同的径向因子 $r^{\lambda - 1}$。自由边和黏结界面的法向均沿环向，故边界上的面力条件只涉及式（10） 的后两行。

3、由应力求位移并检查相容性

对于平面应变问题，材料 $i$ 的面内本构关系可直接用等效弹性模量 ${\overline{E}}_{i}$ 和等效泊松比 ${\overline{\nu}}_{i}$ 表示 $^{[3]}$（这里避免下标混淆，不采用课本上的表达方法）。它们由材料的实际参数 $E_{i},\nu_{i}$ 换算得到。

$$ {\overline{\nu}}_{i} = \frac{\nu_{i}}{1 - \nu_{i}},\quad\quad{\overline{E}}_{i} = \frac{E_{i}}{1 - \nu_{i}^{2}} \tag{11} $$

代入式（11） 的等效参数，径向正应变、环向正应变和剪应变直接写为

$$ \begin{pmatrix} \varepsilon_{r}^{(i)} \\ \varepsilon_{\theta}^{(i)} \\ \gamma_{r\theta}^{(i)} \end{pmatrix} = \frac{1}{{\overline{E}}_{i}}\begin{pmatrix} \sigma_{r}^{(i)} - {\overline{\nu}}_{i}\sigma_{\theta}^{(i)} \\ \sigma_{\theta}^{(i)} - {\overline{\nu}}_{i}\sigma_{r}^{(i)} \\ 2(1 + {\overline{\nu}}_{i})\tau_{r\theta}^{(i)} \end{pmatrix} \tag{12} $$

由式（10） 和式（12），应变也含有公因式 $r^{\lambda - 1}$。由于应变是位移的一阶导数，故可以令径向位移 $u_{r}^{(i)}$、环向位移 $u_{\theta}^{(i)}$ 为

$$ u_{r}^{(i)} = r^{\lambda}P_{i}(\theta),\quad\quad u_{\theta}^{(i)} = r^{\lambda}Q_{i}(\theta) \tag{13} $$

其中 $P_{i},Q_{i}$ 是位移的角函数。将式（13） 代入极坐标几何方程 $^{[3]}$。由左至右依次是位移导数表示的应变和提取径向公因式后的结果。

$$ \begin{pmatrix} \varepsilon_{r}^{(i)} \\ \varepsilon_{\theta}^{(i)} \\ \gamma_{r\theta}^{(i)} \end{pmatrix} = \begin{pmatrix} \frac{\partial u_{r}^{(i)}}{\partial r} \\ \frac{1}{r}\left( u_{r}^{(i)} + \frac{\partial u_{\theta}^{(i)}}{\partial\theta} \right) \\ \frac{1}{r}\frac{\partial u_{r}^{(i)}}{\partial\theta} + \frac{\partial u_{\theta}^{(i)}}{\partial r} - \frac{u_{\theta}^{(i)}}{r} \end{pmatrix} = r^{\lambda - 1}\begin{pmatrix} \lambda P_{i} \\ P_{i} + Q_{i}' \\ P_{i}' + (\lambda - 1)Q_{i} \end{pmatrix} \tag{14} $$

比较式（12） 和式（14） 的第一行，将式（10） 的应力分量代入并约去 $r^{\lambda - 1}$，得到

$$ \lambda P_{i} = \frac{(\lambda + 1)f_{i} + f_{i}^{''} - {\overline{\nu}}_{i}\lambda(\lambda + 1)f_{i}}{{\overline{E}}_{i}} \tag{15} $$

由式（15） 可解出 $P_{i}$。再比较式（12） 和式（14） 的第三行，将 $P_{i}$ 的环向导数和式（10） 的切应力代入，便可求出 $Q_{i}$。整理得

$$ P_{i} = \frac{f_{i}^{''} + (\lambda + 1)(1 - {\overline{\nu}}_{i}\lambda)f_{i}}{\lambda{\overline{E}}_{i}} \tag{16} $$

$$ Q_{i} = \frac{- 2\lambda(1 + {\overline{\nu}}_{i})f_{i}' - {\overline{E}}_{i}P_{i}'}{(\lambda - 1){\overline{E}}_{i}} \tag{17} $$

为简化位移表达式，按式（7） 的两组频率拆开角函数。利用 $f_{i +}^{''} = - k_{+}^{2}f_{i +}$、$f_{i -}^{''} = - k_{-}^{2}f_{i -}$，式（16） 和式（17） 可写成下列形式，其中代换变量 $\kappa_{i}$ 只由等效泊松比决定。

$$ P_{i} = \frac{1 + {\overline{\nu}}_{i}}{{\overline{E}}_{i}}\left\lbrack - k_{+}f_{i +} + (\kappa_{i} - \lambda)f_{i -} \right\rbrack,\quad\quad\kappa_{i} = \frac{3 - {\overline{\nu}}_{i}}{1 + {\overline{\nu}}_{i}} \tag{18} $$

$$ Q_{i} = - \frac{1 + {\overline{\nu}}_{i}}{{\overline{E}}_{i}}\left\lbrack f_{i +}' + \frac{\lambda + \kappa_{i}}{\lambda - 1}f_{i -}' \right\rbrack \tag{19} $$

对 $Q_{i}$ 求导，与 $P_{i}$ 相加，再用两组角函数的二阶导数关系，得

$$ P_{i} + Q_{i}' = \frac{\lambda(1 + {\overline{\nu}}_{i})}{{\overline{E}}_{i}}\left\lbrack (\lambda + 1)f_{i +} + (\lambda + \kappa_{i} - 2)f_{i -} \right\rbrack \tag{20} $$

用式（10） 和式（12） 的第二行直接求环向正应变，约去径向公因式后也得到式（20）。因此由径向应变和剪应变求出的位移同时满足环向应变关系。

4、组装特征矩阵并确定奇异阶次

将第二节的边界条件写成方程。

上自由边位于材料 1 的 $\theta = \theta_{0}$ 处，法向正应力和切应力都为零。

$$ \sigma_{\theta}^{(1)}(r,\theta_{0}) = 0,\quad\quad\tau_{r\theta}^{(1)}(r,\theta_{0}) = 0 \tag{21} $$

下自由边位于材料 2 的 $\theta = - \theta_{0}$ 处，同样要求两项面力为零。

$$ \sigma_{\theta}^{(2)}(r, - \theta_{0}) = 0,\quad\quad\tau_{r\theta}^{(2)}(r, - \theta_{0}) = 0 \tag{22} $$

黏结界面位于 $\theta = 0$。两侧沿同一极坐标方向比较，法向应力与切应力分别相等。

$$ \sigma_{\theta}^{(1)}(r,0) = \sigma_{\theta}^{(2)}(r,0),\quad\quad\tau_{r\theta}^{(1)}(r,0) = \tau_{r\theta}^{(2)}(r,0) \tag{23} $$

完全黏结还要求两侧的径向位移和环向位移分别相等。

$$ u_{r}^{(1)}(r,0) = u_{r}^{(2)}(r,0),\quad\quad u_{\theta}^{(1)}(r,0) = u_{\theta}^{(2)}(r,0) \tag{24} $$

式（21） 至式（24） 一共给出八条条件。把式（10） 的应力和式（13） 的位移代入，分别约去公因式 $r^{\lambda - 1}$ 和 $r^{\lambda}$。以式（23） 中的法向应力连续条件为例，在 $\theta = 0$ 处有 $f_{i}(0) = a_{i} + c_{i}$，因此该条件化为

$$ \lambda(\lambda + 1)(a_{1} + c_{1} - a_{2} - c_{2}) = 0 \tag{25} $$

把八个未知系数按 $(a_{1},b_{1},c_{1},d_{1},a_{2},b_{2},c_{2},d_{2})$ 的顺序排成一列。按这个顺序抄下式（25） 中各未知量前的系数，未出现的未知量填零，得到

$$ \lambda(\lambda + 1)(1,0,1,0, - 1,0, - 1,0) \tag{26} $$

其余七行同样按方程中的系数排列。令 $\boldsymbol{g}(\theta)$ 表示式（7） 中四个三角基函数按原顺序组成的行向量，$\boldsymbol{c}_{i} = (a_{i},b_{i},c_{i},d_{i})^{T}$ 表示材料 $i$ 的系数列向量。这样 $f_{i} = \boldsymbol{g}\boldsymbol{c}_{i}$、$f_{i}' = \boldsymbol{g}'\boldsymbol{c}_{i}$。从式（10） 的后两行约去径向公因子后，两项面力写成 $\boldsymbol{B}(\theta)\boldsymbol{c}_{i}$。矩阵 $\boldsymbol{B}$ 的第一行对应法向应力，第二行对应切应力。

$$ \boldsymbol{B}(\theta) = \begin{pmatrix} \lambda(\lambda + 1)\boldsymbol{g}(\theta) \\ - \lambda\boldsymbol{g}'(\theta) \end{pmatrix} \tag{27} $$

对位移，将 $\theta = 0$ 代入式（18） 和式（19）。此时余弦等于 1，正弦等于 0，正弦的环向导数分别为 $k_{+}$ 和 $k_{-}$。把 $P_{i}(0)$、$Q_{i}(0)$ 的系数各排成一行，就得到 $\boldsymbol{D}_{i}(0)$。

$$ \boldsymbol{D}_{i}(0) = \frac{1 + {\overline{\nu}}_{i}}{{\overline{E}}_{i}}\begin{pmatrix} - k_{+} & 0 & \kappa_{i} - \lambda & 0 \\ 0 & - k_{+} & 0 & - (\lambda + \kappa_{i}) \end{pmatrix} \tag{28} $$

按式（21） 至式（24） 的顺序组装矩阵。自由边条件只含该侧材料的系数，另一侧的四列填零。界面条件写成材料 1 的分量减去材料 2 的分量等于零，因此材料 2 对应的系数块带负号。令 $\boldsymbol{c}$ 为把 $\boldsymbol{c}_{1}$、$\boldsymbol{c}_{2}$ 上下接成的八维列向量，得到

$$ \boldsymbol{Mc} = \begin{pmatrix} \boldsymbol{B}(\theta_{0}) & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{B}( - \theta_{0}) \\ \boldsymbol{B}(0) & - \boldsymbol{B}(0) \\ \boldsymbol{D}_{1}(0) & - \boldsymbol{D}_{2}(0) \end{pmatrix}\begin{pmatrix} \boldsymbol{c}_{1} \\ \boldsymbol{c}_{2} \end{pmatrix} = \boldsymbol{0} \tag{29} $$

式（29） 的第 1、2 行对应上自由边，第 3、4 行对应下自由边，第 5、6 行对应界面应力连续，第 7、8 行对应界面位移连续。矩阵中的零为两行四列的零块，等式右端为八维零列向量。

对任意给定的 $\lambda$，如果 $\boldsymbol{M}$ 可逆，八个系数只能全为零。要存在非零模态，$\lambda$ 就必须使行列式为零。对区间 $0 < \lambda < 1$ 内的实根，定义奇异阶次 $p$。

$$ \det(\boldsymbol{M}) = 0,\quad\quad p = 1 - \lambda \tag{30} $$

载荷激发该模态时，应力按 $r^{- p}$ 增长。$p$ 越大，靠近端点时增长越快，因此若找到多个实根，应取最小的 $\lambda$ 来比较主导奇异阶次。位移的因子是 $r^{\lambda}$，说明位移可以趋于零而其导数和应力同时增大。

四、数值求根与参数分析

1、参数设置与求根检验

取两侧实际泊松比 $\nu_{1} = \nu_{2} = 0.3$，按式（11） 换算后有 ${\overline{\nu}}_{1} = {\overline{\nu}}_{2} \approx 0.428571$。以 ${\overline{E}}_{1}$ 为模量单位。因两侧的换算因子相同，等效模量比与原模量比相等。记模量比为 $m$、材料差异参数为 $t$、辅助角为 $\alpha$。

$$ m = \frac{{\overline{E}}_{2}}{{\overline{E}}_{1}} = \frac{E_{2}}{E_{1}},\quad\quad t = \frac{m - 1}{m + 1},\quad\quad\alpha = \frac{\pi\lambda}{2} \tag{31} $$

在 $- 2 \leq \log_{10}m \leq 2$ 内按 0.05 的步长取 81 个点。利用自由边及面力连续条件消去多余系数，式（29） 可化为标量方程 $^{[4]}$。令 $S = \sin\alpha$、$C = \cos\alpha$，并用等效泊松比定义系数 $b,c$ 进行代换。

$$ b = \frac{2}{1 + {\overline{\nu}}_{1}},\quad\quad c = \frac{1 - {\overline{\nu}}_{1}}{1 + {\overline{\nu}}_{1}} \tag{32} $$

标量特征方程为

$$ F(\lambda) = t^{2}\lbrack(b\lambda)^{2} - (\lambda^{2} + cS^{2})^{2}\rbrack - b^{2}S^{2}C^{2} = 0 \tag{33} $$

在 $0 < \lambda < 1$ 内，根据数值计算方法，可以通过如下方式减小端点值对求根的影响，提高数值稳定性 $^{[6]}$，故求解

$$ \frac{F(\lambda)}{\lambda^{2}(1 - \lambda)} = 0 \tag{34} $$

计算中把平方项之差拆成乘积，并用 $b\lambda - \lambda^{2} - cS^{2} = (1 - \lambda)(\lambda - c) + cC^{2}$ 来替换其中趋于零的部分，避免大数相减导致有效数字丢失。将区间两端各缩进 $10^{- 8}$，用 10001 个等距样点寻找变号区间，再用二分法求根 $^{[4]}$。除 $m = 1$ 外，各取样点均找到一个实根。

对每个候选根，将式（29） 的各行除以各自行范数，记为 $\boldsymbol{N}$。以最小奇异值对应的单位右奇异向量作为模态系数，代回八条条件。对检出的全部实根，$\boldsymbol{N}$ 的最小奇异值及 $\boldsymbol{Nc}$ 的各分量绝对值均小于 $10^{- 14}$。

2、模量比对奇异阶次的影响

由图2 可见，模量比从 1 向任一侧偏离时，奇异阶次增大。比如 $m = 2$、10 和 100 时，$p$ 分别约为 0.0370、0.1985 和 0.2780。相同倍数的模量差异给出相同的阶次。把 $m$ 换成 $m^{- 1}$，式（33） 中 $t$ 的符号反转，$t^{2}$ 保持不变。数值结果中各对倒数模量比的 $p$ 之差均小于 $10^{- 14}$。

当 $m = 1$ 时，两侧是同一种材料，式（33） 在 $0 < \lambda < 1$ 内只剩严格小于零的 $- b^{2}S^{2}C^{2}$，因此无奇异实根。图2 在 $m = 1$ 处用空心点标出邻近参数下 $p \rightarrow 0$ 的极限。模量比在 1 附近变化时，阶次十分接近零，例如模量比为 1.25 时，阶次只有约 0.0043。

图2　奇异阶次随弹性模量比的变化

3、应力径向因子随距离的变化

取参考长度 $L$，用 $\rho$ 表示归一化距离。以 $\rho = 1$ 处为基准，归一化径向因子为 $\rho^{\lambda - 1}$。

$$ \rho = \frac{r}{L} \tag{35} $$

图3 比较 $m = 2$、10 和 100 三种情形。在双对数坐标下，曲线的斜率为 $- p$。当 $\rho$ 从 1 减至 $10^{- 4}$ 时，三个归一化因子分别增至 1.41、6.22 和 12.95。以 $m = 100$ 为例，距离每缩小到原来的十分之一，径向因子约乘以 1.90。

图3　不同弹性模量比下的归一化应力径向因子

4、改变泊松比后的奇异性

在几何与模量比固定的情况下，泊松比也会改变两侧位移如何相合。取 $m = 1.5$、$\nu_{1} = 0.30$，将材料 2 的实际泊松比 $\nu_{2}$ 作为变量。为判断奇异根是否存在，采用 Bogy 对两个黏结直角楔体给出的判据 $^{[4]}$。其中 $E_{i},\nu_{i}$ 均为材料的实际参数，$\mu_{i}$ 为剪切模量。式（18） 中的 $\kappa_{i}$ 也可用实际泊松比计算。

$$ \mu_{i} = \frac{E_{i}}{2(1 + \nu_{i})},\quad\quad\kappa_{i} = 3 - 4\nu_{i},\quad\quad D_{0} = \mu_{1}(\kappa_{2} + 1) + \mu_{2}(\kappa_{1} + 1) \tag{36} $$

用 $\alpha_{D},\beta_{D}$ 表示 Dundurs 的两个无量纲材料参数。令 $\Theta$ 为判断幂律奇异性的量。

$$ \begin{gathered} \alpha_{D} = \frac{\mu_{1}(\kappa_{2} + 1) - \mu_{2}(\kappa_{1} + 1)}{D_{0}} \\ \beta_{D} = \frac{\mu_{1}(\kappa_{2} - 1) - \mu_{2}(\kappa_{1} - 1)}{D_{0}} \\ \Theta = \alpha_{D}(\alpha_{D} - 2\beta_{D}) \end{gathered} \tag{37} $$

当 $\Theta > 0$ 时存在 $0 < \lambda < 1$ 的幂律奇异根。当 $\Theta < 0$ 时，端点处无该类奇异根。$\Theta = 0$ 是两者之间的临界情形。对本组参数，临界线对应 $\alpha_{D} - 2\beta_{D} = 0$。代入式（36） 和式（37），该条件可写为

$$ \mu_{1}\nu_{2} = \mu_{2}\nu_{1},\quad\quad m = \frac{\nu_{2}(1 + \nu_{2})}{\nu_{1}(1 + \nu_{1})} \tag{38} $$

把 $m = 1.5$ 和 $\nu_{1} = 0.30$ 代入式（38），求得临界泊松比约为 0.41378。为计算表1 中的指数，沿用 第四节第1小节的 $S,C$，记 $H = S^{2} - \lambda^{2}$。对泊松比不同的情形，使用 Bogy 的通用特征方程 $^{[4]}$。

$$ \beta_{D}^{2}H^{2} + 2\alpha_{D}\beta_{D}\lambda^{2}H - \alpha_{D}^{2}\lambda^{2}(1 - \lambda^{2}) + S^{2}C^{2} = 0 \tag{39} $$

求根时，将式（39） 左边除以 $1 - \lambda$，并在靠近 $\lambda = 1$ 的区间加密取点。找到变号区间后用二分法求根。将实根代回式（29） 核验，行归一化后的八条条件残差均小于 $10^{- 14}$。

表1　模量比为 1.5 时改变泊松比的特征根

注　特征指数和奇异阶次留空处表示在所考察的区间内无奇异实根。

泊松比提高后，奇异阶次逐渐降低。取 $\nu_{2} = 0.43$ 时，$\Theta < 0$，已越过临界线。下面取 0.30 和 0.43 做有限元对照。

五、COMSOL 有限元对照

1、几何、载荷和网格

在 COMSOL 6.2 的固体力学模块中建立二维平面应变模型。图4 由两个边长均为 $L = 10\, mm$ 的正方形组成。两块材料的公共边为 $y = 0$，其左端点 $O$ 对应图1 的原点，界面上有 $r = x$。几何采用联合体并保留内部边界，使两侧在界面上共用节点。

在上、下外边界分别沿正、负 $y$ 方向施加大小相同的均布面力 $T = 1\, MPa$，左右外边界保持自由。用刚体运动抑制去除整体平移和转动，进行稳态计算。A、B 两组均取 $E_{1} = 2\, GPa$、$E_{2} = 3\, GPa$ 和 $\nu_{1} = 0.30$，只改变 $\nu_{2}$。C 组另取临界泊松比的计算值。D 组取 $E_{2} = 20\, GPa$、$\nu_{2} = 0.30$，用于对照较强的奇异情形。

图4　有限元几何、边界载荷和界面取样位置（A、B、C 组）

采用二次三角形单元，并在 $O$ 点附近加密。粗、细网格分别有 4372 和 19112 个单元。沿界面紧邻 $O$ 的第一段网格边长分别约为 $8.81 \times 10^{- 8}\, m$ 和 $2.13 \times 10^{- 9}\, m$。A、B 两组均在这两种网格上求解。

2、界面正应力与理论指数对照

在黏结界面上，$\sigma_{y}$ 就是前文的法向应力 $\sigma_{\theta}$。沿界面取 281 个按对数分布的点，分别读取两侧的应力并取平均。拟合使用两种网格共有的区间 $3 \times 10^{- 5} \leq \rho \leq 10^{- 3}$，归一化距离 $\rho$ 的定义见式（35）。令 $a$ 为截距、$s_{FE}$ 为拟合斜率，用最小二乘法 $^{[4]}$ 拟合

$$ \ln\left( \frac{\sigma_{y}}{T} \right) \approx a + s_{FE}\ln\rho \tag{40} $$

若奇异项在这段距离内占主导，双对数斜率 $s_{FE}$ 应接近 $- p$。表2 比较理论阶次和计算斜率。图5 中，阴影为拟合区间，曲线和空心符号分别为细网格和粗网格结果。图5(a) 的纵轴为线性刻度，图5(b) 为对数刻度。

图5　界面法向应力随距左端点距离的变化

表2　理论奇异阶次与 COMSOL 界面正应力

注　末列为细网格在 $\rho = 10^{- 4}$ 处的结果。C 组按式（38） 的临界值计算，其泊松比在表中保留五位小数。粗网格栏留空表示该算例仅计算了细网格。

A 组越靠近端点，法向应力越大，细网格斜率约为 $- 0.01356$，接近理论值 $- 0.01367$。B 组接近端点时的法向应力略有降低，符合 第四节第4小节无幂律奇异根的结果。拟合区间内，A、B 两组粗细网格的正应力最大相对差异分别为 0.034% 和 0.017%。图5(b) 的 D 组上升更明显，斜率约为 $- 0.19880$，也接近理论值 $- 0.19847$。

3、临界算例的均匀应力

C 组的应力曲线约为 1 MPa 的水平线。可以用式（12） 解释这一结果。设两侧都处于 $\sigma_{y} = T$、$\sigma_{x} = \tau_{xy} = 0$ 的均匀应力状态。沿界面方向的正应变为

$$ \varepsilon_{x}^{(i)} = - \frac{{\overline{\nu}}_{i}T}{{\overline{E}}_{i}} = - \frac{\nu_{i}(1 + \nu_{i})T}{E_{i}} \tag{41} $$

式（38） 恰好使两种材料的 $\varepsilon_{x}^{(i)}$ 相等，因此两侧沿界面的位移相同。取界面上的法向位移为零，法向位移也连续。这个均匀场也满足上下表面的拉伸载荷和左右自由边条件，所以临界配对在本算例中给出均匀应力。

4、自由边附近的等效应力

图6对照了四个算例在黏结界面左端附近的 von Mises 等效应力。算例 A 在端点附近出现局部高应力区，算例 B 的分布相对平缓。

算例 C 的等效应力在各材料内部近似均匀，材料 1 和材料 2 分别约为 0.889 MPa 和 0.870 MPa。平面应变下，两侧泊松比不同，面外正应力也不同，因此界面法向应力同为约 1 MPa 时，两侧的等效应力仍可不同。算例 D 的高应力集中在端点附近，与图5(b) 中界面法向应力快速上升的趋势相符。

（a）算例 A	（b）算例 B
（c）算例 C，临界泊松比	（d）算例 D，模量比为 10

图6　界面左端附近的 von Mises 等效应力云图，mph文件在 /mph 文件夹中

注　各分图使用独立色标。根据 COMSOL 等效应力取样值绘制，均采取放大处理（无旋转），聚焦原点处。

六、工程应用

芯片与聚合物基板、封装胶与基板的接合面常不覆盖在外表面，温度变化引起的变形差异可使界面端部出现局部剥离应力和剪应力 $^{[5]}$。在接近本文直角黏结的部位，可以把材料在工作温度下的 $E_{i}$ 和 $\nu_{i}$ 代入式（37），用 $\Theta$ 初步筛选候选配对。若调整配方同时改变了弹性模量和泊松比，就用更新后的两组参数重算。

例如，若接头刚度要求使两种材料的模量比约为 1.5，而材料 1 的泊松比约为 0.30，可用 0.41378 作为这一理想模型的筛选参考。可先用供应商数据或材料试验筛选接合材料，再把筛选后的实际参数和温度载荷放入有限元模型。有限元比较可在同一物理距离处读取界面正应力和切应力，便于把不同网格的结果进行对照。对有胶层厚度、端部圆角或温度循环的接头，应按实际几何、热膨胀系数和界面强度评估可靠性。

七、讨论与结论

本算例表明，外表面平直的黏结结构也可以出现端点应力奇异性。在两侧泊松比均为 0.3 的条件下，弹性模量相同时奇异实根消失，模量比越偏离 1，端点应力发散得越快。模量比较大时，继续拉大两侧刚度差异，奇异阶次仍增加，但曲线的变化趋缓。

这种奇异性来自自由边与黏结界面条件在端点相遇。相同面力作用下，刚度不同的材料一般产生不同应变，界面却要求两侧位移相合。应力和位移需要在端点附近同时满足上述条件，位移表达式中的刚度差异由此改变了允许的径向幂次。

八、参考文献

[1] Bogy D B. Edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading[J]. Journal of Applied Mechanics, 1968, 35(3): 460-466. DOI:10.1115/1.3601236.

[2] Bogy D B. Two edge-bonded elastic wedges of different materials and wedge angles under surface tractions[J]. Journal of Applied Mechanics, 1971, 38(2): 377-386. DOI:10.1115/1.3408786.

[3] 吴家龙. 弹性力学[M]. 3版. 北京: 高等教育出版社, 2016.

[4] Bogy D B. On the problem of edge-bonded elastic quarter-planes loaded at the boundary[J]. International Journal of Solids and Structures, 1970, 6(9): 1287-1313. DOI:10.1016/0020-7683(70)90104-6.

[5] Liu X H, Suo Z, Ma Q. Split singularities: stress field near the edge of a silicon die on a polymer substrate[J]. Acta Materialia, 1999, 47(1): 67-76. DOI:10.1016/S1359-6454(98)00345-0.

[6] 丁丽娟, 程杞元. 数值计算方法[M]. 北京: 高等教育出版社, 2011.