古典的な格子生成方法：理論と実装 - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

この記事のポイント

ラプラス方程式を解くことによって滑らかな格子を生成することができます。
座標変換のラプラス方程式を導出しています。
pythonによる数値計算プログラムを載せています。

はじめに
計算式の整理
実装
サンプル
- 参考
参考：SOR v.s. Gauss-Seidel
- 参考
まとめ
GitHub

はじめに

河川流解析でよく用いる境界適合格子の平面二次元計算の格子生成について整理しました。

現在は私はこの方法を使わないのですが、昔のメモが出てきたのでまとめておきます。

計算式の整理

座標変換の準備

$xy$ 座標系から $\xi \eta$ 座標系に変換するため、そのための式を整理しておきます。

$\begin{align} \dfrac{\partial \xi}{\partial x} = \xi_x, \dfrac{\partial \xi}{\partial y} = \xi_y, \dfrac{\partial \eta}{\partial x} = \eta_x, \dfrac{\partial \eta}{\partial y} = \eta_y \\ \dfrac{\partial x}{ \partial \xi} = x_\xi, \dfrac{\partial y}{ \partial \xi} = y_\xi, \dfrac{\partial x}{\partial \eta} = x_\eta, \dfrac{\partial y}{\partial \eta} = y_\eta \end{align}$

$\begin{align} J = \dfrac{1}{x_\xi y_\eta - x_\eta y_\xi} = \eta_y \xi_x - \xi_y \eta_x \end{align}$

$\begin{align} \xi_x &= J \cdot y_\eta \\ \xi_y &= -J \cdot x_\eta \\ \eta_x &= -J \cdot y_\xi \\ \eta_y &= J \cdot x_\xi \end{align}$

$\begin{align} \dfrac{\partial }{\partial x} &= \dfrac{\partial \xi}{\partial x}\dfrac{\partial }{\partial \xi} + \dfrac{\partial \eta}{\partial x}\dfrac{\partial }{\partial \eta} = J\left(y_\eta\dfrac{\partial }{\partial \xi} - y_\xi\dfrac{\partial }{\partial \eta} \right) \\ \dfrac{\partial }{\partial y} &= \dfrac{\partial \xi}{\partial y}\dfrac{\partial }{\partial \xi} + \dfrac{\partial \eta}{\partial y}\dfrac{\partial }{\partial \eta} = J\left( - x_\eta\dfrac{\partial }{\partial \xi} + x_\xi \dfrac{\partial }{\partial \eta} \right) \end{align}$

2階微分についても整理します。後の処理のため、式変形を行っているが、かなり煩雑である。(間違っていたらすいません。)

$\begin{align} \dfrac{\partial^2 }{\partial x^2} &= J^2\left( y_\eta^2 \dfrac{\partial^2 }{\partial \xi^2} - y_\eta y_\xi \dfrac{\partial^2 }{\partial \xi \eta} + y_\xi^2 \dfrac{\partial^2 }{\partial \eta^2} \right) \\ &+J^3 ( y_\eta^2 y_{\xi\xi} -2 y_\xi y_\eta y_{\xi\eta} + y_\xi^2 y_{\eta\eta} ) \left(x_\eta \dfrac{\partial }{\partial \xi} - x_\xi \dfrac{\partial }{\partial \eta} \right) \\ &+J^3 \left( y_\eta^2 x_{\xi\xi} -2 y_\xi y_\eta x_{\xi\eta} + y_\xi^2 x_{\eta\eta} \right) \left(y_\xi \dfrac{\partial }{\partial \eta} - y_\eta \dfrac{\partial }{\partial \xi} \right) \\ % \dfrac{\partial^2 }{\partial y^2} &= J^2\left( x_\eta^2 \dfrac{\partial^2 }{\partial \xi^2} - x_\eta x_\xi \dfrac{\partial^2 }{\partial \xi \eta} + x_\xi^2 \dfrac{\partial^2 }{\partial \eta^2} \right) \\ &+J^3 ( x_\eta^2 y_{\xi\xi} -2 x_\xi x_\eta y_{\xi\eta} + x_\xi^2 y_{\eta\eta} ) \left(x_\eta \dfrac{\partial }{\partial \xi} - x_\xi \dfrac{\partial }{\partial \eta} \right) \\ &+J^3 ( x_\eta^2 x_{\xi\xi} -2 x_\xi x_\eta x_{\xi\eta} + x_\xi^2 x_{\eta\eta} ) \left(y_\xi \dfrac{\partial }{\partial \eta} - y_\eta \dfrac{\partial }{\partial \xi} \right) \\ \end{align}$

