正交各向异性带型中裂纹诱导的温度场分析

Analysis of the Temperature Field Induced by a Crack in an Orthotropic Strip

摘要: 本文研究热荷载下正交各向异性带型中的裂纹问题。利用傅里叶变换，将热边值问题转化为奇异积分方程。采用Lobatto-Chebyshev公式得到线性代数方程系统。给出了裂纹诱导的温度场的数值解，通过数值计算分析了裂缝表面温度差的变化规律，采用图型表明了裂纹位置和裂纹尺寸对温度场的影响。

Abstract: This paper investigates the problem of a crack embedded in an orthotropic strip under a thermal loading. Using the Fourier transform, the thermal boundary value problem is transformed into a singular integral equation. The Lobatto-Chebyshev formula is used to further derive the linear system of algebraic equations. Through the numerical calculation, the temperature jump across the crack surface is obtained. Finally, the influence of crack position and crack size on the temperature field is reflected in the form of a graph.

1. 引言

断裂力学是研究带有裂纹的物体在载荷的作用下裂纹扩展规律的一门学科 [1] 。我们要从晶体学的角度去认识断裂现象的本质。晶体一般由原子(或离子、分子)规则排列而成 [2] 。在一定条件下，比如温度、剪应力、均匀热流，晶体中原子的结合力遭到了破坏，从而诱导了裂纹的产生，当裂纹扩展到一定程度便会发生断裂。固体材料或结构在外力作用下会发生断裂，这不仅在生活中，甚至在工程上都是极为普遍的现象，一个看似简单的现象，其物理本质却难以完全认识透彻 [2] 。断裂现象的研究也因此成为了众多学者们孜孜不倦的追求。复合型材料在力、热、电、磁、光，作用下的断裂现象逐渐引起学者们的关注 [3] 。断裂力学发展至今，学者们已经研究了各种裂纹模式的断裂问题。例如，文献 [4] 中研究了正交各向异性体双线性裂纹的平面热弹性问题，文献 [5] 在开裂纹介质情况下提出了一种热介质裂纹模型，用于研究在热和机械载荷条件下横观各向同性材料的Griffith裂纹问题，文献 [6] 研究了机械、均匀热流载荷条件下正交各向异性体的共线双裂纹问题。再如文献 [7] [8] [9] 研究了单、双裂纹在不同材料、不同影响因子条件下的裂纹问题。

本文侧重研究裂纹位置、裂纹尺寸对正交各向异性带型裂纹温度场的影响，基于热断裂物理边值问题转化为奇异积分方程的思想，通过数值计算求解奇异积分方程，获得了与实际生活的断裂现象较吻合的数值结果，表明了该方法的有效性。本文的研究不仅给出了裂纹位置、裂纹尺寸对温度场的影响，还为热弹性场的分析提供了研究基础，具有一定的实际意义。

2. 研究对象

考虑如图1所示嵌入在正交各向异性带型中的Griffith裂纹问题，该裂纹与带型面平行。建立笛卡尔坐标系，裂纹位于线段 $x\in \left[-a,a\right]$ ，此时 $y=0$ 。假设正交各向异性带型的厚度为 $h_I+h_{II}=2h$ ，上下表面的温度分别为 ${\theta }_0$ 和。研究裂纹的偏移位置、裂纹尺寸对温度场产生的影响，以及裂纹表面温度差的分布情况。正交各向异性带型的温度场满足以下方程 [10] ：

$q_x=-{\lambda }_x\frac{\partial \theta }{\partial x}$ ， $q_y=-{\lambda }_y\frac{\partial \theta }{\partial y},$

其平衡方程为

$\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}=0$

从而得到以下的控制偏微分方程

Figure 1. A crack in an orthotropic strip

图1. 正交异性带型中的裂纹

${\lambda }^2\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}=0,$ (1)

其中 ${\lambda }^2={\lambda }_x/{\lambda }_y$ 。

其热边界条件为：

${\theta }^{\left(I\right)}\left(x,h_I\right)={\theta }_0,{\theta }^{\left(II\right)}\left(x,-h_{II}\right)=-{\theta }_0,x>0,$ (2)

