2017年02月21日
Javaを使った有限体積法による数値計算（２）
第２回目では、さっそくソースをアップする。
今回は、有限体積法の教科書としては、詳しい解説と、例題の丁寧な
回答が掲載されている
『数値流体力学 第2版』
の例題６．２をコード化した。
定常の移流方程式を解いている。
この例題を解いてみたのは、これから解きたい浅水流方程式に
近い形であること。つまり、速度場と密度場が未知であり、
これらを求めるためにSIMPLEアルゴリズムを適用するものであるためである。
また、断面形状も一定でなく、変化しているところも
今後の応用性が高く、参考になると思う。
Javaについては、１か月程度の知識しかなく、
静的メソッドのみの使用と、mainメソッドにほとんどの必要な処理を書いた。
初めて作ったJavaコードなので、至らないところが多い。
勉強しながら、現在はクラス等を用いたコードに修正を加えている。
また、多次元配列の初期化については、
「Arrays#fillメソッドの多次元版の一例」
を利用させていただいた。
連立方程式の解法には、トーマス法を用いている。
格子には等間隔のスタッガード格子を採用している。
また、例題の解答を直接コーディング化ではなく、最後の方に触れられている
不足緩和も係数計算の際に、組み込んでいる。
格子の数や不足緩和の係数を、いろいろ試してみたり、
質量フラックスを計算したりしてみて、有限体積法の面白さを
感じることができると思う。
（注）コードの利用による損害等は一切責任を負いません。
/*
* 有限体積法の計算
* 数値流体力学（第２版）
*
* 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;
}
}
今回は、有限体積法の教科書としては、詳しい解説と、例題の丁寧な
回答が掲載されている
『数値流体力学 第2版』
の例題６．２をコード化した。
定常の移流方程式を解いている。
この例題を解いてみたのは、これから解きたい浅水流方程式に
近い形であること。つまり、速度場と密度場が未知であり、
これらを求めるためにSIMPLEアルゴリズムを適用するものであるためである。
また、断面形状も一定でなく、変化しているところも
今後の応用性が高く、参考になると思う。
Javaについては、１か月程度の知識しかなく、
静的メソッドのみの使用と、mainメソッドにほとんどの必要な処理を書いた。
初めて作ったJavaコードなので、至らないところが多い。
勉強しながら、現在はクラス等を用いたコードに修正を加えている。
また、多次元配列の初期化については、
「Arrays#fillメソッドの多次元版の一例」
を利用させていただいた。
連立方程式の解法には、トーマス法を用いている。
格子には等間隔のスタッガード格子を採用している。
また、例題の解答を直接コーディング化ではなく、最後の方に触れられている
不足緩和も係数計算の際に、組み込んでいる。
格子の数や不足緩和の係数を、いろいろ試してみたり、
質量フラックスを計算したりしてみて、有限体積法の面白さを
感じることができると思う。
（注）コードの利用による損害等は一切責任を負いません。
/*
* 有限体積法の計算
* 数値流体力学（第２版）
*
* 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;
}
}
2017年02月16日
Javaを使った有限体積法による数値計算（１）
ここ数ヶ月、流体の数値計算の勉強を続けている。
２次元の河床変動プログラムを作るのが、当面の最終目標。
河床変動の計算の基礎式は、
１：流体の質量保存則（連続の式）
２；流体の運動量保存則（ナビエストークス方程式）
３；流砂の連続式
４；流砂量式
である。これら４つの方程式を解くことで
流れの場を一般的には求めることができる。
これらの式は複雑に絡み合い、特に２の運動量保存則は、
非線形の方程式（移流方程式）であり、一般に解析的に解くことが
できない。
そこで、数値計算が多用される。
数値計算においては、上記の方程式を離散化して、
解く方法がとられている。
その離散化の方法の違いにより、大きく
１： 格子法
２： 粒子法
に区別される
（他に個別要素法があるが、どう分類されるのかはまだ、勉強不足）
最近は２の粒子法も注目を集めているが、
計算範囲が大きい場合や、時間スケールを長くするには、
計算パワーが必要なことから、まだまだ一般的ではない。
粒子法については
『粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで C/C++ソースコード付』
が詳しくて分かりやすい。
このため、一般には１の格子法が採用される。
格子法は、さらに、
A：差分法
B：有限体積法
C；有限要素法
などに分類される。
差分法は、微分方程式をテイラー展開を利用して離散化する。
また、有限要素法は、重み付きの関数を用いて離散化する。
（有限要素法も勉強不足）
一方で、有限体積法は、コントロールボリュームと呼ばれる
有限な体積の要素に積分することで、その収支を用いて離散化する。
有限体積法では、コントロールボリュームの質量収支等を考慮するので、
質量が保存されやすいという特があり、この点が他の計算方法との違いと
言われている。
連続体としての質量の保存を重視する流体力学の数値計算では、
有限体積法が人気を集めているようである。
他の工学分野（構造力学等）では、有限要素法が一般的なようである。
以上のような理由から、有限体積法を用いて、
数値計算（数値シミュレーション）を進めていくこととした。
なお、使用言語は、
Javaを使っていくことにする（途中で変更するかもしれない）
もともと数値計算でよく用いられているFortranも使っていたが、
Androidアプリ開発のため、Javaを勉強し始めたことと、
オブジェクト指向を理解するためにも、Javaをプログラミング言語として
学んでいくことは、今後に役立つと考えたためである。
有限体積法の勉強に用いた参考文献を以下に示す。
２次元の河床変動プログラムを作るのが、当面の最終目標。
河床変動の計算の基礎式は、
１：流体の質量保存則（連続の式）
２；流体の運動量保存則（ナビエストークス方程式）
３；流砂の連続式
４；流砂量式
である。これら４つの方程式を解くことで
流れの場を一般的には求めることができる。
これらの式は複雑に絡み合い、特に２の運動量保存則は、
非線形の方程式（移流方程式）であり、一般に解析的に解くことが
できない。
そこで、数値計算が多用される。
数値計算においては、上記の方程式を離散化して、
解く方法がとられている。
その離散化の方法の違いにより、大きく
１： 格子法
２： 粒子法
に区別される
（他に個別要素法があるが、どう分類されるのかはまだ、勉強不足）
最近は２の粒子法も注目を集めているが、
計算範囲が大きい場合や、時間スケールを長くするには、
計算パワーが必要なことから、まだまだ一般的ではない。
粒子法については
『粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで C/C++ソースコード付』
が詳しくて分かりやすい。
このため、一般には１の格子法が採用される。
格子法は、さらに、
A：差分法
B：有限体積法
C；有限要素法
などに分類される。
差分法は、微分方程式をテイラー展開を利用して離散化する。
また、有限要素法は、重み付きの関数を用いて離散化する。
（有限要素法も勉強不足）
一方で、有限体積法は、コントロールボリュームと呼ばれる
有限な体積の要素に積分することで、その収支を用いて離散化する。
有限体積法では、コントロールボリュームの質量収支等を考慮するので、
質量が保存されやすいという特があり、この点が他の計算方法との違いと
言われている。
連続体としての質量の保存を重視する流体力学の数値計算では、
有限体積法が人気を集めているようである。
他の工学分野（構造力学等）では、有限要素法が一般的なようである。
以上のような理由から、有限体積法を用いて、
数値計算（数値シミュレーション）を進めていくこととした。
なお、使用言語は、
Javaを使っていくことにする（途中で変更するかもしれない）
もともと数値計算でよく用いられているFortranも使っていたが、
Androidアプリ開発のため、Javaを勉強し始めたことと、
オブジェクト指向を理解するためにも、Javaをプログラミング言語として
学んでいくことは、今後に役立つと考えたためである。
有限体積法の勉強に用いた参考文献を以下に示す。
2016年11月21日
水理学の勉強（１）
働き出して、水理学ってものが
あるのを知っても、なかなか良くわからなかったのですが、
昨年くらいから、腰を据えて勉強をはじめました。
初学者向けに、寒地土木研究所さんの
「現場の水理学」という貴重な資料があるのです。
http://river.ceri.go.jp/contents/tool/suirigaku.html
こちらに、不等流計算など初歩的なところから
のプログラムコードも掲載されています。
言語はfortranです。
でも、それをそのまま打ち込んでも、
動かなかったりしたので、
ほとんどそのままですが、
fortran90の勉強も兼ねて、自分で
ソースコードを作ってみました。
ということで、たまにそのコードを掲載させて
もらいたいと思います。
プログラミングにあたり、参考にした本は、
です。
演習問題３の１階のニュートン法によるもののソース
program newton
implicit none
integer i,j
integer,parameter ::nn=10
double precision q,n,ep,g
double precision aa,bb,cc
double precision fh,f1,dh,hhh
double precision d(nn),dl(nn),z(nn),b(nn),sh(nn),hh(nn)
d=(/ 0d0,500d0,1000d0,1200d0,1800d0,2100d0,2500d0,3000d0,3300d0,3800d0 /)
z=(/ 0d0,.5d0,.9d0,.8d0,2.0d0,2.3d0,3.0d0,3.0d0,3.5d0,4.0d0 /)
b=(/ 300d0,320d0,280d0,250d0,300d0,300d0,320d0,350d0,300d0,250d0 /)
do i=2,nn
dl(i)=d(i)-d(i-1)
enddo
q=1500d0
n=0.025d0
sh(1)=2.5d0
hh(1)=sh(1)+z(1)
ep=1d-3
g=9.8d0
do j=1,nn
if(j==1) then
write(*,*) j,d(j),dl(j),b(j),z(j),sh(j),hh(j)
cycle
endif
hhh=sh(j-1)
aa=q*q/(2*g*b(j)*b(j))
bb=-q*q*n*n*dl(j)/(2*b(j)*b(j))
cc=z(j)-(z(j-1)+sh(j-1)+q*q/(2*g*b(j-1)*b(j-1)*sh(j-1)*sh(j-1))+&
q*q*n*n*dl(j)/(2d0*b(j-1)*b(j-1)*sh(j-1)**(10d0/3d0)))
10fh=hhh+aa/(hhh**2d0)+bb/(hhh**(10d0/3d0))+cc
f1=1-2d0*aa/(hhh**3d0)-10d0*bb/(3d0*hhh**(13d0/3d0))
dh=-fh/f1
if(abs(fh) < ep ) then
sh(j)=hhh
hh(j)=sh(j)+z(j)
write(*,*) j,d(j),dl(j),b(j),z(j),sh(j),hh(j)
cycle
else
hhh=hhh+dh
goto 10
endif
enddo
end program newton
あるのを知っても、なかなか良くわからなかったのですが、
昨年くらいから、腰を据えて勉強をはじめました。
初学者向けに、寒地土木研究所さんの
「現場の水理学」という貴重な資料があるのです。
http://river.ceri.go.jp/contents/tool/suirigaku.html
こちらに、不等流計算など初歩的なところから
のプログラムコードも掲載されています。
言語はfortranです。
でも、それをそのまま打ち込んでも、
動かなかったりしたので、
ほとんどそのままですが、
fortran90の勉強も兼ねて、自分で
ソースコードを作ってみました。
ということで、たまにそのコードを掲載させて
もらいたいと思います。
プログラミングにあたり、参考にした本は、
です。
演習問題３の１階のニュートン法によるもののソース
program newton
implicit none
integer i,j
integer,parameter ::nn=10
double precision q,n,ep,g
double precision aa,bb,cc
double precision fh,f1,dh,hhh
double precision d(nn),dl(nn),z(nn),b(nn),sh(nn),hh(nn)
d=(/ 0d0,500d0,1000d0,1200d0,1800d0,2100d0,2500d0,3000d0,3300d0,3800d0 /)
z=(/ 0d0,.5d0,.9d0,.8d0,2.0d0,2.3d0,3.0d0,3.0d0,3.5d0,4.0d0 /)
b=(/ 300d0,320d0,280d0,250d0,300d0,300d0,320d0,350d0,300d0,250d0 /)
do i=2,nn
dl(i)=d(i)-d(i-1)
enddo
q=1500d0
n=0.025d0
sh(1)=2.5d0
hh(1)=sh(1)+z(1)
ep=1d-3
g=9.8d0
do j=1,nn
if(j==1) then
write(*,*) j,d(j),dl(j),b(j),z(j),sh(j),hh(j)
cycle
endif
hhh=sh(j-1)
aa=q*q/(2*g*b(j)*b(j))
bb=-q*q*n*n*dl(j)/(2*b(j)*b(j))
cc=z(j)-(z(j-1)+sh(j-1)+q*q/(2*g*b(j-1)*b(j-1)*sh(j-1)*sh(j-1))+&
q*q*n*n*dl(j)/(2d0*b(j-1)*b(j-1)*sh(j-1)**(10d0/3d0)))
10fh=hhh+aa/(hhh**2d0)+bb/(hhh**(10d0/3d0))+cc
f1=1-2d0*aa/(hhh**3d0)-10d0*bb/(3d0*hhh**(13d0/3d0))
dh=-fh/f1
if(abs(fh) < ep ) then
sh(j)=hhh
hh(j)=sh(j)+z(j)
write(*,*) j,d(j),dl(j),b(j),z(j),sh(j),hh(j)
cycle
else
hhh=hhh+dh
goto 10
endif
enddo
end program newton
