1.      非圧縮性SPH法とは

SPH(Smoothed Particle Hydrodynamics)法は，連続体の運動を有限個の粒子の運動として，離散化する粒子法の一種である．数値流体解析手法として，粒子法のほかに，空間を格子によって分割し，その格子に物理量を変数として割り当てて計算する格子法が用いられる．格子法では，計算の前処理として格子の作成が必要であるが，粒子法では格子の概念がないため格子作成のための前処理を必要とせず，計算コストを抑えることができる．また，格子法と比較して粒子法の方が，農薬散布や化学物質検知などの複数の相が混ざり合った混相流の計算が容易であるとされている．そのため，本研究では粒子法を用いた．

粒子法は，連続体の運動を有限個の粒子の運動として，離散化する数値解析手法であり，計算点(粒子)が物理量とともに移動するため，支配方程式をラグランジュ形式で表して計算している．その粒子法の一種であるSPH法では，粒子を質点ではなく，各粒子の質量が粒子の位置を代表する点の周りになめらかに分布するとしており，その分布の範囲が粒子の相互作用領域となる．この質量の広がりを通して，粒子系と流体系で物理量の対応関係が付けられ，流体系Navier-Stokes方程式から粒子系の運動方程式が導かれる．SPH法における基礎方程式の導出の過程は格子法と比較し，厳密なものではなく，かなりの程度発見的であるが，得られた方程式はニュートン力学の運動方程式であり，その結果は非常によく流体の運動を再現している(1)．

SPH法の中でも本研究で用いた非圧縮性SPH法では，従来の圧縮性SPH法と異なり，予測子–修正子法にもとづく反復解法によって，粘性や外力による変化があったとしても密度均一化処理を行い，密度一定を保ちつつ，速度場を更新することで，非圧縮性流れを近似的に実現している計算手法である(2)．非圧縮性を満たすために必要である密度均一化については，支配方程式の計算の時間刻みと比較して，十分に小さい時間間隔で密度場に対して密度均一化を行う．

 

 

 

 

2.    計算手法

　SPH法において基礎となる密度場は，カーネル関数によって連続的に質量分配された粒子の質量の足し合わせで与えられる．したがって，ある粒子の位置における密度は，粒子を粒子の周辺粒子とし，粒子間距離をとすると，以下のように表せる．

なお，粒子の質量は空間解像スケールを用いて，図1のような粒子塊を想定し，以下のように定義できる．

 

本研究で用いたカーネル関数は，粒子間距離と空間解像スケールで定義され，標準的に用いられるM4スプライン関数を用いた．空間解像スケールは，粒子1個の影響範囲の半分の長さに等しく，M4スプライン関数を用いたカーネル関数は，式(13)に示す粒子間距離を空間解像スケールで除した無次元化長さを用いて，以下のように表せる(図2)．

 

 

本研究で用いた支配方程式は，3次元Navier-Stokes方程式であり，流体要素の容積変化が無視できる低マッハ数の仮定より，密度，速度，時間，圧力，粘度，外力，ラグランジュ微分，ナブラをそれぞれ，，，，，，，とすると，ラグランジュ形式で以下のように表せる. また，本計算で使用した物性値を表16に示す．

式(1)の解は，時間刻み幅をとすると，以下のような積分方程式を満たす．

　式(6)で示した積分方程式における圧力項および粘性項は，SPHにおいて以下のように表現できる．

 

圧力項

式(6)の圧力項は，以下のように書き換えられる．

ここで，式(7)の各項は，中込の文献よりSPH法では，以下のように表せる(1)．

　したがって，式(8)および式(9)を代入すると式(7)は，以下のように表せる．

この時，圧力は音速の定義式より密度と関連づけられ，以下のように表せる．

　よって，式(11)を代入すると，式(10)は以下のように表せる．

　また，式(13)は中込の文献より流体系と粒子系の整合性を考えることにより，粒子，の目標密度である物性密度をそれぞれ，とすると，以下のように表せる(1)．

 

粘性項

式(2)の粘性項は，中込の文献よりSPH法では，以下のように表せる(1)．

　圧力項と同様に，中込の文献より流体系と粒子系の整合性を考えることにより，以下のように表せる(1)．

 

予測子–修正子法(図3)

先述したように，非圧縮性SPH法では従来の圧縮性SPH法と異なり，予測子–修正子法にもとづく反復解法によって，粘性や外力による変化があったとしても密度均一化処理を行うことで，密度均一つまり非圧縮を実現している．予測子–修正子法では，式(2)で示した3次元Navier-Stokes方程式の解を圧力項のみを残し，粘性項と外力項はそれぞれによる速度変化を加味した予測速度とし，以下のように表している．

ここで，非圧縮性であるため，時刻における密度場は，である．しかし，密度場は予測速度によってそれぞれの粒子座標が変位すると均一な密度場が失われる．この時，予測速度は音速と比較するとであるため，音速と結びつく復元力によって密度勾配は緩和に向かう．したがって，計算の時間刻みと比較して，十分に小さい時間間隔で密度場に対して密度均一化を行うことで，密度場は再び均一になり，が満たされると同時に，新たな速度場が得られる．

