自由边界和刚性边界的反射系数

以 P-SV 为例

Author: Zhu Dengda
Email: zhudengda@mail.iggcas.ac.cn

[1]:

import sympy as sp

动态解

[2]:

# 定义基本变量
k, b, d, a, mu, Omg, kb, ka, w, Delta = sp.symbols(r'k b d, a mu \Omega k_{b} k_{a} \omega \Delta')
lam, rho = sp.symbols(r'\lambda \rho')

[3]:

D = sp.Matrix([
    [k,            b,             k,            -b],
    [a,            k,             -a,           k],
    [2*mu*Omg,     2*k*mu*b,      2*mu*Omg,     -2*k*mu*b],
    [2*k*mu*a,     2*mu*Omg,      -2*k*mu*a,    2*mu*Omg]
])

[4]:

def get_R(D, at_bottom:bool):
    if not at_bottom:
        return - (D[:, 2:])**(-1) * (D[:, :2])
    else:
        return - (D[:, :2])**(-1) * (D[:, 2:])

自由界面（z平面牵引力为0）

[5]:

# z+ 侧， 例如顶层顶界面
R = get_R(D[2:, :], False)
R.simplify()
R

[5]:

$\displaystyle \left[\begin{matrix}\frac{\Omega^{2} + a b k^{2}}{- \Omega^{2} + a b k^{2}} & \frac{2 \Omega b k}{- \Omega^{2} + a b k^{2}}\\\frac{2 \Omega a k}{- \Omega^{2} + a b k^{2}} & \frac{\Omega^{2} + a b k^{2}}{- \Omega^{2} + a b k^{2}}\end{matrix}\right]$

[6]:

# z- 侧， 例如底层底界面
R = get_R(D[2:, :], True)
R.simplify()
R

[6]:

$\displaystyle \left[\begin{matrix}\frac{\Omega^{2} + a b k^{2}}{- \Omega^{2} + a b k^{2}} & \frac{2 \Omega b k}{\Omega^{2} - a b k^{2}}\\\frac{2 \Omega a k}{\Omega^{2} - a b k^{2}} & \frac{\Omega^{2} + a b k^{2}}{- \Omega^{2} + a b k^{2}}\end{matrix}\right]$

刚性界面（位移为0）

[7]:

# z+ 侧
R = get_R(D[:2, :], False)
R.simplify()
R

[7]:

$\displaystyle \left[\begin{matrix}\frac{a b + k^{2}}{a b - k^{2}} & \frac{2 b k}{a b - k^{2}}\\\frac{2 a k}{a b - k^{2}} & \frac{a b + k^{2}}{a b - k^{2}}\end{matrix}\right]$

[8]:

# z- 侧
R = get_R(D[:2, :], True)
R.simplify()
R

[8]:

$\displaystyle \left[\begin{matrix}\frac{a b + k^{2}}{a b - k^{2}} & - \frac{2 b k}{a b - k^{2}}\\- \frac{2 a k}{a b - k^{2}} & \frac{a b + k^{2}}{a b - k^{2}}\end{matrix}\right]$

静态解

[9]:

# 在静态解中，矩阵 D 和相对深度 z-zj 有关, 因此 z+ 侧和 z- 侧需使用两种 D 矩阵
D1 = sp.Matrix([
    [1,         -(1 + 2*k*Delta*d),            1,            -(1 - 2*k*Delta*d)],
    [1,          (1 - 2*k*Delta*d),           -1,            -(1 + 2*k*Delta*d)],
    [2*mu*k,     2*mu*Delta*k*(1 - 2*k*d),   2*mu*k,      2*mu*Delta*k*(1 + 2*k*d)],
    [2*mu*k,    -2*mu*Delta*k*(1 + 2*k*d),  -2*mu*k,      2*mu*Delta*k*(1 - 2*k*d)],
])

z = 0
D2 = sp.Matrix([
    [1,         -(1 + 2*k*Delta*z),            1,            -(1 - 2*k*Delta*z)],
    [1,          (1 - 2*k*Delta*z),           -1,            -(1 + 2*k*Delta*z)],
    [2*mu*k,     2*mu*Delta*k*(1 - 2*k*z),   2*mu*k,      2*mu*Delta*k*(1 + 2*k*z)],
    [2*mu*k,    -2*mu*Delta*k*(1 + 2*k*z),  -2*mu*k,      2*mu*Delta*k*(1 - 2*k*z)],
])

自由界面（z平面牵引力为0）

[10]:

# z+ 侧
R = get_R(D2[2:, :], False)
R.simplify()
R

[10]:

$\displaystyle \left[\begin{matrix}0 & - \Delta\\- \frac{1}{\Delta} & 0\end{matrix}\right]$

[11]:

# z- 侧
R = get_R(D1[2:, :], True)
R.simplify()
R

[11]:

$\displaystyle \left[\begin{matrix}- 2 d k & \Delta \left(- 4 d^{2} k^{2} - 1\right)\\- \frac{1}{\Delta} & - 2 d k\end{matrix}\right]$

刚性界面（位移为0）

[12]:

# z+ 侧
R = get_R(D2[:2, :], False)
R.simplify()
R

[12]:

$\displaystyle \left[\begin{matrix}0 & 1\\1 & 0\end{matrix}\right]$

[13]:

# z- 侧
R = get_R(D1[:2, :], True)
R.simplify()
R

[13]:

$\displaystyle \left[\begin{matrix}2 \Delta d k & 4 \Delta^{2} d^{2} k^{2} + 1\\1 & 2 \Delta d k\end{matrix}\right]$