Skip to content
Snippets Groups Projects
Commit 09dfa160 authored by kutscher's avatar kutscher
Browse files

fix RheologyK17LBMKernel

parent a5cfa683
No related branches found
No related tags found
1 merge request!18Feature/thixotropy
Expand all Hide whitespace changes
Inline Side-by-side
apps/cpu/HerschelBulkleySphere/hbsphere.cpp
+ 11
− 5
View file @ 09dfa160
......@@ -101,6 +101,7 @@ void bflow(string configname)
thix->setOmegaMin(omegaMin);
SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
//noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new HerschelBulkleyModelNoSlipBCAlgorithm()));
//noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelNoSlipBCAlgorithm()));
......@@ -117,9 +118,12 @@ void bflow(string configname)
//fct.DefineConst("H", H);
SPtr<BCAdapter> velocityBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new SimpleVelocityBCAlgorithm()));
//velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
SPtr<BCAdapter> densityBCAdapter(new DensityBCAdapter());
densityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
//densityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonReflectingOutflowBCAlgorithm()));
//BS visitor
BoundaryConditionsBlockVisitor bcVisitor;
......@@ -131,6 +135,8 @@ void bflow(string configname)
SPtr<BCProcessor> bcProc;
bcProc = SPtr<BCProcessor>(new BCProcessor());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CumulantLBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulant4thOrderViscosityLBMKernel());
SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
......@@ -140,8 +146,8 @@ void bflow(string configname)
SPtr<Grid3D> grid(new Grid3D(comm));
grid->setPeriodicX1(false);
grid->setPeriodicX2(false);
grid->setPeriodicX3(false);
grid->setPeriodicX2(true);
grid->setPeriodicX3(true);
grid->setDeltaX(deltax);
grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
......@@ -221,7 +227,7 @@ void bflow(string configname)
//wall interactors
SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, slipBCAdapter, Interactor3D::SOLID));
SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, slipBCAdapter, Interactor3D::SOLID));
SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, slipBCAdapter, Interactor3D::SOLID));
SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, slipBCAdapter, Interactor3D::SOLID));
......@@ -237,10 +243,10 @@ void bflow(string configname)
InteractorsHelper intHelper(grid, metisVisitor);
intHelper.addInteractor(wallXminInt);
intHelper.addInteractor(wallXmaxInt);
intHelper.addInteractor(wallZminInt);
intHelper.addInteractor(wallZmaxInt);
intHelper.addInteractor(wallYminInt);
intHelper.addInteractor(wallYmaxInt);
intHelper.addInteractor(wallZminInt);
intHelper.addInteractor(wallZmaxInt);
intHelper.addInteractor(sphereInt);
intHelper.selectBlocks();
if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
......
src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.cpp
+ 0
− 2
View file @ 09dfa160
......@@ -48,8 +48,6 @@ void SimpleSlipBCAlgorithm::applyBC()
{
//quadratic bounce back
const int invDir = D3Q27System::INVDIR[fdir];
LBMReal q = bcPtr->getQ(invDir);// m+m q=0 stabiler
//vx3=0;
LBMReal velocity = 0.0;
switch (invDir)
{
......
src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp
+ 284
− 228
View file @ 09dfa160
This diff is collapsed.
src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
+ 50
− 26
View file @ 09dfa160
......@@ -107,24 +107,25 @@ void RheologyK17LBMKernel::calculate(int step)
int maxX2 = bcArrayMaxX2-ghostLayerWidth;
int maxX3 = bcArrayMaxX3-ghostLayerWidth;
LBMReal omega = collFactor;
//LBMReal omega = collFactor;
//magic parameter for rheology
LBMReal a = 10;
OxxPyyPzz = c1 / (a * ((c1 / omega) - c1o2) + c1o2);
////magic parameter for rheology
//LBMReal a = 10;
//OxxPyyPzz = c1 / (a * ((c1 / omega) - c1o2) + c1o2);
//OxxPyyPzz = (OxxPyyPzz > c1) ? c1 : OxxPyyPzz;
//LBMReal OxyyPxzz = eight*(-two+omega)*(one+two*omega)/(-eight-fourteen*omega+seven*omega*omega);//one;
//LBMReal OxyyMxzz = eight*(-two+omega)*(-seven+four*omega)/(fiftysix-fifty*omega+nine*omega*omega);//one;
//LBMReal Oxyz = twentyfour*(-two+omega)*(-two-seven*omega+three*omega*omega)/(fourtyeight+c152*omega-c130*omega*omega+twentynine*omega*omega*omega);
LBMReal OxyyPxzz = 8.0*(omega-2.0)*(OxxPyyPzz*(3.0*omega-1.0)-5.0*omega)/(8.0*(5.0-2.0*omega)*omega+OxxPyyPzz*(8.0+omega*(9.0*omega-26.0)));
LBMReal OxyyMxzz = 8.0*(omega-2.0)*(omega+OxxPyyPzz*(3.0*omega-7.0))/(OxxPyyPzz*(56.0-42.0*omega+9.0*omega*omega)-8.0*omega);
LBMReal Oxyz = 24.0*(omega-2.0)*(4.0*omega*omega+omega*OxxPyyPzz*(18.0-13.0*omega)+OxxPyyPzz*OxxPyyPzz*(2.0+omega*(6.0*omega-11.0)))/(16.0*omega*omega*(omega-6.0)-2.0*omega*OxxPyyPzz*(216.0+5.0*omega*(9.0*omega-46.0))+OxxPyyPzz*OxxPyyPzz*(omega*(3.0*omega-10.0)*(15.0*omega-28.0)-48.0));
////LBMReal OxyyPxzz = eight*(-two+omega)*(one+two*omega)/(-eight-fourteen*omega+seven*omega*omega);//one;
////LBMReal OxyyMxzz = eight*(-two+omega)*(-seven+four*omega)/(fiftysix-fifty*omega+nine*omega*omega);//one;
////LBMReal Oxyz = twentyfour*(-two+omega)*(-two-seven*omega+three*omega*omega)/(fourtyeight+c152*omega-c130*omega*omega+twentynine*omega*omega*omega);
//LBMReal OxyyPxzz = 8.0*(omega-2.0)*(OxxPyyPzz*(3.0*omega-1.0)-5.0*omega)/(8.0*(5.0-2.0*omega)*omega+OxxPyyPzz*(8.0+omega*(9.0*omega-26.0)));
//LBMReal OxyyMxzz = 8.0*(omega-2.0)*(omega+OxxPyyPzz*(3.0*omega-7.0))/(OxxPyyPzz*(56.0-42.0*omega+9.0*omega*omega)-8.0*omega);
//LBMReal Oxyz = 24.0*(omega-2.0)*(4.0*omega*omega+omega*OxxPyyPzz*(18.0-13.0*omega)+OxxPyyPzz*OxxPyyPzz*(2.0+omega*(6.0*omega-11.0)))/(16.0*omega*omega*(omega-6.0)-2.0*omega*OxxPyyPzz*(216.0+5.0*omega*(9.0*omega-46.0))+OxxPyyPzz*OxxPyyPzz*(omega*(3.0*omega-10.0)*(15.0*omega-28.0)-48.0));
//LBMReal A = (four + two*omega - three*omega*omega) / (two - seven*omega + five*omega*omega);
//LBMReal B = (four + twentyeight*omega - fourteen*omega*omega) / (six - twentyone*omega + fiveteen*omega*omega);
////LBMReal A = (four + two*omega - three*omega*omega) / (two - seven*omega + five*omega*omega);
////LBMReal B = (four + twentyeight*omega - fourteen*omega*omega) / (six - twentyone*omega + fiveteen*omega*omega);
LBMReal A = (4.0*omega*omega+2.0*omega*OxxPyyPzz*(omega-6.0)+OxxPyyPzz*OxxPyyPzz*(omega*(10.0-3.0*omega)-4.0))/((omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
LBMReal B = (4.0*omega*OxxPyyPzz*(9.0*omega-16.0)-4.0*omega*omega-2.0*OxxPyyPzz*OxxPyyPzz*(2.0+9.0*omega*(omega-2.0)))/(3.0*(omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
//LBMReal A = (4.0*omega*omega+2.0*omega*OxxPyyPzz*(omega-6.0)+OxxPyyPzz*OxxPyyPzz*(omega*(10.0-3.0*omega)-4.0))/((omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
//LBMReal B = (4.0*omega*OxxPyyPzz*(9.0*omega-16.0)-4.0*omega*omega-2.0*OxxPyyPzz*OxxPyyPzz*(2.0+9.0*omega*(omega-2.0)))/(3.0*(omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
for (int x3 = minX3; x3 < maxX3; x3++)
{
......@@ -206,7 +207,7 @@ void RheologyK17LBMKernel::calculate(int step)
(mfbbc-mfbba))/rho;
////////////////////////////////////////////////////////////////////////////////////
LBMReal collFactorF = collFactor;
LBMReal omega = collFactor;
//forcing
///////////////////////////////////////////////////////////////////////////////////////////
......@@ -605,7 +606,7 @@ void RheologyK17LBMKernel::calculate(int step)
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//incl. correction (hat noch nicht so gut funktioniert...Optimierungsbedarf??)
LBMReal dxux = c1o2 * (-omega) *(mxxMyy+mxxMzz)+c1o2 * OxxPyyPzz * (mfaaa-mxxPyyPzz);
LBMReal dxux = c1o2 * (-omega) * (mxxMyy + mxxMzz);// +c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
LBMReal dyuy = dxux+omega * c3o2 * mxxMyy;
LBMReal dzuz = dxux+omega * c3o2 * mxxMzz;
......@@ -615,19 +616,42 @@ void RheologyK17LBMKernel::calculate(int step)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//non Newtonian fluid collision factor
LBMReal shearRate = sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
//collFactorF = getThyxotropyCollFactor(collFactorF, shearRate, rho);
collFactorF = Thixotropy::getHerschelBulkleyCollFactor(collFactorF, shearRate, drho);
LBMReal shearRate = sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (drho + c1);
//omega = getThyxotropyCollFactor(omega, shearRate, rho);
omega = Thixotropy::getHerschelBulkleyCollFactor(omega, shearRate, drho);
//omega = Thixotropy::getBinghamCollFactor(omega, shearRate, drho);
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
mxxMyy += omega * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
mxxMzz += omega * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
mfabb += omega * (-mfabb);
mfbab += omega * (-mfbab);
mfbba += omega * (-mfbba);
if(omega < c1) { omega = c1; } //arbitrary limit (24.09.2020)
//magic parameter for rheology
LBMReal a = 10;
OxxPyyPzz = c1 / (a * ((c1 / omega) - c1o2) + c1o2);
OxxPyyPzz = (OxxPyyPzz > c1) ? c1 : OxxPyyPzz;
LBMReal OxyyPxzz = 8.0 * (omega - 2.0) * (OxxPyyPzz * (3.0 * omega - 1.0) - 5.0 * omega) / (8.0 * (5.0 - 2.0 * omega) * omega + OxxPyyPzz * (8.0 + omega * (9.0 * omega - 26.0)));
LBMReal OxyyMxzz = 8.0 * (omega - 2.0) * (omega + OxxPyyPzz * (3.0 * omega - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * omega + 9.0 * omega * omega) - 8.0 * omega);
LBMReal Oxyz = 24.0 * (omega - 2.0) * (4.0 * omega * omega + omega * OxxPyyPzz * (18.0 - 13.0 * omega) + OxxPyyPzz * OxxPyyPzz * (2.0 + omega * (6.0 * omega - 11.0))) / (16.0 * omega * omega * (omega - 6.0) - 2.0 * omega * OxxPyyPzz * (216.0 + 5.0 * omega * (9.0 * omega - 46.0)) + OxxPyyPzz * OxxPyyPzz * (omega * (3.0 * omega - 10.0) * (15.0 * omega - 28.0) - 48.0));
LBMReal A = (4.0 * omega * omega + 2.0 * omega * OxxPyyPzz * (omega - 6.0) + OxxPyyPzz * OxxPyyPzz * (omega * (10.0 - 3.0 * omega) - 4.0)) / ((omega - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * omega) - 8.0 * omega));
LBMReal B = (4.0 * omega * OxxPyyPzz * (9.0 * omega - 16.0) - 4.0 * omega * omega - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * omega * (omega - 2.0))) / (3.0 * (omega - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * omega) - 8.0 * omega));
//relax
wadjust = OxxPyyPzz+(one-OxxPyyPzz)*fabs((mfaaa-mxxPyyPzz))/(fabs((mfaaa-mxxPyyPzz))+qudricLimitD);
mxxPyyPzz += wadjust*(mfaaa-mxxPyyPzz)-three * (one-c1o2 * OxxPyyPzz) * (vx2 * dxux+vy2 * dyuy+vz2 * dzuz);
//wadjust = OxxPyyPzz+(one-OxxPyyPzz)*fabs((mfaaa-mxxPyyPzz))/(fabs((mfaaa-mxxPyyPzz))+qudricLimitD);
//mxxPyyPzz += wadjust*(mfaaa-mxxPyyPzz)-three * (one-c1o2 * OxxPyyPzz) * (vx2 * dxux+vy2 * dyuy+vz2 * dzuz);
mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) - three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
// mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz)-three * (one-c1o2 * OxxPyyPzz) * (vx2 * dxux+vy2 * dyuy+vz2 * dzuz);//-magicBulk*OxxPyyPzz;
mxxMyy += omega * (-mxxMyy)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vy2 * dyuy);
mxxMzz += omega * (-mxxMzz)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vz2 * dzuz);
//mxxMyy += omega * (-mxxMyy)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vy2 * dyuy);
//mxxMzz += omega * (-mxxMzz)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vz2 * dzuz);
//////////////////////////////////////////////////////////////////////////
//limiter-Scheise Teil 2
......@@ -644,9 +668,9 @@ void RheologyK17LBMKernel::calculate(int step)
//mxxMyy += -(-omega) * (-mxxMyy);
//mxxMzz += -(-omega) * (-mxxMzz);
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
mfabb += omega * (-mfabb);
mfbab += omega * (-mfbab);
mfbba += omega * (-mfbba);
//mfabb += omega * (-mfabb);
//mfbab += omega * (-mfbab);
//mfbba += omega * (-mfbba);
//////////////////////////////////////////////////////////////////////////
//limiter-Scheise Teil 3
......
src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h
+ 1
− 1
View file @ 09dfa160
......@@ -8,7 +8,7 @@
#include "basics/container/CbArray4D.h"
#include "basics/container/CbArray3D.h"
//! \brief compressible cumulant LBM kernel.
//! \brief compressible cumulant LBM kernel with rheological properties of shear and bulk viscosity for non-Newtonian fluids.
//! \details CFD solver that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model
//! \author K. Kutscher, M. Geier
class RheologyK17LBMKernel : public LBMKernel
......
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