弹性力学方程的 Green 函数

基本解和 Green 函数的区别

如果将微分方程写成抽象的算子形式: \[ L[u] = f \text{ on } \Omega \]

特别的, 当函数 \(f\) 为 Dirac 函数 \(\delta\) 时, 微分方程的解称为算子 \(L\) 的基本解: \[ L [G] = \delta \]

注意到上述过程中, 并没有提到方程需要满足什么边界条件, 这也是基本解和 Green 函数之间的区别. Green 函数需要考虑边界的影响, 基本解只需要考虑函数在区域内满足什么微分方程. 在下文中, 我还是按照 Mura 书中的表述, 将以求解基本解方式得到的函数称为 Green 函数. 当微分方程右端项为一般的函数 \(f\) 时, 此时微分方程的解 \(u\) 可以表示为Green \(G\) 与 \(f\) 的卷积:

\[ u = G*f = \int_{-\infty}^{\infty} G(x-x^{\prime}) f(x^{\prime}) \mathrm{d} x^{\prime} \]

张量值函数的卷积定义

对于张量值方程 \[ \boldsymbol{L}[\boldsymbol{u}] = \boldsymbol{f} \] 应用 Fourier 变换, 方程组转化为关于 \(\hat u\) 的代数方程: \[ \hat L_{ij}(\xi) \hat u_{j} = \hat f_i \] 对方程系数矩阵 \(\hat L_{ij}\) 求逆, 再作 Fourier 逆变换, 就得到解 \(u\) 表示为 \[ u_i = \mathscr{F}^{-1} [\hat{L}_{ij}^{-1} \hat{f}_j] \] 定义张量值函数之间的卷积为 \[ A * f: L(\Omega,\mathbb{R}^{n\times n}) \times L(\Omega,\mathbb{R}^{ n}) \mapsto L(\Omega,\mathbb{R}^{n}) \\ A * f \triangleq A_{ij} * f_{j} \] 上式默认了指标求和约定. 这也就是说, 张量的代数运算保存在卷积的积分号内.

张量值方程的基本解定义为 \[ G_{ij} = \mathscr{F}^{-1} [\hat{L}_{ij}^{-1}] \] 因此, 方程解可以表示为 \[ \boldsymbol{u} = \boldsymbol{G} * \boldsymbol{f} \]

含本征应变弹性方程的基本解

对于含本征应变的控制方程: \[ L_{ijkl}~u_{k,lj} - L_{ijkl}~\mu_{kl,j} = 0 \]

我们希望能够找到该方程的 Green 函数 \(G\), 使得方程的解可以表示为 Green 函数与本征应变的卷积. 这需要借助于 Fourier 分析方法. 对方程进行 Fourier 变换, 得到

\[ L_{ijkl}\xi_{l}\xi_{j}~\hat{u}_{k} + i \xi_{j}L_{ijkl}~\hat{\mu}_{kl} = 0 \]

记 \(K_{ik}(\xi) \triangleq L_{ijkl}\xi_{l}\xi_{j}\), \(X_i \triangleq -i \xi_{j}L_{ijkl}~\hat{\mu}_{kl}\), 就得到关于 \(\hat{u}_k\) 的线性方程:

\[ K_{ik} \hat{u}_k = X_i \]

如果 \(L_{ijkl}\) 具有主对称性, 也就是说, \(L_{(ij)[kl]} = L_{[kl](ij)}\), 那么矩阵 \(K_{ik}(\xi)\) 是对称矩阵. 再应用 Fourier 逆变换, 就可以得到解 \(u(x)\):

\[ u_{i}(x) = \mathscr{F}^{-1} [\hat{u}_{i}(\xi)] = \mathscr{F}^{-1}[ -i\xi_{j}~ L_{kjmn}~K_{ik}^{-1}(\xi)\hat{\mu}_{mn}] \]

将三维 Fourier 逆变换的公式代入得到

