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

本記事ではC++を使った2次元ナビエストークス方程式を題材にオブジェクト指向でのコーディングを解説します。
最終的結果は以下のようにParaViewで可視化します。

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

前回の記事で1次元のC++のコードについて解説を行いました。

全体のコード

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

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);


		} // solution ends

	return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

#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);

} // solution ends

return 0;

}

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

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

コードが長いためmain.cpp(メインファイル)だけブログに載せておき、残りのファイルはgithubにアップしています。

ソースコード

1 2 3 4 5 6

全体のコード

【OpenFOAMの理解のためのC++メモ(1)】出力ストリーム、名前空間、関数、配列

【オブジェクト指向C++】1次元非線形偏微分方程式の数値計算

【OpenFOAMの理解のためのC++メモ(3)】関数テンプレート、クラステンプレート

【OpenFOAMの理解のためのC++メモ(2)】クラスの定義、ファイル分割、継承

【オブジェクト指向C++】1次元移流方程式の数値計算

【オブジェクト指向C++】1次元バーガース方程式の数値計算

【オブジェクト指向C++】1次元拡散方程式の数値計算

COMMENT コメントをキャンセル