Skip to content
Snippets Groups Projects
Commit ed091655 authored by AMATERASU\geier's avatar AMATERASU\geier
Browse files

Two phase fields implemented but no energy equation.

parent 254f9664
No related branches found
No related tags found
2 merge requests!171Newest Update,!83Fix MPICommunicator
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1397 additions and 215 deletions
apps/cpu/Multiphase/Multiphase.cpp
+ 1
− 1
View file @ ed091655
......@@ -137,7 +137,7 @@ void run(string configname)
fctF2.SetExpr("vy1");
fctF2.DefineConst("vy1", uLB);
double startTime = 500;
double startTime = 30;
SPtr<BCAdapter> velBCAdapterF1(new MultiphaseVelocityBCAdapter(true, false, false, fctF1, phiH, 0.0, startTime));
SPtr<BCAdapter> velBCAdapterF2(new MultiphaseVelocityBCAdapter(true, false, false, fctF2, phiH, startTime, endTime));
......
apps/cpu/Multiphase/MultiphaseGeier.cfg 0 → 100644
+ 60
− 0
View file @ ed091655
#pathname = E:/Multiphase/HesamCodeWithCumulantsDensRatio
#pathname = E:/Multiphase/HesamCodeWithCumulantsQuartic
#pathname = E:/Multiphase/HesamCode
pathname = E:/Multiphase/HesamCodeCumulantTubeFilter
pathGeo = C:/Users/geier/Documents/VirtualFluids_dev_Kostya/apps/cpu/Multiphase/backup
geoFile=tubeTransformed.stl
#geoFile = JetBreakup2.ASCII.stl
numOfThreads = 4
availMem = 10e9
#Grid
#boundingBox = -1.0 121.0 0.5 629.0 -1.0 121.0 #(Jet Breakup) (Original with inlet length)
#boundingBox = -60.5 60.5 -1.0 -201.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
#blocknx = 22 20 22
#boundingBox = -60.5 60.5 -1.0 -21.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
#boundingBox = -60.5 60.5 -21.0 -1.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
#blocknx = 22 20 22
#dx = 0.5
#boundingBox = 6.0e-3 46.0e-3 -5e-3 5e-3 -5e-3 5e-3
boundingBox = 6.0e-3 16.0e-3 -5e-3 5e-3 -5e-3 5e-3
blocknx = 60 60 60 #20 20 20
dx = 1.66666666667e-4
refineLevel = 0
#Simulation
uLB =0.005# 0.0000005 #inlet velocity
uF2 = 0.0001
Re = 10
nuL =1e-6#1e-2# 1.0e-5 #!1e-2
nuG =1e-6#1e-2# 1.16e-4 #!1e-2
densityRatio = 10000#1000#1000 #30
sigma =0# 1e-4 #4.66e-3 #surface tension 1e-4 ./. 1e-5
interfaceThickness = 5
radius = 615.0 (Jet Breakup)
contactAngle = 110.0
gravity = 0.0
#gravity = -5.04e-6
phi_L = 0.0
phi_H = 1.0
Phase-field Relaxation = 0.6
Mobility = 0.1 #0.02 # 0.01 ./. 0.08, fine correction of Phase-field Relaxation parameter, to activate it need to change in kernel tauH to tauH1
logToFile = false
newStart = true
restartStep = 100000
cpStart = 100000
cpStep = 100000
outTime = 100
endTime = 200000000
\ No newline at end of file
src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
+ 50
− 2
View file @ ed091655
......@@ -155,6 +155,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
datanames.push_back("Vy");
datanames.push_back("Vz");
datanames.push_back("P1");
datanames.push_back("Phi2");
data.resize(datanames.size());
......@@ -162,10 +163,12 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
SPtr<DistributionArray3D> distributionsF = kernel->getDataSet()->getFdistributions();
SPtr<DistributionArray3D> distributionsH = kernel->getDataSet()->getHdistributions();
SPtr<DistributionArray3D> distributionsH2 = kernel->getDataSet()->getH2distributions();
SPtr<PhaseFieldArray3D> divU = kernel->getDataSet()->getPhaseField();
LBMReal f[D3Q27System::ENDF + 1];
LBMReal phi[D3Q27System::ENDF + 1];
LBMReal phi2[D3Q27System::ENDF + 1];
LBMReal vx1, vx2, vx3, rho, p1, beta, kappa;
LBMReal densityRatio = kernel->getDensityRatio();
......@@ -202,6 +205,8 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3, -1);
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField(
new CbArray3D<LBMReal, IndexerX3X2X1>(maxX1, maxX2, maxX3, -999.0));
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField2(
new CbArray3D<LBMReal, IndexerX3X2X1>(maxX1, maxX2, maxX3, -999.0));
for (int ix3 = minX3; ix3 < maxX3; ix3++) {
for (int ix2 = minX2; ix2 < maxX2; ix2++) {
......@@ -211,8 +216,17 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
(*phaseField)(ix1, ix2, ix3) =
((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
(((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
if (distributionsH2) {
distributionsH2->getDistribution(f, ix1, ix2, ix3);
(*phaseField2)(ix1, ix2, ix3) =
((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
(((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
}
else { (*phaseField2)(ix1, ix2, ix3) = 999.0; }
}
}
}
......@@ -244,6 +258,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
float(worldCoordinates[2])));
phi[REST] = (*phaseField)(ix1, ix2, ix3);
phi2[REST] = (*phaseField2)(ix1, ix2, ix3);
if ((ix1 == 0) || (ix2 == 0) || (ix3 == 0)) {
dX1_phi = 0.0;
......@@ -284,6 +299,38 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
dX2_phi = 0.0 * gradX2_phi(phi);
dX3_phi = 0.0 * gradX3_phi(phi);
mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi(phi);
//phi2[E] = (*phaseField2)(ix1 + DX1[E], ix2 + DX2[E], ix3 + DX3[E]);
//phi2[N] = (*phaseField2)(ix1 + DX1[N], ix2 + DX2[N], ix3 + DX3[N]);
//phi2[T] = (*phaseField2)(ix1 + DX1[T], ix2 + DX2[T], ix3 + DX3[T]);
//phi2[W] = (*phaseField2)(ix1 + DX1[W], ix2 + DX2[W], ix3 + DX3[W]);
//phi2[S] = (*phaseField2)(ix1 + DX1[S], ix2 + DX2[S], ix3 + DX3[S]);
//phi2[B] = (*phaseField2)(ix1 + DX1[B], ix2 + DX2[B], ix3 + DX3[B]);
//phi2[NE] = (*phaseField2)(ix1 + DX1[NE], ix2 + DX2[NE], ix3 + DX3[NE]);
//phi2[NW] = (*phaseField2)(ix1 + DX1[NW], ix2 + DX2[NW], ix3 + DX3[NW]);
//phi2[TE] = (*phaseField2)(ix1 + DX1[TE], ix2 + DX2[TE], ix3 + DX3[TE]);
//phi2[TW] = (*phaseField2)(ix1 + DX1[TW], ix2 + DX2[TW], ix3 + DX3[TW]);
//phi2[TN] = (*phaseField2)(ix1 + DX1[TN], ix2 + DX2[TN], ix3 + DX3[TN]);
//phi2[TS] = (*phaseField2)(ix1 + DX1[TS], ix2 + DX2[TS], ix3 + DX3[TS]);
//phi2[SW] = (*phaseField2)(ix1 + DX1[SW], ix2 + DX2[SW], ix3 + DX3[SW]);
//phi2[SE] = (*phaseField2)(ix1 + DX1[SE], ix2 + DX2[SE], ix3 + DX3[SE]);
//phi2[BW] = (*phaseField2)(ix1 + DX1[BW], ix2 + DX2[BW], ix3 + DX3[BW]);
//phi2[BE] = (*phaseField2)(ix1 + DX1[BE], ix2 + DX2[BE], ix3 + DX3[BE]);
//phi2[BS] = (*phaseField2)(ix1 + DX1[BS], ix2 + DX2[BS], ix3 + DX3[BS]);
//phi2[BN] = (*phaseField2)(ix1 + DX1[BN], ix2 + DX2[BN], ix3 + DX3[BN]);
//phi2[BSW] = (*phaseField2)(ix1 + DX1[BSW], ix2 + DX2[BSW], ix3 + DX3[BSW]);
//phi2[BSE] = (*phaseField2)(ix1 + DX1[BSE], ix2 + DX2[BSE], ix3 + DX3[BSE]);
//phi2[BNW] = (*phaseField2)(ix1 + DX1[BNW], ix2 + DX2[BNW], ix3 + DX3[BNW]);
//phi2[BNE] = (*phaseField2)(ix1 + DX1[BNE], ix2 + DX2[BNE], ix3 + DX3[BNE]);
//phi2[TNE] = (*phaseField2)(ix1 + DX1[TNE], ix2 + DX2[TNE], ix3 + DX3[TNE]);
//phi2[TNW] = (*phaseField2)(ix1 + DX1[TNW], ix2 + DX2[TNW], ix3 + DX3[TNW]);
//phi2[TSE] = (*phaseField2)(ix1 + DX1[TSE], ix2 + DX2[TSE], ix3 + DX3[TSE]);
//phi2[TSW] = (*phaseField2)(ix1 + DX1[TSW], ix2 + DX2[TSW], ix3 + DX3[TSW]);
// mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi(phi);
}
distributionsF->getDistribution(f, ix1, ix2, ix3);
......@@ -361,6 +408,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
data[index++].push_back(vx2);
data[index++].push_back(vx3);
data[index++].push_back(p1);
data[index++].push_back(phi2[REST]);
}
}
}
......
src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.h
+ 5
− 0
View file @ ed091655
......@@ -80,18 +80,23 @@ protected:
LBMReal h2[D3Q27System::ENDF + 1];
LBMReal g [D3Q27System::ENDF+1];
LBMReal phi[D3Q27System::ENDF+1];
LBMReal phi2[D3Q27System::ENDF + 1];
LBMReal pr1[D3Q27System::ENDF+1];
LBMReal phi_cutoff[D3Q27System::ENDF+1];
LBMReal gradX1_phi();
LBMReal gradX2_phi();
LBMReal gradX3_phi();
LBMReal gradX1_phi2();
LBMReal gradX2_phi2();
LBMReal gradX3_phi2();
//LBMReal gradX1_pr1();
//LBMReal gradX2_pr1();
//LBMReal gradX3_pr1();
//LBMReal dirgradC_phi(int n, int k);
void computePhasefield();
void findNeighbors(CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr ph /*Phase-Field*/, int x1, int x2, int x3);
void findNeighbors2(CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr ph, int x1, int x2, int x3);
//void findNeighbors(CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr ph /*Phase-Field*/, CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr pf /*Pressure-Field*/, int x1, int x2, int x3);
//void pressureFiltering(CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr pf /*Pressure-Field*/, CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr pf_filtered /*Pressure-Field*/);
......
src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
+ 30
− 1
View file @ ed091655
......@@ -224,7 +224,7 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
feq[dir] = rho*WEIGTH[dir]*(1 + 3*velProd + 4.5*velSq1 - 1.5*(vx1Sq+vx2Sq+vx3Sq));
//geq[dir] = p1*WEIGTH1[dir] + gamma;
geq[dir] = p1*WEIGTH[dir]/(rho*UbMath::c1o3) + gamma;
geq[dir] = p1*WEIGTH[dir]/(rho*UbMath::c1o3) + gamma*rho;
}
......@@ -291,6 +291,35 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
distributionsH->setDistributionInv(f, ix1, ix2, ix3);
if (distributionsH2) {
f[E] = (1.-phi) * feq[E] / rho;
f[W] = (1.-phi) * feq[W] / rho;
f[N] = (1.-phi) * feq[N] / rho;
f[S] = (1.-phi) * feq[S] / rho;
f[T] = (1.-phi) * feq[T] / rho;
f[B] = (1.-phi) * feq[B] / rho;
f[NE] = (1.-phi) * feq[NE] / rho;
f[SW] = (1.-phi) * feq[SW] / rho;
f[SE] = (1.-phi) * feq[SE] / rho;
f[NW] = (1.-phi) * feq[NW] / rho;
f[TE] = (1.-phi) * feq[TE] / rho;
f[BW] = (1.-phi) * feq[BW] / rho;
f[BE] = (1.-phi) * feq[BE] / rho;
f[TW] = (1.-phi) * feq[TW] / rho;
f[TN] = (1.-phi) * feq[TN] / rho;
f[BS] = (1.-phi) * feq[BS] / rho;
f[BN] = (1.-phi) * feq[BN] / rho;
f[TS] = (1.-phi) * feq[TS] / rho;
f[TNE] = (1.-phi) * feq[TNE] / rho;
f[TNW] = (1.-phi) * feq[TNW] / rho;
f[TSE] = (1.-phi) * feq[TSE] / rho;
f[TSW] = (1.-phi) * feq[TSW] / rho;
f[BNE] = (1.-phi) * feq[BNE] / rho;
f[BNW] = (1.-phi) * feq[BNW] / rho;
f[BSE] = (1.-phi) * feq[BSE] / rho;
f[BSW] = (1.-phi) * feq[BSW] / rho;
f[REST] = (1.-phi) * feq[REST] / rho;
distributionsH2->setDistribution(f, ix1, ix2, ix3);
distributionsH2->setDistributionInv(f, ix1, ix2, ix3);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment