【オブジェクト指向C++】2次元ナビエストークス方程式(キャビティ流れ)の数値計算

全体のコード
まとめ
参考書

全体のコード

全体のコードがこちらです。

main.cpp

#include <iostream>
#include <vector>
#include "Mesh.h" //メッシュ設定 
#include "Fields.h" //場の定義 
#include "Diff1d.h" //空間微分
#include "FileWriter.h" //ファイル出力

#include "PrintVector.h" //場のprint文
using std::vector;


int main()
{
	int nx = 20, ny = 20;
	double lengthx = 2.0, lengthy = 2.0;
	Mesh mesh2d_(nx, ny, lengthx, lengthy);//Meshクラスのインスタンス化

	/*
		U,V,P,bの定義
	*/
	#include "creatFields.h"

	// solution
	int ntimestep = 100;
	double dt = 0.001;
	double rho = 1.0;
	double visc = 1e-1;

	// file write
	FileWriter filewriter2d_;
	string name = "2DNS";


	string inlet = "INLET"; string outlet = "OUTLET";
	string right = "RIGHT"; string left = "LEFT";
	int nitmax = 50;
	double dxsq = mesh2d_.dx * mesh2d_.dx;
	double dysq = mesh2d_.dy * mesh2d_.dy;
	double dxsqdysq = dxsq * dysq;
	double deno = 2.0 * (dxsq + dysq);
	double conmined = dxsqdysq / deno;

	double bSpacing = (dxsq * dysq) / (2.0 * (dxsq + dysq));
	double invDelta = 1.0 / dt;
	double invrho = 1.0 / rho;
	double bSpacingRho = bSpacing * rho;
	double two = 2.0;
	double zerov = 0;
	double targetresidual = 1e-03;


	double ulid = 1.0;
	FieldsOperations.fixedValueBoundary(U, inlet, ulid);

	Pnew = P;
	Unew = U;
	Vnew = V;

	print2dfield(U, "U");

	for (int nt = 0; nt < ntimestep; nt++)
	{
		std::cout << "Iteration no : " << nt << std::endl;
		
		U = Unew;
		V = Vnew;
		
		//outcome
		if (nt % 10 == 0)
		{
			filewriter2d_.write2dUVFields(Unew, Vnew, Pnew, nt, name);
		}

		for (int nit = 0; nit < nitmax; nit++)
		{
			P = Pnew;
			b = -1.0 * (fd::gradXCDS(U) && fd::gradXCDS(U))
				- 2.0 * (fd::gradYCDS(U) && fd::gradXCDS(V))
				- (fd::gradYCDS(V) && fd::gradYCDS(V));
			Pnew = (fd::poissonX(P)) + (fd::poissonY(P)) - (bSpacingRho * b);
			// P boundary condition
			double leftboundaryValue = 0.0;
			FieldsOperations.fixedValueBoundary(Pnew, inlet, leftboundaryValue);
			FieldsOperations.zeroGradient(Pnew, outlet);
			FieldsOperations.zeroGradient(Pnew, right);
			FieldsOperations.zeroGradient(Pnew, left);

			double l1normal = FieldsOperations.L1norm(Pnew, P);
			//cout << "L1 norm  Pressure Eqn Residual = "<< nit << "  == " << l1normal << endl;
			if (l1normal < targetresidual)
			{
				std::cout << "pressure eqn is converged : " << nit << std::endl;
				std::cout << "L1 norm  Pressure Eqn Residual = " << l1normal << std::endl;
				break;
			}
		} //  Pressure Eqn end loop
		double l1normal = FieldsOperations.L1norm(Pnew, P);
		std::cout << "L1 norm  Pressure Eqn Residual = " << l1normal << std::endl;


		// du/dt+udu/dx+vdu/dy= -1/ρdp/dx + nu ∇^2 u
		// dv/dt+udv/dx+vdv/dy= -1/ρdp/dy + nu ∇^2 v
		Unew = U
			- dt * (U && fd::gradX(U)) - dt * (V && fd::gradY(U))
			- dt * invrho * fd::gradXCDS(P)
			+ dt * visc * (fd::laplacian2d(U));
		Vnew = V
			- dt * (U && fd::gradX(V)) - dt * (V && fd::gradY(V))
			- dt * invrho * fd::gradYCDS(P)
			+ dt * visc * (fd::laplacian2d(V));


		// UV boundary condition
		FieldsOperations.fixedValueBoundary(Unew, inlet, ulid);
		FieldsOperations.fixedValueBoundary(Unew, outlet, zerov);
		FieldsOperations.fixedValueBoundary(Unew, right, zerov);
		FieldsOperations.fixedValueBoundary(Unew, left, zerov);

		FieldsOperations.fixedValueBoundary(Vnew, inlet, zerov);
		FieldsOperations.fixedValueBoundary(Vnew, outlet, zerov);
		FieldsOperations.fixedValueBoundary(Vnew, right, zerov);
		FieldsOperations.fixedValueBoundary(Vnew, left, zerov);

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

#include <iostream>

#include <vector>

#include "Mesh.h" //メッシュ設定

#include "Fields.h" //場の定義

#include "Diff1d.h" //空間微分

#include "FileWriter.h" //ファイル出力

#include "PrintVector.h" //場のprint文

using std::vector;

int main()

{

int nx = 20, ny = 20;

double lengthx = 2.0, lengthy = 2.0;

Mesh mesh2d_(nx, ny, lengthx, lengthy);//Meshクラスのインスタンス化

U,V,P,bの定義

#include "creatFields.h"

// solution

int ntimestep = 100;

double dt = 0.001;

double rho = 1.0;

double visc = 1e-1;

// file write

FileWriter filewriter2d_;

string name = "2DNS";

string inlet = "INLET"; string outlet = "OUTLET";

string right = "RIGHT"; string left = "LEFT";

int nitmax = 50;

double dxsq = mesh2d_.dx * mesh2d_.dx;

double dysq = mesh2d_.dy * mesh2d_.dy;

double dxsqdysq = dxsq * dysq;

double deno = 2.0 * (dxsq + dysq);

double conmined = dxsqdysq / deno;

double bSpacing = (dxsq * dysq) / (2.0 * (dxsq + dysq));

double invDelta = 1.0 / dt;

double invrho = 1.0 / rho;

double bSpacingRho = bSpacing * rho;

double two = 2.0;

double zerov = 0;

double targetresidual = 1e-03;

double ulid = 1.0;

FieldsOperations.fixedValueBoundary(U, inlet, ulid);

Pnew = P;

Unew = U;

Vnew = V;

print2dfield(U, "U");

for (int nt = 0; nt < ntimestep; nt++)

{

std::cout << "Iteration no : " << nt << std::endl;

U = Unew;

V = Vnew;

//outcome

if (nt % 10 == 0)

{

filewriter2d_.write2dUVFields(Unew, Vnew, Pnew, nt, name);

}

for (int nit = 0; nit < nitmax; nit++)

{

P = Pnew;

b = -1.0 * (fd::gradXCDS(U) && fd::gradXCDS(U))

- 2.0 * (fd::gradYCDS(U) && fd::gradXCDS(V))

- (fd::gradYCDS(V) && fd::gradYCDS(V));

Pnew = (fd::poissonX(P)) + (fd::poissonY(P)) - (bSpacingRho * b);

// P boundary condition

double leftboundaryValue = 0.0;

FieldsOperations.fixedValueBoundary(Pnew, inlet, leftboundaryValue);

FieldsOperations.zeroGradient(Pnew, outlet);

FieldsOperations.zeroGradient(Pnew, right);

FieldsOperations.zeroGradient(Pnew, left);

double l1normal = FieldsOperations.L1norm(Pnew, P);

//cout << "L1 norm Pressure Eqn Residual = "<< nit << " == " << l1normal << endl;

if (l1normal < targetresidual)

{

std::cout << "pressure eqn is converged : " << nit << std::endl;

std::cout << "L1 norm Pressure Eqn Residual = " << l1normal << std::endl;

break;

}

} // Pressure Eqn end loop

double l1normal = FieldsOperations.L1norm(Pnew, P);

std::cout << "L1 norm Pressure Eqn Residual = " << l1normal << std::endl;

// du/dt+udu/dx+vdu/dy= -1/ρdp/dx + nu ∇^2 u

// dv/dt+udv/dx+vdv/dy= -1/ρdp/dy + nu ∇^2 v

Unew = U

- dt * (U && fd::gradX(U)) - dt * (V && fd::gradY(U))

- dt * invrho * fd::gradXCDS(P)

+ dt * visc * (fd::laplacian2d(U));

Vnew = V

- dt * (U && fd::gradX(V)) - dt * (V && fd::gradY(V))

- dt * invrho * fd::gradYCDS(P)

+ dt * visc * (fd::laplacian2d(V));

// UV boundary condition

FieldsOperations.fixedValueBoundary(Unew, inlet, ulid);

FieldsOperations.fixedValueBoundary(Unew, outlet, zerov);

FieldsOperations.fixedValueBoundary(Unew, right, zerov);

FieldsOperations.fixedValueBoundary(Unew, left, zerov);

FieldsOperations.fixedValueBoundary(Vnew, inlet, zerov);

FieldsOperations.fixedValueBoundary(Vnew, outlet, zerov);

FieldsOperations.fixedValueBoundary(Vnew, right, zerov);

FieldsOperations.fixedValueBoundary(Vnew, left, zerov);

main.cpp（メインファイル)以外のコードはgithubにアップしています。

その他のヘッダーファイルは以下となります。

#include “Mesh.h” //メッシュ設定
#include “Fields.h” //場の定義
#include “Diff1d.h” //空間微分
#include “FileWriter.h” //ファイル出力
#include “creatFields.h” //場の定義
#include “Diff1d.h” //微分の離散化定義

ソースコード

結果はParaViewで確認すると以下のようになります。

キャビティ流れ pic.twitter.com/BMWxa9fjXa
— カマキリ🐲Python頑張る昆虫 (@t_kun_kamakiri) July 29, 2022

まとめ

C++を使ってオブジェクト指向を意識してナビエストークス方程式を解くプログラムを組みました。

今回は「MAC（Marker And Cell）法」と呼ばれる非圧縮流体の解析手法のでナビエストークス方程式、およびナビエストークス方程式の発散と連続の式から得られる圧力方程式を連立して計算を行いました。

プログラムのメインは

		// du/dt+udu/dx+vdu/dy= -1/ρdp/dx + nu ∇^2 u
		// dv/dt+udv/dx+vdv/dy= -1/ρdp/dy + nu ∇^2 v
		Unew = U
			- dt * (U && fd::gradX(U)) - dt * (V && fd::gradY(U))
			- dt * invrho * fd::gradXCDS(P)
			+ dt * visc * (fd::laplacian2d(U));
		Vnew = V
			- dt * (U && fd::gradX(V)) - dt * (V && fd::gradY(V))
			- dt * invrho * fd::gradYCDS(P)
			+ dt * visc * (fd::laplacian2d(V));

// du/dt+udu/dx+vdu/dy= -1/ρdp/dx + nu ∇^2 u

// dv/dt+udv/dx+vdv/dy= -1/ρdp/dy + nu ∇^2 v

Unew = U