Lagrangian subgrid bubble model

docs/documentation/case.md

@@ -799,6 +799,58 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 
 - `bc_y%beg = -2` to enable reflective boundary conditions
 
+### 15. Lagrangian subgrid bubble model
+
+| Parameter                 | Type    | Description                                               |
+| ---:                      | :---:   | :---                                                      |
+| `lag_bubbles`             | Logical | Lagrangian subgrid bubble model switch                    |
+| `lag_solver_approach`     | Integer | 1: One-way coupling, 2: two-way coupling                  |
+| `lag_bubble_model`        | Integer | Bubble dynamics model. 1: Keller-Miksis equation          |
+| `lag_cluster_type`        | Integer | Method to find p_inf                                      |
+| `lag_pressure_corrector`  | Logical | Cell pressure correction term                             |
+| `lag_adap_dt`             | Logical | Activates the adaptive rkck time stepping algorithm       |
+| `lag_smooth_type`         | Integer | Smoothing function. 1: Gaussian, 2:Delta 3x3              |
+| `lag_heatTransfer_model`  | Logical | Activates the interface heat transfer model               |
+| `lag_massTransfer_model`  | Logical | Activates the interface mass transfer model               |
+| `lag_write_bubbles`       | Logical | Write files to track the bubble evolution each time step  |
+| `lag_write_bubble_stats`  | Logical | Write the maximum and minimum radius of each bubble       |
+| `lag_nBubs_glb`           | Integer | Global number of bubbles                                  |
+| `lag_epsilonb`            | Real    | Standard deviation scaling for the gaussian function      |
+| `lag_rkck_tolerance`      | Real    | Adaptive RKCK tolerance                                   |
+| `lag_charwidth`           | Real    | Domain virtual depth (z direction, for 2D simulations)    |
+| `lag_valmaxvoid`          | Real    | Maximum void fraction permitted                           |
+| `csonhost`                | Real    | Liquid speed of sound                                     |
+| `vischost`                | Real    | Liquid viscosity                                          |
+| `Thost`                   | Real    | Liquid temperature                                        |
+| `gammagas`                | Real    | Heat capacity ratio of the gas                            |
+| `gammavapor`              | Real    | Heat capacity ratio of the vapor                          |
+| `pvap`                    | Real    | Vapor pressure at the reference temperature               |
+| `cpgas`                   | Real    | Specific heat capacity of the gas                         |
+| `cpvapor`                 | Real    | Specific heat capacity of the vapor                       |
+| `kgas`                    | Real    | Thermal conductivity of the gas                           |
+| `kvapor`                  | Real    | Thermal conductivity of the vapor                         |
+| `Rgas`                    | Real    | Specific gas constant od the gas                          |
+| `Rvap`                    | Real    | Specific gas constant od the vapor                        |
+| `diffcoefvap`             | Real    | Vapor diffusivity in the gas                              |
+| `sigmabubble`             | Real    | Surface tension                                           |
+
+- `lag_solver_approach` Specifies the Euler-Lagrange coupling method: [1] enables a one-way coupling approach, where the bubbles do not influence the Eulerian field. [2] activates the two-way coupling approach based on [Maeda and Colonius (2018)](references.md#Maeda18), where the effect of the bubbles is added in the Eulerian field as source terms.
+
+- `lag_bubble_model` [1] Activates the Keller-Miksis equation to model the bubble dynamics.
+
+- `lag_cluster_type` Specifies method to find p_inf (pressure that drives the bubble dynamics): [1] activates the bilinear interpolation of the pressure field, while [2] enables the bubble dynamic closure based on [Maeda and Colonius (2018)](references.md#Maeda18), the full model is obtained when 'lag_pressure_corrector' is true.
+
+- `lag_adap_dt` Activates the adaptive 4th/5th order Runge—Kutta–Cash–Karp (RKCK) time-stepping algorithm. A maximum error between the 4th and 5th order Runge-Kutta-Cash-Karp solutions for the same time step size is calculated. If the error is smaller than a tolerance ('lag_rkck_tolerance'), then the algorithm employs the 5th order solution, while if not, both eulerian/lagrangian variables are re-calculated with a smaller time step size.
+
+- `lag_smooth_type` Specifies the smoothening method of projecting the lagrangian bubbles in the Eulerian field: [1] activates the gaussian kernel function described in  [Maeda and Colonius (2018)](references.md#Maeda18), while [2] activates the delta kernel function where the effect of the bubble is only seen in the specific bubble location cell.
+
+- `lag_heatTransfer_model` Activates the heat transfer model at the bubble's interface based on ([Preston et al., 2007](references.md#Preston07)).
+
+- `lag_massTransfer_model` Activates the mass transfer model at the bubble's interface based on ([Preston et al., 2007](references.md#Preston07)).
+
+- `lag_nBubs_glb` Total number of bubbles. Their initial conditions need to be specified in the ./input/lag_bubbles.dat file. See the example cases for additional information.
+
+
 ## Enumerations
 
 ### Boundary conditions

docs/documentation/references.md

@@ -30,6 +30,8 @@
 
 - <a id="Maeda17">Maeda, K. and Colonius, T. (2017). A source term approach for generation of one-way acoustic waves in the euler and navier–stokes equations. Wave Motion, 75:36–49.</a>
 
+- <a id="Maeda18">Maeda, K. and Colonius, T. (2018). Eulerian–lagrangian method for simulation of cloud cavitation. Journal of computational physics, 371:994–1017.</a>
+
 - <a id="Meng16">Meng, J. C. C. (2016). Numerical simulations of droplet aerobreakup. PhD thesis, California Institute of Technology.</a>
 
 - <a id="Pirozzoli13">Pirozzoli, S., and Colonius, T. (2013). Generalized characteristic relaxation boundary conditions for unsteady compressible flow simulations. Journal of Computational Physics, 248:109-126.</a>

examples/3D_lag_bubbles_bubblescreen/case.py

@@ -0,0 +1,187 @@
+#!/usr/bin/env python3
+
+import math
+import json
+
+# Bubble screen
+# Description: A planar acoustic wave interacts with a bubble cloud 
+# in water. The background field is modeled in using an Eulerian framework, 
+# while the bubbles are tracked using a Lagrangian framework.
+
+# Reference values for nondimensionalization 
+x0 = 1.e-03         # length - m
+rho0 = 1.e+03       # density - kg/m3
+c0 = 1475.          # speed of sound - m/s
+p0 = rho0*c0*c0     # pressure - Pa
+T0 = 298            # temperature - K
+
+# Host properties (water)
+gamma_host  = 2.7466     # Specific heat ratio
+pi_inf_host = 792.02e+06 # Stiffness - Pa
+mu_host = 1e-3           # Dynamic viscosity - Pa.s
+c_host = 1475.           # speed of sound - m/s
+rho_host = 1000          # density kg/m3
+T_host = 298             # temperature K
+
+# Lagrangian bubbles' properties
+R_uni = 8314        # Universal gas constant - J/kmol/K
+MW_g = 28.0         # Molar weight of the gas - kg/kmol
+MW_v = 18.0         # Molar weight of the vapor - kg/kmol
+gamma_g = 1.4       # Specific heat ratio of the gas
+gamma_v = 1.333     # Specific heat ratio of the vapor
+pv = 2350           # Vapor pressure of the host - Pa
+cp_g = 1.e3         # Specific heat of the gas - J/kg/K
+cp_v = 2.1e3        # Specific heat of the vapor - J/kg/K
+k_g = 0.025         # Thermal conductivity of the gas - W/m/K
+k_v = 0.02          # Thermal conductivity of the vapor - W/m/K
+diffVapor = 2.5e-5  # Diffusivity coefficient of the vapor - m2/s
+sigBubble = 0.069   # Surface tension of the bubble - N/m
+
+# Acoustic source properties
+patm = 101325.      # Atmospheric pressure - Pa
+pamp = 1.e5         # Amplitude of the acoustic source - Pa
+freq = 300e+03      # Source frequency - Hz
+wlen = c_host/freq  # Wavelength - m
+
+# Domain and time set up
+
+xb = -12.e-3    # Domain boundaries - m (x direction)
+xe = 12.e-3
+yb = -2.5e-3    # Domain boundaries - m (y direction)
+ye = 2.5e-3
+zb = -2.5e-3    # Domain boundaries - m (z direction)
+ze = 2.5e-3
+
+Nx = 240        # number of elements into x direction
+Ny = 50         # number of elements into y direction
+Nz = 50         # number of elements into z direction
+
+dt = 7.5e-9     # constant time-step - sec
+
+# ==============================================================================
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================
+    'run_time_info'                : 'T',
+    # ==========================================================
+
+    # Computational Domain Parameters ==========================
+    'x_domain%beg'                 : xb/x0,
+    'x_domain%end'                 : xe/x0,
+    'y_domain%beg'                 : yb/x0,
+    'y_domain%end'                 : ye/x0,
+    'z_domain%beg'                 : zb/x0,
+    'z_domain%end'                 : ze/x0,
+    'stretch_z'                    : 'F',
+    'stretch_y'                    : 'F',
+    'stretch_x'                    : 'F',
+    'm'                            : Nx,
+    'n'                            : Ny,
+    'p'                            : Nz,
+    'dt'                           : dt*(c0/x0),
+    't_step_start'                 : 0,
+    't_step_stop'                  : 3000,
+    't_step_save'                  : 500,
+    # ==========================================================
+
+    # Simulation Algorithm Parameters ==========================
+    'model_eqns'                   : 2,
+    'num_fluids'                   : 1,
+    'num_patches'                  : 1,
+    'viscous'                      : 'T',
+    'mpp_lim'                      : 'F',
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1.0E-16,
+    'mapped_weno'                  :'T',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     :-6,
+    'bc_x%end'                     :-6,
+    'bc_y%beg'                     :-1,
+    'bc_y%end'                     :-1,
+    'bc_z%beg'                     :-1,
+    'bc_z%end'                     :-1,
+    # ==========================================================
+
+    # Acoustic source ==========================================
+    'acoustic_source'              : 'T',
+    'num_source'                   : 1,
+    'acoustic(1)%support'          : 3,
+    'acoustic(1)%pulse'            : 1,
+    'acoustic(1)%npulse'           : 1,
+    'acoustic(1)%mag'              : pamp/p0,
+    'acoustic(1)%wavelength'       : wlen/x0,
+    'acoustic(1)%length'           : 2*(ze-zb)/x0,
+    'acoustic(1)%height'           : 2*(ye-yb)/x0,
+    'acoustic(1)%loc(1)'           : -7.e-03/x0,
+    'acoustic(1)%loc(2)'           : 0.,
+    'acoustic(1)%loc(3)'           : 0.,
+    'acoustic(1)%dir'              : 0.,
+    'acoustic(1)%delay'            : 0.,
+    # ==========================================================
+
+    # Formatted Database Files Structure Parameters ============
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'parallel_io'                  :'T',
+    # ==========================================================
+
+    # Patch 1: Water (left) ====================================
+    'patch_icpp(1)%geometry'       : 9,
+    'patch_icpp(1)%x_centroid'     : 0.,
+    'patch_icpp(1)%y_centroid'     : 0.,
+    'patch_icpp(1)%z_centroid'     : 0.,
+    'patch_icpp(1)%length_x'       : 2*(xe-xb)/x0,
+    'patch_icpp(1)%length_y'       : 2*(ye-yb)/x0,
+    'patch_icpp(1)%length_z'       : 2*(ze-zb)/x0,
+    'patch_icpp(1)%vel(1)'         : 0.,
+    'patch_icpp(1)%vel(2)'         : 0.,
+    'patch_icpp(1)%vel(3)'         : 0.,
+    'patch_icpp(1)%pres'           : patm/p0,
+    'patch_icpp(1)%alpha_rho(1)'   : rho_host/rho0,
+    'patch_icpp(1)%alpha(1)'       : 1.,
+    # ==========================================================
+
+    # Lagrangian Bubbles ===========================
+     'lag_bubbles'                 : 'T',
+     'lag_nBubs_glb'               : 1194, # Number of bubbles
+     'lag_adap_dt'                 : 'F',
+     'lag_solver_approach'         : 2,
+     'lag_cluster_type'            : 2,
+     'lag_pressure_corrector'      : 'T',
+     'lag_smooth_type'             : 1,
+     'lag_bubble_model'            : 1,
+     'lag_heatTransfer_model'      : 'T',
+     'lag_massTransfer_model'      : 'T',
+     'lag_epsilonb'                : 1.0,
+     'lag_rkck_tolerance'          : 1.0e-05,
+     'lag_valmaxvoid'              : 0.9,
+     'lag_write_bubbles'           : 'F',
+     'lag_write_bubble_stats'      : 'F',
+     'csonhost'                    : c_host/c0,
+     'vischost'                    : mu_host/(rho0*x0*c0), 
+     'Thost'                       : T_host/T0,
+     'gammagas'                    : gamma_g,
+     'gammavapor'                  : gamma_v,
+     'pvap'                        : pv/p0,
+     'cpgas'                       : cp_g*(T0/(c0*c0)),
+     'cpvapor'                     : cp_v*(T0/(c0*c0)),
+     'kgas'                        : k_g*(T0/(x0*rho0*c0*c0*c0)),
+     'kvapor'                      : k_v*(T0/(x0*rho0*c0*c0*c0)),
+     'Rgas'                        : (R_uni/MW_g)*(T0/(c0*c0)),
+     'Rvap'                        : (R_uni/MW_v)*(T0/(c0*c0)),
+     'diffcoefvap'                 : diffVapor/(x0*c0),
+     'sigmabubble'                 : sigBubble/(rho0*x0*c0*c0),
+    # ==========================================================
+
+    # Fluids Physical Parameters ===============================
+    'fluid_pp(1)%gamma'            : 1.0/(gamma_host-1.0),
+    'fluid_pp(1)%pi_inf'           : gamma_host*(pi_inf_host/p0)/(gamma_host-1.0),
+    'fluid_pp(1)%Re(1)'            : 1.0/(mu_host/(rho0*c0*x0)),
+    # ==========================================================
+ }))
+
+# ==============================================================================

examples/3D_lag_bubbles_bubblescreen/input/README.txt

@@ -0,0 +1,5 @@
+
+The user input file 'input/lag_bubbles.dat' contains the initial conditions of the lagrangian bubbles.
+Each row represents the initial state of one specific bubble, which are:
+
+xPosition  yPosition  zPosition  xVel   yVel   zVel   radius  interfaceVelocity