sandbox/ecipriano/run/expansion.c

    Thermal Expansion of a Liquid Droplet

    In this test case, we simulate the thermal expansion of an initially cold n-heptane droplet in a hot isothermal environment, neglecting the phase change. The aim of this simulation is to test the convergence of the volumetric source term due to density changes, using the multicomponent.h solver with variable properties.

    Evolution of the temperature field

    Phase Change Setup

    We define the number of gas and liquid species, and we set the initial composition of the two phases, which does not change in time. By setting the thermodynamic equilibrium Keq to 0 we suppress the phase change.

    #define NGS 2
    #define NLS 1
    
    char* gas_species[NGS] = {"NC7H16", "N2"};
    char* liq_species[NLS] = {"NC7H16"};
    char* inert_species[1] = {"N2"};
    double gas_start[NGS] = {0., 1.};
    double liq_start[NLS] = {1.};
    double inKeq[NLS] = {0.};

    We declare all the properties necessary for the multicomponent phase change model. The properties are set to null values because they are overwritten by the variable-properties formulation, which computes all the physical properties as a function of the thermodynamic state of the mixture.

    double inDmix1[NLS] = {0.};
    double inDmix2[NGS] = {0., 0.};
    double lambda1 = 0.;
    double lambda2 = 0.;
    double dhev = 0.;
    double cp1 = 0.;
    double cp2 = 0.;

    We set the initial temperature of the liquid and of the gas phase. The gas phase temperature can be changed by the user through a compilation variable, for example by setting -DTEMPERATURE=400

    #ifndef TEMPERATURE
    # define TEMPERATURE 350
    #endif
    double TL0 = 300.;
    double TG0 = TEMPERATURE;

    We solve the temperature field, and we reduce the tolerance for the calculation of the variable properties. The interfacial temperature is not computed from the jump conditon, it is just set to the gas phase temperature value, in order to keep the gas phase temperature constant. The advection term of the momentum equation is solved in non-conservative form by setting the variable NO_ADVECTION_DIV to 1.

    #define SOLVE_TEMPERATURE
    #define NO_ADVECTION_DIV 1
    #define T_PROP 1.e-6
    #define FIXED_INTERFACE_TEMPERATURE TG0

    Simulation Setup

    In this simulation, the choice of the Navier-Stokes equations solver is not important, because there is no velocity jump due to the phase change. We can choose between the centered solver with expansion source term, and the centered solver with velocity jump. The two-phase.h solver is extended to variable physical properties by including a policy for the calculation of such properties. In this case, we use correlations implemented in OpenSMOKE++ (refer to Cipriano et al., 2024 for details). The solution of the temperature field is performed by the multicomponent phase change model.

    #include "axi.h"
    #if JUMP
    # include "navier-stokes/velocity-jump.h"
    #else
    # include "navier-stokes/centered-evaporation.h"
    # include "navier-stokes/centered-doubled.h"
    #endif
    #include "opensmoke-properties.h"
    #include "two-phase.h"
    #include "tension.h"
    #include "evaporation.h"
    #include "multicomponent-varprop.h"
    #include "view.h"

    Boundary conditions

    We initialize the droplet at the lower left corner of the domain, and we exploit axial symmetry. We set symmetry boundary conditions on the sides in contact with the droplet (left, bottom), and outflow boundary conditions on the other sides of the domain.

    #if JUMP
    u1.n[top] = neumann (0.);
    u1.t[top] = neumann (0.);
    u2.n[top] = neumann (0.);
    u2.t[top] = neumann (0.);
    p[top] = dirichlet (0.);
    ps[top] = dirichlet (0.);
    pg[top] = dirichlet (0.);
    
    u1.n[right] = neumann (0.);
    u1.t[right] = neumann (0.);
    u2.n[right] = neumann (0.);
    u2.t[right] = neumann (0.);
    p[right] = dirichlet (0.);
    ps[right] = dirichlet (0.);
    pg[right] = dirichlet (0.);
    #else
    u.n[top] = neumann (0.);
    u.t[top] = neumann (0.);
    p[top] = dirichlet (0.);
    uext.n[top] = neumann (0.);
    uext.t[top] = neumann (0.);
    pext[top] = dirichlet (0.);
    
    u.n[right] = neumann (0.);
    u.t[right] = neumann (0.);
    p[right] = dirichlet (0.);
    uext.n[right] = neumann (0.);
    uext.t[right] = neumann (0.);
    pext[right] = dirichlet (0.);
    #endif

    Simulation Data

    We declare the maximum and minimum levels of refinement, the initial diameter (1 mm), and additional data for post-processing.

    int maxlevel, minlevel = 2;
    double D0 = 1e-3, effective_radius0, d_over_d02 = 1.;
    double volume0 = 0., rho0 = 0., rhof = 0.;
    
    int main (void) {

    We set the kinetics folder, which defines the species of the simulation, and it is used by OpenSMOKE++ for the calculation of the thermodynamic and transport properties. The path is relative to OpenSMOKEppInterface/kinetics.

      kinfolder = "evaporation/n-heptane-in-nitrogen";

    We set additional simulation properties. Be careful, these are “dummy” properties which are overwritten by the variable-properties formulation already at the first iteration. To start the simulation without any problem of division by zero, the viscosity value should be non-null.

      rho1 = 1.; rho2 = 1.;
      mu1 = 1.; mu2 = 1.;

    The thermodynamic pressure is set to 10 atm, we consider a Pref value in SI units (Pa).

      Pref = 10*101325.;

    We change the dimensions of the domain as a function of the initial diameter of the droplet.

      L0 = 1.5*D0;

    We use a constant surface tension value.

      f.sigma = 0.03;

    We run the simulation at different maximum levels of refinement.

      for (maxlevel = 6; maxlevel <= 6; maxlevel++) {
        init_grid (1 << maxlevel);
        run();
      }
    }
    
    #define circle(x,y,R)(sq(R) - sq(x) - sq(y))

    We initialize the volume fraction field according to D0.

    event init (i = 0) {
      fraction (f, circle (x, y, 0.5*D0));

    We compute initial variables useful for post-processing.

      effective_radius0 = pow (3.*statsf(f).sum, 1./3.);
      volume0 = 4./3.*pi*pow (effective_radius0, 3.);
      
      foreach (reduction(+:mLiq0))
        mLiq0 += rho1v[]*f[]*dv();

    We set the molecular weights of the chemial species involved in the simulation (by default inMW=1).

      for (int jj=0; jj<NGS; jj++)
        inMW[jj] = OpenSMOKE_MW (jj);

    We compute variables for the analytical steady-state solution.

      ThermoState ts0;
      ts0.T = TL0;
      ts0.P = Pref;
      ts0.x = (double[]){1.};
    
      ThermoState tsf;
      tsf.T = TG0;
      tsf.P = Pref;
      tsf.x = (double[]){1.};
    
      rho0 = tp1.rhov (&ts0);
      rhof = tp1.rhov (&tsf);
    }

    On the top and right sides we set Dirichlet boundary conditions for the temperature and mass fraction fields.

    event bcs (i = 0) {
      scalar NC7H16 = YGList[OpenSMOKE_IndexOfSpecies ("NC7H16")];
      scalar N2    = YGList[OpenSMOKE_IndexOfSpecies ("N2")];
    
      NC7H16[top] = dirichlet (0.);
      NC7H16[right] = dirichlet (0.);
    
      N2[top] = dirichlet (1.);
      N2[right] = dirichlet (1.);
    
      TG[top] = dirichlet (TG0);
      TG[right] = dirichlet (TG0);
    }

    We adapt the grid according to the temperature and the velocity fields.

    #if TREE
    event adapt (i++) {
      adapt_wavelet_leave_interface ({T,u.x,u.y}, {f},
          (double[]){1.e-2,1.e-1,1.e-1}, maxlevel, minlevel, 1);
    }
    #endif

    Post-Processing

    The following lines of code are for post-processing purposes.

    Output Files

    We write on a file squared diameter, the exact droplet diameter at steady-state, and the ratio between the current and the initial liquid mass values.

    event output_data (i++) {
      char name[80];
      sprintf (name, "OutputData-%d", maxlevel);
      static FILE * fp = fopen (name, "w");
    
      double effective_radius = pow (3.*statsf(f).sum, 1./3.);
      d_over_d02 = sq (effective_radius / effective_radius0);
    
      double mLiq = 0.;
      foreach(reduction(+:mLiq))
        mLiq += rho1v[]*f[]*dv();
    
      double exact_volume = volume0*rho0/rhof;
      double exact_radius = pow (3./4./pi*exact_volume, 1./3.);
      double exact_d_over_d02 = sq (exact_radius/effective_radius0);
      double relerr = fabs (effective_radius - exact_radius)/effective_radius;
    
      fprintf (fp, "%g %g %g %g %g %g %g\n", t, t/sq(D0*1e3), effective_radius,
          d_over_d02, mLiq/mLiq0, relerr, exact_d_over_d02);
    }

    Movie

    We write images with the temperature and the divergence of the liquid velocity at time 0.1 s.

    event pictures (t = 0.1) {
      clear();
      view (ty = -0.5);
      draw_vof ("f", filled = -1, fc = {1.,1.,1.});
      squares ("T", min = TL0, max = TG0, linear = true);
      isoline ("T", n = 12., min = TL0, max = TG0);
      mirror ({1.,0.}) {
        draw_vof ("f", filled = -1, fc = {1.,1.,1.});
        squares ("T", min = TL0, max = TG0, linear = true);
        isoline ("T", n = 12., min = TL0, max = TG0);
      }
      save ("temperature.png");
    
      scalar drhodtplot[];
      foreach()
        drhodtplot[] = drhodtext[]*f[];
    
      clear();
      view (ty = -0.5);
      draw_vof ("f", filled = -1, fc = {1.,1.,1.});
      squares ("drhodtplot", spread = -1, linear = true);
      isoline ("drhodtplot", n = 12.,
          min = statsf(drhodtplot).min, max = statsf(drhodtplot).max);
      mirror ({1.,0.}) {
        draw_vof ("f", filled = -1, fc = {1.,1.,1.});
        squares ("drhodtplot", spread = -1, linear = true);
        isoline ("drhodtplot", n = 12.,
            min = statsf(drhodtplot).min, max = statsf(drhodtplot).max);
      }
      save ("divergence.png");
    }

    We write a video with the evolution of the temperature field.

    event movie (t += 0.01; t <= 3) {
      if (maxlevel == 6) {
        clear();
        box();
        view (tx = -0.5, ty = -0.5);
        draw_vof ("f");
        squares ("T", min = TL0, max = statsf(T).max, linear = true);
        save ("movie.mp4");
      }
    }

    Results

    set border linewidth 2
    set grid
    set key bottom right
    set xlabel "t/D_0^2 [s/mm^2]"
    set ylabel "(D/D_0)^2 [-]"
    set xrange[0:3]
    set yrange[1:1.12]
    set size square
    
    basedir350 = "../expansion-T350/"
    basedir375 = "../expansion-T375/"
    basedir400 = "../expansion-T400/"
    
    set label "∆T = 50 K"  at 2.45,1.05  left font ",11" tc rgb "black"
    set label "∆T = 75 K"  at 2.45,1.078 left font ",11" tc rgb "black"
    set label "∆T = 100 K" at 2.45,1.11  left font ",11" tc rgb "black"
    
    plot basedir350."OutputData-7" u 2:7 every 50000 w p ps 0.8 lc 1 pt 6 dt 2 t "Exact", \
         basedir350."OutputData-5" u 2:4 w l lw 2 lc 1 dt 4 t "LEVEL 5", \
         basedir350."OutputData-6" u 2:4 w l lw 2 lc 1 dt 3 t "LEVEL 6", \
         basedir350."OutputData-7" u 2:4 w l lw 2 lc 1 dt 2 t "LEVEL 7", \
         basedir350."OutputData-8" u 2:4 w l lw 2 lc 1 dt 1 t "LEVEL 8", \
         basedir375."OutputData-7" u 2:7 every 50000 w p ps 0.8 lc 2 pt 6 dt 2 notitle, \
         basedir375."OutputData-5" u 2:4 w l lw 2 lc 2 dt 4 notitle, \
         basedir375."OutputData-6" u 2:4 w l lw 2 lc 2 dt 3 notitle, \
         basedir375."OutputData-7" u 2:4 w l lw 2 lc 2 dt 2 notitle, \
         basedir375."OutputData-8" u 2:4 w l lw 2 lc 2 dt 1 notitle, \
         basedir400."OutputData-7" u 2:7 every 50000 w p ps 0.8 lc 3 pt 6 dt 2 notitle, \
         basedir400."OutputData-5" u 2:4 w l lw 2 lc 3 dt 4 notitle, \
         basedir400."OutputData-6" u 2:4 w l lw 2 lc 3 dt 3 notitle, \
         basedir400."OutputData-7" u 2:4 w l lw 2 lc 3 dt 2 notitle, \
         basedir400."OutputData-8" u 2:4 w l lw 2 lc 3 dt 1 notitle
    Expansion of the Square Diameter (script)

    Expansion of the Square Diameter (script)

    reset
    set border linewidth 2
    
    basedir = "../expansion-T350/"
    stats basedir."OutputData-5" using (last5=$6) nooutput
    stats basedir."OutputData-6" using (last6=$6) nooutput
    stats basedir."OutputData-7" using (last7=$6) nooutput
    stats basedir."OutputData-8" using (last8=$6) nooutput
    
    set print "Errors-350.csv"
    print sprintf ("%d %.12f", 2**5, last5)
    print sprintf ("%d %.12f", 2**6, last6)
    print sprintf ("%d %.12f", 2**7, last7)
    print sprintf ("%d %.12f", 2**8, last8)
    unset print
    
    basedir = "../expansion-T375/"
    stats basedir."OutputData-5" using (last5=$6) nooutput
    stats basedir."OutputData-6" using (last6=$6) nooutput
    stats basedir."OutputData-7" using (last7=$6) nooutput
    stats basedir."OutputData-8" using (last8=$6) nooutput
    
    set print "Errors-375.csv"
    print sprintf ("%d %.12f", 2**5, last5)
    print sprintf ("%d %.12f", 2**6, last6)
    print sprintf ("%d %.12f", 2**7, last7)
    print sprintf ("%d %.12f", 2**8, last8)
    unset print
    
    basedir = "../expansion-T400/"
    stats basedir."OutputData-5" using (last5=$6) nooutput
    stats basedir."OutputData-6" using (last6=$6) nooutput
    stats basedir."OutputData-7" using (last7=$6) nooutput
    stats basedir."OutputData-8" using (last8=$6) nooutput
    
    set print "Errors-400.csv"
    print sprintf ("%d %.12f", 2**5, last5)
    print sprintf ("%d %.12f", 2**6, last6)
    print sprintf ("%d %.12f", 2**7, last7)
    print sprintf ("%d %.12f", 2**8, last8)
    unset print
    
    set xlabel "N"
    set ylabel "Relative Error"
    
    set logscale x 2
    set logscale y
    
    set xr[2**4:2**9]
    
    set size square
    set grid
    
    f1(x) = a1*x**(-b1)
    f2(x) = a2*x**(-b2)
    f3(x) = a3*x**(-b3)
    
    fit f1(x) "Errors-350.csv" via a1,b1
    fit f2(x) "Errors-375.csv" via a2,b2
    fit f3(x) "Errors-400.csv" via a3,b3
    
    ftitle(a,b) = sprintf("%.2fx^{-%.3f}", a, b)
    
    plot "Errors-350.csv" w p pt 1 ps 2 lc 1 t "∆T =  50 K", \
         "Errors-375.csv" w p pt 2 ps 2 lc 2 t "∆T =  75 K", \
         "Errors-400.csv" w p pt 3 ps 2 lc 3 t "∆T = 100 K", \
          f1(x) w l lw 2 lc 1 t ftitle (a1, b1), \
          f2(x) w l lw 2 lc 2 t ftitle (a2, b2), \
          f3(x) w l lw 2 lc 3 t ftitle (a3, b3)
    Relative Errors (script)

    Relative Errors (script)