From d7672e5f4ab420819ad353d7b084141537d0982e Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Tue, 29 Dec 2020 13:49:44 +0100
Subject: [PATCH] add velocity BC with interpolation / add new K17 model for
Bingham fluid
---
apps/cpu/HerschelBulkleyModel/hbf1.cfg | 8 ++--
apps/cpu/HerschelBulkleyModel/hbflow.cpp | 25 ++++++++--
apps/cpu/rheometer/rheometer.cfg | 8 ++--
apps/cpu/rheometer/rheometer.cpp | 31 +++++++-----
src/cpu/VirtualFluids.h | 1 +
.../BoundaryConditions/BCAlgorithm.h | 1 +
.../BinghamModelVelocityBCAlgorithm.h | 27 +++++++++++
.../RheologicalVelocityBCAlgorithm.cpp | 47 +++++++++++++++++++
.../RheologicalVelocityBCAlgorithm.h | 22 +++++++++
.../ThixotropyNoSlipBCAlgorithm.cpp | 6 +--
.../LBM/BinghamModelLBMKernel.h | 4 +-
.../LBM/RheologyK17LBMKernel.cpp | 6 ++-
src/cpu/VirtualFluidsCore/LBM/Thixotropy.h | 20 ++++----
.../LBM/ThixotropyModelLBMKernel2.cpp | 25 +++++++---
14 files changed, 184 insertions(+), 47 deletions(-)
create mode 100644 src/cpu/VirtualFluidsCore/BoundaryConditions/BinghamModelVelocityBCAlgorithm.h
create mode 100644 src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
create mode 100644 src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h
diff --git a/apps/cpu/HerschelBulkleyModel/hbf1.cfg b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
index 2cf4375f0..b5f8ef94d 100644
--- a/apps/cpu/HerschelBulkleyModel/hbf1.cfg
+++ b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
@@ -1,4 +1,4 @@
-pathname = d:/temp/BinghamNewOmegaTest
+pathname = d:/temp/BinghamNewConvergenceAcousticScaling_it1/128
numOfThreads = 4
availMem = 8e9
@@ -19,7 +19,7 @@ boundingBox = 2 32 2
deltax = 1
#nuLB = 0.05
-forcing = 8e-7
+forcing = 8e-6
#n = 0.8
#tau0 = 3e-6
@@ -28,5 +28,5 @@ n = 0.4
Re = 1
Bn = 0.01
-outTime = 1000
-endTime = 30000
\ No newline at end of file
+outTime = 10000
+endTime = 300000
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleyModel/hbflow.cpp b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
index bd4879c96..e0a24aba9 100644
--- a/apps/cpu/HerschelBulkleyModel/hbflow.cpp
+++ b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
@@ -93,8 +93,25 @@ void bflow(string configname)
double U = velocity;
double Gamma = U / d;
- double k = 0.05; // (U * d) / (Re * std::pow(Gamma, n - 1));
- double tau0 = 3e-6;// Bn* k* std::pow(Gamma, n);
+ double scaleFactor = 4.0;
+
+ // Diffusive Scaling
+
+ //double k = 0.005; // (U * d) / (Re * std::pow(Gamma, n - 1));
+ //double tau0 = 3e-5 / (scaleFactor * scaleFactor);// *4.0);//*4.0);// Bn* k* std::pow(Gamma, n);
+ //forcing /= scaleFactor * scaleFactor * scaleFactor;
+ //endTime *= scaleFactor * scaleFactor;
+ //deltax /= scaleFactor;
+
+ // Acoustic Scaling
+
+ double k = 0.005 * scaleFactor;
+ double tau0 = 3e-5;
+ forcing /= scaleFactor;
+ endTime *= scaleFactor;
+ deltax /= scaleFactor;
+
+ outTime = endTime;
double beta = 14;
double c = 10; // 1.0 / 6.0;
@@ -123,8 +140,8 @@ void bflow(string configname)
bcProc = SPtr<BCProcessor>(new BCProcessor());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new PowellEyringModelLBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
- //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
- SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
+ SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
kernel->setForcingX1(forcing);
kernel->setWithForcing(true);
kernel->setBCProcessor(bcProc);
diff --git a/apps/cpu/rheometer/rheometer.cfg b/apps/cpu/rheometer/rheometer.cfg
index 0330c9297..c65d95096 100644
--- a/apps/cpu/rheometer/rheometer.cfg
+++ b/apps/cpu/rheometer/rheometer.cfg
@@ -1,4 +1,4 @@
-outputPath = d:/temp/rheometer/rheometerBinghamNewEquilibriumSBB/rheometerBingham_nu_1.5e-3_nnn
+outputPath = d:/temp/rheometer/rheometerBinghamqQBB/rheometerBingham_nu_1.5e-3_newBC_K17
viscosityPath = d:/Projects/VirtualFluidsCombined/apps/cpu/rheometer
@@ -16,10 +16,10 @@ deltax = 1
refineLevel = 0
OmegaLB = 1e-4
-tau0 = 1e-6
+tau0 = 5e-7
-resolution = 32
-scaleFactor = 1
+resolution = 128
+scaleFactor = 4
newStart = true
restartStep = 100000
diff --git a/apps/cpu/rheometer/rheometer.cpp b/apps/cpu/rheometer/rheometer.cpp
index c7bc4f3e6..47851a944 100644
--- a/apps/cpu/rheometer/rheometer.cpp
+++ b/apps/cpu/rheometer/rheometer.cpp
@@ -63,13 +63,21 @@ void bflow(string configname)
LBMReal rhoLB = 0.0;
- OmegaLB /= scaleFactor;
- tau0 /= scaleFactor;
-
- endTime *= scaleFactor;
- outTime = endTime;
- cpStart = endTime;
- cpStep = endTime;
+ //akoustic
+ OmegaLB /= scaleFactor;
+ nuLB *=scaleFactor;
+ endTime *= scaleFactor;
+ outTime = endTime;
+ cpStart = endTime;
+ cpStep = endTime;
+
+//diffusive
+ //OmegaLB /= scaleFactor * scaleFactor;
+ //tau0 /= scaleFactor * scaleFactor;
+ //endTime *= scaleFactor * scaleFactor;
+ //outTime = endTime;
+ //cpStart = endTime;
+ //cpStep = endTime;
SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
// double uWorld = (N * PI) / 30.0; //0.0037699111843
@@ -149,8 +157,9 @@ void bflow(string configname)
SPtr<BCAdapter> velocityBCAdapter(new VelocityBCAdapter(true, true, true, fctVx, fctVy, fctVz, 0, BCFunction::INFCONST));
//velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
- velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new SimpleVelocityBCAlgorithm()));
+ //velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new SimpleVelocityBCAlgorithm()));
//velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
+ velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelVelocityBCAlgorithm()));
//SPtr<BCAdapter> densityBCAdapter(new DensityBCAdapter());
//densityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
@@ -159,7 +168,7 @@ void bflow(string configname)
//BS visitor
BoundaryConditionsBlockVisitor bcVisitor;
- //bcVisitor.addBC(noSlipBCAdapter);
+ bcVisitor.addBC(noSlipBCAdapter);
//bcVisitor.addBC(slipBCAdapter);
bcVisitor.addBC(velocityBCAdapter);
//bcVisitor.addBC(densityBCAdapter);
@@ -169,9 +178,9 @@ void bflow(string configname)
//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 RheologyK17LBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
- SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
kernel->setBCProcessor(bcProc);
//kernel->setForcingX1(forcing);
//kernel->setWithForcing(true);
diff --git a/src/cpu/VirtualFluids.h b/src/cpu/VirtualFluids.h
index e02a5cdf1..e033e94a9 100644
--- a/src/cpu/VirtualFluids.h
+++ b/src/cpu/VirtualFluids.h
@@ -132,6 +132,7 @@
#include <BoundaryConditions/HerschelBulkleyModelNoSlipBCAlgorithm.h>
#include <BoundaryConditions/SimpleSlipBCAlgorithm.h>
#include <BoundaryConditions/PowellEyringModelNoSlipBCAlgorithm.h>
+#include <BoundaryConditions/BinghamModelVelocityBCAlgorithm.h>
#include <Connectors/Block3DConnector.h>
#include <Connectors/Block3DConnectorFactory.h>
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
index f9427d41a..3adbed448 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
@@ -66,6 +66,7 @@ public:
static const char SimpleVelocityBCAlgorithm = 16;
static const char SimpleSlipBCAlgorithm = 17;
static const char PowellEyringModelNoSlipBCAlgorithm = 18;
+ static const char BinghamModelVelocityBCAlgorithm = 19;
public:
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BinghamModelVelocityBCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/BinghamModelVelocityBCAlgorithm.h
new file mode 100644
index 000000000..0155babfc
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BinghamModelVelocityBCAlgorithm.h
@@ -0,0 +1,27 @@
+#ifndef BinghamModelVelocityBCAlgorithm_h__
+#define BinghamModelVelocityBCAlgorithm_h__
+
+#include "RheologicalVelocityBCAlgorithm.h"
+#include "Thixotropy.h"
+
+class BinghamModelVelocityBCAlgorithm : public RheologicalVelocityBCAlgorithm
+{
+public:
+ BinghamModelVelocityBCAlgorithm()
+ {
+ BCAlgorithm::type = BCAlgorithm::BinghamModelVelocityBCAlgorithm;
+ BCAlgorithm::preCollision = false;
+ }
+ ~BinghamModelVelocityBCAlgorithm() {}
+ SPtr<BCAlgorithm> clone() override
+ {
+ SPtr<BCAlgorithm> bc(new BinghamModelVelocityBCAlgorithm());
+ return bc;
+ }
+protected:
+ LBMReal getThyxotropyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho) const override
+ {
+ return Thixotropy::getBinghamCollFactor(omegaInf, shearRate, drho);
+ }
+};
+#endif // BinghamModelVelocityBCAlgorithm_h__
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
new file mode 100644
index 000000000..2614c703b
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
@@ -0,0 +1,47 @@
+#include "RheologicalVelocityBCAlgorithm.h"
+#include "DistributionArray3D.h"
+#include "BoundaryConditions.h"
+
+RheologicalVelocityBCAlgorithm::RheologicalVelocityBCAlgorithm()
+{
+ //BCAlgorithm::type = BCAlgorithm::RheologicalVelocityBCAlgorithm;
+ //BCAlgorithm::preCollision = false;
+}
+//////////////////////////////////////////////////////////////////////////
+RheologicalVelocityBCAlgorithm::~RheologicalVelocityBCAlgorithm()
+{
+}
+//////////////////////////////////////////////////////////////////////////
+void RheologicalVelocityBCAlgorithm::addDistributions(SPtr<DistributionArray3D> distributions)
+{
+ this->distributions = distributions;
+}
+//////////////////////////////////////////////////////////////////////////
+void RheologicalVelocityBCAlgorithm::applyBC()
+{
+ LBMReal f[D3Q27System::ENDF+1];
+ LBMReal feq[D3Q27System::ENDF+1];
+ distributions->getDistributionInv(f, x1, x2, x3);
+ LBMReal rho, vx1, vx2, vx3, drho;
+ calcMacrosFct(f, drho, vx1, vx2, vx3);
+ calcFeqFct(feq, drho, vx1, vx2, vx3);
+
+ LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
+ LBMReal collFactorF = getThyxotropyCollFactor(collFactor, shearRate, drho);
+
+ rho = 1.0+drho*compressibleFactor;
+
+ for (int fdir = D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++)
+ {
+ if (bcPtr->hasVelocityBoundaryFlag(fdir))
+ {
+ const int invDir = D3Q27System::INVDIR[fdir];
+ LBMReal q = bcPtr->getQ(invDir);// m+m q=0 stabiler
+ LBMReal velocity = bcPtr->getBoundaryVelocity(invDir);
+ LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity*rho)/(1.0+q));
+ distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
+ }
+ }
+
+}
+
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h
new file mode 100644
index 000000000..60237e3db
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h
@@ -0,0 +1,22 @@
+#ifndef RheologicalVelocityBCAlgorithm_h__
+#define RheologicalVelocityBCAlgorithm_h__
+
+#include "BCAlgorithm.h"
+#include <PointerDefinitions.h>
+
+class DistributionArray3D;
+
+class RheologicalVelocityBCAlgorithm : public BCAlgorithm
+{
+public:
+ RheologicalVelocityBCAlgorithm();
+ ~RheologicalVelocityBCAlgorithm();
+ virtual SPtr<BCAlgorithm> clone() override { UB_THROW(UbException("LBMReal clone() - belongs in the derived class")); }
+ void addDistributions(SPtr<DistributionArray3D> distributions);
+ void applyBC() override;
+protected:
+ virtual LBMReal getThyxotropyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho) const = 0; // { UB_THROW(UbException("LBMReal getThyxotropyCollFactor() - belongs in the derived class")); }
+};
+
+#endif // VelocityBoundaryCondition_h__
+
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp
index 083b8703e..0ee6d281d 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp
@@ -24,13 +24,13 @@ void ThixotropyNoSlipBCAlgorithm::applyBC()
calcMacrosFct(f, rho, vx1, vx2, vx3);
calcFeqFct(feq, rho, vx1, vx2, vx3);
+ LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
+ LBMReal collFactorF = getThyxotropyCollFactor(collFactor, shearRate, rho);
+
for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++)
{
if (bcPtr->hasNoSlipBoundaryFlag(fDir))
{
- LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
- LBMReal collFactorF = getThyxotropyCollFactor(collFactor, shearRate, rho);
-
//quadratic bounce back
const int invDir = D3Q27System::INVDIR[fDir];
LBMReal q = bcPtr->getQ(invDir);
diff --git a/src/cpu/VirtualFluidsCore/LBM/BinghamModelLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/BinghamModelLBMKernel.h
index 03f6a7c56..8a3f1e58f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/BinghamModelLBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/BinghamModelLBMKernel.h
@@ -1,12 +1,12 @@
#ifndef BinghamModelLBMKernel_H
#define BinghamModelLBMKernel_H
-#include "ThixotropyModelLBMKernel2.h"
+#include "ThixotropyModelLBMKernel.h"
#include "Thixotropy.h"
//! \brief Cumulant LBM kernel + Bingham plastic model
//! \author K. Kutscher, M. Geier
-class BinghamModelLBMKernel : public ThixotropyModelLBMKernel2
+class BinghamModelLBMKernel : public ThixotropyModelLBMKernel
{
public:
BinghamModelLBMKernel() {};
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
index c7bd879cb..9df968e12 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
@@ -618,8 +618,8 @@ 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) / (drho + c1);
//omega = getThyxotropyCollFactor(omega, shearRate, rho);
- omega = Thixotropy::getHerschelBulkleyCollFactor(omega, shearRate, drho);
- //omega = Thixotropy::getBinghamCollFactor(omega, shearRate, drho);
+ //omega = Thixotropy::getHerschelBulkleyCollFactor(omega, shearRate, drho);
+ omega = Thixotropy::getBinghamCollFactor(omega, shearRate, drho);
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
dxux = c1o2 * (-omega) * (mxxMyy + mxxMzz);// +c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
@@ -639,6 +639,8 @@ void RheologyK17LBMKernel::calculate(int step)
if(omega < c1) { omega = c1; } //arbitrary limit (24.09.2020)
+ omega = collFactor;
+
//magic parameter for rheology
LBMReal a = 10;
OxxPyyPzz = c1 / (a * ((c1 / omega) - c1o2) + c1o2);
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
index 50caabbe8..0b2aee946 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
@@ -66,17 +66,17 @@ inline LBMReal Thixotropy::getBinghamCollFactor(LBMReal omegaInf, LBMReal shearR
//LBMReal omega = omegaInf / (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a))))))))));
//20 iterations
- LBMReal omega = omegaInf / (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a))))))))))))))))))));
+ //LBMReal omega = omegaInf / (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a))))))))))))))))))));
- //LBMReal omega = omegaInf*cs2 * shearRate * rho / (cs2 * shearRate * rho + omegaInf * tau0);
- //LBMReal shearRateNew = shearRate * (omega / omegaInf);
-
- //for (int i = 0; i < 20; i++)
- //{
- // omega = omegaInf * cs2 * shearRateNew * rho / (cs2 * shearRateNew * rho + omegaInf * tau0);
- // shearRateNew = shearRate * (omega / omegaInf);
- //}
- //omega = omegaInf * cs2 * shearRateNew * rho / (cs2 * shearRateNew * rho + omegaInf * tau0);
+ LBMReal omega = omegaInf*cs2 * shearRate * rho / (cs2 * shearRate * rho + omegaInf * tau0);
+ LBMReal shearRateNew = shearRate * (omega / omegaInf);
+
+ for (int i = 0; i < 3; i++)
+ {
+ omega = omegaInf * cs2 * shearRateNew * rho / (cs2 * shearRateNew * rho + omegaInf * tau0);
+ shearRateNew = shearRate * (omega / omegaInf);
+ }
+ omega = omegaInf * cs2 * shearRateNew * rho / (cs2 * shearRateNew * rho + omegaInf * tau0);
//if (omega < 0.2)
// omega = 0.2;
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
index 0b416730d..8839818a7 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
@@ -481,7 +481,7 @@ void ThixotropyModelLBMKernel2::calculate(int step)
LBMReal Dyz = -three * collFactorF * mfabb;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//non Newtonian fluid collision factor
- //LBMReal shearRate = sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+ LBMReal shearRate = sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
LBMReal shearFactor = sqrt(c1o2 * ((mfcaa - mfaaa * c1o3) * (mfcaa - mfaaa * c1o3) + (mfaca - mfaaa * c1o3) * (mfaca - mfaaa * c1o3) + (mfaac - mfaaa * c1o3) * (mfaac - mfaaa * c1o3)) + mfbba * mfbba + mfbab * mfbab + mfabb * mfabb) + UbMath::Epsilon<LBMReal>::val();
@@ -504,13 +504,24 @@ void ThixotropyModelLBMKernel2::calculate(int step)
SPtr<Thixotropy> thix = Thixotropy::getInstance();
LBMReal tau0 = thix->getYieldStress();
- mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz + ((mxxPyyPzz-mfaaa)/shearFactor*tau0)) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
- mxxMyy += collFactorF * (-mxxMyy + mxxMyy/shearFactor*tau0) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vy2 * dyuy);
- mxxMzz += collFactorF * (-mxxMzz + mxxMzz/shearFactor*tau0) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vz2 * dzuz);
+ mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz /*+ ((mxxPyyPzz-mfaaa)/shearFactor*tau0)*/) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+ //mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+ //mxxMyy += collFactorF * (-mxxMyy + mxxMyy/shearFactor*tau0) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vy2 * dyuy);
+ //mxxMzz += collFactorF * (-mxxMzz + mxxMzz/shearFactor*tau0) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vz2 * dzuz);
+
+ //mfabb += collFactorF * (-mfabb + mfabb/shearFactor*tau0);
+ //mfbab += collFactorF * (-mfbab + mfbab/shearFactor*tau0);
+ //mfbba += collFactorF * (-mfbba + mfbba/shearFactor*tau0);
+
+ collFactorF = collFactor * (c1 - tau0 / shearFactor);
+
+ mxxMyy += collFactorF * (-mxxMyy/* + mxxMyy / shearFactor * tau0*/) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vy2 * dyuy);
+ mxxMzz += collFactorF * (-mxxMzz/* + mxxMzz / shearFactor * tau0*/) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vz2 * dzuz);
+
+ mfabb += collFactorF * (-mfabb/* + mfabb / shearFactor * tau0*/);
+ mfbab += collFactorF * (-mfbab/* + mfbab / shearFactor * tau0*/);
+ mfbba += collFactorF * (-mfbba/* + mfbba / shearFactor * tau0*/);
- mfabb += collFactorF * (-mfabb + mfabb/shearFactor*tau0);
- mfbab += collFactorF * (-mfbab + mfbab/shearFactor*tau0);
- mfbba += collFactorF * (-mfbba + mfbba/shearFactor*tau0);
// linear combinations back
mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
--
GitLab