2017年02月21日

Javaを使った有限体積法による数値計算(2)

第2回目では、さっそくソースをアップする。

今回は、有限体積法の教科書としては、詳しい解説と、例題の丁寧な
回答が掲載されている
数値流体力学 第2版
の例題6.2をコード化した。

定常の移流方程式を解いている。
この例題を解いてみたのは、これから解きたい浅水流方程式に
近い形であること。つまり、速度場と密度場が未知であり、
これらを求めるためにSIMPLEアルゴリズムを適用するものであるためである。
また、断面形状も一定でなく、変化しているところも
今後の応用性が高く、参考になると思う。

Javaについては、1か月程度の知識しかなく、
静的メソッドのみの使用と、mainメソッドにほとんどの必要な処理を書いた。
初めて作ったJavaコードなので、至らないところが多い。

勉強しながら、現在はクラス等を用いたコードに修正を加えている。
また、多次元配列の初期化については、
Arrays#fillメソッドの多次元版の一例
を利用させていただいた。

連立方程式の解法には、トーマス法を用いている。
格子には等間隔のスタッガード格子を採用している。
また、例題の解答を直接コーディング化ではなく、最後の方に触れられている
不足緩和も係数計算の際に、組み込んでいる。

格子の数や不足緩和の係数を、いろいろ試してみたり、
質量フラックスを計算したりしてみて、有限体積法の面白さを
感じることができると思う。

(注)コードの利用による損害等は一切責任を負いません。

/*
* 有限体積法の計算
* 数値流体力学(第2版)
*
* 2017/2/11
* KZ
*
*/

package FVM;

import java.util.Arrays;