具体的な密度均一化手法としては，以下のとおりである．まず，時間刻みに対して十分に密度均一化時間刻み幅()において粘性および外力から算出した速度(以降，速度予測子)に基づいて粒子の座標を変位させる．この際に生じる圧力勾配によって，その速度予測子を修正する速度修正子を算出し正しい変位を実施する予測子–修正子法に基づく反復解法を用いて，式(18)を解く．

 

この時，圧力勾配は式(11)のように密度勾配と関連付けられるため，密度均一化の時間刻み幅は式(11)を用いて，以下のように表せる．

ここで，予測子–修正子法で用いる速度予測子を，速度修正子を，速度修正子に基づいて粒子を移動させた場合における変位を，速度修正子に基づいて粒子を移動させた場合における圧力の変化量を，速度修正子に基づいて粒子を移動させた場合における密度をとし，式(19)で示した密度均一化の時間刻み幅を用いると，式(18)は以下のように書き換えられる．

　ここで，式(13)を代入すると式(23)は，計算の対象とする粒子の速度をとすると，以下のように表せる．

 

本研究では，中込の文献にある安定した計算を行うための方法を採用し，粒子の計算密度がその粒子の目標密度（物性密度）を下回ったらその項を0として計算しないこととした．密度が実際の値より大きくなった粒子の影響は受けることができるので，実質的な密度均一化操作に支障は出ない(1)．伊藤の研究によると，この密度場補正において，反復回数を増やしすぎると数値拡散による計算上の粘性が増大する傾向があり，逆に反復回数を減らしすぎると密度場が物性値から乖離するため，カルマン渦列を再現した予備実験より，反復回数を1回とした(2)．また，式(23)で得られる速度に式(1)の粘性項と外力項を加えたものが，次時刻の速度修正子の初期値，すなわち予測子となり，以下のように表せる．

 

これにより，1つの計算ステップ終了し，次の計算ステップに移る．この計算の流れを図4にまとめた．

 

散布図 が含まれている画像 自動的に生成された説明

図1　粒子塊

 

図2　カーネル関数

 

ダイアグラム, 概略図 自動的に生成された説明

図3　予測子-修正子法

 

グラフィカル ユーザー インターフェイス 自動的に生成された説明

図4　非圧縮SPH法における計算フローチャート

 

 

簡単にSPH法について解説したが，以下の論文を参考にした,一読をおすすめする．※ts-kikuchilab内に保存済み

(1)       中込照明, “粒子法（SPH法）とその計算プログラムの実装についての解説”, Scientific and Educational Reports of the Faculty of Science and Technology, Kochi University, 2019 , 2 巻, No.1

(2)       伊藤真澄, “非圧縮性SPH法を用いたアーク溶接プロセス中流動現象の数値シミュレーション”, 東北大学大学院工学研究科機械システムデザイン工学専攻 2015年度博士学位論文, 2016

(3)       Manuel Hirschler, Ulrich Nieken, ‘Study of implicit time-integration in truly incompressible SPH’, SPHERIC, 2017

 

 

 

 

3.      計算コード

・定義

#define NMAX                  (300)

#define DT                         (5.0e-4f)

#define P_SPACING        (5.0e-4f)

#define H                            P_SPACING

#define A0                                    0.004f

#define B0                                    0.002f

#define AB                                    0.0045f

#define RHO_WATER                1000.0f

#define M_WATER         (RHO_WATER * P_SPACING * P_SPACING * P_SPACING)

#define C_WATER           (1.453e3f)  

#define NU_WATER                   (1.0038e-6f)

#define DAMPING_STEP                            (P_SPACING / C_WATER)  //3.441156*10^-7

#define DAMPLMT                  5

#define BG                                                      0.994f

#define SURFACE_MACRO          72.75e-3f

#define GAMMA                                           ( SURFACE_MACRO * P_SPACING / M_WATER )

 

NMAX

計算回数

DT

時間間隔

P_SPACING

隣り合う粒子との距離(中心と中心)

H

空間解像スケール

A0,B0

粒子を配置する楕円の大きさ

AB

計算における原点と表示の原点の調整

RHO_WATER

水の密度

M_WATER

水粒子1つ分の質量

C_WATER

水の音速

NU_WATER

水の動粘度

DAMPING_STEP

音速cと空間解像スケールhに依存する密度均一化時間刻み幅Δt’=h/c(<<Δt)

DAMPLMT

密度均一化のための反復回数

BG

目標密度と計算密度が同一にならないための補正係数

SURFACE_MACRO

20℃の水の表面張力 [N/m]

GAMMA 

粒子が持つ表面張力 [m/s^2]

EQT

表面張力均一化のための反復回数

 

・粒子数の計算(void set_particle_number)

