OpenFOAMでシュレーディンガー方程式の検討
オープンCAEのアドベントカレンダー2023の投稿記事です。
はじめにこちらの記事は自分のメモなので、全く有益な情報はないかもしれません。
それにうまくいっていません(笑)
というか言い訳をすると、ちょっと考えて、ちょっとやってみたという程度です。
ゆえに、進捗程度にお聞きください。
やりたいこと
やりたいことは本記事のタイトル通りです。
OpenFOAMは流体解析のオープンソースとして有名です。
流体解析がメインですが、偏微分方程式解法がクラスとしてまとめられていたりするので、原理的にはどういった偏微分方程式もカスタマイズして使用することができるはずです。
そこで、量子力学の基礎の方程式となるシュレディンガー方程式をOpenFOAMでできないかを考えます。
シュレディンガー方程式を書くと以下のような、複素関数の時間発展の式です。
シュレディンガー方程式は全く流体と関係がないわけではなく、似た形でグロスピタエフスキー方程式と呼ばれる多体系のボース粒子が従う方程式があります。
グロス＝ピタエフスキー方程式（グロス＝ピタエフスキーほうていしき、英: Gross–Pitaevskii equation; GPE）は、ボソン間相互作用が擬ポテンシャルとして表される理想的なボソン多体系の、ハートリー＝フォック近似の下での基底状態を記述するモデルである。
上記の方程式はボース＝アインシュタイン凝縮体の一粒子波動関数に対するモデル方程式となっている。グロス＝ピタエフスキー方程式はギンツブルグ＝ランダウ方程式と似た形をしており、また非線形シュレーディンガー方程式として言及されることも多い。
説明は端折り気味ですが、ここから量子流体という量子力学的な効果を持った流体現象を引き越します。
シュレディンガー方程式方程式は1粒子の確率密度に関係した関数が従う方程式ですが、グロスピタエフスキーは多体系の数密度が従う方程式です。
ですが、以下に示すように粒子間の相互作用が加わっただけでほとんど式はシュレディンガー方程式です。
量子力学的効果と言われているのは、上記のグロスピタエフスキー方程式から、ナビエストークス方程式と似た式を導くことができることから来ており、量子力学的な流体現象を示唆しています。
詳しくは、書籍でどうぞ_(._.)_
いきなり難しいことはできないので、まずはシュレディンガー方程式ができないか考えてみます。
※先に言っておきますが、現状できていません！
が、少しだけ考えたので良い方法があれば是非ともトライして共有してください_(._.)_
OpenFOAMでシュレディンガー方程式を考える
似たような例がないか探すと・・・・ありました。
さらに、実装方法はこちらを参考にできそうです。
まず、シュレディンガー方程式を方程式は複素関数を扱う式であり、式の中に虚数が入っています。
OpenFOAMには
volScalarField：圧力、温度、エネルギーなど
volVectorField：流速など
のようなクラスが存在しており、基本的に実関数を扱うクラスです。
volComplexScalarFieldみたいなものがあればいいですよね。
仮にあったとしても、OpenFOAMで用意されている行列式を解くことができるのか怪しいです。
そこで、実部と虚部に分けることを考えます。
これで、解くべき方程式は実関数のみになりました。
laplacianFoamをカスタマイズ
ほとんどlaplacianFoamと同じ式なので、適当にカスタマイズします。
laplacianFoamをMylaplacianFoamにコピーして編集していきます。
Make/files
MySchroedingerFoam.C
EXE = $(FOAM_USER_APPBIN)/MySchroedingerFoamMylaplacianFoam.C
#include "fvCFD.H"
#include "fvOptions.H"
#include "simpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"Laplace equation solver for a scalar quantity."
);
#include "postProcess.H"
#include "addCheckCaseOptions.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
simpleControl simple(mesh);
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating temperature distribution\n" << endl;
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
while (simple.correctNonOrthogonal())
{
// solve(fvm::ddt(T) - fvm::laplacian(DT, T));
fvScalarMatrix PsiRealEqn
(
fvm::ddt(PsiImg) - fvc::laplacian(DThbar/(2*mass), PsiReal)
//==
// fvOptions(PsiImg)
);
//add
/*
fvcにして陽解法にしている
空間微分にクランクニコルソン法を適用する場合は、fvc::laplacian(DT, Ti)ではなく
0.5*( fvm::laplacian(DT, Ti) + fvc::laplacian(DT, Ti) )
のように、TnとToの平均を計算するようにする
*/
fvScalarMatrix PsiImgEqn
(
fvm::ddt(PsiReal) + fvc::laplacian(DThbar/(2*mass), PsiImg)
//==
// fvOptions(PsiReal)
);
//fvOptions.constrain(PsiRealEqn);
PsiRealEqn.solve();
//fvOptions.correct(PsiReal);
//add
//fvOptions.constrain(PsiImgEqn);
PsiImgEqn.solve();
//fvOptions.correct(PsiImg);
//add end
}
forAll(PsiReal,i)
{
absPsi[i] = sqr(pow(PsiReal[i],2) + pow(PsiImg[i],2) );
}
//Info << "absT : " << absT << endl;
//#include "write.H"
runTime.write(); //add
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
createFields.H
Info<< "Reading field PsiReal\n" << endl;
volScalarField PsiReal
(
IOobject
(
"PsiReal",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field PsiImg\n" << endl;
volScalarField PsiImg
(
IOobject
(
"PsiImg",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField absPsi
(
IOobject
(
"absPsi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading diffusivity DT\n" << endl;
volScalarField DT
(
IOobject
(
"DT",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimViscosity, Zero)
);
Info<< "Reading diffusivity DThbar\n" << endl;
// m/s
dimensionSet dimhbar(1, 2, -1, 0, 0, 0, 0); // [kg m s K mol A Cd]
volScalarField DThbar
(
IOobject
(
"DThbar",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimhbar, Zero)
);
Info<< "Reading mass\n" << endl;
volScalarField mass
(
IOobject
(
"mass",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimMass, Zero)
);
Info << "dimhbar : " << dimhbar << endl;
if (!DThbar.headerOk())
{
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Info << "dimhbar : " << dimhbar << endl;
DT = dimensionedScalar("DT", dimViscosity, transportProperties);
DThbar = dimensionedScalar("DThbar", dimhbar , transportProperties);
mass = dimensionedScalar("mass", dimMass , transportProperties);
Info << "mass : " << mass << endl;
// Info << "DThbar/(2*mass) : " << DThbar.value()/(2*mass.value()) << endl;
}
#include "createFvOptions.H"
PsiReal：実部
PsiImg：虚部
absPsi：絶対値
DT：デフォルトであった定数（今回は使っていない）
DThbar：$${\hbar}$$（次元を合わせる必要がある）
mass：質量（次元を合わせる必要がある）
以上適当にカスタマイズ終了。
ケースファイルをコピー
とりあえず計算できるのかを確認します。
sysem/blockMeshDict
scale 1.0;
xmax 1;
ymax 0.1;
zmax 0.1;
vertices
(
( 0.0 0.0 0.0)
( $xmax 0.0 0.0)
( $xmax 0.0 $zmax)
( 0.0 0.0 $zmax)
( 0.0 $ymax 0.0)
( $xmax $ymax 0.0)
( $xmax $ymax $zmax)
( 0.0 $ymax $zmax)
);
blocks
(
hex (0 1 5 4 3 2 6 7) (100 1 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
XMin
{
type patch;
faces
(
(0 4 7 3)
);
}
XMax
{
type patch;
faces
(
(1 5 6 2)
);
}
ZMin
{
type empty;
faces
(
(0 1 5 4)
);
}
ZMax
{
type empty;
faces
(
(3 2 6 7)
);
}
YMin
{
type wall;
faces
(
(0 1 2 3)
);
}
YMax
{
type wall;
faces
(
(4 5 6 7)
);
}
);
mergePatchPairs
(
);メッシュは適当にこんな感じです。
system/fvSchems
ddtSchemes
{
default Euler;
}
gradSchemes
{
default none;
grad(T) Gauss linear;
grad(Ti) Gauss linear;
grad(PsiReal) Gauss linear;
grad(PsiImg) Gauss linear;
}
divSchemes
{
default none;
//div(phi,T) Gauss upwind grad(T);
}
laplacianSchemes
{
default none;
// laplacian(DT,T) Gauss linear corrected;
// laplacian(DT,Ti) Gauss linear corrected;
// laplacian(DThbar,T) Gauss linear corrected;
// laplacian(DThbar,Ti) Gauss linear corrected;
laplacian((DThbar|(2*mass)),PsiReal) Gauss linear corrected;
laplacian((DThbar|(2*mass)),PsiImg) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}system/fvSolution
solvers
{
PsiReal
{
solver PBiCGStab;
preconditioner DIC;//DILU;
tolerance 1e-06;
relTol 0;
}
PsiImg
{
$PsiReal
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;//2;
}
constant/transportProperties/
/DT 0.001;
DThbar 0.1;
mass 1.0;0.orig/absPsi
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object absPsi;
}
dimensions [0 0 0 0 0 0 0];//[0 0 0 1 0 0 0];
internalField uniform 0.0;
boundaryField
{
XMin
{
type zeroGradient;
}
XMax
{
type zeroGradient;
}
ZMin
{
type empty;
}
ZMax
{
type empty;
}
YMin
{
type slip;
}
YMax
{
type slip;
}
}
0.orig/PsiReal
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object PsiReal;
}
dimensions [0 0 0 0 0 0 0];//[0 0 0 1 0 0 0];
internalField uniform 0.0;
boundaryField
{
XMin
{
type zeroGradient;
}
XMax
{
type zeroGradient;
}
ZMin
{
type empty;
}
ZMax
{
type empty;
}
YMin
{
type slip;
}
YMax
{
type slip;
}
}
0.orig/PsiImg
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object PsiImg;
}
dimensions [0 0 0 0 0 0 0];//[0 0 0 1 0 0 0];
internalField uniform 0.0;
boundaryField
{
XMin
{
type zeroGradient;
}
XMax
{
type zeroGradient;
}
ZMin
{
type empty;
}
ZMax
{
type empty;
}
YMin
{
type slip;
}
YMax
{
type slip;
}
}
system/controlDict
application MySchroedingerFoam;//scalarTransportFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 0.1;
deltaT 1e-6;
writeControl runTime;//timeStep;
writeInterval 0.001;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
初期状態を作るためにsetFieldsDictを使用します。
system/setFieldsDict
defaultFieldValues
(
volScalarFieldValue PsiReal 0.0
//volVectorFieldValue U ( 0.1 0 0 )
);
regions
(
boxToCell
{
box ( 0.3 -1 -1 ) ( 0.7 1 1 );
fieldValues
(
volScalarFieldValue PsiReal 1
);
}
);適当な範囲を$${\psi_{rel}=1}$$にしました。
以下のコマンドで初期状態を作成します。
setFields
では、計算をさせます。
今回作った新たなソルバを実行します。
MySchroedingerFoamなにこれ？って感じですが、エラー無く計算できたので良かったです（良くはない・・・）
まずは、ひな形を作るのが大事ですよね・・・・
理論解があるシュレディンガー方程式の問題を設定して、計算の確からしさを検証する必要があります。
色々と今後手直しできればと思います。
まとめ
OpenFOAMでシュレディンガー方程式は解けるのか、少し考えてみました。
laplacianFoamをカスタマイズすれば可能であることがわかりましたが、複素関数を扱うシュレディンガー方程式は少々工夫が必要です。
とりあえず、適当にカスタマイズしてみましたが、「今回は計算ができた」程度のお粗末な内容ですが、今後手直ししていければと思います。
どなたか、興味持った方もやっていただければOpenFOAMも量子力学を扱えるようになる・・・・のかな？と思います。
X➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）
Youtube➡カマキリ理系ブログ
この記事が気に入ったらサポートをしてみませんか？