Masassiah Blog

Masassiah のブログです。主に読書で得た気づきをまとめています。

高次精度移流方程式解法 CIP 法

2021年3月7日更新

はじめに

流体力学分野において,東工大の矢部らによって提案された高次精度移流方程式解法である Constrained Interpolation Profile Method(CIP 法)を試してみた。
まずは,一次元の移流方程式に適用してみます。
方向分離解法を使えば,多次元化することもできます(M 型 CIP 法や C 型 CIP 法と呼ばれている)(M型CIP法やC型CIP法)。
今後,電磁波や音波の一次元伝搬解析に CIP 法を適用させ,さらには,多次元化にも挑戦してもたいと思います。

CIP 法の特徴

CIP 法は,関数の値とその空間微分値を用いて,空間格子間を補間し,更新する手法である。微分値を用いることで,波形を崩さずに更新できるという特長がある。具体的には,関数値と微分値が,同じ速度で移流することを利用して,格子点間で 3 次関数補間(Hermite 補間)曲線を作成し,値を移流させる。

CIP法とJavaによるCGシミュレーション

CIP法とJavaによるCGシミュレーション

 

 

MATLAB ブログラム

下図に示すガウシアンパルスが速度 1 で伝搬していく様子を CIP 法によりシミュレーションする。

ガウシアンパルス

ガウシアンパルス
%===Parameter===%
NX=101;NT=100;
c=1;
dx=1;
CFL=0.5;
dt=CFL*dx/c;
x=0:dx:dx*(NX-1);
N=2;
sigma=5;
x0=50;
%===Memory Occupation===%
F=zeros(N,NX);
G=zeros(N,NX);
FN=zeros(N,NX);
GN=zeros(N,NX);
%===Initial Condition===%
F(1,:)=0.5*exp(-(x-x0).^2/2/sigma^2);
F(2,:)=0.5*exp(-(x-x0).^2/2/sigma^2);
G(1,:)=-(x-x0).*F(1,:)/sigma^2;
G(2,:)=-(x-x0).*F(2,:)/sigma^2;
%===Time Step===%
disp('Star Time Step !') 
for T=0:NT
    %===figure===%
    plot(x,F(1,:)+F(2,:),'b-');
    axis([0,max(x),-0.1,1.1]);
    xlabel('x');ylabel('F');
    grid on; drawnow;
    %===Step===%
    for n=1:N
        if n==1;
            D=-dx;xx=-c*dt;iup0=-1;
            for m=2:NX iup=m+iup0;
                a=(G(n,m)+G(n,iup))./D^2+2*(F(n,m)-F(n,iup))./D^3;
                b=3*(F(n,iup)-F(n,m))./D^2-(2*G(n,m)+G(n,iup))./D;
                FN(n,m)=a*xx^3+b*xx^2+G(n,m)*xx+F(n,m);
                GN(n,m)=3*a*xx^2+2*b*xx+G(n,m);
            end
        end
        if n==2;
            D= dx;xx= c*dt;iup0= 1;
            for m=NX-1:-1:1 iup=m+iup0;
                a=(G(n,m)+G(n,iup))./D^2+2*(F(n,m)-F(n,iup))./D^3;
                b=3*(F(n,iup)-F(n,m))./D^2-(2*G(n,m)+G(n,iup))./D;
                FN(n,m)=a*xx^3+b*xx^2+G(n,m)*xx+F(n,m);
                GN(n,m)=3*a*xx^2+2*b*xx+G(n,m);
            end
        end
    end
    %===Reset===%
    F=FN;G=GN;
end