add K17-kernel implementation with rheological properties and add the output of effective viscosity

apps/cpu/HerschelBulkleySphere/hbsphere.cfg

-outputPath = d:/temp/Herschel-BulkleySphere
+outputPath = d:/temp/Herschel-Bulkley_sphere_MP_test
 
 numOfThreads = 4
 availMem = 8e9
 logToFile = false
 
-blocknx = 50 50 50
-boundingBox = 150 50 50  #30*20=600**3=216000000
+blocknx = 25 25 25
+boundingBox = 250 50 50  #30*20=600**3=216000000
 deltax = 1
 
 refineLevel = 0
 
 radius = 5
+sphereCenter = 25 25 25
+
 velocity = 1e-3
 n = 0.9
 Re = 1

apps/cpu/HerschelBulkleySphere/hbsphere.cpp

@@ -131,7 +131,8 @@ void bflow(string configname)
       SPtr<BCProcessor> bcProc;
       bcProc = SPtr<BCProcessor>(new BCProcessor());
 
-      SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
+      SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
+      //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
       //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
       kernel->setBCProcessor(bcProc);
       //kernel->setForcingX1(forcing);
@@ -329,11 +330,14 @@ void bflow(string configname)
       SPtr<CalculateForcesCoProcessor> fp = make_shared<CalculateForcesCoProcessor>(grid, forceSch, outputPath + "/forces/forces.txt", comm, velocity, area);
       fp->addInteractor(sphereInt);
 
+      SPtr<WriteThixotropyQuantitiesCoProcessor> writeThixotropicMQCoProcessor(new WriteThixotropyQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+
       SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
       SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
       calculator->addCoProcessor(npr);
       calculator->addCoProcessor(fp);
       calculator->addCoProcessor(writeMQCoProcessor);
+      calculator->addCoProcessor(writeThixotropicMQCoProcessor);
       calculator->addCoProcessor(restartCoProcessor);
 
       if (myid == 0) UBLOG(logINFO, "Simulation-start");

src/cpu/VirtualFluids.h

@@ -223,6 +223,7 @@
 #include <LBM/BinghamModelLBMKernel.h>
 #include <LBM/HerschelBulkleyModelLBMKernel.h>
 #include <LBM/ThixotropyInterpolationProcessor.h>
+#include <LBM/RheologyK17LBMKernel.h>
 
 
 #include <geometry3d/CoordinateTransformation3D.h>

src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp

@@ -11,6 +11,7 @@
 #include <numeric>
 #include "basics/writer/WbWriterVtkXmlASCII.h"
 #include "ThixotropyExpLBMKernel.h"
+#include "Thixotropy.h"
 
 using namespace std;
 
@@ -119,7 +120,8 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 
 	//Diese Daten werden geschrieben:
 	datanames.resize(0);
-	datanames.push_back("lambda");
+	datanames.push_back("viscosity");
+	//datanames.push_back("lambda");
 	//datanames.push_back("ShearRate");
 	//datanames.push_back("collFactor");
 	//datanames.push_back("Fluxx");
@@ -139,23 +141,23 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
    SPtr<ILBMKernel> kernel = block->getKernel();
    SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();          
    SPtr<DistributionArray3D> distributionsF = kernel->getDataSet()->getFdistributions(); 
-	SPtr<DistributionArray3D> distributionsH = kernel->getDataSet()->getHdistributions();
+	//SPtr<DistributionArray3D> distributionsH = kernel->getDataSet()->getHdistributions();
 	//LBMReal collFactorF = staticPointerCast<ThixotropyExpLBMKernel>(kernel)->getCollisionFactorF();
 	LBMReal collFactor = kernel->getCollisionFactor();
 	LBMReal f[D3Q27System::ENDF + 1];
-	LBMReal h[D3Q27System::ENDF + 1];
-	LBMReal lambda, gammaDot;
+	//LBMReal h[D3Q27System::ENDF + 1];
+	//LBMReal viscosity=0; // lambda, gammaDot;
 	
 	//knotennummerierung faengt immer bei 0 an!
-	unsigned int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+	int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
 
 	int minX1 = 0;
 	int minX2 = 0;
 	int minX3 = 0;
 
-	int maxX1 = (int)(distributionsH->getNX1());
-	int maxX2 = (int)(distributionsH->getNX2());
-	int maxX3 = (int)(distributionsH->getNX3());
+	int maxX1 = (int)(distributionsF->getNX1());
+	int maxX2 = (int)(distributionsF->getNX2());
+	int maxX3 = (int)(distributionsF->getNX3());
 
 	//int minX1 = 1;
 	//int minX2 = 1;
@@ -187,25 +189,33 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 						float(val<2>(org) - val<2>(nodeOffset) + ix2*dx),
 						float(val<3>(org) - val<3>(nodeOffset) + ix3*dx)));
 