${\theta }^{\left(I\right)}\left(x,0\right)={\theta }^{\left(II\right)}\left(x,0\right),x>a,$ (3)

$q_y^{\left(I\right)}\left(x,0\right)=q_y^{\left(II\right)}\left(x,0\right),x>a,$ (4)

$q_y^{\left(I\right)}\left(x,0\right)=q_y^{\left(II\right)}\left(x,0\right)=q_c,0<x<a.$ (5)

是热流密度， $q_x,q_y$ 分别是 $x$ 轴和 $y$ 轴方向的热流分量， $\theta$ 是温度，而 $\lambda ,{\lambda }_x,{\lambda }_y$ 是方程热传导系数，它们具有重要的物理意义。

我们将研究偏微分方程(1)在以上热物理边值条件下的数值解。综合考虑带型边界初始温度、裂纹位置、裂纹尺寸对温度场的影响。

3. 求解过程

傅里叶变换在求解微分方程中有着重要的应用。其解题基本思想是：通过傅里叶变换将复杂的方程化成简单的代数方程，然后求解代数方程，再对代数方程的解作反变换，就得到了原方程的解 [11] 。本文通过傅里叶变换，用物理边值条件(2)~(5)来求解偏微分方程(1)的解，即本文重点研究的温度场。再采用Lobatto-Chebyshev积分公式 [12] 求解温度场的数值解。

在求解之前，我们先引进傅里叶余弦变换 [11]

$\bar{f}\left(\xi \right)=\frac{2}{\text{π}}{\displaystyle {\int }_0^{+\infty }f\left(x\right)\cos \left(\xi x\right)}\text{d}x,$

$f\left(x\right)={\displaystyle {\int }_0^{+\infty }\bar{f}}\left(\xi \right)\cos \left(\xi x\right)\text{d}\xi .$

以及傅里叶变换的性质之一

$f^{\left(n\right)}\left(x\right)={\left(-i\xi \right)}^n\bar{f}\left(\xi \right).$

同时，介绍如下求解常微分方程的特征值法。求解常微分方程

$\frac{{\text{d}}^{2}f}{\text{d}{y}^2}+p\frac{\text{d}f}{\text{d}y}+qf=0,$ (6)

其特征方程为

$m^2+pm+q=0.$

当特征方程存在两个不同实数根，设为 $m_1$ ， $m_2$ ，则该微分方程的通解可表为

$f=c_1{e}^{m_{1}y}+c_2{e}^{m_{2}y},$ (7)

当特征方程存在两个相同的实数根，即，则该微分方程的通解可表为

$f=c_1{e}^{m_{1}y}+c_{2}y{e}^{m_{1}y},$ (8)

其中 $c_1$ ， $c_2$ 为任意常数。又由于指数函数与双曲函数之间的转换关系

$e^x=\sinh \left(x\right)+\cosh \left(x\right),$

所以(7)、(8)可分别表示为

$f={\displaystyle \underset{i=1}{\overset{2}{\sum }}{c}_i\left[\cosh \left(m_{i}y\right)+\sinh \left(m_{i}y\right)\right]},$ (9)

$f=\left(c_1+c_{2}y\right)\left[\cosh \left(m_{1}y\right)+\sinh \left(m_{1}y\right)\right].$ (10)

微分方程(6)的通解(9)或(10)将运用在温度场的求解过程中。

3.1. 温度场的求解

假设 $I$ 和 $II$ 分别表示在坐标系中 $y>0$ 和 $y<0$ 的区域。由于该裂纹关于 $y$ 轴对称，所以本文只取 $y$ 轴右侧的带型进行研究。为了满足温度场的物理边值条件，温度场的解由两部分组成 [13] ，一部分为与 $x$ 和 $y$ 都有关的 ${\theta }_1\left(x,y\right)$ ，另一部分为只与 $$ 有关的 ${\theta }_2\left(y\right)$ ，即 $\theta \left(x,y\right)={\theta }_1\left(x,y\right)+{\theta }_2\left(y\right)$ 。其中 ${\theta }_1\left(x,y\right)$ 满足方程(1)，对方程(1)两边取关于变量 $x$ 的傅里叶余弦变换，求解后再作反变换，获得温度场方程(1)的解 ${\theta }_1\left(x,y\right)$ ，可表示为：

