Coupled level set-VOF方法
-
各位大佬们好，我是研究波浪与结构物相互作用的，之前看到了下面这篇文章 https://www.sciencedirect.com/science/article/pii/S0301932206001832?via%3Dihub ，因为用VOF捕捉自由液面会存在钝化的问题，特别是想研究破碎波的时候，水气交界面的网格需要非常细，于是就想用一下这个方法。
根据这篇文章网上有公开的一部分代码，我就照着植入到了interFoam里面，取了个名字叫clsinterFoam，下面几张是我完成后对比两个求解器气泡上升的结果，可以看到clsinterFoam确实对交界面的捕捉更到位（左边是interFoam，右边是clsinterFoam）。
但是当我换到溃坝的算例的时候，出现了奇怪的问题，后处理结果来看，水体好像变得不连续了，最后会扩散到整个计算域。看上去像是没有了重力，但是算例里面我是加了的，想请教一下各位大佬，看看问题出在哪儿。
以下是在interFoam上改的部分代码：#include "mappingPsi.H" #include "solveLSFunction.H" #include "calcNewCurvature.H" turbulence->validate(); if (!LTS) { #include "CourantNo.H" #include "setInitialDeltaT.H" } Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { #include "readDyMControls.H" if (LTS) { #include "setRDeltaT.H" } else { #include "CourantNo.H" #include "alphaCourantNo.H" #include "setDeltaT.H" } runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) { if (pimple.firstIter() || moveMeshOuterCorrectors) { mesh.update(); if (mesh.changing()) { // Do not apply previous time-step mesh compression flux // if the mesh topology changed if (mesh.topoChanging()) { talphaPhi1Corr0.clear(); } gh = (g & mesh.C()) - ghRef; ghf = (g & mesh.Cf()) - ghRef; MRF.update(); if (correctPhi) { // Calculate absolute flux // from the mapped surface velocity phi = mesh.Sf() & Uf(); #include "correctPhi.H" // Make the flux relative to the mesh motion fvc::makeRelative(phi, U); mixture.correct(); } if (checkMeshCourantNo) { #include "meshCourantNo.H" } } } #include "alphaControls.H" #include "alphaEqnSubCycle.H" #include "mappingPsi.H" #include "solveLSFunction.H" #include "calcNewCurvature.H" #include "updateFlux.H" mixture.correct(); if (pimple.frozenFlow()) { continue; } #include "UEqn.H" // --- Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" } if (pimple.turbCorr()) { turbulence->correct(); } } // reInitialise the alpha equation if (runTime.outputTime()) { Info << "Overwriting alpha" << nl << endl; alpha1 = H; volScalarField& alpha10 = const_cast<volScalarField&>(alpha1.oldTime()); alpha10 = H.oldTime(); //const_cast<volScalarField&>(alpha1.storeOldTime()()) = H.oldTime(); } runTime.write(); runTime.printExecutionTime(Info); } Info<< "End\n" << endl; return 0; //mappingPsi.H // mapping alpha value to psi0 psi0 == (double(2.0)*alpha1-double(1.0))*gamma; //solvevelSFunction.H // solve Level-Set function as the re-initialization equation Info<< "solve the reinitialization equation" << nl << endl; psi == psi0; for (int corr=0; corr<int(epsilon.value()/deltaTau.value()); corr++) { psi = psi + psi0/mag(psi0)*(double(1)-mag(fvc::grad(psi)*dimChange))*deltaTau; psi.correctBoundaryConditions(); } // update Dirac function forAll(mesh.cells(),celli) { if(mag(psi[celli]) > epsilon.value()) delta[celli] = double(0); else delta[celli] = double(1.0)/(double(2.0)*epsilon.value())*(double(1.0)+Foam::cos(M_PI*psi[celli]/epsilon.value())); }; // update Heaviside function forAll(mesh.cells(),celli) { if(psi[celli] < -epsilon.value()) H[celli] = double(0); else if(epsilon.value() < psi[celli]) H [celli] = double(1); else H[celli] = double(1.0)/double(2.0)*(double(1.0)+psi[celli]/epsilon.value()+Foam::sin(M_PI*psi[celli]/epsilon.value())/M_PI); }; //calcNewCurvature.H // calculate normal vector volVectorField gradPsi(fvc::grad(psi)); surfaceVectorField gradPsif(fvc::interpolate(gradPsi)); surfaceVectorField nVecfv(gradPsif/(mag(gradPsif)+scalar(1.0e-6)/dimChange)); surfaceScalarField nVecf(nVecfv & mesh.Sf()); // calculate new curvature based on psi (LS function) C == -fvc::div(nVecf); //updataFlux.H { word alphaScheme("div(phi,alpha)"); word alpharScheme("div(phirb,alpha)"); // Standard face-flux compression coefficient surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf())); surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef(); //surfaceScalarField phic(mag(phi/mesh.magSf())); // Add the optional isotropic compression contribution if (icAlpha > 0) { phic *= (1.0 - icAlpha); phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U)); } // Do not compress interface at non-coupled boundary faces // (inlets, outlets etc.) forAll(phic.boundaryField(), patchi) { fvsPatchScalarField& phicp = phicBf[patchi]; if (!phicp.coupled()) { phicp == 0; } } surfaceScalarField phir(phic*nVecf); surfaceScalarField phiAlpha ( fvc::flux ( phi, alpha1, alphaScheme ) + fvc::flux ( -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), alpha1, alpharScheme ) ); MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2; } -
@Zhujh 您好，也采用了与您基本相同的CLSVOF方法进行计算，然而对于毛细张力主导的问题，这个方法表现出了更强的虚假流动（寄生流动）的问题，如无重力状态下的水滴。请问您遇到过类似的问题么？同时我也参考fluent里面的两种抑制虚假流动方法（密度和H函数），但是收效甚微，想问下您有什么建议吗
-
@Shihang-Chen 十分感谢，我去查一下
-
@Shihang-Chen 关于spurious currents/velocities，我们此前在OpenFOAM里植入了一个非常简单但是可行的方法，代码量与植入难度不大，具体请参考：https://www.sciencedirect.com/science/article/pii/S1877750323002557
-
@capillaryFix 太好了，感谢感谢，我会好好看看。我目前一直使用的消除寄生流的方法是基于sharpen surface model来的，这个效果也是很好。
-
@capillaryFix 这个流弊，老铁
-
@capillaryFix ，您好。我仔细研究了您的文章，有两点疑惑，向您请教：
- 最后一步是关于capillary flux的过滤是基于原始的$F_{st,f}$，过滤结束的$\phi^{filt}_{cf}$是放在这里吗？
solve ( UEqn == fvc::reconstruct ( ( - ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh) ) * mesh.magSf()+phi_filt_cf ), mesh.solver(U.select(finalIter)) );和
surfaceScalarField phig ( ( - ghf*fvc::snGrad(rho) )*rAUf*mesh.magSf()+ phi_filt_cf*rAUf );- 这步$F^{f}{st,f} = F{st,f} − F^{filt}_{st,f}$不知道怎么植入进方程，是不是关于capillary flux的过滤应该基于这步过滤的结果
谢谢、谢谢
-
@capillaryFix 是的，我求解了Eqn.30 -> Eqn.34 ->Eqn.33，但是Eqn.33的结果用在哪里了？