void set_particle_number() {

              float a0 = A0, b0 = B0, caldl = 0.0f;

              int is = (int) (a0 * 2.0f / (float) H) + 1, js = (int) (b0 * 2.0f / (float) H) + 1, ks = (int) (b0 * 2.0f / (float) H) + 1;

              float inix = -a0 + (float) H * 0.5f, iniy = -b0 + (float) H * 0.5f, iniz = -b0 + (float) H * 0.5f;

              float x0, y0, z0;

              for (int k = 0; k < ks; k++) {

                            z0 = iniz + (float) k * (float) H;

                            for (int j = 0; j < js; j++) {

                                          y0 = iniy + (float) j * (float) H;

                                          for (int i = 0; i < is; i++) {

                                                        x0 = inix + (float) i * (float) H;

                                                        caldl = x0 * x0 / a0 / a0 + y0 * y0 / b0 / b0 + z0 * z0 / b0 / b0;

                                                        if (caldl < 1.0f) {

                                                                      np++;

                                                        }

                                          }

                            }

              }

              cout << "The number of water particles : " << np << endl;

}

 

ここでは，指定した範囲内(本計算ではA0・B0を定数とした楕円内)に隣り合う粒子との距離(粒子の大きさ)P_SPACING(H)[mm]で、いくつ配置できるか計算する．

 

・粒子の配置(void set_particle)

void set_particle(float xp[], float yp[], float zp[], float up[], float vp[], float wp[], float dp[]) {

              float a0 = A0, b0 = B0, caldl = 0.0f;

              int is = (int) (a0 * 2.0f / (float) H) + 1, js = (int) (b0 * 2.0f / (float) H) + 1, ks = (int) (b0 * 2.0f / (float) H) + 1;

              float inix = -a0 + (float) H * 0.5f, iniy = -b0 + (float) H * 0.5f, iniz = -b0 + (float) H * 0.5f;

              float x0, y0, z0;

              np = 0;

              for (int k = 0; k < ks; k++) {

                            z0 = iniz + (float) k * (float) H;

                            for (int j = 0; j < js; j++) {

                                          y0 = iniy + (float) j * (float) H;

                                          for (int i = 0; i < is; i++) {

                                                        x0 = inix + (float) i * (float) H;

                                                        caldl = x0 * x0 / a0 / a0 + y0 * y0 / b0 / b0 + z0 * z0 / b0 / b0;

                                                        if (caldl < 1.0f) {

                                                                      xp[np] = x0 + AB;

                                                                      yp[np] = y0 + AB;

                                                                      zp[np] = z0 + AB;

                                                                      up[np] = 0.0f;

                                                                      vp[np] = 0.0f;

                                                                      wp[np] = 0.0f;

                                                                      dp[np] = (float) RHO_WATER;

                                                                      np++;

                                                        }

                                          }

                            }

              }

}

 

先ほど計算した範囲内に配置できる粒子数を最大値として，粒子を実際に配置し，それぞれに番号を振る．また，ABによって計算における原点と表示の原点の調整している．

 

void gnuplot_locationは省略

 

・粒子の移動 (void realmove)

void realmove(float xp[], float yp[], float zp[], float xpp[], float ypp[], float zpp[], float up[], float vp[], float wp[], unsigned int np0) {

              for (unsigned int i = 0; i < np0; i++) {

                            xp[i] = xpp[i] + up[i] * DT;

                            yp[i] = ypp[i] + vp[i] * DT;

                            zp[i] = zpp[i] + wp[i] * DT;

              }

 

}

 

粒子の座標を速度と時間変化によって移動させている()．初期値については[void set_particle] にてそれぞれ0となっている．

 

・密度の計算(void getDensity)

void getDensity(float dp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, kernel, densityTemp, divH = 1.0f / H, c1 = divH * divH * divH / M_PI, mass = M_WATER;

              for (unsigned int i = 0; i < np0; i++) {

                            densityTemp = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          disx = xp[j] - xp[i];

                                          disy = yp[j] - yp[i];

                                          disz = zp[j] - zp[i];

                                          dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                          q = dis * divH;

                                          if (q < 2.0f) {

                                                        if (q < 1.0f)

                                                                      kernel = c1 * (1.0f - 1.5f * q * q + 0.75f * q * q * q);

                                                        else

                                                                      kernel = c1 * (2.0f - 3.0f * q + 1.5f * q * q - 0.25f * q * q * q);

 

                                                        densityTemp += mass * kernel;

                                          }

                            }

                            dp[i] = densityTemp;

              }

}

　

　ここでは，対象粒子の密度を他の粒子との影響によって計算している．

 

・加速度の計算(void getAcceleration)

