codedsource源项不收敛
-
我在fvoptions里面想设置一个区域性的阻力源项 ，但是因为需要实现的结构比较特殊，所以没法用toposet去进行捕捉获取，然后知乎看见一个大佬用codedsource写的一个特定区域的动量源项合计自己试试去写一个这种，为了尝试设置的阻力源项可不可用我就toposet了一个区域合计试试，但是无论我的源项是纯数字还是仿darcy-forchheimer公式类型的他用interFoam跑都不收敛，找了好几天论坛视频啥的也找不到原因，vscode报错也只是给代码报错，偶尔代码看起来没错误但是跑起来就是不收敛。现在目前写的是如下这样的，因为回头要耦合dem，所以of版本是5.x。
momentumSource { type vectorCodedSource; active yes; name sourceTime; vectorCodedSourceCoeffs { selectionMode cellZone; // all; cellZone pZone; fields (U); codeInclude #{ #}; codeCorrect #{ // Pout<< "**codeCorrect**" << endl; #}; codeAddSup #{ // Pout<< "**codeAddSup**" << endl; // const vectorField& C = mesh_.C(); const scalarField& V = mesh_.V(); vectorField& Usource = eqn.source(); const vectorField& U = mesh().lookupObject<volVectorField>("U"); const scalarField& Rho = mesh().lookupObject<volScalarField>("rho"); const scalarField& nu = mesh().lookupObject<volScalarField>("nu"); const scalarField& magU = mag(U); scalar A = 10; scalar B = 10; forAll(V,i) { Usource -= (A * nu * Rho + B * Rho * magU / 2) * U ; } // Info << "***codeAddSup***" << nl; #}; codeSetValue #{ // Pout<< "**codeSetValue**" << endl; #}; // Dummy entry. Make dependent on above to trigger recompilation code #{ $codeInclude $codeCorrect $codeAddSup $codeSetValue #}; } sourceTimeCoeffs { $vectorCodedSourceCoeffs; } }请各位大佬各位巨佬帮着指点一下哪里出现了问题。
-
给源项公式改了一下，但是发现一个非常难受的事情，本身这个东西是仿照darcy-forchheimer写的，但是openfoam原生的模型的D和F可以提到1e10这么高的数量级，但是我这个写的感觉只能到1e2这样子，1e3开始就会浮点溢出，但是1e2这个水平只是能看出来这个特定区域确实有阻力存在，并不能实现极大多孔介质阻力并形成类固体区域的状态，所以就很好奇这个是什么情况，新改的代码是这样：
momentumSource { type vectorCodedSource; active yes; name sourceTime; vectorCodedSourceCoeffs { selectionMode all; // cellZone pZone; fields (U); codeInclude #{ #}; codeCorrect #{ // Pout<< "**codeCorrect**" << endl; #}; codeAddSup #{ // Pout<< "**codeAddSup**" << endl; // const vectorField& C = mesh_.C(); const scalarField& V = mesh_.V(); vectorField& Usource = eqn.source(); const vectorField& U = mesh().lookupObject<volVectorField>("U"); const scalarField& Rho = mesh().lookupObject<volScalarField>("rho"); const scalarField& nu = mesh().lookupObject<volScalarField>("nu"); // const scalarField& magU = mag(U); scalar A = 1e3; scalar B = 1e3; // vector C(0,1e4,0); forAll(V,i) { const scalar x = mesh_.C()[i][0]; const scalar y = mesh_.C()[i][1]; if(x < 0.5 && x > 0 && y < 0.5 && y > 0.45) { Usource[i] += (nu[i] * A + mag(U[i])* B * 0.5 ) * Rho[i] * U[i]* V[i]; // Usource = (A * U[i] + B * mag(U[i]) * U[i]) * V[i]; // Usource[i] += - C * V[i]; } } // Info << "***codeAddSup***" << nl; #}; codeSetValue #{ // Pout<< "**codeSetValue**" << endl; #}; // Dummy entry. Make dependent on above to trigger recompilation code #{ $codeInclude $codeCorrect $codeAddSup $codeSetValue #}; } sourceTimeCoeffs { $vectorCodedSourceCoeffs; } } -
然后又去pimpleFoam跑了一下，代码改成这样的：
momentumSource { type vectorCodedSource; active yes; name sourceTime; vectorCodedSourceCoeffs { selectionMode all; // cellZone pZone; fields (U); codeInclude #{ #}; codeCorrect #{ // Pout<< "**codeCorrect**" << endl; #}; codeAddSup #{ // Pout<< "**codeAddSup**" << endl; // const vectorField& C = mesh_.C(); const scalarField& V = mesh_.V(); vectorField& Usource = eqn.source(); const vectorField& U = mesh().lookupObject<volVectorField>("U"); // const scalarField& Rho = mesh().lookupObject<volScalarField>("rho"); const scalarField& nu = mesh().lookupObject<volScalarField>("nu"); // const scalarField& magU = mag(U); scalar A = 1e2; scalar B = 1e2; // vector C(0,1e4,0); forAll(V,i) { const scalar x = mesh_.C()[i][0]; const scalar y = mesh_.C()[i][1]; if(x < 0.5 && x > 0 && y < 0.5 && y > 0.45) { Usource[i] += (1e-5 * A + mag(U[i])* B * 0.5 ) * U[i]* V[i]; // Usource = (A * U[i] + B * mag(U[i]) * U[i]) * V[i]; // Usource[i] += - C * V[i]; } } // Info << "***codeAddSup***" << nl; #}; codeSetValue #{ // Pout<< "**codeSetValue**" << endl; #}; // Dummy entry. Make dependent on above to trigger recompilation code #{ $codeInclude $codeCorrect $codeAddSup $codeSetValue #}; } sourceTimeCoeffs { $vectorCodedSourceCoeffs; } }跑完的速度场是这样：
就完全处于A和B只能在1e2这个数量级，但凡再大一点就直接浮点溢出了。。。完全不知道是咋回事。。