public class FvmTest007 {

public static void main(String args[]){

int nn = 5; // 圧力ポイントの数
int it_max = 200; // 最大回数
int count = 0; // 計算回数の記憶用
double r_limit = 10e-5; // 残差の制限


double m_0 = 1.0; // 初期仮定の質量フラックス
double p_0 = 10.0; // 上流端の圧力
double rho = 1.0; // 密度
double alpha = 0.7; // 不足緩和の係数 (nn,alpha)=(10,0.7),(20,0.5)
double residual = 0.0; // 残差

//double[] A_p = {0.5 ,0.4 ,0.3 ,0.2 ,0.1}; //圧力ポイントの面積
double[] A_p = new double[nn];
for(int i=0; i < nn; i++){
A_p[i] = 0.5 - (0.5-0.1)/(nn-1) * i;
}
//double[] A_u = {0.45,0.35,0.25,0.15};
double[] A_u = new double[nn-1];
for(int i=0; i < nn-1; i++){
A_u[i] = ( A_p[i] + A_p[i+1] ) / 2.0;
}

double[] x_u = new double[nn-1]; // 速度場
double[] u_star = new double[nn-1]; // 補正後の速度
double[] u_old = new double[nn-1]; // 速度の古いもの保存(残差計算用)
//double[] x_p = new double[nn]; // 圧力場
//double[] x_p = {10.0, 7.5, 5.0, 2.5, 0.0}; // 圧力場と仮定初期値
double[] x_p = new double[nn];
for(int i = 0; i < nn; i++){
x_p[i] = 10.0 - (10.0 - 0.0)/(nn-1) * i;
}
double[] p_prime = new double[nn]; // 圧力の補正値
double[] p_star = new double[nn]; // 補正後の圧力
double[] d = new double[nn-1]; // 圧力補正のパラメータ

double u_in = 2.0; // 上流端の流速(境界条件)

//方程式の係数
double F_w = 0.0;
double F_e = 0.0;
double a_W = 0.0;
double a_E = 0.0;
double a_P = 0.0;
double S_u = 0.0;

double b_prime = 0.0;

//連立方程式の係数行列
double[][] matrixA = new double[nn-1][nn-1];
double[] matrixB = new double[nn-1];


double[][] matrixC = new double[nn][nn];
double[] matrixD = new double[nn];


//配列の初期値
Arrays.fill(x_u,0.0);
Arrays.fill(u_star, 0.0);
Arrays.fill(p_prime, 0.0);
Arrays.fill(p_star, 0.0);
Arrays.fill(d, 0.0);
ArraysfillEx.fill(matrixA, 0.0); //多次元配列の初期値
Arrays.fill(matrixB, 0.0);
ArraysfillEx.fill(matrixC, 0.0); //多次元配列の初期値
Arrays.fill(matrixD, 0.0);



//速度場の初期仮定の設定
for(int i = 0; i < nn-1 ; i++){
x_u[i] = m_0 / ( rho * A_u[i] );
}

System.out.println("x_u : "+Arrays.toString(x_u)); // 行列の標準出力
System.out.println("x_p : "+Arrays.toString(x_p)); // 行列の標準出力
System.out.println("A_u : "+Arrays.toString(A_u));
System.out.println("A_p : "+Arrays.toString(A_p));



//計算の開始
for(int it = 0; it < it_max ; it++){

System.out.println("Count : " + it);

// 速度の更新(残差計算用)
for(int i = 0; i < nn-1 ; i++){
u_old[i] = x_u[i];
}
residual = 0.0; // 残差の初期化


//運動量保存則 始まり
//速度の格子点//
//格子点 1(0)
u_in = x_u[0] * (A_u[0] / A_p[0]); System.out.println("u_in : " + u_in);
F_w = rho * u_in * A_p[0]; //System.out.println(F_w);
F_e = rho * (x_u[0] + x_u[1])/2 * A_p[1]; //System.out.println(F_e);

a_W = 0.0;
a_E = 0.0;
a_P = F_e + F_w * 0.5 * Math.pow((A_u[0]/A_p[0]),2.);


S_u = (x_p[0]-x_p[1])*A_u[0] + F_w*(A_u[0]/A_p[0])*x_u[0] + (1-alpha)*a_P/alpha * x_u[0];
a_P = a_P/alpha; //不足緩和

matrixA[0][0] = a_P;
matrixA[0][1] = -1.0 * a_E;

matrixB[0] = S_u;

d[0] = A_u[0]/ a_P * alpha;//不足緩和


//格子点 2,3(1 〜 nn-3)
for(int i=1 ; i < nn-2;i++){
F_w = rho * (x_u[i-1] + x_u[i])/2 * A_p[i];
F_e = rho * (x_u[i] + x_u[i+1])/2 * A_p[i+1];
a_W = Max(F_w,0.0);
a_E = Max(0.0, -1.0*F_e);
a_P = a_W + a_E + (F_e - F_w);

S_u = (x_p[i]-x_p[i+1])*A_u[i] + (1-alpha)*a_P/alpha * x_u[i];//不足緩和
a_P = a_P/alpha; //不足緩和

matrixA[i][i-1] = -1.0 * a_W;
matrixA[i][i] = a_P;
matrixA[i][i+1] = -1.0 * a_E;

matrixB[i] = S_u;

d[i] = A_u[i]/ a_P * alpha;
}


//格子点 4(nn-2)
F_w = rho * (x_u[nn-3] + x_u[nn-2])/2 * A_p[nn-2]; //System.out.println(F_w);
//F_e = rho * (x_u[nn-2] + x_u[1])/2 * A_p[1]; //System.out.println(F_e);
if(it == 0) {
F_e = m_0; // 最初の計算だけ 要検討!!!
}else{
F_e = rho * x_u[nn-2] * A_u[nn-2];
}

a_W = F_w;
a_E = 0.0;
a_P = a_W + a_E + (F_e - F_w);

S_u = (x_p[nn-2]- x_p[nn-1])*A_u[nn-2] + (1-alpha)*a_P/alpha * x_u[nn-2];
a_P = a_P/alpha; //不足緩和

matrixA[nn-2][nn-3] = -1.0 * a_W;
matrixA[nn-2][nn-2] = a_P;

matrixB[nn-2] = S_u;

d[nn-2] = A_u[nn-2]/ a_P * alpha;



//System.out.println("matA: "+Arrays.deepToString(matrixA));
//System.out.println("matB: "+Arrays.toString(matrixB));
//System.out.println(" d : "+Arrays.toString(d));

//連立方程式を解く

u_star = thomas(nn-1,matrixA,matrixB);

System.out.println("u_star: " + Arrays.toString(u_star));

//運動量保存則 終わり

//圧力の補正式 始まり
//圧力の格子点 1,A (0)
matrixC[0][0] = -1.0;

//圧力の格子点 2から4,A (1 - nn-2)
for(int i = 1; i < nn-1 ; i++){
a_W = rho * d[i-1] * A_u[i-1];
a_E = rho * d[i] * A_u[i];
F_w = rho * u_star[i-1] * A_u[i-1];
F_e = rho * u_star[i] * A_u[i];
a_P = a_W + a_E;
b_prime = F_w - F_e;

matrixC[i][i-1] = -1.0 * a_W;
matrixC[i][i] = a_P;
matrixC[i][i+1] = -1.0 * a_E;

matrixD[i] = b_prime;
}

//圧力の格子点 5,E (nn-1)
matrixC[nn-1][nn-1] = -1.0;


//System.out.println("matC: "+Arrays.deepToString(matrixC));
//System.out.println("matD: "+Arrays.toString(matrixD));

//連立方程式を解く
p_prime = thomas(nn,matrixC,matrixD);

System.out.println("p_prime: " + Arrays.toString(p_prime));

//圧力の補正式 終わり

//圧力の補正 始り
for(int i=0;i < nn;i++){
x_p[i] = x_p[i] + alpha * p_prime[i];
}
System.out.println("calculated x_p : " + Arrays.toString(x_p));
//圧力の補正 終わり

//速度の補正 始まり
for(int i=0;i < nn-1;i++){
x_u[i] = u_star[i] + d[i] * (p_prime[i] - p_prime[i+1]);
}
System.out.println("calculated x_u : " + Arrays.toString(x_u));
//速度の補正 終わり
//格子点Aの圧力の補正
x_p[0] = p_0 - 0.5 * rho * Math.pow(x_u[0], 2.0) * Math.pow(A_u[0]/A_p[0],2.0);
System.out.println("calculated2 x_p : " + Arrays.toString(x_p));
// 終わり

// 残差の計算
for(int i = 0; i < nn-1 ; i++){
residual += Math.pow(u_old[i]-x_u[i],2.0);
}
System.out.println("residual : " + residual);

if(residual < r_limit) {
count = it;
break;
}

}

System.out.println("¥n Count: " + count);
System.out.println("x_u : " + Arrays.toString(x_u));
System.out.println("x_p : " + Arrays.toString(x_p));
System.out.println("A_u : "+Arrays.toString(A_u));
System.out.println("A_p : "+Arrays.toString(A_p));



}



//以下、メソッド

//大きい数字を返す
public static double Max(double a,double b){
if(a>b){
return a;
} else {
return b;
}
}


//トーマス法
public static double[] thomas(int n,double[][] A,double[] B){

double[] x = new double[n];

double[] g = new double[n];
double[] s = new double[n];

//配列の初期値
Arrays.fill(x,0.0);
Arrays.fill(g, 0.0);
Arrays.fill(s, 0.0);

g[0] = A[0][0];
s[0] = B[0];

for(int i = 1 ; i < n ; i++){
s[i] = B[i] - A[i][i-1]*s[i-1]/g[i-1];
g[i] = A[i][i] - A[i][i-1]*A[i-1][i]/g[i-1];
}

x[n-1] = s[n-1]/g[n-1];

for(int i = n-2 ; i >= 0; i--){
x[i] = (s[i]-A[i][i+1]*x[i+1])/g[i];
}

return x;
}
}

  


Posted by 和歌山サイエンスカフェインフィニティ at 23:53Comments(0)気になる水理学