二次元偏微分方程式数値解 - 重力・空気抵抗あり膜の３Dシミュレーション -

重力・空気抵抗項ありの二次元波動方程式の数値計算＆３D可視化

#弦の振動の運動方程式から波動方程式導出(1次元)

```math
\begin{eqnarray}
ρΔx\frac{∂^2u(x,t)}{∂^2t}&=&(Tsinθ_{x+Δx}-Tsinθ_x)-ρΔxg-kΔx\frac{∂u(x,t)}{∂t}\\
\end{eqnarray}
```
$ρ$ : 線密度
$T$ : 張力
$g$ : 重力加速度
$k$ : 抵抗係数
$θ_{〇}$ : $x=〇$での弦と$x$軸のなす角
右辺第1・2項は両端に働いてる張力の$u$方向成分、第3項は重力、第4項は速度に比例した抵抗
$θ$が小さいときは、$sinθ≒tanθ=\frac{∂u(x,t)}{∂x}$より

```math
\begin{eqnarray}
ρΔx\frac{∂^2u(x,t)}{∂^2t}&≒&T(\frac{∂u(x+Δx,t)}{∂x}-\frac{∂u(x,t)}{∂x})-ρΔxg-kΔx\frac{∂u(x,t)}{∂t}\\
&≒&TΔx\frac{∂^2u(x,t)}{∂^2x}-ρΔxg-kΔx\frac{∂u(x,t)}{∂t}\\
\frac{∂^2u(x,t)}{∂^2t}&≒&\frac{T}{ρ}\frac{∂^2u(x,t)}{∂^2x}-g-\frac{k}{ρ}\frac{∂u(x,t)}{∂t}
\end{eqnarray}
```

#2次元

```math
\begin{eqnarray}
\frac{∂^2u(x,y,t)}{∂^2t}&≒&\frac{T}{ρ}(\frac{∂^2u(x,y,t)}{∂^2x}+\frac{∂^2u(x,y,t)}{∂^2y})-g-\frac{k}{ρ}\frac{∂u(x,y,t)}{∂t}
\end{eqnarray}
```

#離散化

```math
\begin{eqnarray}
\frac{1}{Δt^2}(u(t+Δt,x,y)+u(t-Δt,x,y)-2u(t,x,y))&=&\frac{1}{Δ^2}(u(t,x-Δ,y)+u(t,x+Δ,y)+u(t,x,y-Δ)+u(t,x,y+Δ)-4u(t,x,y))-g-\frac{k}{ρ}\frac{(u(t,x,y)-u(t-Δt,x,y))}{Δt}\\
u(t+Δt,x,y)&=&2u(t,x,y)-u(t-Δt,x,y)+\frac{Δt^2}{Δ^2}(u(t,x-Δ,y)+u(t,x+Δ,y)+u(t,x,y-Δ)+u(t,x,y+Δ)-4u(t,x,y))-Δt^2g-\frac{kΔt}{ρ}(u(t,x,y)-u(t-Δt,x,y))…①\\
(Δ=Δx=Δy)
\end{eqnarray}
```

#今回の条件
・正方形の膜
・端は固定
・いくつかの場所に上下から振動しながら突き抜ける棒を設置

#可視化
![membrane3.PNG](https://qiita-image-store.s3.amazonaws.com/0/38382/9cfca687-faad-258e-2336-201638cea424.png)
![membrane2.PNG](https://qiita-image-store.s3.amazonaws.com/0/38382/ee2c3458-c1e4-f914-134e-1383b863bb7e.png)

・動画
<blockquote class="twitter-tweet" data-lang="ja"><p lang="und" dir="ltr"><a href="https://t.co/9NwPOwwIsX">https://t.co/9NwPOwwIsX</a></p>&mdash; ぽんでりおん (@popondeli) <a href="https://twitter.com/popondeli/status/878336566027878400">2017年6月23日</a></blockquote>
<script async src="//platform.twitter.com/widgets.js" charset="utf-8"></script>

#アルゴリズム

```cpp
…
    //uの初期化
…
    //u(t,x,y)の更新
    for(int i = 1; i < MESH_NX-1; ++i) //端以外
    {
        for(int j = 1; j < MESH_NY-1; ++j) //端以外
        {
            u[t_index][i][j] = calcHeight(i, j, dt, mtime, t_index);
        }
    }
    //棒が突き抜けている位置は固定
    if (u[t_index][25][25] < 6.0f*sin(mtime))
        u[t_index][25][25] = 6.0f*sin(mtime);
    if (u[t_index][75][75] < -8.0f*cos(0.7f*mtime))
        u[t_index][75][75] = -8.0f*cos(0.7f*mtime);
    if (u[t_index][55][45] < -8.0f*cos(0.2f*mtime))
        u[t_index][55][45] = -8.0f*cos(0.2f*mtime);
    if (u[t_index][48][52] > 10.0f*cos(0.6f*mtime))
        u[t_index][48][52] = 10.0f*cos(0.6f*mtime);
        
    t_index++;
    if (t_index > 2)
        t_index = 0;
…

float calcHeight(int i, int j, float dt, float time, int t_index)
{
	float height;
	// height[t][x][y]、時間は2つ前(2Δt)のデータまで使用するため、t_indexを回して使いまわす。
	switch(t_index)
	{
	case 0:
		height = 2.0f*u[2][i][j] - u[1][i][j] + (TX*(u[2][i+1][j]-2.0f*u[2][i][j]+u[2][i-1][j]) + TY*(u[2][i][j+1]-2.0f*u[2][i][j]+u[2][i][j-1]))*dt*dt*0.2f + (FPerMass(i,j,time)*dt - K*(u[2][i][j]-u[1][i][j]))*dt;//①の計算
		break;
	case 1:
		height = 2.0f*u[0][i][j] - u[2][i][j] + (TX*(u[0][i+1][j]-2.0f*u[0][i][j]+u[0][i-1][j]) + TY*(u[0][i][j+1]-2.0f*u[0][i][j]+u[0][i][j-1]))*dt*dt*0.2f + (FPerMass(i,j,time)*dt - K*(u[0][i][j]-u[2][i][j]))*dt;//①の計算
		break;
	case 2:
		height = 2.0f*u[1][i][j] - u[0][i][j] + (TX*(u[1][i+1][j]-2.0f*u[1][i][j]+u[1][i-1][j]) + TY*(u[1][i][j+1]-2.0f*u[1][i][j]+u[1][i][j-1]))*dt*dt*0.2f + (FPerMass(i,j,time)*dt - K*(u[1][i][j]-u[0][i][j]))*dt;//①の計算
		break;
	}
	return height;
}

float FPerMass(int i, int j, float t)
{
	return -0.02f;//重力
}
```