void getAcceleration(float axp[], float ayp[], float azp[], float dp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, der_k, accelX, accelY, accelZ, ap, divDis;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI, mc2 = M_WATER * C_WATER * C_WATER, bgt = BG * RHO_WATER;

              for (unsigned int i = 0; i < np0; i++) {

                            accelX = 0.0f;

                            accelY = 0.0f;

                            accelZ = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        divDis = 1.0f / dis;

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    der_k = c1 * (-3.0f * q + 2.25f * q * q) * divH;

                                                                      else

                                                                                    der_k = c1 * (-3.0f + 3.0f * q - 0.75f * q * q) * divH;

 

                                                                      ap = mc2 * ((dp[i] - bgt) / (dp[i] * dp[i]) + (dp[j] - bgt) / (dp[j] * dp[j])) * divDis * der_k;

 

                                                                      accelX += ap * disx;

                                                                      accelY += ap * disy;

                                                                      accelZ += ap * disz;

                                                        }

                                          }

                            }

                            axp[i] = accelX;

                            ayp[i] = accelY;

                            azp[i] = accelZ;

              }

}

 

　ここでは，対象粒子の加速度をNavier-Stokes方程式の圧力項によって計算している．

 

・仮の粒子座標の計算(void imgMove)

void imgMove(float axp[], float ayp[], float azp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float timeStep = DAMPING_STEP;

              float temp;

              temp = 0.5f * timeStep * timeStep;

              for (unsigned int i = 0; i < np0; i++) {

                            xp[i] += axp[i] * temp;

                            yp[i] += ayp[i] * temp;

                            zp[i] += azp[i] * temp;

              }

}

 

　ここでは，対象粒子が持つ加速度に対する移動距離である。密度場補正のため反復するので仮の移動である．音速と空間解像スケールに依存する密度均一化時間刻み幅が時間刻みで，が時間刻みでないことに注意．

「getDensity、getAcceleration、imgMove」を繰り返す．

 

・表面張力(void surfaceAttraction)

void surfaceAttraction(float axsp[], float aysp[], float azsp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, kernel, accelX, accelY, accelZ, temp;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI, gamma = GAMMA;

              for (unsigned int i = 0; i < np0; i++) {

                            accelX = 0.0f;

                            accelY = 0.0f;

                            accelZ = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    kernel = q;

                                                                      else

                                                                                    kernel = 2.0f - q;

 

                                                                      temp = gamma * kernel / dis;

                                                                      accelX += temp * disx;

                                                                      accelY += temp * disy;

                                                                      accelZ += temp * disz;

                                                        }

                                          }

                            }

                            axsp[i] = accelX;

                            aysp[i] = accelY;

                            azsp[i] = accelZ;

              }

}

　

　ここでは，計算粒子が持つ表面張力を人工引力として計算している。

 

・表面張力の平均化(void aspSmoothing)

void aspSmoothing(float axsp[], float aysp[], float azsp[], float axp[], float ayp[], float azp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, kernel, accelX, accelY, accelZ, weight, divWeight;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI;

              for (unsigned int i = 0; i < np0; i++) {

                            accelX = axsp[i] * c1;

                            accelY = aysp[i] * c1;

                            accelZ = azsp[i] * c1;

                            weight = c1;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    kernel = c1 * (1.0f - 1.5f * q * q + 0.75f * q * q * q);

                                                                      else

                                                                                    kernel = c1 * (2.0f - 3.0f * q + 1.5f * q * q - 0.25f * q * q * q);

 

                                                                      accelX += axsp[j] * kernel;

                                                                      accelY += aysp[j] * kernel;

                                                                      accelZ += azsp[j] * kernel;

                                                                      weight += kernel;

                                                        }

                                          }

                            }

                            divWeight = 1.0f / weight;

                            axp[i] = accelX * divWeight;

                            ayp[i] = accelY * divWeight;

                            azp[i] = accelZ * divWeight;

              }

}

 

ここでは，重み関数(カーネル関数)を用いて計算粒子とその周辺粒子が持つ表面張力(人工引力)を平均化している．

 

・表面張力による移動(void surfaceMove)

void surfaceMove(float axsp[], float aysp[], float azsp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float temp;

              temp = 0.5f * DT * DT;

              for (unsigned int i = 0; i < np0; i++) {

                            xp[i] += axsp[i] * temp;

                            yp[i] += aysp[i] * temp;

                            zp[i] += azsp[i] * temp;

              }

}

 

　ここでは，平均化された表面張力(人工引力)より計算粒子を移動させている．

 

・速度の計算(void velocityStep)

void velocityStep(float xp[], float yp[], float zp[], float xpp[], float ypp[], float zpp[], float up[], float vp[], float wp[], unsigned int np0) {

              float divDT;

              divDT = 1.0f / DT;

              for (unsigned int i = 0; i < np0; i++) {

                            up[i] = (xp[i] - xpp[i]) * divDT;

                            vp[i] = (yp[i] - ypp[i]) * divDT; // - G_DIRECT*DT;

                            wp[i] = (zp[i] - zpp[i]) * divDT;

              }

}

 

　ここでは，移動した距離をもとに計算粒子の速度を計算している．

 

・粘性の計算(void viscosity)

