java 超定方程求解_增加一个标量方程求解器

Bertha 。 2023-01-10 13:14 196阅读 0赞

![0bb704fad742007c7d271550859d41b1.png][]

有了网格，就可以来写CFD程序了。

之所以用java，是想利用java的面向对象特性，把CFD代码变得直观一点，以方便自己学习。

当然，相比C/C++，Java程序的速度会慢一些。

先来求解一个标量方程看看。CFD求解的输运方程通常为：

![554b0add21c548518d2c0e43b3e09f79.png][]

从左到右分别是瞬态项，对流项，扩散项和源项。

使用cell-centered有限体积法，求解的步骤为：

*  将方程在单元上离散
 *  求解代数方程组

**1.离散：**

离散的目的是要得到单元间的代数关系，形成如下的代数方程（NB表示C的相邻单元）：

![equation?tex=a\_%7BC%7D%5Cphi\_%7BC%7D%3D%5Csum\_%7B1%5Csim+n%7D%5E%7B%7D%7Ba\_%7BNB%7D%5Cphi\_%7BNB%7D%7D%2Bs\_%7BC%7D][equation_tex_a_7BC_7D_5Cphi_7BC_7D_3D_5Csum_7B1_5Csim_n_7D_5E_7B_7D_7Ba_7BNB_7D_5Cphi_7BNB_7D_7D_2Bs_7BC_7D]

![4995bc7f1e67ba75d18543e9ed265dee.png][]

忽略瞬态项，将输运方程在单元上积分，方程变为：

![equation?tex=%5Csum\_%7Bf%7D%5E%7B%7D%7B%5Crho+VS%5Cphi\_%7Bf%7D%3D%5Csum\_%7Bf%7D%5E%7B%7D%7B%5CGamma+grad%5Cphi++%5C+S%7D%7D%2BSource][equation_tex_5Csum_7Bf_7D_5E_7B_7D_7B_5Crho_VS_5Cphi_7Bf_7D_3D_5Csum_7Bf_7D_5E_7B_7D_7B_5CGamma_grad_5Cphi_5C_S_7D_7D_2BSource]

意思是：

*边界上流走的Φ=边界上扩散进来的Φ+单元自身产生的Φ*

使用cell-centered有限体积法时，基本变量都是位于单元中心的。因此上面方程里界面变量Φf要写成中心变量ΦC和ΦNB的表达式。

假设变量线性分布，使用最简单的中心差分，界面上的对流和扩散量可以写成：

对流通量：

![equation?tex=%5Crho+VS%5Cphi\_%7Bf%7D%3D%5Crho+VS%5Cphi\_%7BC%7D%281-%5Clambda%29%2B%5Crho+VS%5Cphi\_%7BNB%7D%5Clambda][equation_tex_5Crho_VS_5Cphi_7Bf_7D_3D_5Crho_VS_5Cphi_7BC_7D_281-_5Clambda_29_2B_5Crho_VS_5Cphi_7BNB_7D_5Clambda]

其中lambda为插值系数，f距离C和NB一样长时lambda=0.5

扩散通量：

![equation?tex=%5CGamma+grad%5Cphi+%5C+S%3D%5CGamma+%5Cfrac%7Bn%7D%7B%5Cleft%7C+N+%5Cright%7C%7DS%28%5Cphi\_%7BNB%7D-%5Cphi\_%7BC%7D%29][equation_tex_5CGamma_grad_5Cphi_5C_S_3D_5CGamma_5Cfrac_7Bn_7D_7B_5Cleft_7C_N_5Cright_7C_7DS_28_5Cphi_7BNB_7D-_5Cphi_7BC_7D_29]

其中n为N的单位向量 。

对单元C的每个面的通量进行求和，且做一个小忽略，方程最后变成

![equation?tex=%5Csum\_f%28-%5Crho+VS%5Clambda%2B%5CGamma+%5Cfrac%7Bn%7D%7B%5Cleft%7C+N+%5Cright%7C%7DS%29%5Cphi\_%7BC%7D%3D%5Csum\_f%28-%5Crho+VS%5Clambda%2B%5CGamma+%5Cfrac%7Bn%7D%7B%5Cleft%7C+N+%5Cright%7C%7DS%29%5Cphi\_%7BNB%7D%2Bsrc][equation_tex_5Csum_f_28-_5Crho_VS_5Clambda_2B_5CGamma_5Cfrac_7Bn_7D_7B_5Cleft_7C_N_5Cright_7C_7DS_29_5Cphi_7BC_7D_3D_5Csum_f_28-_5Crho_VS_5Clambda_2B_5CGamma_5Cfrac_7Bn_7D_7B_5Cleft_7C_N_5Cright_7C_7DS_29_5Cphi_7BNB_7D_2Bsrc]

