[OpenFOAM 소스코드 파해치기] 구배(Gradient)

일본어 원문링크

작성일 : 2009년 12월   8일
번역일 : 2016년  2월   10일

크롬 브라우저로 보시는 것을 권장해 드립니다.

 ------ OpenFOAM 소스코드 파해치기 시리즈 ------

OpenFOAM 소스코드 파해치기 목차로 이동

 ------------------------------------------------------------

 ----- OpenFOAM 소스코드를 다루는 문법 기본 -----

OpenFOAM 소스코드를 다루는 문법 기본으로 이동

 ------------------------------------------------------------

 

들어가기

유동장의 구배(Gradient)를 계산해 보자.

사용 버전

OpenFOAM 1.6

 

파일

이번에 사용파는 파일이다.(09-grad.tar.gz).

• Make/files
• Make/options
• system/controlDict
• system/fvSchemes
• system/fvSolution
• constant/polyMesh/blockMeshDict
• 0/T
• 0/gradT
• 0/gradT2
• heat.C

이전의 열전도 프로그램을 사용하여 온도장의 구배를 계산해보도록 하자. 구배를 벡터필드로 정의하고 2가지의 방법으로 계산하여 각각을 비교한다. 각각의 계산을 위해 경계조건 파일로 0/gradT, 0/gradT2를 준비한다. 구배계산에는 함수 grad를 사용하므로 이에 대한 설정을 fvSchemes에 적어놓아야 한다.

프로그램

 

설정파일의 설명은 생략하며, 프로그램 설명만 하겠다.

heat.C (일부)

        ...
	volScalarField T(
		IOobject("T", runTime.timeName(), mesh,
			IOobject::MUST_READ, IOobject::AUTO_WRITE),
		mesh
	);

	/* grad T */
	volVectorField gradT(
		IOobject("gradT", runTime.timeName(), mesh,
			IOobject::MUST_READ, IOobject::AUTO_WRITE),
		mesh
	);

	volVectorField gradT2(
		IOobject("gradT2", runTime.timeName(), mesh,
			IOobject::MUST_READ, IOobject::AUTO_WRITE),
		mesh
	);
	
	...
	
	/* calc grad T */
	gradT = fvc::grad(T);
	gradT.write();

	gradT2 *= 0.;
	surfaceScalarField Ts = fvc::interpolate(T);
	forAll(mesh.owner(), i){
		vector phi = Ts[i]*mesh.Sf()[i];
		gradT2[mesh.owner()[i]] += phi;
		gradT2[mesh.neighbour()[i]] -= phi;
	}

	forAll(mesh.boundary(), i){
		forAll(mesh.boundary()[i], j){
			gradT2[mesh.boundary()[i].faceCells()[j]]
			  += T.boundaryField()[i][j]
			    * mesh.Sf().boundaryField()[i][j];
		}
	}

	Field &fGradT2 = gradT2;
	fGradT2 /= mesh.V();

	gradT2.correctBoundaryConditions();

	gradT2.write();
	...

구배를 2가지 방법으로 계싼하도록 한다.  gradT는 fvc::grad(T) 함수로 계산하고 있다.  gradT2는 직접 계산하고 있다. OpenFOAM은 대부분 유한체적법의 코드이며 이번에 사용하는것도 유한체적법에 해당하므로 gradT2의 계산에는 유한체적법의 개념으로 계산하도록 한다. 즉, 구배의 검사체적(Control Volume)적분에 가우스 발산정리를 사용하여 얻어진 식을 셀체적으로 나누어 이산화하여 계산한다. 즉, 셀 하나하나에 대하여 각 셀의 모든 면에 대해 면의 면적벡터(크기가 면의 면적이고 방향이 셀의 바깥쪽으로 면에 대한 법선방향)와 면의 온도를 곱한 것을 합하는 조작을 수행한다. (이번에 사용된 코드는 src/finiteVolume/finiteVolume/gradSchemes/GaussGrad/gaussGrad.C 를 참고로 하고 있습니다).

자세히 보면,

gradT = fvc::grad(T);

grad(T)가 T의 구배 그 자체가 된다.

gradT2 *= 0.;

단순히 초기화를 수행하는 코드에 해당한다.

surfaceScalarField Ts = fvc::interpolate(T);

T의 셀의 면에 대한 값의 집합에 해당한다.

        forAll(mesh.owner(), i){
		vector phi = Ts[i]*mesh.Sf()[i];
		gradT2[mesh.owner()[i]] += phi;
		gradT2[mesh.neighbour()[i]] -= phi;
	}

