diff --git a/apps/cpu/rheometer/rheometer.cfg b/apps/cpu/rheometer/rheometer.cfg
index 530ece5ff031f9feadf34ef7e86c6caf99a3073a..0330c92976b4e67a72114f7c4eb556fe7d40f5e7 100644
--- a/apps/cpu/rheometer/rheometer.cfg
+++ b/apps/cpu/rheometer/rheometer.cfg
@@ -1,11 +1,13 @@
-outputPath = d:/temp/rheometer/rheometerBingham128_8_test_New_Omega
+outputPath = d:/temp/rheometer/rheometerBinghamNewEquilibriumSBB/rheometerBingham_nu_1.5e-3_nnn
+
+viscosityPath = d:/Projects/VirtualFluidsCombined/apps/cpu/rheometer
numOfThreads = 4
availMem = 8e9
logToFile = false
blocknx = 8 8 1
-boundingBox = 128 128 1
+#boundingBox = 32 32 1
deltax = 1
#boundingBox = 0.02 0.02 0.00125
@@ -13,26 +15,17 @@ deltax = 1
refineLevel = 0
-radius = 5
-sphereCenter = 25 25 25
-
-#R0=0.02 #m
-#Ri=0.0075 #m
-
-uLB = 8e-4
-N = 20 #rpm
-
-
-n = 1
-Re = 10
-Bn = 0.01
+OmegaLB = 1e-4
+tau0 = 1e-6
+resolution = 32
+scaleFactor = 1
newStart = true
restartStep = 100000
-cpStart = 100000
-cpStep = 100000
+cpStart = 10000
+cpStep = 10000
outTime = 10000
-endTime = 1600000
\ No newline at end of file
+endTime = 100000
\ No newline at end of file
diff --git a/apps/cpu/rheometer/rheometer.cpp b/apps/cpu/rheometer/rheometer.cpp
index cc67f85af4408e3107d9796d3a77f9f8c63bae2e..c7bc4f3e6dac523bf4b30918ee434068ece9a1bd 100644
--- a/apps/cpu/rheometer/rheometer.cpp
+++ b/apps/cpu/rheometer/rheometer.cpp
@@ -14,10 +14,11 @@ void bflow(string configname)
config.load(configname);
string outputPath = config.getValue<string>("outputPath");
+ string viscosityPath = config.getValue<string>("viscosityPath");
int numOfThreads = config.getValue<int>("numOfThreads");
vector<int> blocknx = config.getVector<int>("blocknx");
- vector<double> boundingBox = config.getVector<double>("boundingBox");
- double nuLB = 1.5e-3;//config.getValue<double>("nuLB");
+ //vector<double> boundingBox = config.getVector<double>("boundingBox");
+ //double nuLB = 1.5e-3;//config.getValue<double>("nuLB");
double endTime = config.getValue<double>("endTime");
double outTime = config.getValue<double>("outTime");
double availMem = config.getValue<double>("availMem");
@@ -25,18 +26,19 @@ void bflow(string configname)
bool logToFile = config.getValue<bool>("logToFile");
double restartStep = config.getValue<double>("restartStep");
double deltax = config.getValue<double>("deltax");
- // double radius = config.getValue<double>("radius");
double cpStep = config.getValue<double>("cpStep");
double cpStart = config.getValue<double>("cpStart");
bool newStart = config.getValue<bool>("newStart");
- double uLB = config.getValue<double>("uLB");
- //double nuLB = config.getValue<double>("nuLB");
+ double OmegaLB = config.getValue<double>("OmegaLB");
double tau0 = config.getValue<double>("tau0");
double scaleFactor = config.getValue<double>("scaleFactor");
- // double N = config.getValue<double>("N");
- //vector<double> sphereCenter = config.getVector<double>("sphereCenter");
+ double resolution = config.getValue<double>("resolution");
- outputPath = outputPath + "3";
+ ConfigurationFile viscosity;
+ viscosity.load(viscosityPath + "/viscosity.cfg");
+ double nuLB = viscosity.getValue<double>("nuLB");
+
+ outputPath = outputPath + "/rheometerBingham_" + config.getValue<string>("resolution") + "_" + config.getValue<string>("OmegaLB");
SPtr<Communicator> comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
@@ -61,7 +63,7 @@ void bflow(string configname)
LBMReal rhoLB = 0.0;
- uLB /= scaleFactor;
+ OmegaLB /= scaleFactor;
tau0 /= scaleFactor;
endTime *= scaleFactor;
@@ -83,9 +85,9 @@ void bflow(string configname)
double g_minX2 = 0;
double g_minX3 = 0;
- double g_maxX1 = boundingBox[0];
- double g_maxX2 = boundingBox[1];
- double g_maxX3 = boundingBox[2];
+ double g_maxX1 = resolution;// boundingBox[0];
+ double g_maxX2 = resolution;// boundingBox[1];
+ double g_maxX3 = 1.0; // boundingBox[2];
//double g_minX1 = -boundingBox[0]/2.0;
//double g_minX2 = -boundingBox[1] / 2.0;
@@ -132,13 +134,13 @@ void bflow(string configname)
mu::Parser fctVx;
//fctVx.SetExpr("omega*(r-x2)");
fctVx.SetExpr("-Omega*(x2-r)");
- fctVx.DefineConst("Omega", uLB);
+ fctVx.DefineConst("Omega", OmegaLB);
//fctVx.DefineConst("r", R0);
fctVx.DefineConst("r", g_maxX1*0.5);
mu::Parser fctVy;
fctVy.SetExpr("Omega*(x1-r)");
- fctVy.DefineConst("Omega", uLB);
+ fctVy.DefineConst("Omega", OmegaLB);
//fctVy.DefineConst("r", R0);
fctVy.DefineConst("r", g_maxX2 * 0.5);
@@ -157,7 +159,7 @@ void bflow(string configname)
//BS visitor
BoundaryConditionsBlockVisitor bcVisitor;
- bcVisitor.addBC(noSlipBCAdapter);
+ //bcVisitor.addBC(noSlipBCAdapter);
//bcVisitor.addBC(slipBCAdapter);
bcVisitor.addBC(velocityBCAdapter);
//bcVisitor.addBC(densityBCAdapter);
@@ -184,11 +186,6 @@ void bflow(string configname)
SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), outputPath + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
- //sphere
- //SPtr<GbObject3D> sphere(new GbSphere3D(sphereCenter[0], sphereCenter[1], sphereCenter[2], radius));
- //GbSystem3D::writeGeoObject(sphere.get(), outputPath + "/geo/sphere", WbWriterVtkXmlBinary::getInstance());
- //SPtr<D3Q27Interactor> sphereInt(new D3Q27Interactor(sphere, grid, noSlipBCAdapter, Interactor3D::SOLID));
-
//////////////////////////////////////////////////////////////////////////
//restart
SPtr<UbScheduler> mSch(new UbScheduler(cpStep, cpStart));
@@ -215,7 +212,7 @@ void bflow(string configname)
UBLOG(logINFO, "Parameters:");
//UBLOG(logINFO, "forcing = " << forcing);
UBLOG(logINFO, "rho = " << rhoLB);
- UBLOG(logINFO, "uLB = " << uLB);
+ UBLOG(logINFO, "uLB = " << OmegaLB);
UBLOG(logINFO, "nuLB = " << nuLB);
// UBLOG(logINFO, "Re = " << (U * d) / (k * std::pow(Gamma, n - 1)));
// UBLOG(logINFO, "Bn = " << tau0 /(k * std::pow(Gamma, n)));
@@ -228,7 +225,8 @@ void bflow(string configname)
UBLOG(logINFO, "number of threads = " << numOfThreads);
UBLOG(logINFO, "number of processes = " << comm->getNumberOfProcesses());
UBLOG(logINFO, "blocknx = " << blocknx[0] << " " << blocknx[1] << " " << blocknx[2]);
- UBLOG(logINFO, "boundingBox = " << boundingBox[0] << " " << boundingBox[1] << " " << boundingBox[2]);
+ UBLOG(logINFO, "resolution = " << resolution);
+ // UBLOG(logINFO, "boundingBox = " << boundingBox[0] << " " << boundingBox[1] << " " << boundingBox[2]);
// UBLOG(logINFO, "sphereCenter = " << sphereCenter[0] << " " << sphereCenter[1] << " " << sphereCenter[2]);
UBLOG(logINFO, "output path = " << outputPath);
UBLOG(logINFO, "Preprozess - start");
@@ -254,33 +252,15 @@ void bflow(string configname)
//walls
- GbCuboid3DPtr wallZmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_minX3));
- if (myid == 0) GbSystem3D::writeGeoObject(wallZmin.get(), outputPath + "/geo/wallZmin", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr wallZmax(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_maxX3, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(wallZmax.get(), outputPath + "/geo/wallZmax", WbWriterVtkXmlASCII::getInstance());
-
- //GbCuboid3DPtr wallYmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_minX2, g_maxX3 + blockLength));
- //if (myid == 0) GbSystem3D::writeGeoObject(wallYmin.get(), outputPath + "/geo/wallYmin", WbWriterVtkXmlASCII::getInstance());
-
- //GbCuboid3DPtr wallYmax(new GbCuboid3D(g_minX1 - blockLength, g_maxX2, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
- //if (myid == 0) GbSystem3D::writeGeoObject(wallYmax.get(), outputPath + "/geo/wallYmax", WbWriterVtkXmlASCII::getInstance());
-
- //GbCuboid3DPtr wallXmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_minX1, g_maxX2 + blockLength, g_maxX3 + blockLength));
- //if (myid == 0) GbSystem3D::writeGeoObject(wallXmin.get(), outputPath + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
-
- //GbCuboid3DPtr wallXmax(new GbCuboid3D(g_maxX1, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
- //if (myid == 0) GbSystem3D::writeGeoObject(wallXmax.get(), outputPath + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
+ //GbCuboid3DPtr wallZmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_minX3));
+ //if (myid == 0) GbSystem3D::writeGeoObject(wallZmin.get(), outputPath + "/geo/wallZmin", WbWriterVtkXmlASCII::getInstance());
- //wall interactors
- SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
-
- //SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, slipBCAdapter, Interactor3D::SOLID));
- //SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, slipBCAdapter, Interactor3D::SOLID));
+ //GbCuboid3DPtr wallZmax(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_maxX3, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ //if (myid == 0) GbSystem3D::writeGeoObject(wallZmax.get(), outputPath + "/geo/wallZmax", WbWriterVtkXmlASCII::getInstance());
- //SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, velocityBCAdapter, Interactor3D::SOLID));
- //SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, densityBCAdapter, Interactor3D::SOLID));
+ ////wall interactors
+ //SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ //SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
////////////////////////////////////////////
//METIS
@@ -289,10 +269,6 @@ void bflow(string configname)
/////delete solid blocks
if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - start");
InteractorsHelper intHelper(grid, metisVisitor);
- //intHelper.addInteractor(wallXminInt);
- //intHelper.addInteractor(wallXmaxInt);
- //intHelper.addInteractor(wallYminInt);
- //intHelper.addInteractor(wallYmaxInt);
//intHelper.addInteractor(wallZminInt);
//intHelper.addInteractor(wallZmaxInt);
intHelper.addInteractor(statorInt);
diff --git a/apps/cpu/rheometer/viscosity.cfg b/apps/cpu/rheometer/viscosity.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..bf2822c7b1d6dd42fdcbc513a4e6b39264ab4180
--- /dev/null
+++ b/apps/cpu/rheometer/viscosity.cfg
@@ -0,0 +1 @@
+nuLB = 1.5e-3
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/LBM/BinghamModelLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/BinghamModelLBMKernel.h
index 8a3f1e58fd8209c96bbc82f71cfbf3316fa2cf6c..03f6a7c56e7dabfa88b867a9ebdca9975ac80c07 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 "ThixotropyModelLBMKernel.h"
+#include "ThixotropyModelLBMKernel2.h"
#include "Thixotropy.h"
//! \brief Cumulant LBM kernel + Bingham plastic model
//! \author K. Kutscher, M. Geier
-class BinghamModelLBMKernel : public ThixotropyModelLBMKernel
+class BinghamModelLBMKernel : public ThixotropyModelLBMKernel2
{
public:
BinghamModelLBMKernel() {};
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
index da74d7a67710f8e980eaeadd5d36e7dc40596c3f..0b416730df883856682efd0d2a1eaa3ad73a802d 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
@@ -5,6 +5,7 @@
#include <math.h>
#include "DataSet3D.h"
#include "LBMKernel.h"
+#include "Thixotropy.h"
#define PROOF_CORRECTNESS
@@ -480,8 +481,11 @@ 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);
- collFactorF = getThyxotropyCollFactor(collFactorF, shearRate, rho);
+ //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();
+
+ //collFactorF = getThyxotropyCollFactor(collFactorF, shearRate, rho);
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//relax
@@ -497,13 +501,16 @@ void ThixotropyModelLBMKernel2::calculate(int step)
//mfbab += getThyxotropyCollFactor(collFactorF, std::abs(Dxz) / (rho + one), rho) * (-mfbab);
//mfbba += getThyxotropyCollFactor(collFactorF, std::abs(Dxy) / (rho + one), rho) * (-mfbba);
- mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
- mxxMyy += collFactorF * (-mxxMyy) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vy2 * dyuy);
- mxxMzz += collFactorF * (-mxxMzz) - 3. * (1. - c1o2 * collFactorF) * (vx2 * dxux - vz2 * dzuz);
+ 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);
- mfabb += collFactorF * (-mfabb);
- mfbab += collFactorF * (-mfbab);
- mfbba += collFactorF * (-mfbba);
+ 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);