-					distributionsH->getDistribution(h, ix1, ix2, ix3);
+					//distributionsH->getDistribution(h, ix1, ix2, ix3);
 
-					lambda = D3Q27System::getDensity(h);
+					//lambda = D3Q27System::getDensity(h);
 
-					if (UbMath::isNaN(lambda) || UbMath::isInfinity(lambda) )
-						UB_THROW(UbException(UB_EXARGS, "lambda is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" + block->toString() +
-							", node=" + UbSystem::toString(ix1) + "," + UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+					//if (UbMath::isNaN(lambda) || UbMath::isInfinity(lambda) )
+					//	UB_THROW(UbException(UB_EXARGS, "lambda is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" + block->toString() +
+					//		", node=" + UbSystem::toString(ix1) + "," + UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
 
-					distributionsF->getDistribution(f, ix1, ix2, ix3);
-					LBMReal rho = D3Q27System::getDensity(f);
+					//distributionsF->getDistribution(f, ix1, ix2, ix3);
+					//LBMReal rho = D3Q27System::getDensity(f);
 					
 					//gammaDot = BinghamModel::getShearRate(f, collFactor);
 
 					//LBMReal collFactorF = collFactor - 1e-6 / (gammaDot + one * 1e-9);
 					//collFactorF = (collFactorF < 0.5) ? 0.5 : collFactorF;
 
-					data[index++].push_back(lambda);
+					//data[index++].push_back(lambda);
 					//data[index++].push_back(gammaDot);
 					//data[index++].push_back(collFactorF);
+
+					distributionsF->getDistribution(f, ix1, ix2, ix3);
+					LBMReal rho = D3Q27System::getDensity(f);
+					LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
+					LBMReal omega = Thixotropy::getHerschelBulkleyCollFactor(collFactor, shearRate, rho);
+					LBMReal viscosity = (omega == 0) ? 0 : c1o3 * (c1/omega-c1o2);
+
+					data[index++].push_back(viscosity);
 				}
 			}
 		}
@@ -229,7 +239,7 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 					&& (NET = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0
 					&& (NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0)
 				{
-					cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
+					cells.push_back(makeUbTuple((unsigned int)SWB, (unsigned int)SEB, (unsigned int)NEB, (unsigned int)NWB, (unsigned int)SWT, (unsigned int)SET, (unsigned int)NET, (unsigned int)NWT));
 				}
 			}
 		}

src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h

+#ifndef RheologyK17LBMKernel_h__
+#define RheologyK17LBMKernel_h__
+
+#include "LBMKernel.h"
+#include "BCProcessor.h"
+#include "D3Q27System.h"
+#include "basics/utilities/UbTiming.h"
+#include "basics/container/CbArray4D.h"
+#include "basics/container/CbArray3D.h"
+
+//! \brief   compressible cumulant LBM kernel. 
+//! \details CFD solver that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model
+//! \author  K. Kutscher, M. Geier
+class RheologyK17LBMKernel :  public LBMKernel
+{
+public:
+   //! This option set relaxation parameter: NORMAL  
+   enum Parameter{NORMAL, MAGIC};
+public:
+   RheologyK17LBMKernel();
+   virtual ~RheologyK17LBMKernel(void);
+   virtual void calculate(int step);
+   virtual SPtr<LBMKernel> clone();
+   double getCalculationTime() override;
+   //! The value should not be equal to a shear viscosity
+   void setBulkViscosity(LBMReal value);
+protected:
+   virtual void initDataSet();
+   LBMReal f[D3Q27System::ENDF+1];
+
+   UbTimer timer;
+
+   CbArray4D<LBMReal,IndexerX4X3X2X1>::CbArray4DPtr localDistributions;
+   CbArray4D<LBMReal,IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributions;
+   CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr   zeroDistributions;
+
+   mu::value_type muX1,muX2,muX3;
+   mu::value_type muDeltaT;
+   mu::value_type muNu;
+   LBMReal forcingX1;
+   LBMReal forcingX2;
+   LBMReal forcingX3;
+   
+   // bulk viscosity
+   LBMReal OxxPyyPzz; //omega2 (bulk viscosity)
+   LBMReal bulkViscosity;
+
+};
+#endif // RheologyK17LBMKernel_h__
+
+