${\theta }_1^{\left(j\right)}\left(x,y\right)=\frac{2}{\text{π}}{\displaystyle {\int }_0^{\infty }\left[A^{\left(j\right)}\left(\xi \right)\cosh \left(\lambda \xi y\right)+B^{\left(j\right)}\left(\xi \right)\sinh \left(\lambda \xi y\right)\right]}\cos \left(\xi x\right)\text{d}\xi ,$

其中 $j=I,II$ 。此外， ${\theta }_2\left(y\right)$ 满足如下3个方程

$\frac{{\partial }^2{\theta }_2\left(y\right)}{\partial y}=0,$ (11)

${\theta }_2\left(h_I\right)={\theta }_0,$ (12)

${\theta }_2\left(-h_{II}\right)=-\theta$ (13)

用特征方程法求解常微分方程(11)，获得 ${\theta }_2\left(y\right)$ 的通解

${\theta }_2\left(y\right)=A_{0}y+B_0,$

又根据(12)、(13)，求得

$A_0=\frac{2{\theta }_0}{h_I+h_{II}},B_0=\frac{{\theta }_0\left(h_{II}-h_I\right)}{h_I+h_{II}}.$

所以，温度场的解为

${\theta }^{\left(j\right)}\left(x,y\right)=\frac{2}{\text{π}}{\displaystyle {\int }_0^{\infty }\left[A^{\left(j\right)}\left(\xi \right)\cosh \left(\lambda \xi y\right)+B^{\left(j\right)}\left(\xi \right)\sinh \left(\lambda \xi y\right)\right]}\cos \left(\xi x\right)\text{d}\xi +A_{0}y+B_0.$ (14)

其中 $A^{\left(j\right)}\left(\xi \right)$ ， $B^{\left(j\right)}\left(\xi \right)$ 为待求函数。将(14)应用于热边界条件(2)，得 $A^{\left(j\right)}\left(\xi \right)$ 与 $B^{\left(j\right)}\left(\xi \right)$ 的关系

$B^{\left(j\right)}\left(\xi \right)={\delta }^{\left\{j\right\}}\coth \left(\lambda \xi h_i\right)A^{\left(j\right)}\left(\xi \right),$

将上式代入(14)中，温度场的解可表示为

${\theta }^{\left(j\right)}\left(x,y\right)=\frac{2}{\text{π}}{\displaystyle {\int }_0^{\infty }\frac{\sinh \left[\lambda \xi \left(h_j+{\delta }^{\left\{j\right\}}y\right)\right]}{\sinh \left(\lambda \xi h_j\right)}{A}^{\left(j\right)}\left(\xi \right)}\cos \left(\xi x\right)\text{d}\xi +A_{0}y+B_0,$ (15)

综上，

$q_y^{\left(j\right)}\left(x,y\right)=-{\lambda }_y{\delta }^{\left\{j\right\}}\frac{2}{\text{π}}{\displaystyle {\int }_0^{\infty }\frac{\cosh \left[\lambda \xi \left(h_j+{\delta }^{\left\{j\right\}}y\right)\right]}{\sinh \left(\lambda \xi h_j\right)}}{A}^{\left(j\right)}\left(\xi \right)\cos \left(\xi x\right)\text{d}\xi -{\lambda }_y{A}_0.$ (16)

其中， ${\delta }^{\left\{I\right\}}=-1$ ， ${\delta }^{\left\{II\right\}}=1$ 。

为了求解 $A^{\left(j\right)}\left(\xi \right)$ ，首先引进辅助函数

${\phi }_1\left(x\right)=\frac{\partial \left[{\theta }^{\left(I\right)}\left(x,0\right)-{\theta }^{\left(II\right)}\left(x,0\right)\right]}{\partial x},$ (17)