mesh.owner()로 격자파일의 owner로부터 읽어들인 정보를 참조하고 있다. 즉, 셀의 내부면에 대한 루프이다. mesh.Sf()는 셀의 면의 면적벡터 리스트이다. 면의 온도와 면적 벡터를 곱하고 그것을 해당하는 셀의 구배에 더하고 있다. 유한체적법에 의한 직관적인 기술방법으로는 셀의 루프 내에서 셀의 면에 대한 루프를 돌려 면의 값을 면적 벡터와 곱해 합산해 셀의 구배를 계산하여야 하나, 여기서는 전체의 면의 루프를 돌려, 각각의 면에 관련된 셀에 대해 면의 값과 면적벡터의 곱을 더하는 기술방식으로 되어있다. mesh.neighbour()는 격자파일 neighbour의 정보를 출력한다. mesh.owner()[i]와 mesh.neighbour()[i]로 더할때 phi의 부호가 다른것은 옆 셀에 대해 대상 면에 대한 면적벡터의 방향이 반대(셀의 내부방향)이기 때문이다.

 	forAll(mesh.boundary(), i){
		forAll(mesh.boundary()[i], j){
			gradT2[mesh.boundary()[i].faceCells()[j]]
			  += T.boundaryField()[i][j]
			    * mesh.Sf().boundaryField()[i][j];
		}
	}
 

경계에 대한 처리이다.

Field &fGradT2 = gradT2;
fGradT2 /= mesh.V();

마지막으로 셀의 체적으로 나누어준다. 셀 하나씩 하는것은 코드가 지저분해지므로 위와 같은 방식으로 계산한다.

실행

컴파일하여 실행해 본다.

$ wmake
$ blockMesh
$ ./heat

paraFoam에서 결과를 확인해보자.

추가정보 1

 

하는김에 라플라시안도 계산해 보도록 하자(09-grad2.tar.gz).

heat.C (일부)

	...
	/* laplacian T */
	volScalarField laplacianT(
		IOobject("laplacianT", runTime.timeName(), mesh,
			IOobject::MUST_READ, IOobject::AUTO_WRITE),
		mesh
	);

	volScalarField laplacianT2(
		IOobject("laplacianT2", runTime.timeName(), mesh,
			IOobject::MUST_READ, IOobject::AUTO_WRITE),
		mesh
	);
	
	...
	
	/* calc laplacian T */
	laplacianT = fvc::laplacian(DT, T);
	laplacianT.write();

	laplacianT2 *= 0.;
	surfaceVectorField gradTs = fvc::interpolate(gradT2);
	forAll(mesh.owner(), i){
		scalar phi = DT.value()*dot(gradTs[i], mesh.Sf()[i]);
		laplacianT2[mesh.owner()[i]] += phi;
		laplacianT2[mesh.neighbour()[i]] -= phi;
	}

	forAll(mesh.boundary(), i){
		forAll(mesh.boundary()[i], j){
			laplacianT2[mesh.boundary()[i].faceCells()[j]]
			  += dot(gradT2.boundaryField()[i][j],
			    mesh.Sf().boundaryField()[i][j]);
		}
	}

	Field &fLaplacianT2 = laplacianT2;
	fLaplacianT2 /= mesh.V();
	...

구배의 계산과 같다. T의 라플라시안은 T가 방정식을 만족하면 0이 될 것이다.

추가정보 2

 

구배계산의 프로그램을 조금 수정해 보자(09-grad3.tar.gz).

heat.C (일부)

	gradT2 *= 0.;
	surfaceVectorField phi = fvc::interpolate(T) * mesh.Sf();
	forAll(mesh.owner(), i){
		gradT2[mesh.owner()[i]] += phi[i];
		gradT2[mesh.neighbour()[i]] -= phi[i];
	}

phi를 벡터가 아닌 벡터필드로 한번에 계산하도록 해 보았다.

추가정보 3

 

라플라시안계산의 프로그램도 같은 방법으로 수정해 보자(09-grad4.tar.gz).

heat.C (일부)

	laplacianT2 *= 0.;
	{
		surfaceScalarField phi = fvc::interpolate(gradT2) & mesh.Sf();
		phi *= DT.value();
		forAll(mesh.owner(), i){
			laplacianT2[mesh.owner()[i]] += phi[i];
			laplacianT2[mesh.neighbour()[i]] -= phi[i];
		}
	}

변수명이 충돌하지 않도록 고려하고 있으며,

SurfaceScalarField phi = fvc::interpolate(gradT2) & mesh.Sf();

벡터필드간의 내적을 취할때에는 &를 사용한다.

OpenFOAM 소스코드 파해치기 목차로 이동