\[ \begin{aligned} u_{i}(x) &= - \frac{1}{(2\pi)^3} \int i\xi_{j} L_{kjmn} K_{ik}^{-1}(\xi) e^{i\xi \cdot x} \mathrm{d} \xi \int \mu_{mn}(x^{\prime}) e^{-i\xi \cdot x^{\prime}} \mathrm{d} x^{\prime}\\ &= - \int L_{kjmn} \frac{\partial}{\partial x_j} \left( \frac{1}{(2\pi)^3} \int K_{ik}^{-1}(\xi) e^{i\xi \cdot (x-x^{\prime})} \mathrm{d} \xi \right) \mu_{mn}(x^{\prime}) \mathrm{d} x^{\prime} \end{aligned} \]

定义 Green 函数为弹性方程逆矩阵的 Fourier 逆变换:

\[ G_{ik}(x) \triangleq \mathscr{F}^{-1}[K_{ik}^{-1}](x) =\frac{1}{(2\pi)^3} \int K_{ik}^{-1}(\xi) e^{i\xi \cdot x} \mathrm{d} \xi \]

那么位移 \(u_i\) 可以表示为 Green 函数与本征应变的卷积

\[ u_{i}(x) = -\int L_{jkmn} G_{ik,j}(x-x^{\prime}) \mu_{mn}(x^{\prime}) \mathrm{d} x^{\prime} \]

式中, Green 函数下标 \(\square_{,j}\) 表示对变量 \(x_j\) 求偏导 (注意区别变量 \(x\) 与 \(x^{\prime}\)). 或者将上式记为卷积的形式: \[ u_i = G_{ik,j} * \left( -L_{jkmn} \mu_{mn} \right) \] ## 验证 Green 函数所满足的微分方程

对 Green 函数求二阶导数, 得到 \[ G_{ij,kl}(x) = \frac{1}{(2\pi)^3} \int -K_{ij}^{-1}(\xi) \xi_k \xi_l e^{i\xi \cdot x} \mathrm{d} \xi \]

方程两端乘 \(L_{jknl}\), 等式依旧成立, 得到

\[ \begin{align} L_{jknl}G_{ij,kl}(x) &= -\frac{1}{(2\pi)^3} \int K_{ij}^{-1}(\xi) L_{jknl} \xi_k \xi_l e^{i\xi \cdot (x-x^{\prime})} \mathrm{d} \xi \end{align} \]

注意到 \(K_{jn}(\xi) := L_{jknl}\xi_{k}\xi_{l}\), 所以有 \[ L_{jknl}G_{ij,kl}(x) = -\delta_{in} \frac{1}{(2\pi)^3}\int e^{i\xi \cdot x} \mathrm{d} \xi = -\delta_{in} \delta(x) \] 整理一下指标, 得到 \[ L_{ijkl}G_{nk,lj}(x) + \delta_{in} \delta(x) = 0 \]

展开指标 \(i\), \(n\), 就得到 9 个方程: \[ \begin{pmatrix} L_{1jkl}G_{1k,lj} + \delta(x) & L_{1jkl}G_{2k,lj} & L_{1jkl}G_{3k,lj} \\ L_{2jkl}G_{1k,lj} & L_{2jkl}G_{2k,lj} + \delta(x) & L_{2jkl}G_{3k,lj} \\ L_{3jkl}G_{1k,lj} & L_{3jkl}G_{2k,lj} & L_{3jkl}G_{3k,lj} + \delta(x) \end{pmatrix} = \mathbf{0}_{3\times 3} \] 方程的每一列代表了空间中三个方向上作用在 \(x=0\)​ 处的单位力, Green 函数 \(G_{nk}(x)\) 表示在 \(n\) 方向的单位集中力作用下, 弹性方程 \(k\) 方向的位移响应. 如果 \(L_{ijkl}\) 具有主对称性, 那么矩阵 \(K\) \(N\) 是对称矩阵, 所以 \[ G_{nk} = G_{kn} \] 所以, \(G_{nk}(x)\) 同样也可以解释为: 在 \(k\) 方向的单位集中力作用下, 弹性方程 \(n\) 方向的位移响应. 这恰好与互易原理的结论相同.

各项同性材料的 Green 函数表示