根据温度场的解(15)得

${\theta }^{\left(I\right)}\left(x,0\right)=\frac{2}{\text{π}}{\displaystyle {\int }_0^{\infty }{A}^{\left(I\right)}\left(\xi \right)\cos \left(\xi x\right)\text{d}\xi +B_0},$ (18)

${\theta }^{\left(II\right)}\left(x,0\right)=\frac{2}{\text{π}}{\displaystyle {\int }_0^{\infty }{A}^{\left(II\right)}\left(\xi \right)\cos \left(\xi x\right)\text{d}\xi +B_0},$ (19)

将方程(18)、(19)代入(17)且由边值条件(3)，可得

$A^{\left(I\right)}\left(\xi \right)-A^{\left(II\right)}\left(\xi \right)=-\frac{1}{\xi }{\displaystyle {\int }_0^a{\phi }_1\left(s\right)\sin \left(\xi s\right)\text{d}s}.$ (20)

然后，由条件(4)得

$A^{\left(II\right)}\left(\xi \right)=-\frac{\coth \left(\lambda \xi h_I\right)}{\coth \left(\lambda \xi h_{II}\right)}{A}^{\left(I\right)}\left(\xi \right),$ (21)

联立(20)、(21)有

$A^{\left(I\right)}\left(\xi \right)=-\frac{1}{\xi }{\displaystyle {\int }_0^a\frac{\coth \left(\lambda \xi h_{II}\right)}{\coth \left(\lambda \xi h_I\right)+\coth \left(\lambda \xi h_{II}\right)}}{\phi }_1\left(s\right)\sin \left(\xi s\right)\text{d}s,$ (22)

$A^{\left(II\right)}\left(\xi \right)=\frac{1}{\xi }{\displaystyle {\int }_0^a\frac{\coth \left(\lambda \xi h_I\right)}{\coth \left(\lambda \xi h_I\right)+\coth \left(\lambda \xi h_{II}\right)}}{\phi }_1\left(s\right)\sin \left(\xi s\right)\text{d}s,$ (23)

将方程(22)、(23)分别代入方程(18)、(19)，可得裂纹面的温度分布公式

${\theta }^{\left(I\right)}\left(x,0\right)=-\frac{\text{2}}{\text{π}}{\displaystyle {\int }_0^a\left[{\displaystyle {\int }_0^{\infty }\frac{1}{\xi }\frac{\coth \left(\lambda \xi h_{II}\right)}{\coth \left(\lambda \xi h_I\right)+\coth \left(\lambda \xi h_{II}\right)}\sin \left(\xi s\right)\cos \left(\xi x\right)\text{d}\xi }\right]}{\phi }_1\left(s\right)\text{d}s+B_0,$

将(22)代入方程(16)，并根据边界条件(5)得到以下积分方程

$\frac{2}{\text{π}}{\displaystyle {\int }_0^a\left[{\displaystyle {\int }_0^{\infty }\frac{1}{\xi }\frac{{\lambda }_y\lambda \xi \coth \left(\lambda \xi h_I\right)\coth \left(\lambda \xi h_{II}\right)}{\coth \left(\lambda \xi h_I\right)+\coth \left(\lambda \xi h_{II}\right)}\sin \left(\xi s\right)\cos \left(\xi x\right)\text{d}\xi }\right]}{\phi }_1\left(s\right)\text{d}s=-\left(q_c+{\lambda }_y{A}_0\right).$ (24)

3.2. 奇异积分方程的求解

奇异积分方程在断裂力学方面有着重要的应用。对奇异积分方程数值解的研究，不少学者们已经取得了相对成熟的成果。例如，配位法、基于Fourier展开的Galekrin-Petrvo方法、样条方法等 [14] 。而本文采用的是Lobatto-Chebyshev积分公式。

为了方便求解奇异积分方程，假设

