From 46abab8cd3e7d31304d1c914a82d28a46364cd21 Mon Sep 17 00:00:00 2001
From: niikonst <niikonst@blogin2.usr.hlrn.de>
Date: Thu, 29 Jun 2023 13:48:30 +0200
Subject: [PATCH] change Sharp Interface implementation
---
apps/cpu/JetBreakup/JetBreakup.cpp | 148 ++++++++----------
apps/cpu/RisingBubble2D/RisingBubble2D.cpp | 51 +++---
.../MultiphaseScaleDistributionLBMKernel.cpp | 63 +++++---
.../LBM/MultiphaseSharpInterfaceLBMKernel.cpp | 106 +++++++++----
4 files changed, 213 insertions(+), 155 deletions(-)
diff --git a/apps/cpu/JetBreakup/JetBreakup.cpp b/apps/cpu/JetBreakup/JetBreakup.cpp
index f4b74ca37..235987e83 100644
--- a/apps/cpu/JetBreakup/JetBreakup.cpp
+++ b/apps/cpu/JetBreakup/JetBreakup.cpp
@@ -177,6 +177,19 @@ void run(string configname)
UBLOG(logINFO, "path = " << pathname);
}
+ // bounding box
+ real g_minX1 = 0;
+ real g_minX2 = 0;
+ real g_minX3 = 0;
+
+ //double g_maxX1 = 8.0*D;
+ //double g_maxX2 = 2.5*D;
+ //double g_maxX3 = 2.5*D;
+
+ real g_maxX1 = 8.0 * D;
+ real g_maxX2 = 5 * D;
+ real g_maxX3 = 5 * D;
+
SPtr<LBMUnitConverter> conv(new LBMUnitConverter());
// const int baseLevel = 0;
@@ -189,7 +202,8 @@ void run(string configname)
// kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel());
// kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
//kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
- kernel = SPtr<LBMKernel>(new MultiphaseSimpleVelocityBaseExternalPressureLBMKernel());
+ //kernel = SPtr<LBMKernel>(new MultiphaseSimpleVelocityBaseExternalPressureLBMKernel());
+ kernel = make_shared<MultiphaseScaleDistributionLBMKernel>();
kernel->setWithForcing(true);
kernel->setForcingX1(0.0);
@@ -209,6 +223,7 @@ void run(string configname)
kernel->setContactAngle(theta);
kernel->setInterfaceWidth(interfaceWidth);
//dynamicPointerCast<MultiphasePressureFilterLBMKernel>(kernel)->setPhaseFieldBC(0.0);
+ kernel->setSigma(sigma);
SPtr<BCSet> bcProc(new BCSet());
// BCSetPtr bcProc(new ThinWallBCSet());
@@ -240,25 +255,26 @@ void run(string configname)
// BC Adapter
//////////////////////////////////////////////////////////////////////////////
mu::Parser fctF1;
- // fctF1.SetExpr("vy1*(1-((x1-x0)^2+(x3-z0)^2)/(R^2))");
- // fctF1.SetExpr("vy1*(1-(sqrt((x1-x0)^2+(x3-z0)^2)/R))^0.1");
- fctF1.SetExpr("vy1");
- fctF1.DefineConst("vy1", 0.0);
- fctF1.DefineConst("R", 8.0);
- fctF1.DefineConst("x0", 0.0);
- fctF1.DefineConst("z0", 0.0);
- // SPtr<BC> velBCF1(
- // new MultiphaseVelocityBC(false, true, false, fctF1, phiH, 0.0, BCFunction::INFCONST));
+ //u_actual(r) = U_max * (4 / (D^2) * (R^2 - r^2))
+ //fctF1.SetExpr("vx1*(1-((x2-y0)^2+(x3-z0)^2)/(R^2))");
+ //fctF1.SetExpr("vx1*(1-(sqrt((x2-y0)^2+(x3-z0)^2))/(R))");
+ //fctF1.SetExpr("vy1*(1-(sqrt((x2-x0)^2+(x3-z0)^2)/R))^0.1");
+ fctF1.SetExpr("vx1");
+ //fctF1.DefineConst("vx1", 0);
+ fctF1.DefineConst("vx1", U_LB);
+ fctF1.DefineConst("R", 0.5*D);
+ fctF1.DefineConst("y0", (g_minX2+g_maxX2)/2);
+ fctF1.DefineConst("z0", (g_minX3+g_maxX3)/2);
+ SPtr<BC> velBCF1(new MultiphaseVelocityBC(true, false, false, fctF1, phiH, 0.0, BCFunction::INFCONST));
mu::Parser fctF2;
- fctF2.SetExpr("vy1");
- fctF2.DefineConst("vy1", U_LB);
+ fctF2.SetExpr("vx1");
+ fctF2.DefineConst("vx1", U_LB);
- real startTime = 1;
- SPtr<BC> velBCF1(
- new MultiphaseVelocityBC(true, false, false, fctF1, phiH, 0.0, startTime));
- SPtr<BC> velBCF2(
- new MultiphaseVelocityBC(true, false, false, fctF2, phiH, startTime, endTime));
+ // real startTime = 1;
+ // SPtr<BC> velBCF1(new MultiphaseVelocityBC(true, false, false, fctF1, phiH, 0.0, startTime));
+ //SPtr<BC> velBCF2(new MultiphaseVelocityBC(true, false, false, fctF2, phiH, startTime, endTime));
+ SPtr<BC> velBCF2(new MultiphaseVelocityBC(true, false, false, fctF2, phiL, 0.0, BCFunction::INFCONST));
SPtr<BC> noSlipBC(new NoSlipBC());
noSlipBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseNoSlipBCStrategy()));
@@ -283,7 +299,7 @@ void run(string configname)
// BC visitor
MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
bcVisitor.addBC(noSlipBC);
- bcVisitor.addBC(denBC); // Ohne das BB?
+ //bcVisitor.addBC(denBC); // Ohne das BB?
bcVisitor.addBC(velBCF1);
//SPtr<D3Q27Interactor> inflowF1Int;
@@ -295,18 +311,7 @@ void run(string configname)
// if (newStart) {
- // bounding box
- real g_minX1 = 0;
- real g_minX2 = 0;
- real g_minX3 = 0;
- //double g_maxX1 = 8.0*D;
- //double g_maxX2 = 2.5*D;
- //double g_maxX3 = 2.5*D;
-
- real g_maxX1 = 1.0 * D; // 8.0 * D;
- real g_maxX2 = 2.0 * D;
- real g_maxX3 = 2.0 * D;
// geometry
SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
@@ -330,18 +335,11 @@ void run(string configname)
// GbSystem3D::writeGeoObject(geoInflowF1.get(), pathname + "/geo/geoInflowF1",
// WbWriterVtkXmlASCII::getInstance());
- GbCylinder3DPtr geoInflow(new GbCylinder3D(g_minX1 - 2.0*dx, g_maxX2 / 2.0, g_maxX3 / 2.0, g_minX1,
- g_maxX2 / 2.0,
- g_maxX3 / 2.0, D / 2.0));
- if (myid == 0)
- GbSystem3D::writeGeoObject(geoInflow.get(), pathname + "/geo/geoInflow",
- WbWriterVtkXmlASCII::getInstance());
+ GbCylinder3DPtr geoInflow(new GbCylinder3D(g_minX1 - 2.0*dx, g_maxX2 / 2.0, g_maxX3 / 2.0, g_minX1, g_maxX2 / 2.0, g_maxX3 / 2.0, D / 2.0));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), pathname + "/geo/geoInflow", WbWriterVtkXmlASCII::getInstance());
- GbCylinder3DPtr geoSolid(new GbCylinder3D(g_minX1 - 2.0 * dx, g_maxX2 / 2.0, g_maxX3 / 2.0, g_minX1-dx,
- g_maxX2 / 2.0, g_maxX3 / 2.0, 1.5*D / 2.0));
- if (myid == 0)
- GbSystem3D::writeGeoObject(geoSolid.get(), pathname + "/geo/geoSolid",
- WbWriterVtkXmlASCII::getInstance());
+ GbCylinder3DPtr geoSolid(new GbCylinder3D(g_minX1 - 2.0 * dx, g_maxX2 / 2.0, g_maxX3 / 2.0, g_minX1-dx, g_maxX2 / 2.0, g_maxX3 / 2.0, 1.5*D / 2.0));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoSolid.get(), pathname + "/geo/geoSolid", WbWriterVtkXmlASCII::getInstance());
// GbCylinder3DPtr cylinder2(
@@ -384,38 +382,26 @@ void run(string configname)
SPtr<D3Q27Interactor> outflowInt(new D3Q27Interactor(geoOutflow, grid, denBC, Interactor3D::SOLID));
// Create boundary conditions geometry
- GbCuboid3DPtr wallXmin(
- new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_minX1, g_maxX2 + 2.0*dx, g_maxX3));
+ GbCuboid3DPtr wallXmin(new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_minX1, g_maxX2 + 2.0*dx, g_maxX3));
GbSystem3D::writeGeoObject(wallXmin.get(), pathname + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
- GbCuboid3DPtr wallXmax(
- new GbCuboid3D(g_maxX1, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3));
+ GbCuboid3DPtr wallXmax(new GbCuboid3D(g_maxX1, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3));
GbSystem3D::writeGeoObject(wallXmax.get(), pathname + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
- GbCuboid3DPtr wallZmin(
- new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_minX3));
+ GbCuboid3DPtr wallZmin(new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_minX3));
GbSystem3D::writeGeoObject(wallZmin.get(), pathname + "/geo/wallZmin", WbWriterVtkXmlASCII::getInstance());
- GbCuboid3DPtr wallZmax(
- new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_maxX3, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3 + 2.0*dx));
+ GbCuboid3DPtr wallZmax(new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_maxX3, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3 + 2.0*dx));
GbSystem3D::writeGeoObject(wallZmax.get(), pathname + "/geo/wallZmax", WbWriterVtkXmlASCII::getInstance());
- GbCuboid3DPtr wallYmin(
- new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_minX2, g_maxX3));
+ GbCuboid3DPtr wallYmin(new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_minX2, g_maxX3));
GbSystem3D::writeGeoObject(wallYmin.get(), pathname + "/geo/wallYmin", WbWriterVtkXmlASCII::getInstance());
- GbCuboid3DPtr wallYmax(
- new GbCuboid3D(g_minX1 - 2.0*dx, g_maxX2, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3));
+ GbCuboid3DPtr wallYmax(new GbCuboid3D(g_minX1 - 2.0*dx, g_maxX2, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3));
GbSystem3D::writeGeoObject(wallYmax.get(), pathname + "/geo/wallYmax", WbWriterVtkXmlASCII::getInstance());
// Add boundary conditions to grid generator
- SPtr<D3Q27Interactor> wallXminInt(
- new D3Q27Interactor(wallXmin, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallXmaxInt(
- new D3Q27Interactor(wallXmax, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallZminInt(
- new D3Q27Interactor(wallZmin, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallZmaxInt(
- new D3Q27Interactor(wallZmax, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallYminInt(
- new D3Q27Interactor(wallYmin, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallYmaxInt(
- new D3Q27Interactor(wallYmax, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBC, Interactor3D::SOLID));
//cylInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(cylinder1, grid, velBCF1, Interactor3D::SOLID));
//cylInt->addBC(velBCF2);
@@ -423,19 +409,22 @@ void run(string configname)
// Interactor3D::SOLID));
inflowInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, velBCF1, Interactor3D::SOLID));
- inflowInt->addBC(velBCF2);
+ //inflowInt->addBC(velBCF2);
+
+ GbCylinder3DPtr geoAirInflow(new GbCylinder3D(g_minX1 - 2.0*dx, g_maxX2 / 2.0, g_maxX3 / 2.0, g_minX1, g_maxX2 / 2.0, g_maxX3 / 2.0, D * 0.9 / 2.0));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoAirInflow.get(), pathname + "/geo/geoAirInflow", WbWriterVtkXmlASCII::getInstance());
+ SPtr<D3Q27Interactor> inflowAirInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoAirInflow, grid, velBCF2, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> solidInt =
- SPtr<D3Q27Interactor>(new D3Q27Interactor(geoSolid, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> solidInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoSolid, grid, noSlipBC, Interactor3D::SOLID));
InteractorsHelper intHelper(grid, metisVisitor, true);
//intHelper.addInteractor(cylInt);
//intHelper.addInteractor(tubes);
- intHelper.addInteractor(outflowInt);
+ //intHelper.addInteractor(outflowInt);
// intHelper.addInteractor(cyl2Int);
intHelper.addInteractor(wallXminInt);
- //intHelper.addInteractor(wallXmaxInt);
+ intHelper.addInteractor(wallXmaxInt);
intHelper.addInteractor(wallZminInt);
intHelper.addInteractor(wallZmaxInt);
intHelper.addInteractor(wallYminInt);
@@ -490,8 +479,8 @@ void run(string configname)
//fct1.SetExpr("phiL");
//fct1.DefineConst("phiL", phiL);
real x1c = 0; // (g_maxX1 - g_minX1-1)/2; //
- real x2c = (g_maxX2 - g_minX2)/2;
- real x3c = (g_maxX3 - g_minX3)/2;
+ real x2c = (g_minX2 + g_maxX2)/2;
+ real x3c = (g_minX3 + g_maxX3)/2;
mu::Parser fct1;
fct1.SetExpr("0.5-0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
@@ -572,8 +561,8 @@ void run(string configname)
grid->accept(bcVisitor);
- ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
- //TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ //ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
grid->accept(setConnsVisitor);
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
@@ -581,16 +570,17 @@ void run(string configname)
t_ast = 7.19;
t = (int)(t_ast/(U_LB/(D_LB)));
visSch->addSchedule(t,t,t); //t=7.19
- SPtr<WriteMultiphaseQuantitiesSimulationObserver> pp(new WriteMultiphaseQuantitiesSimulationObserver(
- grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
+ // SPtr<WriteMultiphaseQuantitiesSimulationObserver> pp(new WriteMultiphaseQuantitiesSimulationObserver(
+ // grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
+ SPtr<WriteSharpInterfaceQuantitiesSimulationObserver> pp(new WriteSharpInterfaceQuantitiesSimulationObserver(grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
pp->update(0);
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
SPtr<NUPSCounterSimulationObserver> npr(new NUPSCounterSimulationObserver(grid, nupsSch, numOfThreads, comm));
- SPtr<UbScheduler> timeBCSch(new UbScheduler(1, startTime, startTime));
- auto timeDepBC = make_shared<TimeDependentBCSimulationObserver>(TimeDependentBCSimulationObserver(grid, timeBCSch));
- timeDepBC->addInteractor(inflowInt);
+ // SPtr<UbScheduler> timeBCSch(new UbScheduler(1, startTime, startTime));
+ // auto timeDepBC = make_shared<TimeDependentBCSimulationObserver>(TimeDependentBCSimulationObserver(grid, timeBCSch));
+ // timeDepBC->addInteractor(inflowInt);
#ifdef _OPENMP
omp_set_num_threads(numOfThreads);
@@ -600,7 +590,7 @@ void run(string configname)
SPtr<Simulation> simulation(new Simulation(grid, stepGhostLayer, endTime));
simulation->addSimulationObserver(npr);
simulation->addSimulationObserver(pp);
- simulation->addSimulationObserver(timeDepBC);
+ // simulation->addSimulationObserver(timeDepBC);
simulation->addSimulationObserver(rcp);
if (myid == 0)
diff --git a/apps/cpu/RisingBubble2D/RisingBubble2D.cpp b/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
index c9a28efc4..0fb06bc25 100644
--- a/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
+++ b/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
@@ -58,14 +58,15 @@ void run(string configname)
if (myid == 0) UBLOG(logINFO, "2D Rising Bubble: Start!");
if (logToFile) {
-#if defined(__unix__)
- if (myid == 0) {
- const char *str = pathname.c_str();
- mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- }
-#endif
+// #if defined(__unix__)
+// if (myid == 0) {
+// const char *str = pathname.c_str();
+// mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+// }
+// #endif
if (myid == 0) {
+ UbSystem::makeDirectory(pathname);
stringstream logFilename;
logFilename << pathname + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt";
UbLog::output_policy::setStream(logFilename.str());
@@ -147,8 +148,8 @@ void run(string configname)
// kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
// kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
//kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
- //kernel = make_shared<MultiphaseScaleDistributionLBMKernel>();
- kernel = make_shared<MultiphaseSharpInterfaceLBMKernel>();
+ kernel = make_shared<MultiphaseScaleDistributionLBMKernel>();
+ //kernel = make_shared<MultiphaseSharpInterfaceLBMKernel>();
mu::Parser fgr;
fgr.SetExpr("-rho*g_y");
fgr.DefineConst("g_y", g_y);
@@ -187,7 +188,7 @@ void run(string configname)
SPtr<Grid3D> grid(new Grid3D(comm));
grid->setDeltaX(dx);
grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
- grid->setPeriodicX1(false);
+ grid->setPeriodicX1(true);
grid->setPeriodicX2(false);
grid->setPeriodicX3(true);
grid->setGhostLayerWidth(2);
@@ -246,8 +247,8 @@ void run(string configname)
SPtr<WriteBlocksSimulationObserver> ppblocks(new WriteBlocksSimulationObserver(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
InteractorsHelper intHelper(grid, metisVisitor, true);
- intHelper.addInteractor(wallXminInt);
- intHelper.addInteractor(wallXmaxInt);
+ //intHelper.addInteractor(wallXminInt);
+ //intHelper.addInteractor(wallXmaxInt);
intHelper.addInteractor(wallYminInt);
intHelper.addInteractor(wallYmaxInt);
intHelper.selectBlocks();
@@ -360,13 +361,24 @@ void run(string configname)
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
// visSch->addSchedule(307200,307200,307200); //t=2
// visSch->addSchedule(1228185,1228185,1228185);
- double t_ast, t;
- t_ast = 2;
- t = (int)(t_ast / std::sqrt(g_y / D));
+
+ //Tlb = (np.sqrt(2 * dLB / gLB)) / (np.sqrt(2 * dLT / gLT))
+ double gLT = 0.98;
+ double dLT = 0.5;
+ double Tlb = (sqrt(2 * D / g_y)) / (sqrt(2. * dLT / gLT));
+ UBLOG(logINFO, "Tlb = " << Tlb);
+ double t = Tlb * 3.;
visSch->addSchedule(t, t, t); // t=2
- t_ast = 3;
- t = (int)(t_ast / std::sqrt(g_y / D));
- visSch->addSchedule(t, t, t); // t=3
+ UBLOG(logINFO, "T3 = " << t);
+
+
+ // double t_ast, t;
+ // t_ast = 2;
+ // t = (int)(t_ast / std::sqrt(g_y / D));
+ // visSch->addSchedule(t, t, t); // t=2
+ // t_ast = 3;
+ // t = (int)(t_ast / std::sqrt(g_y / D));
+ // visSch->addSchedule(t, t, t); // t=3
// t_ast = 4;
// t = (int)(t_ast/std::sqrt(g_y/D));
// visSch->addSchedule(t,t,t); //t=4
@@ -391,7 +403,10 @@ void run(string configname)
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
SPtr<NUPSCounterSimulationObserver> npr(new NUPSCounterSimulationObserver(grid, nupsSch, numOfThreads, comm));
- // omp_set_num_threads(numOfThreads);
+ omp_set_num_threads(numOfThreads);
+
+ endTime = t + 1000;
+ UBLOG(logINFO, "endTime = " << endTime);
SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
SPtr<Simulation> simulation(new Simulation(grid, stepGhostLayer, endTime));
diff --git a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
index bae43ebf0..169ea1593 100644
--- a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
+++ b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
@@ -427,6 +427,17 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
// (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
// (mfbbc - mfbba));
D3Q27System::calcIncompMacroscopicValues(ff, rhoG, vx, vy, vz);
+ if (withForcing) {
+
+ real forcingX1 = muForcingX1.Eval();
+ real forcingX2 = muForcingX2.Eval();
+ real forcingX3 = muForcingX3.Eval();
+
+ vx += (forcingX1)*deltaT * c1o2;
+ vy += (forcingX2)*deltaT * c1o2;
+ vz += (forcingX3)*deltaT * c1o2;
+ }
+
//if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3)) { }
//else { rhoG = 0.0; vx = 0.0; vy = 0.0; vz = 0.0; }
//// very bad save the world procedure!!!!
@@ -1916,34 +1927,34 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
this->swapDistributions();
- for (int x3 = minX3; x3 < maxX3; x3++) {
- for (int x2 = minX2; x2 < maxX2; x2++) {
- for (int x1 = minX1; x1 < maxX1; x1++) {
+// for (int x3 = minX3; x3 < maxX3; x3++) {
+// for (int x2 = minX2; x2 < maxX2; x2++) {
+// for (int x1 = minX1; x1 < maxX1; x1++) {
- int x1p = x1 + 1;
- int x2p = x2 + 1;
- int x3p = x3 + 1;
- findNeighbors(phaseFieldOld, x1, x2, x3);
-
- //if (((*phaseField)(x1, x2, x3) > c1o2) && (((*phaseFieldOld)(x1, x2, x3) <= c1o2)))
- {//Refill liquid
- real vx;
- real vy;
- real vz;
-
-
- distribution->getDistribution(ff, x1, x2, x3);
- real rhoL;
- D3Q27System::calcIncompMacroscopicValues(ff, rhoL, vx, vy, vz);
- //if (vz != 0) {
-
- // std::cout << "precol: rhoL=" << rhoL << " vx=" << vx << " vy=" << vy << " vz=" << vz << "ffRest=" << ff[DIR_000] << " x=" << x1 << " y=" << x2 << " z=" << x3 << "\n";
- //}
- }
+// int x1p = x1 + 1;
+// int x2p = x2 + 1;
+// int x3p = x3 + 1;
+// findNeighbors(phaseFieldOld, x1, x2, x3);
+
+// //if (((*phaseField)(x1, x2, x3) > c1o2) && (((*phaseFieldOld)(x1, x2, x3) <= c1o2)))
+// {//Refill liquid
+// real vx;
+// real vy;
+// real vz;
+
+
+// distribution->getDistribution(ff, x1, x2, x3);
+// real rhoL;
+// D3Q27System::calcIncompMacroscopicValues(ff, rhoL, vx, vy, vz);
+// //if (vz != 0) {
+
+// // std::cout << "precol: rhoL=" << rhoL << " vx=" << vx << " vy=" << vy << " vz=" << vz << "ffRest=" << ff[DIR_000] << " x=" << x1 << " y=" << x2 << " z=" << x3 << "\n";
+// //}
+// }
- }
- }
-}
+// }
+// }
+// }
diff --git a/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp b/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp
index af5ba7d63..e77ee2043 100644
--- a/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp
+++ b/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp
@@ -350,24 +350,32 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
rhoG = sumRho / sumWeight;// uncheck excpetion: what if there is no adequate neighbor?
for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
if ((phi[fdir] > c1o2) ) {
+ // real fL = distribution->getDistributionInvForDirection(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir], D3Q27System::INVDIR[fdir]);
+
+ if ((phi[D3Q27System::INVDIR[fdir]] > c1o2)) {
real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
- vBC = (vBC + vDir) / (c2o1 + vBC - vDir);
- real fL = distribution->getDistributionInvForDirection(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir], D3Q27System::INVDIR[fdir]);
-
- if ((phi[D3Q27System::INVDIR[fdir]] > c1o2)) {
- ///here we need reconstruction from scrach
- real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
- real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoG, (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+ vBC = (vBC + vDir) / (c2o1 + vBC - vDir);
+
+ //27.04.23
+ real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ //real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
+ //real dvDir = (vBC - vIDir)*c1o2;
+ real dvDir = (vBC - vDir);
+
+ LBMReal fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
+ LBMReal fBC = (fGEQ - c3o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
//real fBC = (fL - feqOLD) * (c1o1 / collFactorG - c1o1) / (c1o1 / collFactorL - c1o1) + feqNew;
//real fG = distribution->getDistributionInvForDirection(x1, x2, x3, fdir);
//real fGEQOld = D3Q27System::getIncompFeqForDirection(fdir, (*rhoNode)(x1, x2, x3), vx, vy, vz);
//real fGEQNew = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
- real fBC = (fL - feqOLD) * (c1o1 / collFactorG - c1o1) / (c1o1 / collFactorL - c1o1) + feqNew;// fL -feqOLD + feqNew;
+ //real fBC = (fL - feqOLD) * (c1o1 / collFactorG - c1o1) / (c1o1 / collFactorL - c1o1) + feqNew;// fL -feqOLD + feqNew;
//real fBC = fGG - c6o1 * WEIGTH[fdir] * (vBC);
distribution->setDistributionForDirection(fBC, x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir], D3Q27System::INVDIR[fdir]);
///// other possibility is tor replace the node itself instead of the neighbor (only c1o1 of them is allowed!)
@@ -401,13 +409,21 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
//real dvDir = vBC - vDir;
+ //27.04.23
+ real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ //real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
+ //real dvDir = (vBC - vIDir)*c1o2;
+ real dvDir = (vBC - vDir);
+
vBC = (vBC + vDir) / (c2o1 + vBC - vDir);
real fL = distribution->getDistributionInvForDirection(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir], D3Q27System::INVDIR[fdir]);
real fG = distribution->getDistributionInvForDirection(x1, x2, x3, fdir);
- real fBC = fG - c6o1 * WEIGTH[fdir] * (vBC);
- //real fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
+ //real fBC = fG - c6o1 * WEIGTH[fdir] * (vBC);
+ real fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
// real fBC = (-fGInv + fGInvEQ + fGEQ - c6o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
- //real fBC = (fGEQ - c3o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
+ real fBC = (fGEQ - c3o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
//real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
//real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoG, (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
@@ -482,14 +498,22 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
//real dvDir = vBC - vDir;
+ //27.04.23
+ real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
+ //real dvDir = (vBC - vIDir)*c1o2;
+ real dvDir = (vBC - vDir);
+
vBC = (vBC + vDir) / (c2o1 + vBC - vDir);
real fL = distribution->getDistributionInvForDirection(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir], D3Q27System::INVDIR[fdir]);
real fG = distribution->getDistributionInvForDirection(x1, x2, x3, fdir);
- real fBC = fG - c6o1 * WEIGTH[fdir] * (vBC);
+ //real fBC = fG - c6o1 * WEIGTH[fdir] * (vBC);
- //real fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
+ real fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
// real fBC = (-fGInv + fGInvEQ + fGEQ - c6o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
- //real fBC = (fGEQ - c3o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
+ real fBC = (fGEQ - c3o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
//real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
//real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoG, (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
@@ -577,9 +601,27 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
if (((phi[fdir] <= c1o2) ))//&& (((*phaseField)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir])) > c1o2)))
{
- real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoG, vx, vy, vz);
- real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoL, vx, vy, vz);
- ff[D3Q27System::INVDIR[fdir]] = (ff[D3Q27System::INVDIR[fdir]] - feqOLD) * (c1o1 / collFactorL - c1o1) / (c1o1 / collFactorG - c1o1) + feqNew;
+ real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+ real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+ real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+ real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
+ real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
+ vBC = (vBC + vDir) / (c2o1 + vBC - vDir);
+
+ //27.04.23
+ real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+ real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
+ //real dvDir = (vBC - vIDir)*c1o2;
+ real dvDir = (vBC - vDir);
+
+
+ real fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
+ ff[D3Q27System::INVDIR[fdir]] = (fGEQ - c3o1 * WEIGTH[fdir] * dvDir * (c1o1 / collFactorL - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
+ // real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoG, vx, vy, vz);
+ // real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoL, vx, vy, vz);
+ // ff[D3Q27System::INVDIR[fdir]] = (ff[D3Q27System::INVDIR[fdir]] - feqOLD) * (c1o1 / collFactorL - c1o1) / (c1o1 / collFactorG - c1o1) + feqNew;
distribution->setDistributionForDirection(ff[D3Q27System::INVDIR[fdir]], x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir], D3Q27System::INVDIR[fdir]);
}
}
@@ -918,21 +960,21 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
real Dxz = -c3o1 * collFactorM * mfbab;
real Dyz = -c3o1 * collFactorM * mfabb;
- if (phi[DIR_000] > c1o2) {
- /// QR eddyviscosity:
- real eddyR = -(Dxy * Dxy + Dxz * Dxz + c1o3 * dxux * dxux) * (dxux) - (Dxy * Dxy + Dyz * Dyz + c1o3 * dyuy * dyuy) * dyuy - (Dxz * Dxz + Dyz * Dyz + c1o3 * dzuz * dzuz) * dzuz - c2o1 * Dxy * Dxz * Dyz;
- real eddyQ = Dxy * Dxz + Dxy * Dyz + Dxz * Dyz + c1o2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz);
- real nuEddy = 5.0e1 * (eddyR / (eddyQ + 1e-100)) * (dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi);
- nuEddy = (nuEddy < c1o1 / collFactorM) ? c1o1 / collFactorM : nuEddy;
- collFactorM = c1o1 / nuEddy;
- // collFactorM = c1o1 / (c1o1 / collFactorM +1.e2*nuEddy*(dX1_phi*dX1_phi+dX2_phi*dX2_phi+dX3_phi*dX3_phi));
- collFactorM = (collFactorM < 1.8) ? 1.8 : collFactorM;
- OxyyPxzz = 8.0 * (collFactorM - 2.0) * (OxxPyyPzz * (3.0 * collFactorM - 1.0) - 5.0 * collFactorM) / (8.0 * (5.0 - 2.0 * collFactorM) * collFactorM + OxxPyyPzz * (8.0 + collFactorM * (9.0 * collFactorM - 26.0)));
- OxyyMxzz = 8.0 * (collFactorM - 2.0) * (collFactorM + OxxPyyPzz * (3.0 * collFactorM - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * collFactorM + 9.0 * collFactorM * collFactorM) - 8.0 * collFactorM);
- Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
- A = (4.0 * collFactorM * collFactorM + 2.0 * collFactorM * OxxPyyPzz * (collFactorM - 6.0) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (10.0 - 3.0 * collFactorM) - 4.0)) / ((collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
- BB = (4.0 * collFactorM * OxxPyyPzz * (9.0 * collFactorM - 16.0) - 4.0 * collFactorM * collFactorM - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * collFactorM * (collFactorM - 2.0))) / (3.0 * (collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
- }
+ // if (phi[DIR_000] > c1o2) {
+ // /// QR eddyviscosity:
+ // real eddyR = -(Dxy * Dxy + Dxz * Dxz + c1o3 * dxux * dxux) * (dxux) - (Dxy * Dxy + Dyz * Dyz + c1o3 * dyuy * dyuy) * dyuy - (Dxz * Dxz + Dyz * Dyz + c1o3 * dzuz * dzuz) * dzuz - c2o1 * Dxy * Dxz * Dyz;
+ // real eddyQ = Dxy * Dxz + Dxy * Dyz + Dxz * Dyz + c1o2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz);
+ // real nuEddy = 5.0e1 * (eddyR / (eddyQ + 1e-100)) * (dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi);
+ // nuEddy = (nuEddy < c1o1 / collFactorM) ? c1o1 / collFactorM : nuEddy;
+ // collFactorM = c1o1 / nuEddy;
+ // // collFactorM = c1o1 / (c1o1 / collFactorM +1.e2*nuEddy*(dX1_phi*dX1_phi+dX2_phi*dX2_phi+dX3_phi*dX3_phi));
+ // collFactorM = (collFactorM < 1.8) ? 1.8 : collFactorM;
+ // OxyyPxzz = 8.0 * (collFactorM - 2.0) * (OxxPyyPzz * (3.0 * collFactorM - 1.0) - 5.0 * collFactorM) / (8.0 * (5.0 - 2.0 * collFactorM) * collFactorM + OxxPyyPzz * (8.0 + collFactorM * (9.0 * collFactorM - 26.0)));
+ // OxyyMxzz = 8.0 * (collFactorM - 2.0) * (collFactorM + OxxPyyPzz * (3.0 * collFactorM - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * collFactorM + 9.0 * collFactorM * collFactorM) - 8.0 * collFactorM);
+ // Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
+ // A = (4.0 * collFactorM * collFactorM + 2.0 * collFactorM * OxxPyyPzz * (collFactorM - 6.0) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (10.0 - 3.0 * collFactorM) - 4.0)) / ((collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
+ // BB = (4.0 * collFactorM * OxxPyyPzz * (9.0 * collFactorM - 16.0) - 4.0 * collFactorM * collFactorM - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * collFactorM * (collFactorM - 2.0))) / (3.0 * (collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
+ // }
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// non Newtonian fluid collision factor
--
GitLab