void viscosity(float xp[], float yp[], float zp[], float up[], float vp[], float wp[], float axp[], float ayp[], float azp[], float dp[], unsigned int np0) {

              float disx, disy, disz, dis, q, der_k, ut1, vt1, wt1, temp;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI, nu = NU_WATER;

              for (unsigned int i = 0; i < np0; i++) {

                            ut1 = 0.0f;

                            vt1 = 0.0f;

                            wt1 = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    der_k = c1 * (-3.0f * q + 2.25f * q * q) * divH;

                                                                      else

                                                                                    der_k = c1 * (-3.0f + 3.0f * q - 0.75f * q * q) * divH;

 

                                                                      temp = (dp[i] + dp[j]) * der_k / (dp[i] * dp[j] * dis);

                                                                      ut1 -= temp * (up[j] - up[i]);

                                                                      vt1 -= temp * (vp[j] - vp[i]);

                                                                      wt1 -= temp * (wp[j] - wp[i]);

                                                        }

                                          }

                            }

                            temp = M_WATER * nu * DT;

                            axp[i] = up[i] + temp * ut1;

                            ayp[i] = vp[i] + temp * vt1;

                            azp[i] = wp[i] + temp * wt1;

              }

}

 

　ここでは，粘性を計算し，計算粒子が持つ加速度に加算している．

 

 

 

・プログラムの全容

（注）windowsの設定上「＼(半角)」が「\」に置き換わっている場合がある。

#include <iostream>

#include <fstream>

#include <cmath>

#include <time.h>

#include <string.h>

#include <stdlib.h>

#include <stdio.h>

#include <sys/stat.h>

#include <sys/types.h>

using namespace std;

 

#define NMAX                    (300)

#define DT                    (5.0e-4f)

#define P_SPACING            (5.0e-4f)

#define H                            P_SPACING

#define A0                                         0.004f

#define B0                                         0.002f

#define AB                                        0.0045f

#define RHO_WATER          1000.0f

#define M_WATER              (RHO_WATER * P_SPACING * P_SPACING * P_SPACING) /* mass for water particle */

#define C_WATER               (1.453e3f)   // sonic speed in water

#define NU_WATER            (1.0038e-6f) // kinematic viscosity

#define DAMPING_STEP    (P_SPACING / C_WATER)

#define DAMPLMT                           5

#define BG                                        0.994f

#define SURFACE_MACRO               72.75e-3f//（20℃水 72.75e-3[N/m]）

#define GAMMA                 ( SURFACE_MACRO * P_SPACING / M_WATER ) //[m/s^2]

 

unsigned int n = 0, np = 0;

 

void set_particle_number() {

              float a0 = A0, b0 = B0, caldl = 0.0f;

              int is = (int)(a0 * 2.0f / (float)H) + 1, js = (int)(b0 * 2.0f / (float)H) + 1, ks = (int)(b0 * 2.0f / (float)H) + 1;

              float inix = -a0 + (float)H * 0.5f, iniy = -b0 + (float)H * 0.5f, iniz = -b0 + (float)H * 0.5f;

              float x0, y0, z0;

              for (int k = 0; k < ks; k++) {

                            z0 = iniz + (float)k * (float)H;

                            for (int j = 0; j < js; j++) {

                                          y0 = iniy + (float)j * (float)H;

                                          for (int i = 0; i < is; i++) {

                                                        x0 = inix + (float)i * (float)H;

                                                        caldl = x0 * x0 / a0 / a0 + y0 * y0 / b0 / b0 + z0 * z0 / b0 / b0;

                                                        if (caldl < 1.0f) {

                                                                      np++;

                                                        }

                                          }

                            }

              }

              cout << "The number of water particles : " << np << endl;

              cout << "The number of calculation : " << NMAX << endl;

}

void set_particle(float xp[], float yp[], float zp[], float up[], float vp[], float wp[], float dp[]) {

              float a0 = A0, b0 = B0, caldl = 0.0f;

              int is = (int)(a0 * 2.0f / (float)H) + 1, js = (int)(b0 * 2.0f / (float)H) + 1, ks = (int)(b0 * 2.0f / (float)H) + 1;

              float inix = -a0 + (float)H * 0.5f, iniy = -b0 + (float)H * 0.5f, iniz = -b0 + (float)H * 0.5f;

              float x0, y0, z0;

              np = 0;

              for (int k = 0; k < ks; k++) {

                            z0 = iniz + (float)k * (float)H;

                            for (int j = 0; j < js; j++) {

                                          y0 = iniy + (float)j * (float)H;

                                          for (int i = 0; i < is; i++) {

                                                        x0 = inix + (float)i * (float)H;

                                                        caldl = x0 * x0 / a0 / a0 + y0 * y0 / b0 / b0 + z0 * z0 / b0 / b0;

                                                        if (caldl < 1.0f) {

                                                                      xp[np] = x0 + AB;

                                                                      yp[np] = y0 + AB;

                                                                      zp[np] = z0 + AB;

                                                                      up[np] = 0.0f;

                                                                      vp[np] = 0.0f;

                                                                      wp[np] = 0.0f;

                                                                      dp[np] = (float)RHO_WATER;

                                                                      np++;

                                                        }

                                          }

                            }

              }

}

