今回は一次元定常移流拡散方程式を有限要素法で解いていきます。C++コード付きです。一次元定常移流拡散方程式は、前に書いた記事「有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか?」で紹介した、有限要素法を学ぶ場合に二番目に解くべき方程式です。
一次元定常移流拡散方程式とは
のような二階の常微分方程式のことをいいます。ここで、
前回の一次元Poisson方程式に移流行列を加えるだけなので簡単です。今回も二次元配列を使うコードと一次元配列だけで処理するコードの両方を公開します。境界条件は前回と同じく3パターン用意しました。すなわち、
Case 1, と
でDirichlet条件を課す
( と
の値を指定)
Case 2, でNeumann条件、
でDirichlet条件を課す
( と
の値を指定)
Case 3, でDirichlet条件、
でNeumann条件を課す
( と
の値を指定)
です。それぞれコード上ではbc=1、bc=2、bc=3に対応しています。連立一次方程式の解法はGauss-Seidel法を用いています。
まずは二次元配列を使う方です。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <vector> using namespace std; //inner product// inline double ip(vector<double> &a, vector<double> &b, int &Node) { double value=0.; int i; for(i=0;i<Node;i++) value+=a[i]*b[i]; return value; } //Gauss-Seidel inline void GS(vector<vector<double> > &A,vector<double> &x,vector<double> &b, int &Node, double &eps) { int i,k; double in,sum,max; for(k=0;k<100000;k++) { max=0.0; for(i=0;i<Node;i++) { in=x[i]; x[i]=(b[i]-(ip(A[i],x,Node)-A[i][i]*x[i]))/A[i][i]; if(max<abs(in-x[i])) max=abs(in-x[i]); } if(eps>max) break; } } inline void mat(vector<vector<double> > &A,vector<double> &b, vector<double> &f, int &Ele, double &dx,double &D,double &V) { int i; for(i=0;i<Ele;i++) { //Diffusion// A[i][i]+=D/dx;A[i][i+1]+=-D/dx; A[i+1][i]+=-D/dx;A[i+1][i+1]+=D/dx; //Advection// A[i][i]+=-V/2.0;A[i][i+1]+=V/2.0; A[i+1][i]+=-V/2.0;A[i+1][i+1]+=V/2.0; //Source// b[i]+=(2.0*f[i]+1.0*f[i+1])*dx/6.0; b[i+1]+=(1.0*f[i]+2.0*f[i+1])*dx/6.0; } } inline void boundary(int &bc,double &D0,double &D1,double &N0,double &N1,vector<vector<double> > &A,vector<double> &b,int &Node) { int i; if(bc==1) { for(i=0;i<Node;i++) { A[0][i]=0.0; A[Node-1][i]=0.0; } A[0][0]=1.0;b[0]=D0; A[Node-1][Node-1]=1.0;b[Node-1]=D1; } if(bc==2) { for(i=0;i<Node;i++) { A[Node-1][i]=0.0; } A[Node-1][Node-1]=1.0;b[Node-1]=D1; b[0]+=-N0; } if(bc==3) { for(i=0;i<Node;i++) { A[0][i]=0.0; } A[0][0]=1.0;b[0]=D0; b[Node-1]+=N1; } } int main() { int i; int Ele=100; int Node=Ele+1; double L=1.0; double dx=L/Ele; double D=0.01; double V=0.1; double eps=pow(2.0,-50); int bc=1;//bc=1 both Diriclet bc=2 left Neumann bc=3 right Neumann double D0=0.0; double D1=1.0; double N0=0.0; double N1=2.0; vector<double> b(Node,0.); vector<double> x(Node,0.); vector<double> f(Node,0.); vector<vector<double> > A(Node,(vector<double>(Node,0.))); ofstream fk; fk.open("answer.txt"); mat(A,b,f,Ele,dx,D,V); boundary(bc,D0,D1,N0,N1,A,b,Node); GS(A,x,b,Node,eps); for(i=0;i<Node;i++) cout<<x[i]<<endl; for(i=0;i<Node;i++) fk<<dx*i<<" "<<x[i]<<endl; return 0; }
次に一次元配列だけで処理する方です。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <vector> using namespace std; //Gauss-Seidel inline void GS(vector<double> &D,vector<double> &L,vector<double> &R,vector<double> &x,vector<double> &b, int &Node, double &eps) { int i,k; double in,sum,max; for(k=0;k<100000;k++) { max=0.0; in=x[0]; x[0]=(b[0]-R[0]*x[1])/D[0]; if(max<abs(in-x[0])) max=abs(in-x[0]); for(i=1;i<Node-1;i++) { in=x[i]; x[i]=(b[i]-(L[i]*x[i-1]+R[i]*x[i+1]))/D[i]; if(max<abs(in-x[i])) max=abs(in-x[i]); } in=x[Node-1]; x[Node-1]=(b[Node-1]-L[Node-1]*x[Node-2])/D[Node-1]; if(max<abs(in-x[Node-1])) max=abs(in-x[Node-1]); if(eps>max) break; } } inline void mat(vector<double> &D,vector<double> &L,vector<double> &R,vector<double> &b, vector<double> &f, int &Ele, double &dx, double &V, double &Diff) { int i; for(i=0;i<Ele;i++) { //Diffusion// D[i]+=Diff/dx;R[i]+=-Diff/dx; L[i+1]+=-Diff/dx;D[i+1]+=Diff/dx; //Advection// D[i]+=-V/2.0;R[i]+=V/2.0; L[i+1]+=-V/2.0;D[i+1]+=V/2.0; //Source// b[i]+=(2.0*f[i]+1.0*f[i+1])*dx/6.0; b[i+1]+=(1.0*f[i]+2.0*f[i+1])*dx/6.0; } } inline void boundary(int &bc,double &D0,double &D1,double &N0,double &N1,vector<double> &D,vector<double> &L,vector<double> &R,vector<double> &b,int &Node) { if(bc==1) { D[0]=1.0;R[0]=0.0;b[0]=D0; L[Node-1]=0.0;D[Node-1]=1.0;b[Node-1]=D1; } if(bc==2) { L[Node-1]=0.0;D[Node-1]=1.0;b[Node-1]=D1; b[0]+=-N0; } if(bc==3) { D[0]=1.0;R[0]=0.0;b[0]=D0; b[Node-1]+=N1; } } int main() { int i; int Ele=100; int Node=Ele+1; double LL=1.0; double dx=LL/Ele; double eps=pow(2.0,-50); double Diff=0.01; double V=0.1; int bc=1;//bc=1 both Diriclet bc=2 left Neumann bc=3 right Neumann double D0=0.0; double D1=1.0; double N0=0.0; double N1=2.0; vector<double> b(Node,0.); vector<double> x(Node,0.); vector<double> f(Node,0.); vector<double> D(Node,0.); vector<double> L(Node,0.); vector<double> R(Node,0.); ofstream fk; fk.open("answer.txt"); mat(D,L,R,b,f,Ele,dx,V,Diff); boundary(bc,D0,D1,N0,N1,D,L,R,b,Node); GS(D,L,R,x,b,Node,eps); for(i=0;i<Node;i++) cout<<x[i]<<endl; for(i=0;i<Node;i++) fk<<dx*i<<" "<<x[i]<<endl; return 0; }
こちらの離散化も詳しく書く予定です。Galerkin法では移流が卓越する場合をきちんと扱えないので、移流方向を考慮した離散化「Petrov-Galerkin法」を使う必要があります。それも後程。書きたいこと、書かねばならぬことが山のようにあります。ひとつずつ着実に消化していきたいです。