平面地震波斜め入射解析の定常応答解析の定式化

地表面が自由境界の一様な半無限地盤にSV波の定常波が入射した際の変位と応力を計算します。式展開は以下の論文を参考にします:point_down:
[地震動のSV波動特性の研究](https://www.jstage.jst.go.jp/article/jscej1969/1979/289/1979_289_1/_pdf)(1979 金子)
![SV波入射.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44756/5f650b5a-7c37-8722-b316-58c73deb08bf.png)
SV波が入射した際にはSV波とP波が反射されることが知られています。P波の変位ポテンシャル$\phi$とS波の変位ポテンシャル$\psi$を用いると変位は次式のように表されます:writing_hand:

```math
\begin{align}
u&=\dfrac{\partial\phi}{\partial x}+\dfrac{\partial\psi}{\partial z}\\
w&=\dfrac{\partial\phi}{\partial z}-\dfrac{\partial\psi}{\partial x}\\
\end{align}
```



```math
\begin{align}
\psi_{0}=I_{S}\exp\left(i\left(fx+sz-pt\right)\right) \\
\phi_{1}=A_{1}\exp\left(i\left(fx-rz-pt\right)\right) \\
\psi_{2}=B_{2}\exp\left(i\left(fx-sz-pt\right)\right) \\
\end{align}
```
ここで、$I_S$、$A_{1}$、$B_{2}$はそれぞれ、変位ポテンシャルの振幅です。また、座標のZ軸の方向が論文は下向きになっているのに対し、本計算では上向きとしているため、zに関する符号が逆向きになっていることに注意してください。これは、物理的には上昇波と下降波の式が逆になることを意味します:arrow_up_down:ここで、$f$と$s$は次の式で表されます。

```math
\begin{align}
f&=\left(\dfrac{2\pi}{L_0}\right)\sin\theta_0 \\
s&=\left(\dfrac{2\pi}{L_0}\right)\cos\theta_0 \\
L_0&=V_ST_0
\end{align}
```

$V_S$は地盤のS波速度、$T_0$は定常波の波長を示します。この変位ポテンシャルを変位の式に代入します。

```math
\begin{align}
u&=\left.\dfrac{\partial\phi}{\partial x}+\dfrac{\partial\psi}{\partial z}\right|_{\begin{array}{c}
\phi=\phi_{1}\\
\psi=\psi_{0}+\psi_{2}\\
\\
\end{array}}\\
&=\dfrac{\partial\phi_1}{\partial x}+\dfrac{\partial \psi_0}{\partial z}+\dfrac{\partial \psi_2}{\partial z}\\
&=if\phi_1+is\psi_0-is\psi_2
\end{align}
```

```math
\begin{align}
w&=\left.\dfrac{\partial\phi}{\partial z}-\dfrac{\partial\psi}{\partial x}\right|_{\begin{array}{c}
\phi=\phi_{1}\\
\psi=\psi_{0}+\psi_{2}
\end{array}} \\
&=\dfrac{\partial\phi_1}{\partial z}-\dfrac{\partial \psi_0}{\partial x}-\dfrac{\partial \psi_2}{\partial x} \\
&=-ir\phi_1-if\psi_0-if\psi_2
\end{align}
```

また、変位を空間方向に対して微分すると応力を計算することができます。

```math
\begin{align}
\sigma_z&= \lambda\dfrac{\partial u}{\partial x}+\left(\lambda+2\mu\right) \dfrac{\partial w}{\partial z}\\
&=\lambda \left(-f^2\phi_1-fs\psi_0+fs\psi_2\right)
+\left(\lambda+2\mu\right)\left(-r^2\phi_1+fs\psi_0-fs\psi_2\right)\\
&=2\mu fs\psi_0+\left(-\lambda f^2-\lambda r^2-2\mu r^2 \right)\phi_1+2\mu fs\psi_2\\
\tau_{zx}&=\mu \left(\dfrac{\partial u}{\partial z}+\dfrac{\partial w}{\partial x}\right)\\
&=\mu \left( fr\phi_1-s^2\psi_0-s^2\psi_2+fr\phi_1 +fs\psi_0-fs\psi_2\right)\\
&=\left(-\mu s^2 + \mu fs\right)\psi_0 + 2\mu fr\phi_1 -\mu \left(s^2+fs\right) \psi_2
\end{align}
```

ただし、$\lambda$と$\mu$は地盤の物性値を示します。Lameの弾性定数です。**以上で、本問題における変位と応力を計算することが可能とな・・・っていません**:frowning2:未知の変数があります。変位ポテンシャルで出てきた振幅$A_1$、$B_2$です。この変数を求めることができれば、本問題における変位と応力を計算することができます。$I_S$もありますが、ここでは入射波の振幅は既知のものとしそれにより、$A_1$、$B_2$がどのように計算できるかという問題とします。

本問題における境界条件は地表面が自由表面であることのみですので、これに絡めて解きましょう:bulb:自由表面($z=L_Z$)では応力は$0$になります。

```math
\begin{align}
\left.\sigma_z\right|_{z=L_Z}=0 \\
\left.\tau_{zx}\right|_{z=L_Z}=0
\end{align}
```

先程展開した応力の式を代入すると次の式が得られます:relieved:

```math
\begin{align}
\left(-\lambda f^2-\lambda r^2-2\mu r^2 \right)\phi_1+2\mu fs\psi_2 &= -2\mu fs\psi_0 \\
2\mu fr\phi_1 -\left(s^2+fs\right) \psi_2 &= \left(\mu s^2 - \mu fs\right)\psi_0
\end{align}
```

これを連立方程式とすると次式になります。

```math
\begin{align}
\begin{bmatrix}
-\lambda f^2-\lambda r^2-2\mu r^2 & 2\mu fs \\
2\mu fr & \mu\left(-s^2-fs\right) \\
\end{bmatrix}
\left\{
\begin{array}{r} \phi_1 \\ \psi_2\end{array}
\right\}
= \left\{
\begin{array}{r}
-2\mu fs \psi_0 \\ \left(\mu s^2 - \mu fs\right) \psi_0
\end{array}
\right\}
\end{align}
```

複素数は扱いにくいため、実部のみの比較を考えます:bulb:

```math
\begin{align}
\begin{bmatrix}
\left(-\lambda f^2-\lambda r^2-2\mu r^2\right)\cos\left(-rL_Z\right) & \left(2\mu fs\right)\cos\left(-rL_Z\right) \\
\left(2\mu fr\right) \cos\left(-rL_Z\right) & \left(-\mu s^2-\mu fs\right) \cos\left(-rL_Z\right) \\
\end{bmatrix}
\left\{
\begin{array}{r}
A_{1} \\
B_{2}
\end{array}
\right\}
= \left\{
\begin{array}{r}
-2\mu fs \cos\left(sL_Z\right) I_S \\
\left(\mu s^2 - \mu fs\right) \cos\left(sL_Z\right) I_S
\end{array}
\right\}
\end{align}
```

この式を解けば、SV波の斜め入射時の$A_1$と$B_2$を計算することができ。変位と応力の理論解を得ることができます:robot: