sandbox/Antoonvh/smoke.c

An atmospheric boundary-layer filled with a smoke cloud.

We simulate a boudary layer that is filled with a so-called smoke cloud. This cloud cools at the top due to the emission of longwave radiation. The case was introduced as an intercomparison case by Bretherton et al. (1999). This page shows the implementation of the case in Basilisk using an LES formulation. However, on this page we run the case in 2 dimensions to that we can view the results.

Set-up

In order to invoke the LES formulation we include the Navier-stokes solver and the SGS headerfile.

//#include "grid/octree.h" // <- uncomment for 3D as was used in Van Hooft et al. (2018)
#include "navier-stokes/centered.h"
#include "SGS.h"

The maximum and minimum levels of refinment are set. We also update the radiation tendencies once every 5 timesteps. This way we limit the effect of the fact that the current implementation of the tendency term due to radiation is not very efficient and should be improved in the future. Also we set some other global variables that according to the case description by Bretherton et al. (1999).

int maxlevel=8,minlevel=5;
int Radupdate = 5;
double T0 = 291.5, F0 = 60., p0=1.1436, Ka=0.02, Cp=1004.;

double temp = 675*sq(dimension+1)/(4-dimension); // ~1 hour in 2D or 3 hour in 3D

We declare fields for the smoke cloud fraction (S), the (virtual) potential temperature (T), the radiative flux field (F), and the field that stores the temperature tendency (co). We also tell the solver that we want to advect S and T with the flow and diffuse it according to the sub-grid scale model. Finally, an helper acceleration field av is declared that helps with the implementation of the ‘gravity’.

scalar S[],T[],F[],co[];
scalar * tracers = {T,S}; // Temperature and smoke fraction get both advected and diffused according to the SGS model.
face vector av[];

int main(){
  F0=F0/sq((4-dimension));
  L0 = 3200;
  a=av; 
  X0=Y0=Z0=0.;
  init_grid(1<<(minlevel));
  co.refine=refine_linear;
  run();
}

Initialization

The fields T and S field are initialized with an inversion at z=700m (i.e. y=700 in our simulation) with a width of 25 metres.

event init (t = 0){
  DT=0.01;
  TOLERANCE=1e-8; 
  for (scalar s in tracers)
    s.gradient=minmod2;
  foreach(){
    S[] = 0.5+0.5*tanh((700-y)/12.5);
    T[] = T0+3.5*tanh((y-700)/12.5)+(0.1*noise()*(y<700));
    T[] += 0.001*(y-712.5)*(y>712.5);  
  }
  boundary(all);
  while(adapt_wavelet({S},(double[]){0.05}, maxlevel,minlevel).nf>10){
    foreach(){
      S[] = 0.5+0.5*tanh((700-y)/12.5);
      T[] = T0+3.5*tanh((y-700)/12.5)+(0.01*noise()*(y<700));
      T[] += 0.001*(y-712.5)*(y>712.5);  
    }
    boundary(all);
  }
}

Gravity is implemented via a buoyancy formulation based on the virtual potential temperature.

event acceleration(i++){
  foreach_face(y) 
    av.y[]=9.81*((T[]+T[0,-1]-T0))/(2*T0);
}

Radiation

Bretherton et al. (1999) discribed how the radiation tendencies should be evaluated in a staggered-grid manner. We neglect that and do it in our own way so that we do not need to take proper care of face fluxes etc.

Based on the smoke cloud fraction field (S) we can evaluate the readiative fluxes (F), according to:

\displaystyle F = F_0 e^{-K_a S_l},

where F_0 and K_a are specified constants, and S_l is defined for a point at height z as follows:

\displaystyle S_l = \int_z^{z_{top}} \rho\ S\ dz,

where \rho is taken as the reference density \rho_0, and with our grid orientation and set-up; z_{top}=Y0+L0=L0.

The tendencies due to the radiation are applied each timestep. However the tendency field co, is only updated each Radupdate (i.e. 5) timesteps.

event radiation(i++){
  if (i%Radupdate==0){
    boundary({S});
    foreach(){  //Calculate for each grid cell, there is room to improve this!
      double sl = 0;
      double xp = x;
      double yp = y;
      double zp = z;
      while (yp<L0){ // This is the not smart part of evaluating the radiation scheme. 
        Point point = locate(xp,yp,zp); 
        yp=y;// jump to the found cell centered height.
        sl+= p0*interpolate(S,xp,yp,zp)*Delta;
        yp+= Delta/1.5;
      }
      F[] = F0*exp(-Ka*sl);
    }
    boundary({F});

We have found the flux, now we take the vertical derivative (1D radiative divergence) to obtain the tendency.

    foreach(){
      co[] =-(F[0,1]-F[0,-1])/(2*p0*Cp*Delta);
    }
  }

We copy the tendency into a helper field (cooling). We do this because the cooling field will be modified by the solver such that we cannot use it in case we do not update it next time step nor in the output, see here.

  boundary({co});
  scalar cooling[]; 
  foreach()
    cooling[]=co[];
  boundary({cooling});
  diffusion(T,dt,zerof,beta = cooling);// Do cooling, not diffusion!
}

Adaptation

We use the well-known wavelet-based technique for error estimation and perform the corresonding grid adaptation.

event adapt(i++){
  double eS = 0.1/(dimension-1);
  double eT = 0.7/(dimension-1);    
  adapt_wavelet({S,T},(double[]){eS,eT}, maxlevel,minlevel,all);
}

This back_to_defaults event is added so that we can start with a very small timestep (DT) and tolerance for the Poisson problems (TOLERANCE), so that the iterative solver has some low-impact iterations before it finds the not-consistently-initialized pressure field.

event back_to_defaults(i=50){
  DT=1.;
  CFL = 0.7;
  TOLERANCE=1e-3; 
}

2D output

In 2D we output three movies.

event mov(t+=10;t<=temp){
  if (dimension==2){
    double mc = 4./(3600*sq(4-dimension)); //Cooling scale in K/sec
    scalar lev[];
    foreach()
      lev[]=level;
    boundary({S,T}); //Linear interp for movie
    static FILE * fpg = popen("ppm2mp4 S.mp4","w"); 
    output_ppm(S,fpg,linear=true,n=pow(2,maxlevel),min=0,max=1,map=gray);
    
    static FILE * fpc = popen("ppm2mp4 c.mp4","w"); 
    output_ppm(co,fpc,linear=true,n=pow(2,maxlevel),min=-mc,max=mc);
    
    static FILE * fpl = popen("ppm2mp4 l.mp4","w"); 
    output_ppm(lev,fpl,n=pow(2,maxlevel),min=minlevel,max=maxlevel);
  }
  else if (dimension==3){
    //no movie, no nothing
    
  }
  fprintf(ferr,"%d %g\n",i,t);
}

results

We can view the movies. First we choose to show the smoke-cloud fraction field.

We can observe the emergence of subsiding shells in the boundary layer. Next we have a look at the process that drives the flow; the radiative cooling:

Looks OK. Notice that the depth of the cooling layer is controlled by the “optical thickness” parameter Ka that was tuned to help the participating models by supressing the scale separation. Finally we look at the evolution of the grid’s level of refinement:

Well done adapt_wavelet function!

Reference:

Bretherton, Christopher S., et al. “An intercomparison of radiatively driven entrainment and turbulence in a smoke cloud, as simulated by different numerical models.” Quarterly Journal of the Royal Meteorological Society 125.554 (1999): 391-423.

Read,print,copy and obtain it via this link.