2式よりラプラシアン $\nabla \cdot \nabla$ は次のとおりになります。

$\begin{align} \nabla \cdot \nabla &= \dfrac{\partial^2 }{\partial x^2} + \dfrac{\partial^2 }{\partial y^2} \\ &=J^2\left( \alpha \dfrac{\partial^2 }{\partial \xi^2} - 2\beta \dfrac{\partial^2 }{\partial \xi \eta} + \gamma \dfrac{\partial^2 }{\partial \eta^2} \right) \\ &+J^3 ( \alpha x_{\xi\xi} -2 \beta x_{\xi\eta} + \gamma x_{\eta\eta} ) \left(y_\xi \dfrac{\partial }{\partial \eta} - y_\eta \dfrac{\partial }{\partial \xi} \right) \\ &+J^3 ( \alpha y_{\xi\xi} -2 \beta y_{\xi\eta} + \gamma y_{\eta\eta} ) \left(x_\eta \dfrac{\partial }{\partial \xi} - x_\xi \dfrac{\partial }{\partial \eta} \right) \\ \end{align}$

ここに、

$\begin{align} \alpha &= x_\eta^2 + y_\eta^2 \\ \beta &= x_\eta x_\xi + y_\eta y_\xi \\ \gamma &= x_\xi^2 + y_\xi^2 \end{align}$

ラプラス方程式による座標変換

座標変換では、ラプラス方程式（またはポアソン方程式）を解くことによって滑らかな格子を生成することができます。

$\begin{align} \dfrac{\partial^2 \xi}{\partial x^2} + \dfrac{\partial^2 \xi}{\partial y^2} &= 0 \\ \dfrac{\partial^2 \eta}{\partial x^2} + \dfrac{\partial^2 \eta}{\partial y^2} &= 0 \end{align}$

これら2式について、前述の座標変換を用いて、 $\xi,\eta$ による微分の形に変形します。2式を連立して式変形をすると下式が導出されます。

$\begin{align} \alpha x_{\xi\xi} -2 \beta x_{\xi\eta} + \gamma x_{\eta\eta} &=0 \\ \alpha y_{\xi\xi} -2 \beta y_{\xi\eta} + \gamma y_{\eta\eta} &=0 \end{align}$

この2式を解くことにより、滑らかな計算格子を生成することできます。

詳細は以下の参考書籍をご確認いただければと思います。

参考書籍

リンク

実装

前述のポアソン方程式による座標変換を実装します。数値解法にはSORを使用しています。言語はpythonです。

def relaxation(xi,yi,omega=1.9, nit=300):
    dxi, deta = float(1), float(1)
    x = xi.copy()
    y = yi.copy()
    outerr = []
    for it in range(nit):
        xn = x.copy()
        yn = y.copy()
        err = []
        for i in range(1, x.shape[0]-1):
            for j in range(1, x.shape[1]-1):
                alfa = ((x[i,j+1]-x[i,j-1])/2/deta)**2 + ((y[i,j+1]-y[i,j-1])/2/deta)**2
                beta = (x[i+1,j]-x[i-1,j])/2/dxi*(x[i,j+1]-x[i,j-1])/2/deta + (y[i+1,j]-y[i-1,j])/2./dxi*(y[i,j+1]-y[i,j-1])/2/deta
                gamma = ((x[i+1,j]-x[i-1,j])/2/dxi)**2 + ((y[i+1,j]-y[i-1,j])/2/dxi)**2.
                
                # SOR
                dd = (alfa*(x[i+1,j]+x[i-1,j])/dxi**2 \
                          - 2*beta*(x[i+1,j+1]-x[i-1,j+1]-(x[i+1,j-1]-x[i-1,j-1]))/4/dxi/deta \
                          + gamma*(x[i,j+1]+x[i,j-1])/deta**2)/(2*alfa/dxi**2+2*gamma/deta**2) 
                
                err.append( dd - x[i,j] )
                x[i,j] += omega*(dd - x[i,j])
                
                dd = (alfa*(y[i+1,j]+y[i-1,j])/dxi**2 \
                          - 2*beta*(y[i+1,j+1]-y[i-1,j+1]-(y[i+1,j-1]-y[i-1,j-1]))/4/dxi/deta \
                          + gamma*(y[i,j+1]+y[i,j-1])/deta**2)/(2*alfa/dxi**2+2*gamma/deta**2)
                
                err.append( dd - y[i,j] )
                y[i,j] += omega*(dd - y[i,j])
        
        serr = np.sqrt( np.sum(np.array(err)**2) )
        print(serr)
        outerr.append(serr)
        
    return x,y,outerr