void gnuplot_location(float xp[], float yp[], float zp[]) {

              char str2[20];

              sprintf(str2, "./img");

              mkdir(str2, 0755);

              sprintf(str2, "./img/location");

              mkdir(str2, 0755);

              FILE* gp;

              gp = popen("gnuplot", "w");

              fprintf(gp, "set terminal png font \"Times,45\" size 1920,1080\n");

              fprintf(gp, "set border lw 5.0\n");

              fprintf(gp, "set ticscale 0.8\n");

              fprintf(gp, "set ylabel \"{/Times-Italic y} [mm]\" rotate by 90 offset 0.0,0\n");

              fprintf(gp, "set xlabel \"{/Times-Italic x} [mm]\" offset 0,0.0\n");

              fprintf(gp, "set xrange[-0.0045:0.0045]\n");

              fprintf(gp, "set yrange[-0.0045:0.0045]\n");

              fprintf(gp, "set zrange[-0.0045:0.0045]\n");

              fprintf(gp, "set xtics (\"-4\" -0.004, \"-3\" -0.003, \"-2\" -0.002, \"-1\" -0.001, \"0\" 0, \"4\" 0.004, \"3\" 0.003, \"2\" 0.002, \"1\" 0.001,)\n");

              fprintf(gp, "set ytics (\"-4\" -0.004, \"-3\" -0.003, \"-2\" -0.002, \"-1\" -0.001, \"0\" 0, \"4\" 0.004, \"3\" 0.003, \"2\" 0.002, \"1\" 0.001,)\n");

              fprintf(gp, "set ztics (\"-4\" -0.004, \"-3\" -0.003, \"-2\" -0.002, \"-1\" -0.001, \"0\" 0, \"4\" 0.004, \"3\" 0.003, \"2\" 0.002, \"1\" 0.001,)\n");

              fprintf(gp, "set size ratio 1.0\n");

              fprintf(gp, "set nokey\n");

              fprintf(gp, "set view map\n");

              fprintf(gp, "set output \"%s/png%d.png\"\n", str2, n);

              fprintf(gp, "splot \"-\" u 1:2:3 with points pt 7 ps 3 lc rgb \"blue\"\n");

              for (int i = 0; i < np; i++) {

                            fprintf(gp, "%lf, %lf, %lf\n", xp[i] - AB, yp[i] - AB, zp[i] - AB);

              }

              fprintf(gp, "e\n");

              fflush(gp);

              pclose(gp);

}

void realmove(float xp[], float yp[], float zp[], float xpp[], float ypp[], float zpp[], float up[], float vp[], float wp[], unsigned int np0) {

              for (unsigned int i = 0; i < np0; i++) {

                            xp[i] = xpp[i] + up[i] * DT;

                            yp[i] = ypp[i] + vp[i] * DT;

                            zp[i] = zpp[i] + wp[i] * DT;

              }

 

}

void getDensity(float dp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, kernel, densityTemp, divH = 1.0f / H, c1 = divH * divH * divH / M_PI, mass = M_WATER;

              for (unsigned int i = 0; i < np0; i++) {

                            densityTemp = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          disx = xp[j] - xp[i];

                                          disy = yp[j] - yp[i];

                                          disz = zp[j] - zp[i];

                                          dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                          q = dis * divH;

                                          if (q < 2.0f) {

                                                        if (q < 1.0f)

                                                                      kernel = c1 * (1.0f - 1.5f * q * q + 0.75f * q * q * q);

                                                        else

                                                                      kernel = c1 * (2.0f - 3.0f * q + 1.5f * q * q - 0.25f * q * q * q);

 

                                                        densityTemp += mass * kernel;

                                          }

                            }

                            dp[i] = densityTemp;

              }

}

void getAcceleration(float axp[], float ayp[], float azp[], float dp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, der_k, accelX, accelY, accelZ, ap, divDis;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI, mc2 = M_WATER * C_WATER * C_WATER, bgt = BG * RHO_WATER;

              for (unsigned int i = 0; i < np0; i++) {

                            accelX = 0.0f;

                            accelY = 0.0f;

                            accelZ = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        divDis = 1.0f / dis;

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    der_k = c1 * (-3.0f * q + 2.25f * q * q) * divH;

                                                                      else

                                                                                    der_k = c1 * (-3.0f + 3.0f * q - 0.75f * q * q) * divH;

 

                                                                      ap = mc2 * ((dp[i] - bgt) / (dp[i] * dp[i]) + (dp[j] - bgt) / (dp[j] * dp[j])) * divDis * der_k;

 

                                                                      accelX += ap * disx;

                                                                      accelY += ap * disy;

                                                                      accelZ += ap * disz;

                                                        }

                                          }

                            }

                            axp[i] = accelX;

                            ayp[i] = accelY;

                            azp[i] = accelZ;

              }

}