对于各向同性材料, 方程左端的系数矩阵为 \[ \hat K \triangleq \mu \xi^{T} \xi ~\mathbf{I} + (\lambda + \mu)\xi \xi^{T} = \mu |\xi|^2 \mathbf{I} + (\lambda + \mu) \begin{pmatrix} \xi_1 \xi_1 & \xi_1 \xi_2 & \xi_1 \xi_3 \\ \xi_2 \xi_1 & \xi_2 \xi_2 & \xi_2 \xi_3 \\ \xi_3 \xi_1 & \xi_3 \xi_2 & \xi_3 \xi_3 \end{pmatrix} \] 式中, \(\mathbf{I}\) 为 3 阶单位矩阵, \(\xi = (\xi_1, \xi_2, \xi_3)^{\top}\), \(|\xi|^2 = \xi_1^2 + \xi_2^2 + \xi_3^2\). 对于形式为 \(I + \alpha vv^T\) 的矩阵, 其逆矩阵的形式依旧为 \(I + \beta vv^T\), 因此矩阵 \(\hat{K}\) 的逆为 \[ \hat{K}^{-1} = \frac{1}{\mu} \frac{1}{|\xi|^2} \mathbf{I} - \frac{\lambda + \mu}{\mu (\lambda + 2\mu)} \frac{1}{|\xi|^4} \xi \xi^{\top}, \quad \det(\hat{K}) = (\lambda + 2\mu)\mu^2 |\xi|^6 \] 各项同性材料的 Green 函数为 \[ G(x) = \frac{1}{(2\pi)^3} \frac{1}{\mu} \int_{\mathbb{R}^3} \frac{1}{|\xi|^2} \mathbf{I} ~\mathrm{d} \xi - \frac{1}{(2\pi)^3} \frac{\lambda + \mu}{\mu (\lambda + 2\mu)} \int_{\mathbb{R}^3} \frac{1}{|\xi|^4} \xi \xi^{\top} \mathrm{d} \xi \] 将积分分成两部分进行计算 \[ I_0(x) = \int_{\mathbb{R}^3} \frac{1}{|\xi|^2} e^{i\xi \cdot x} ~\mathrm{d} \xi, \quad I_{ij}(x) = - \frac{\partial^2}{\partial x_i \partial x_j} \int_{\mathbb{R}^3} \frac{1}{|\xi|^4} e^{i\xi \cdot x} \mathrm{d} \xi \] 一波操作得到 \[ I_0(x) = 2\pi^2 \frac{1}{|x|}, \quad I_{ij} = \pi^2 \frac{\partial^2 |x|}{\partial x_i \partial x_j} \] 代入方程中得到 \[ G_{ij}(x) = \frac{1}{4\pi \mu} \frac{1}{|x|} \delta_{ij} - \frac{\lambda + \mu}{8\pi \mu (\lambda + 2\mu)} \frac{\partial^2}{\partial x_i \partial x_j} |x| \] > 对于径向对称幂函数 \(1/|\xi|^a\), Fourier 逆变换依旧得到径向对称幂函数 \(C_a/|x|^{d-a}\), 这时需要确定常数 \(C_a\), 当 \(a < d\) 时, 有 > \[ > C_a=(2 \pi)^{\frac{n}{2}} \frac{2^{\frac{n-a}{2}} \Gamma\left(\frac{n-a}{2}\right)}{2^{\frac{a}{2}} \Gamma\left(\frac{a}{2}\right)} . > \] > Radial functions and the Fourier transform > > 对于 \(a > d\) 的情况, 积分是不收敛的, 但注意到 \(I_{ij}\) 还是关于 \(x_i\) 的二阶导数, 这就抵消了 2 幂次, 积分是收敛的. 尝试计算 \(I_{11}(x)\), 并令 \(x=(0,0,1)\), 最终得到 > \[ > C_{4} = 4\pi \int_{0}^{\infty} \frac{1}{r^3} \sin r - \frac{1}{r^2} \cos r ~\mathrm{d}r > \] > 有意思的是, 如果直接通过计算积分 \(\int_{\mathbb{R}^3} \frac{1}{|\xi|^4} e^{i\xi \cdot x} \mathrm{d} \xi\) 求解常数 \(C_4\), 最终会得到上式积分项中关于 \(\sin r\) 部分 > \[ > C_4 = 4\pi \int_{0}^{\infty} \frac{1}{r^3} \sin r ~\mathrm{d}r > \] > 而这个积分是不收敛的