这样就得到了我们想要的代数方程：

最后形成的系数矩阵A如下图

![92b8c07f5133aad7a9ed71431589f898.png][]

**代码实现：**

在face2d.class里，写一个calEqnCoef方法，来计算该面的系数，然后直接保存在矩阵A和s里：

/*in face2d class:   */
       void calScalarEqnCoef() {
         
            //mass flux 
            flux = rho * (u * Sx + v * Sy); 
            //diffusion coef
            diff = gamma * (Sx * Nx + Sy * Ny) / (Nx * Nx + Ny * Ny);
            //here use center difference
            a = diff - lamda * flux;
            //equation coef for cell C: a_C
            Eqn.A[C.id][C.id] += a;
    
            if (NB != null) {
         
                //equation coef for cell C: a_NB
                Eqn.A[C.id][NB.id] = -a;
    
                //equation coef for cell NB
                //for cell NB, diffusion coef is same, but flux is reverse.
                double a1 = diff + (1 - lamda) * flux; 
                Eqn.A[NB.id][NB.id] += a1;
                Eqn.A[NB.id][C.id] = -a1;
            }
            if (isInlet || isWall) {
         
                //move a*phi to source item.
                Eqn.s[C.id] += a * t;
            } else if (isOutlet) {
         
                //zero gradient condition, remove the contribution from a_C 
                Eqn.A[C.id][C.id] += -a;
            }
        }

这样只要对所有面进行一遍循环，就建立了系数矩阵A,s.

// in solver.class
        void calEqnCoef() {
         
            for (Face2d f : model.faceList) {
         
                f.calScalarEqnCoef();
            }
        }

**2.求解代数方程：**

网上能找到不少关于代数方程求解的java库，也可以自己写一个。这里就参考维基百科词条Bicgstab，写了一个共轭梯度法, 具体步骤见wiki词条或文末的代码。

// in solver.class
     void solveEqnf() {
         
         x = Matrix.solve(Eqn.A, Eqn.s, 50);
     }

求解了一个(静止)温度场，并且和FLUENT的结果进行了对比：

java T: min=50.04, max=79.25

Fluent T min=50.07, max=79.25

![96702b69dc363a095bca444ebdc1dd80.png][]

![07dde5bf9ae248d08cec089940e20fbd.png][]

附：BICGSTAB 求解代数方程：

/*in Matrix.class */
    /*BICGSTAB method to solve Ax=b */
    public static double[] solve(double[][] A, double[] b, int NLOOP) {
         
            int n=b.length;
            double[] X0 = new double[n];
            double[] R0 = sub(b,(multiply(A,X0)));
            double[] RR = copy(R0);
            double rho_0 = 1;
            double alpha = 1;
            double omega_0 = 1;
            double[] V0 = new double[n];
            double[] P0 = new double[n];
    
            double[]  P;
            double[]  X=new double[n] ;
            double[]  V;
            double[]  S;
            double[]  T;
            double[]  R;
    
            int nloop = 0;
            while (nloop < NLOOP) {
         
                //main loop
                double rho = dot(RR, R0);
                //System.out.println("rho is "+rho);
                double beta = rho / rho_0 * alpha / omega_0;
                P = add(R0,scale(sub(P0,scale(V0,omega_0)),beta));
                V = multiply(A,P);
                double temp = dot(RR, V);
                //check if it is 0, to avoid divide by zero
                if (temp == 0) {
         
                    X = new double[n];
                    return X;
                }
                alpha = rho / temp;
    
                S = sub(R0,scale(V,alpha));
                T = multiply(A,S);
                double omega = dot(T, S) / dot(T, T);
                X = add(add(X0,scale(P,alpha)),scale(S,omega));
                R = sub(S,scale(T,omega));
                nloop++;
                //update value
                R0 = R;
                P0 = P;
                X0 = X;
                V0 = V;
                rho_0 = rho;
                omega_0 = omega;
            }
            return X;
        }

