Merge pull request #503 from TravisMitchell/overlay/thermo

Move thermocapillary model to TCLB/overlay framework

Author: llaniewski

README.md

@@ -106,6 +106,29 @@ The DEM codes that TCLB can be integrated with are:
 
 Refer to the documentation for instructions on compilation and coupling.
 
+## Models
+For users looking to apply existing LBM methods, common/supported models are below. Note extensions to these models exist using the [TCLB's overlay](https://github.com/CFD-GO/TCLB_overlay) framework and TCLB optional compile flags.
+
+**Two-Dimensional**
+- d2q9: MRT LBM for single-phase flow.
+- d2q9_les: MRT LBM with Smagorinski LES turbulence model.
+- d2q9q9_cm_cht: thermal LBM with Boussinesq approx for coupling and cumulant or cascaded relaxation kernels.
+- d2q9_pf_velocity: multiphase LBM based on the phase field model and incompressible LBM.
+
+**Three-Dimensional**
+- d3q27_cumulant: cumulant LBM with options for:
+     * Interpolated bounceback.
+     * Smagorinski LES turbulence model.
+- d3q27q27_cm_cht: thermal LBM with Boussinesq approx for coupling and cumulant or cascaded collision relaxation kernels.
+- d3q27_pf_velocity: multiphase LBM based on the phase field model and incompressible LBM.
+     * Options for various contact angle implementations (surface energy or geometric)
+     * [Thermocapillary flow extension](https://github.com/TravisMitchell/thermocapillary)
+
+**Particle (DEM) Coupled**
+- d3q27_PSM: Applies the partially saturated model for coupling particles in single phase flow.
+     * Options for Two-Relaxation-Time kernel
+     * Option for Non-Equilibiurm Bounce-Back and Superposition for the DEM-LBM coupling.
+
 ## About
 
 ### Authors

models/multiphase/d3q27_pf_velocity/Dynamics.R

@@ -108,14 +108,6 @@ if (Options$geometric){
     AddField("PhaseF",stencil3d=1, group="PF")
 }
 
-if (Options$thermo){
-    source("thermocapillary.R")
-
-    save_initial_PF = c(save_initial_PF,"Thermal")
-    save_iteration  = c(save_iteration, "Thermal")
-    load_iteration  = c(load_iteration, "Thermal")
-}
-
 ######################
 ########STAGES########
 ######################
@@ -138,26 +130,11 @@ if (Options$thermo){
 		AddStage("calcWall" , "calcWallPhase", save=Fields$name=="PhaseF", load=DensityAll$group %in% c("nw", "solid_boundary", extra_fields_to_load_for_bc))
 		AddStage("calcWallPhase_correction", "calcWallPhase_correction", save=Fields$name=="PhaseF", load=DensityAll$group %in% c("nw", "solid_boundary"))
 	}
-	if (Options$thermo){
-		AddStage("CopyDistributions", "TempCopy",  save=Fields$group %in% c("g","h","Vel","nw", "PF","Thermal"))
-		AddStage("CopyThermal","ThermalCopy", save=Fields$name %in% c("Temp","Cond","SurfaceTension"), load=DensityAll$name %in% c("Temp","Cond","SurfaceTension"))
-		AddStage("RK_1", "TempUpdate1", save=Fields$name=="RK1", load=DensityAll$name %in% c("U","V","W","Cond","PhaseF","Temp"))
-		AddStage("RK_2", "TempUpdate2", save=Fields$name=="RK2", load=DensityAll$name %in% c("U","V","W","RK1","Cond","PhaseF","Temp"))
-		AddStage("RK_3", "TempUpdate3", save=Fields$name=="RK3", load=DensityAll$name %in% c("U","V","W","RK1","RK2","Cond","PhaseF","Temp"))
-		AddStage("RK_4", "TempUpdate4", save=Fields$name %in% c("Temp","SurfaceTension"), load=DensityAll$name %in% c("U","V","W","RK1","RK2","RK3","Cond","PhaseF","Temp"))
-
-		AddStage("NonLocalTemp","BoundUpdate", save=Fields$name %in% c("Temp","SurfaceTension"), load=DensityAll$name %in% c("Temp"))
-	}
 
 #######################
 ########ACTIONS########
 #######################
-	if (Options$thermo){	
-		AddAction("TempToSteadyState", c("CopyDistributions","RK_1", "RK_2", "RK_3", "RK_4","NonLocalTemp"))
-		AddAction("Iteration", c("BaseIter", "calcPhase", "calcWall","RK_1", "RK_2", "RK_3", "RK_4","NonLocalTemp"))
-		AddAction("IterationConstantTemp", c("BaseIter", "calcPhase", "calcWall","CopyThermal"))
-		AddAction("Init"     , c("PhaseInit","WallInit" , "calcWall","BaseInit"))
-	} else if (Options$geometric) {
+	if (Options$geometric) {
         calcGrad <- if (Options$isograd)  "calcPhaseGrad" else "calcPhaseGrad_init"
         AddAction("Iteration", c("BaseIter", "calcPhase",  calcGrad, "calcWall_CA", "calcWallPhase_correction"))
 	    AddAction("Init"     , c("PhaseInit","WallInit_CA" , "calcPhaseGrad_init"  , "calcWall_CA", "calcWallPhase_correction", "BaseInit"))

models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt

@@ -42,6 +42,7 @@
 //		xx/xx/2022: Extension of phase field to d3q27 option to move away from crappy q15 lattice.
 //		xx/xx/2022: Development of geometric wetting B.C. option as well as existing surface energy.
 //		xx/xx/2023: Code extension for wetting on curved boundaries. 
+// 		30/01/2024: Thermocapillary support moved to TCLB overlay: https://github.com/TravisMitchell/thermocapillary
 
 #include <math.h>
 #define PI 3.14159265
@@ -217,13 +218,8 @@ CudaDeviceFunction real_t calcMu(real_t C)
 				myLaplace('lpPhi', 'PhaseF')
 			?>
 	#endif
-	#ifdef OPTIONS_thermo
-		mu = 4.0*(12.0*SurfaceTension(0,0,0)/IntWidth)*(C-PhaseField_l)*(C-PhaseField_h)*(C-pfavg)
-			- (1.5 *SurfaceTension(0,0,0)*IntWidth) * lpPhi;
-	#else
-		mu = 4.0*(12.0*sigma/IntWidth)*(C-PhaseField_l)*(C-PhaseField_h)*(C-pfavg)
-			- (1.5 *sigma*IntWidth) * lpPhi;
-	#endif
+	mu = 4.0*(12.0*sigma/IntWidth)*(C-PhaseField_l)*(C-PhaseField_h)*(C-pfavg)
+		 - (1.5 *sigma*IntWidth) * lpPhi;
 	return mu;
 }
 
@@ -288,35 +284,7 @@ CudaDeviceFunction void InitFromFieldsStage()
 }
 
 CudaDeviceFunction void specialCases_Init()
-{
-	#ifdef OPTIONS_thermo
-		Temp   = T_init;
-		if (fabs(dT) > 0){
-			Temp = T_init + dT*Y;
-		}
-		if (fabs(dTx) > 0){
-			Temp = T_init + dTx*X;
-		}
-		#ifdef OPTIONS_planarBenchmark
-			if ( (NodeType & NODE_ADDITIONALS) == NODE_BWall) { //bottom wall
-				real_t x, omega;
-				x = (X-0.5) - myL;
-				omega = 3.1415926535897 / myL;
-				Temp = T_h + T_0 * cos(omega * x);
-				printf("y,x=%.4lf,%.4lf\n", Y,x);
-			} else if ( (NodeType & NODE_ADDITIONALS) == NODE_TWall) {
-				Temp = T_c;
-				printf("y,x=%.4lf,%.4lf\n", Y, X);
-			}
-			PhaseF = 0.5 + PLUSMINUS * (0.5) * tanh( (Y - MIDPOINT)/(IntWidth/2) );
-		#endif
-		if (surfPower > 1) {
-			SurfaceTension = sigma + sigma_TT*pow((Temp(0,0,0) - T_ref),surfPower) * (1.0/surfPower);
-		} else {
-			SurfaceTension = sigma + sigma_T*(Temp(0,0,0) - T_ref);
-		}
-		Cond = interp(PhaseF, k_h, k_l);
-	#endif											
+{								
 	// Pre-defined Initialisation patterns:
 	// Diffuse interface sphere
 	   // BubbleType = -1 refers to light fluid bubble.
@@ -418,10 +386,6 @@ CudaDeviceFunction void Init_distributions()
 	?>
         pnorm  = <?R C(sum(g)) ?>;
 		PhaseF = <?R C(sum(h)) ?>;
-	#ifdef OPTIONS_thermo
-		Temp   = Temp(0,0,0);
-		Cond   = interp(PhaseF, k_h, k_l);
-	#endif
 	#ifdef OPTIONS_OutFlow
 			if ((NodeType & NODE_BOUNDARY) == NODE_EConvect){
 					<?R if (Options$OutFlow){
@@ -526,38 +490,9 @@ CudaDeviceFunction void calc_Fb(real_t *fx, real_t *fy, real_t *fz, real_t rho){
 }
 
 CudaDeviceFunction void calc_Fs(real_t *fx, real_t *fy, real_t *fz, real_t mu, vector_t gPhi){
-	#ifdef OPTIONS_thermo
-		Temp = Temp(0,0,0);
-		SurfaceTension = SurfaceTension(0,0,0);
-		Cond = Cond(0,0,0);
-		real_t tmpSig, delta_s, dotTMP, magnPhi, magnPhi2;
-		magnPhi = sqrt(gPhi.x*gPhi.x + gPhi.y*gPhi.y + gPhi.z*gPhi.z);
-		magnPhi2 = magnPhi*magnPhi;
-		vector_t gradT;
-		<?R
-			IsotropicGrad('gradT','Temp')
-		?>
-		dotTMP = dotProduct(gradT,gPhi);
-		if (surfPower < 2) {
-			delta_s = 1.5*IntWidth*sigma_T;
-			*fx = mu * gPhi.x + delta_s*( magnPhi2*gradT.x - dotTMP*gPhi.x );
-			*fy = mu * gPhi.y + delta_s*( magnPhi2*gradT.y - dotTMP*gPhi.y );
-			*fz = mu * gPhi.z + delta_s*( magnPhi2*gradT.z - dotTMP*gPhi.z );
-		} else {
-			vector_t gradSig;
-			<?R
-				IsotropicGrad('gradSig','SurfaceTension')
-			?>
-			delta_s = 1.5*IntWidth;
-			*fx = mu * gPhi.x + delta_s*gradSig.x*( magnPhi2*gradT.x - dotTMP*gPhi.x );
-			*fy = mu * gPhi.y + delta_s*gradSig.y*( magnPhi2*gradT.y - dotTMP*gPhi.y );
-			*fz = mu * gPhi.z + delta_s*gradSig.z*( magnPhi2*gradT.z - dotTMP*gPhi.z );			
-		}
-	#else
-		*fx = mu * gPhi.x;
-		*fy = mu * gPhi.y;
-		*fz = mu * gPhi.z;
-	#endif
+	*fx = mu * gPhi.x;
+	*fy = mu * gPhi.y;
+	*fz = mu * gPhi.z;
 }
 
 #ifndef OPTIONS_BGK
@@ -836,13 +771,6 @@ CudaDeviceFunction void updateTrackers(real_t C){
 			}
 	}
 }
-//#############//
-//######THERMOCAPILLARY UPDATE######//
-#ifdef OPTIONS_thermo
-	<?RT models/multiphase/d3q27_pf_velocity/thermo.c.Rt ?>
-#endif
-//#############//
-
 
 #ifdef OPTIONS_BGK
 CudaDeviceFunction void CollisionBGK()

models/multiphase/d3q27_pf_velocity/conf.mk

@@ -1,6 +1,6 @@
 ADJOINT=0
 TEST=FALSE
-OPT="(q27 + OutFlow  + BGK + thermo*planarBenchmark)*autosym*geometric*staircaseimp*isograd*tprec"
+OPT="(q27 + OutFlow  + BGK)*autosym*geometric*staircaseimp*isograd*tprec"
 # q27 - Q27 lattice structure for phasefield
 #
 # OutFlow - include extra velocity stencil for outflowing boundaries