$M\left(\xi \right)=\frac{{\lambda }_y\lambda \xi \coth \left(\lambda \xi h_I\right)\coth \left(\lambda \xi h_{II}\right)}{\coth \left(\lambda \xi h_I\right)+\coth \left(\lambda \xi h_{II}\right)}$ ，

所以，方程(24)可以写为

$\frac{2}{\text{π}}{\displaystyle {\int }_0^a\left[{\displaystyle {\int }_0^{\infty }\frac{M\left(\xi \right)}{\xi }\sin \left(\xi s\right)}\cos \left(\xi x\right)\text{d}\xi \right]}{\phi }_1\left(s\right)\text{d}s={\tau }_0.$ (25)

根据文献 [15] ，有

$\frac{2}{\text{π}}{\displaystyle {\int }_0^{\infty }\sin \left(\xi s\right)\cos \left(\xi x\right)\text{d}\xi =\frac{1}{\text{π}}}\left(\frac{1}{s-x}+\frac{1}{s+x}\right),$

且假设

$\Lambda =\underset{x\to +\infty }{\lim }\frac{M\left(\xi \right)}{\xi },$

$\Lambda$ 是由材料性质决定的常数，可以通过数值计算解得。所以方程(25)可化为以下奇异积分方程

$\frac{1}{\text{π}}{\displaystyle {\int }_{-a}^a\frac{\Lambda }{s-x}}{\phi }_1\left(s\right)\text{d}s+\frac{1}{\text{π}}{\displaystyle {\int }_{-a}^{a}k\left(s,x\right){\phi }_1}\left(s\right)\text{d}s={\tau }_0,$ (26)

其中，

$k\left(s,x\right)={\displaystyle {\int }_0^{\infty }\left[\frac{M\left(\xi \right)}{\xi }-\Lambda \right]}\sin \left(\xi s\right)\cos \left(\xi x\right)\text{d}\xi .$

引进以下变量变换

$x=a\bar{x}$ ， $s=a\bar{s}$ ， ${\phi }_1\left(s\right)=\frac{Y\left(\bar{s}\right)}{\sqrt{1-{\bar{s}}^2}},$

方程(26)可表示为

$\frac{1}{\text{π}}{\displaystyle {\int }_{-1}^1\frac{\Lambda }{\bar{s}-\bar{x}}}\frac{Y\left(\bar{s}\right)}{\sqrt{1-{\bar{s}}^2}}\text{d}\bar{s}+\frac{1}{\text{π}}{\displaystyle {\int }_{-1}^1\bar{k}\left(\bar{s},\bar{x}\right)\frac{Y\left(\bar{s}\right)}{\sqrt{1-{\bar{s}}^2}}}\text{d}\bar{s}={\tau }_0,$ (27)

其中，

$\bar{k}\left(\bar{s},\bar{x}\right)=a{\displaystyle {\int }_0^{\infty }\left[\frac{M\left(\xi \right)}{\xi }-\Lambda \right]}\sin \left(\xi a\bar{s}\right)\cos \left(\xi a\bar{x}\right)\text{d}\xi .$

由单值性条件(3)得

${\displaystyle {\int }_{-1}^1{\phi }_1}\left(x\right)\text{d}x=0.$ (28)

一般的，奇异积分方程(27)可以通过数值计算方法求解。本文对方程(27)、(28)应用

Lobatto-Chebyshev积分公式，获得以下代数方程组

$\frac{1}{L}{\displaystyle \underset{l=0}{\overset{L}{\sum }}\frac{{\lambda }_l\Lambda }{{\bar{s}}_l-{\bar{x}}_m}}Y\left({\bar{s}}_l\right)+\frac{1}{L}{\displaystyle \underset{l=0}{\overset{L}{\sum }}{\lambda }_l}\bar{k}\left(\bar{s},\bar{x}\right)Y\left({\bar{s}}_l\right)={\tau }_0,$ (29)

${\displaystyle \underset{l=0}{\overset{L}{\sum }}{\lambda }_l}Y\left({\bar{s}}_l\right)=0,$ (30)

其中

${\bar{x}}_m=\cos \left[\left(2m-1\right)\text{π}/\left(2L\right)\right]$ ， $m=1,2,\cdots ,L,$