[0bb704fad742007c7d271550859d41b1.png]: /images/20221119/b1da3cd7133f4eb49e1ea0cbdaac6d23.png
[554b0add21c548518d2c0e43b3e09f79.png]: /images/20221119/9c9506bd2aab4071a8c6d78d306518f7.png
[equation_tex_a_7BC_7D_5Cphi_7BC_7D_3D_5Csum_7B1_5Csim_n_7D_5E_7B_7D_7Ba_7BNB_7D_5Cphi_7BNB_7D_7D_2Bs_7BC_7D]: https://www.zhihu.com/equation?tex=a_%7BC%7D%5Cphi_%7BC%7D%3D%5Csum_%7B1%5Csim+n%7D%5E%7B%7D%7Ba_%7BNB%7D%5Cphi_%7BNB%7D%7D%2Bs_%7BC%7D
[4995bc7f1e67ba75d18543e9ed265dee.png]: https://img-blog.csdnimg.cn/img_convert/4995bc7f1e67ba75d18543e9ed265dee.png
[equation_tex_5Csum_7Bf_7D_5E_7B_7D_7B_5Crho_VS_5Cphi_7Bf_7D_3D_5Csum_7Bf_7D_5E_7B_7D_7B_5CGamma_grad_5Cphi_5C_S_7D_7D_2BSource]: https://www.zhihu.com/equation?tex=%5Csum_%7Bf%7D%5E%7B%7D%7B%5Crho+VS%5Cphi_%7Bf%7D%3D%5Csum_%7Bf%7D%5E%7B%7D%7B%5CGamma+grad%5Cphi++%5C+S%7D%7D%2BSource
[equation_tex_5Crho_VS_5Cphi_7Bf_7D_3D_5Crho_VS_5Cphi_7BC_7D_281-_5Clambda_29_2B_5Crho_VS_5Cphi_7BNB_7D_5Clambda]: https://www.zhihu.com/equation?tex=%5Crho+VS%5Cphi_%7Bf%7D%3D%5Crho+VS%5Cphi_%7BC%7D%281-%5Clambda%29%2B%5Crho+VS%5Cphi_%7BNB%7D%5Clambda
[equation_tex_5CGamma_grad_5Cphi_5C_S_3D_5CGamma_5Cfrac_7Bn_7D_7B_5Cleft_7C_N_5Cright_7C_7DS_28_5Cphi_7BNB_7D-_5Cphi_7BC_7D_29]: https://www.zhihu.com/equation?tex=%5CGamma+grad%5Cphi+%5C+S%3D%5CGamma+%5Cfrac%7Bn%7D%7B%5Cleft%7C+N+%5Cright%7C%7DS%28%5Cphi_%7BNB%7D-%5Cphi_%7BC%7D%29
[equation_tex_5Csum_f_28-_5Crho_VS_5Clambda_2B_5CGamma_5Cfrac_7Bn_7D_7B_5Cleft_7C_N_5Cright_7C_7DS_29_5Cphi_7BC_7D_3D_5Csum_f_28-_5Crho_VS_5Clambda_2B_5CGamma_5Cfrac_7Bn_7D_7B_5Cleft_7C_N_5Cright_7C_7DS_29_5Cphi_7BNB_7D_2Bsrc]: https://www.zhihu.com/equation?tex=%5Csum_f%28-%5Crho+VS%5Clambda%2B%5CGamma+%5Cfrac%7Bn%7D%7B%5Cleft%7C+N+%5Cright%7C%7DS%29%5Cphi_%7BC%7D%3D%5Csum_f%28-%5Crho+VS%5Clambda%2B%5CGamma+%5Cfrac%7Bn%7D%7B%5Cleft%7C+N+%5Cright%7C%7DS%29%5Cphi_%7BNB%7D%2Bsrc
[92b8c07f5133aad7a9ed71431589f898.png]: https://img-blog.csdnimg.cn/img_convert/92b8c07f5133aad7a9ed71431589f898.png
[96702b69dc363a095bca444ebdc1dd80.png]: https://img-blog.csdnimg.cn/img_convert/96702b69dc363a095bca444ebdc1dd80.png
[07dde5bf9ae248d08cec089940e20fbd.png]: https://img-blog.csdnimg.cn/img_convert/07dde5bf9ae248d08cec089940e20fbd.png