Write a new force in Navier-Stokes solver

We give the main stages to add a new force in the Navier-Stokes solver. You must modify three files: LBMScheme_NS_AC_Comp.h, Index_NS_AC_Comp.h, Models_NS_AC_Comp.h

File LBMScheme_NS_AC_Comp.h

Declare the force with RVect<dim> ForceMA ;

RVect<dim> ForceMA ;
     if (Model.Marangoni == 1) {
          ForceMA      = Model.force_MA<dim>(lbmState);
     }
     else if(Model.Marangoni == 0) {
             ForceMA[IX] = 0.0 ;
             ForceMA[IY] = 0.0 ;
     }

where Marangoni is an integer which is read in the input file and defined in Model. Next each component of that new force is added in the total force ForceTot

ForceTot[IX] = (1.0 - psi)*(ForceG[IX]+ForceP[IX]+ForceTS[IX]+ForceV[IX]+ForceMA[IX]) + ForcePsolid[IX] + ForcePenal[IX];
     ForceTot[IY] = (1.0 - psi)*(ForceG[IY]+ForceP[IY]+ForceTS[IY]+ForceV[IY]+ForceMA[IY]) + ForcePsolid[IY] + ForcePenal[IY];

Those components can be saved for visualization with paraview

this->set_lbm_val(IJK, IFMARAX, ForceMA[IX]);
this->set_lbm_val(IJK, IFMARAY, ForceMA[IY]);

where

• IJK is the current node

• IFMARAX and IFMARAY are indices

• ForceMA[IX] and ForceMA[IY] are the values of each component of the Marangoni force.

File Index_NS_AC_Comp.h

Those indices must be declared in the enumeration of file Index_NS_AC_Comp.h:

enum ComponentIndex {
...
   IFMARAX, /*!< Force Marangoni X */
   IFMARAY, /*!< Force Marangoni Y */
   ...
}

File Models_NS_AC_Comp.h

1. Read and declare parameters of Marangoni and surfactant

   Read the parameters in the input file

   // Read parameters for Marangoni force
   Marangoni       = configMap.getInteger("params_marangoni", "marangoni_force", 0  );
        Closure_Model   = configMap.getInteger("params_marangoni", "Closure_model"  , 0  ); //1 for logarithm, 0 for linear, default linear
        //param for linear law
   sigma_marangoni = configMap.getFloat  ("params_marangoni", "sigma_marangoni", 0.0);
   dsigmadcomp     = configMap.getFloat  ("params_marangoni", "dsigmadcomp"    , 0.0);
        //param for log law
   beta_log        = configMap.getFloat  ("params_marangoni","beta_log"        , 0.0);
   C_infinity      = configMap.getFloat  ("params_marangoni","C_infinity"      , 1.0); //maximum concentration
   

   Declare the variables

   //! Input parameters for surfactant dynamics
   real_t beta_surf,k_surf,eps_surf;
   //! Input Parameters for Marangoni force
   int Marangoni,Closure_Model;
   real_t sigma_marangoni, dsigmadcomp, beta_log, C_infinity;
   

2. Define function force_MA<dim>(lbmState);

   template <int dim>
   KOKKOS_INLINE_FUNCTION RVect<dim> force_MA(const LBMState &lbmState) const {
   
        //according to Hanyang Mo et al. (2023)
   
        if (dim == 2) {
                const real_t grad_phi_sqr = (SQR(lbmState[IDPHIDX])+SQR(lbmState[IDPHIDY])); //grad phi squared
                const real_t grad_sigmax = sigma_prime(lbmState[IC])*lbmState[IDCDX]; //d(sigma)/dx = d(sigma)/dC * d(C)/dx
                const real_t grad_sigmay = sigma_prime(lbmState[IC])*lbmState[IDCDY]; //d(sigma)/dy = d(sigma)/dC * d(C)/dy
   
                const real_t dot_prod = lbmState[IDPHIDX]*grad_sigmax + lbmState[IDPHIDY]*grad_sigmay ;// grad sigma . grad phi
   
                RVect<dim> term;
                term[IX] = 1.5*W*(grad_phi_sqr*grad_sigmax - dot_prod*lbmState[IDPHIDX]);
                term[IY] = 1.5*W*(grad_phi_sqr*grad_sigmay - dot_prod*lbmState[IDPHIDY]);
   
                return term;
           }
   }
   

3. Define \(d\sigma/dc\)

   KOKKOS_INLINE_FUNCTION real_t sigma_func(real_t comp ) const {
        if(Closure_Model == 1){
                     return sigma_marangoni*(1+beta_log*log(1.0-comp/C_infinity));
             }
   
        else{
                     return sigma_marangoni + comp*dsigmadcomp;
             }
   }
   

   KOKKOS_INLINE_FUNCTION real_t sigma_prime(real_t comp) const {
   if(Closure_Model == 1){
                     return -sigma_marangoni*beta_log/(C_infinity-comp);
             }
   else{
                return dsigmadcomp;
             }
   }
   

Section author: Alain Cartalade