${\bar{s}}_l=\cos \left(l\text{π}/L\right)$ ， $l=0,1,2,\cdots ,L,$

${\lambda }_0={\lambda }_L=1/2$ ， ${\lambda }_1=\cdots ={\lambda }_{L-1}=1.$

在上式中的 $L$ 取偶数。

只要代数方程(29)、(30)得以求解，那么裂纹上下表面的温度差就可以确定了。比如，

$\begin{array}{c}{\displaystyle {\int }_x^a{\phi }_1}\left(s\right)\text{d}s={\displaystyle {\int }_x^a\frac{\partial \left[{\theta }^{\left(I\right)}\left(s,0\right)-{\theta }^{\left(II\right)}\left(s,0\right)\right]}{\partial s}}\text{d}s\\ =\left[{\theta }^{\left(I\right)}\left(s,0\right)-{\theta }^{\left(II\right)}\left(s,0\right)\right]|\begin{array}{c}a\\ x\end{array}\\ =-\left[{\theta }^{\left(I\right)}\left(x,0\right)-{\theta }^{\left(II\right)}\left(x,0\right)\right],\end{array}$

所以，

${\theta }^{\left(I\right)}\left(x,0\right)-{\theta }^{\left(II\right)}\left(x,0\right)=-a{\displaystyle {\int }_{\bar{x}}^1\frac{Y\left(\bar{s}\right)}{\sqrt{1-{\bar{s}}^2}}}\text{d}\bar{s}.$

特别地，当 $\bar{x}=0$ ，

${\theta }^{\left(I\right)}\left(x,0\right)-{\theta }^{\left(II\right)}\left(x,0\right)\approx -a{\displaystyle \underset{l=0}{\overset{L/2}{\sum }}{\lambda }_l}Y\left({\bar{s}}_l\right).$ (31)

综上，我们便获得了裂纹上下表面温度差的数值结果。

4. 数值实例

例：为了体现本文方法的有效性，取陶瓷材料Tyrannohex [6] 的参数进行数值计算，有 ${\lambda }_x=3.08$ ， ${\lambda }_y=2.81$ 。取 $L=10$ ， $q_c=0$ 。其中，通过数值计算得 $\Lambda =1.471$ 。

根据公式(29)、(30)、(31)，我们可以得到以下数值结果。当我们令裂纹面的热流密度 $q_c=0$ ，则该带型裂纹就变成热不渗透裂纹，此时只考虑裂纹位置、裂纹尺寸、带型边界初始温度对温度场的影响。图2中的 $a/h$ 可表示裂纹相对于带型的尺寸， $h_I/h$ 表示裂纹相对于带型中心的位置，其中 $\Delta \theta ={\theta }^{\left(I\right)}\left(x,0\right)-{\theta }^{\left(II\right)}\left(x,0\right)$ 。通过观察图2，我们不难发现，若 $h_I/h$ 固定，当裂纹尺寸相对于带型宽度的比值越大，裂纹上下表面的温度差越大；若 $a/h$ 固定，当裂纹逐渐偏离带型边界，越向带型中间位置靠拢，裂纹上下表面的温度差越大。此外，裂纹上下表面的温度差也与带型边界初始温度成正比。

5. 结论

在本文中，通过傅里叶余弦变换和数值计算，给出了裂纹面的温度分布公式以及裂纹上下表面温度

Figure 2. The influences of crack location and size on $\Delta \theta /{\theta }_0$

图2. 裂纹位置以及裂纹尺寸对 $\Delta \theta /{\theta }_0$ 的影响曲线

差的数值结果。并且通过数值结果验证了裂纹上下表面的温度差与裂纹位置、裂纹尺寸、带型边界初始温度有关。本文不仅为裂纹温度场的求解提供了新的思路，采用的方法能够广泛的应用于各种裂纹问题的求解，所得到的结果还为热弹性场的分析奠定了研究基础。

基金项目

广西自然科学基金(2016GXNSFAA380261)。