void imgMove(float axp[], float ayp[], float azp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float timeStep = DAMPING_STEP;

              float temp;

              temp = 0.5f * timeStep * timeStep;

              for (unsigned int i = 0; i < np0; i++) {

                            xp[i] += axp[i] * temp;

                            yp[i] += ayp[i] * temp;

                            zp[i] += azp[i] * temp;

              }

}

void surfaceAttraction(float axsp[], float aysp[], float azsp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, kernel, accelX, accelY, accelZ, temp;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI, gamma = GAMMA;

              for (unsigned int i = 0; i < np0; i++) {

                            accelX = 0.0f;

                            accelY = 0.0f;

                            accelZ = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    kernel = q;

                                                                      else

                                                                                    kernel = 2.0f - q;

 

                                                                      temp = gamma * kernel / dis;

                                                                      accelX += temp * disx;

                                                                      accelY += temp * disy;

                                                                      accelZ += temp * disz;

                                                        }

                                          }

                            }

                            axsp[i] = accelX;

                            aysp[i] = accelY;

                            azsp[i] = accelZ;

              }

}

void aspSmoothing(float axsp[], float aysp[], float azsp[], float axp[], float ayp[], float azp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float disx, disy, disz, dis, q, kernel, accelX, accelY, accelZ, weight, divWeight;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI;

              for (unsigned int i = 0; i < np0; i++) {

                            accelX = axsp[i] * c1;

                            accelY = aysp[i] * c1;

                            accelZ = azsp[i] * c1;

                            weight = c1;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    kernel = c1 * (1.0f - 1.5f * q * q + 0.75f * q * q * q);

                                                                      else

                                                                                    kernel = c1 * (2.0f - 3.0f * q + 1.5f * q * q - 0.25f * q * q * q);

 

                                                                      accelX += axsp[j] * kernel;

                                                                      accelY += aysp[j] * kernel;

                                                                      accelZ += azsp[j] * kernel;

                                                                      weight += kernel;

                                                        }

                                          }

                            }

                            divWeight = 1.0f / weight;

                            axp[i] = accelX * divWeight;

                            ayp[i] = accelY * divWeight;

                            azp[i] = accelZ * divWeight;

              }

}

void surfaceMove(float axsp[], float aysp[], float azsp[], float xp[], float yp[], float zp[], unsigned int np0) {

              float temp;

              temp = 0.5f * DT * DT;

              for (unsigned int i = 0; i < np0; i++) {

                            xp[i] += axsp[i] * temp;

                            yp[i] += aysp[i] * temp;

                            zp[i] += azsp[i] * temp;

              }

}

void velocityStep(float xp[], float yp[], float zp[], float xpp[], float ypp[], float zpp[], float up[], float vp[], float wp[], unsigned int np0) {

              float divDT;

              divDT = 1.0f / DT;

              for (unsigned int i = 0; i < np0; i++) {

                            up[i] = (xp[i] - xpp[i]) * divDT;

                            vp[i] = (yp[i] - ypp[i]) * divDT; // - G_DIRECT*DT;

                            wp[i] = (zp[i] - zpp[i]) * divDT;

              }

}

void viscosity(float xp[], float yp[], float zp[], float up[], float vp[], float wp[], float axp[], float ayp[], float azp[], float dp[], unsigned int np0) {

              float disx, disy, disz, dis, q, der_k, ut1, vt1, wt1, temp;

              float divH = 1.0f / H, c1 = divH * divH * divH / M_PI, nu = NU_WATER;

              for (unsigned int i = 0; i < np0; i++) {

                            ut1 = 0.0f;

                            vt1 = 0.0f;

                            wt1 = 0.0f;

                            for (unsigned int j = 0; j < np0; j++) {

                                          if (j != i) {

                                                        disx = xp[j] - xp[i];

                                                        disy = yp[j] - yp[i];

                                                        disz = zp[j] - zp[i];

                                                        dis = sqrtf(disx * disx + disy * disy + disz * disz);

                                                        q = dis * divH;

                                                        if (q < 2.0f) {

                                                                      if (q < 1.0f)

                                                                                    der_k = c1 * (-3.0f * q + 2.25f * q * q) * divH;

                                                                      else

                                                                                    der_k = c1 * (-3.0f + 3.0f * q - 0.75f * q * q) * divH;

 

                                                                      temp = (dp[i] + dp[j]) * der_k / (dp[i] * dp[j] * dis);

                                                                      ut1 -= temp * (up[j] - up[i]);

                                                                      vt1 -= temp * (vp[j] - vp[i]);

                                                                      wt1 -= temp * (wp[j] - wp[i]);

                                                        }

                                          }

                            }

                            temp = M_WATER * nu * DT;

                            axp[i] = up[i] + temp * ut1;

                            ayp[i] = vp[i] + temp * vt1;

                            azp[i] = wp[i] + temp * wt1;

              }

}

int main() {

              set_particle_number();

 

              float* h_xp, * h_yp, * h_zp, * h_up, * h_vp, * h_wp, * h_dp, * h_xpp, * h_ypp, * h_zpp, * h_axp, * h_ayp, * h_azp, * h_axsp, * h_aysp, * h_azsp;

              h_xp = (float*)malloc(np * sizeof(float));

              h_yp = (float*)malloc(np * sizeof(float));

              h_zp = (float*)malloc(np * sizeof(float));

              h_up = (float*)malloc(np * sizeof(float));

              h_vp = (float*)malloc(np * sizeof(float));

              h_wp = (float*)malloc(np * sizeof(float));

              h_dp = (float*)malloc(np * sizeof(float));

              h_xpp = (float*)malloc(np * sizeof(float));

              h_ypp = (float*)malloc(np * sizeof(float));

              h_zpp = (float*)malloc(np * sizeof(float));

              h_axp = (float*)malloc(np * sizeof(float));

              h_ayp = (float*)malloc(np * sizeof(float));

              h_azp = (float*)malloc(np * sizeof(float));

              h_axsp = (float*)malloc(np * sizeof(float));

              h_aysp = (float*)malloc(np * sizeof(float));

              h_azsp = (float*)malloc(np * sizeof(float));

              memset(h_xp, 0, np * sizeof(float));

              memset(h_yp, 0, np * sizeof(float));

              memset(h_zp, 0, np * sizeof(float));

              memset(h_up, 0, np * sizeof(float));

              memset(h_vp, 0, np * sizeof(float));

              memset(h_wp, 0, np * sizeof(float));

              memset(h_dp, 0, np * sizeof(float));

              memset(h_xpp, 0, np * sizeof(float));

              memset(h_ypp, 0, np * sizeof(float));

              memset(h_zpp, 0, np * sizeof(float));

              memset(h_axp, 0, np * sizeof(float));

              memset(h_ayp, 0, np * sizeof(float));

              memset(h_azp, 0, np * sizeof(float));

              memset(h_axsp, 0, np * sizeof(float));

              memset(h_aysp, 0, np * sizeof(float));

              memset(h_azsp, 0, np * sizeof(float));

 

              set_particle(h_xp, h_yp, h_zp, h_up, h_vp, h_wp, h_dp);

              gnuplot_location(h_xp, h_yp, h_zp);

 

              for (n = 1; n < NMAX + 1; n++) {

                            memcpy(h_xpp, h_xp, np * sizeof(float));

                            memcpy(h_ypp, h_yp, np * sizeof(float));

                            memcpy(h_zpp, h_zp, np * sizeof(float));

                            realmove(h_xp, h_yp, h_zp, h_xpp, h_ypp, h_zpp, h_up, h_vp, h_wp, np);

 

                            for (int loop0 = 0; loop0 < DAMPLMT; loop0++) {

                                          getDensity(h_dp, h_xp, h_yp, h_zp, np);

                                          getAcceleration(h_axp, h_ayp, h_azp, h_dp, h_xp, h_yp, h_zp, np);

                                          imgMove(h_axp, h_ayp, h_azp, h_xp, h_yp, h_zp, np);

                            }

 

                            surfaceAttraction(h_axsp, h_aysp, h_azsp, h_xp, h_yp, h_zp, np);

                            for (int loop0 = 0; loop0 < 3; loop0++) {

                                          aspSmoothing(h_axsp, h_aysp, h_azsp, h_axp, h_ayp, h_azp, h_xp, h_yp, h_zp, np);

                                          memcpy(h_axsp, h_axp, np * sizeof(float));

                                          memcpy(h_aysp, h_ayp, np * sizeof(float));

                                          memcpy(h_azsp, h_azp, np * sizeof(float));

                            }

                            surfaceMove(h_axsp, h_aysp, h_azsp, h_xp, h_yp, h_zp, np);

 

                            gnuplot_location(h_xp, h_yp, h_zp);

 

                            velocityStep(h_xp, h_yp, h_zp, h_xpp, h_ypp, h_zpp, h_up, h_vp, h_wp, np);

                            viscosity(h_xp, h_yp, h_zp, h_up, h_vp, h_wp, h_axp, h_ayp, h_azp, h_dp, np);

                            memcpy(h_up, h_axp, np * sizeof(float));

                            memcpy(h_vp, h_ayp, np * sizeof(float));

                            memcpy(h_wp, h_azp, np * sizeof(float));

 

                            cout << n << endl;

              }

 

              free(h_xp);

              free(h_yp);

              free(h_zp);

              free(h_up);

              free(h_vp);

              free(h_wp);

              free(h_dp);

              free(h_xpp);

              free(h_ypp);

              free(h_zpp);

              free(h_axp);

              free(h_ayp);

              free(h_azp);

              free(h_axsp);

              free(h_aysp);

              free(h_azsp);

 

              cout << "end" << endl;

 

              return 0;

}