Lagrangian subgrid bubble model

.gitignore

@@ -56,6 +56,8 @@ examples/*/workloads/
 examples/*/run-*/
 examples/*/logs/
 examples/**/*.f90
+!examples/3D_lag_bubbles_shbubcollapse/input/lag_bubbles.dat
+!examples/3D_lag_bubbles_bubblescreen/input/lag_bubbles.dat
 workloads/
 
 benchmarks/*batch/*/

benchmarks/viscous_weno5_sgb_acoustic/case.py

@@ -5,7 +5,7 @@
 # - viscosity enabled
 # - weno_Re_flux : T
 # - weno_order : 5
-# - bubbles : T
+# - bubbles_euler : T
 # - bubble_model : 3
 # - acoustic_source : T
 
@@ -216,7 +216,7 @@
     # ==========================================================
 
     # Bubbles ==================================================
-    'bubbles'                      : 'T',
+    'bubbles_euler'                : 'T',
     'bubble_model'                 : 3,
     'polytropic'                   : 'T',
     'polydisperse'                 : 'F',

docs/documentation/case.md

@@ -630,56 +630,73 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon
 
 - `%%bb_lowest_freq` specifies the lower frequency bound of the broadband acoustic wave. The upper frequency bound will be calculated as `%%bb_lowest_freq + %%bb_num_freq * %%bb_bandwidth`. The wave is no longer broadband below the lower bound and above the upper bound.
 
-### 9. Ensemble-Averaged Bubble Model
+### 9. Sub-grid Bubble Models
+
+| Parameter         | Type    | Description                                    |
+| ---:              | :----:  |          :---                                  |
+| `bubbles_euler`   | Logical	| Ensemble-averaged bubble modeling	|
+| `bubbles_lagrange`| Logical	| Volume-averaged bubble modeling	|
+| `bubble_model`    | Integer	| [1] Gilmore; [2] Keller--Miksis |
+| `mu_l0` *	        | Real 		|	Liquid viscosity (only specify in liquid phase)  |
+| `ss` *		        | Real 		|	Surface tension (only specify in liquid phase) |
+| `pv` *		        | Real 		|	Vapor pressure (only specify in liquid phase) |
+| `gamma_v` † 	    | Real 	  |	Specific heat ratio |
+| `M_v` †     	    | Real 		| Molecular weight |
+| `mu_v` †	        | Real 		|	Viscosity |
+| `k_v` †	          | Real 		|	Thermal conductivity |
+| `cp_v` †	        | Real 		|	Specific heat capacity |
+
+These options work only for gas-liquid two-component flows.
+Component indexes are required to be 1 for liquid and 2 for gas.
+
+- \* These parameters should be prepended with patch index $1$ that is filled with liquid: `fluid_pp(1)%`.
+- †  These parameters should be prepended with patch indexes filled with liquid and gas: `fluid_pp(1)%` and `fluid_pp(2)%`.
+
+This table lists the sub-grid bubble model parameters, which can be utilized in both the ensemble-averaged and volume-averaged bubble models.
+
+- `bubbles_euler` activates the ensemble-averaged bubble model.
+
+- `bubbles_lagrange` activates the volume-averaged bubble model.
+
+- `bubble_model` specified a model for spherical bubble dynamics by an integer of 1 and 2.
+`bubble_model = 1`, `2`, and `3` correspond to the Gilmore, Keller-Miksis, and Rayleigh-Plesset models.
+
+- `mu_l0`, `ss`, and `pv`, `gamma_v`, `M_v`, `mu_v`, `k_v`, and `cp_v` specify simulation parameters for the non-polytropic gas compression model.
+`mu_l0`, `ss`, and `pv` correspond to the liquid viscosity, surface tension, and vapor pressure, respectively.
+`gamma_v`, `M_v`, `mu_v`, `k_v`, and `cp_v` specify the specific heat ratio, molecular weight, viscosity, thermal conductivity and specific heat capacity of a chosen component (`cp_v` only for ``bubbles_lagrange = 'T'``).
+Implementation of the parameters into the model follow [Ando (2010)](references.md#Ando10).
+
+#### 9.1 Ensemble-Averaged Bubble Model
 
 | Parameter      | Type    | Description                                    |
 | ---:           | :----:  |          :---                                  |
-| `bubbles` 		 | Logical	| Ensemble-averaged bubble modeling	 |
-| `bubble_model` | Integer	| [1] Gilmore; [2] Keller--Miksis |
-| `polytropic`   | Logical	| Polytropic gas compression |
-| `thermal` 		 | Integer	| Thermal model: [1] Adiabatic; [2] Isothermal; [3] Transfer |
+| `bubbles_euler`| Logical	| Ensemble-averaged bubble modeling	|
+| `polytropic`      | Logical	| Polytropic gas compression |
+| `thermal` 		    | Integer	| Thermal model: [1] Adiabatic; [2] Isothermal; [3] Transfer |
 | `R0ref` 			 | Real		  | Reference bubble radius |
-| `polydisperse`   | Logical	| Polydispersity in equilibrium bubble radius R0|
+| `polydisperse`   | Logical	| Polydispersity in equilibrium bubble radius R0 |
 | `nb` 			     | Integer	| Number of bins: [1] Monodisperse; [$>1$] Polydisperse |
 | `poly_sigma` 	       | Real 		|	Standard deviation for probability density function of polydisperse bubble populations |
 | `R0_type` 	       | Integer 		|	Quadrature rule for probability density function of polydisperse bubble populations |
 | `Ca` 			     | Real		  | Cavitation number |
 | `Web` 			   | Real		  | Weber number |
 | `Re_inv` 		   | Real		  | Inverse Reynolds number |
-| `mu_l0` *	     | Real 		|	Liquid viscosity (only specify in liquid phase)  |
-| `ss` *		     | Real 		|	Surface tension (only specify in liquid phase) |
-| `pv` *		     | Real 		|	Vapor pressure (only specify in liquid phase) |
-| `gamma_v` † 	 | Real 	  |	Specific heat ratio |
-| `M_v` †     	 | Real 		| Molecular weight |
-| `mu_v` †	     | Real 		|	Viscosity |
-| `k_v` †	       | Real 		|	Thermal conductivity |
-| `qbmm` 	       | Logical 		|	Quadrature by  method of moments|
-| `dist_type` 	       | Integer 		|	Joint probability density function for bubble radius and velocity (only for ``qbmm = 'T'``)|
-| `sigR` 	       | Real 		|	Standard deviation for the probability density function of bubble radius (only for ``qbmm = 'T'``) |
-| `sigV` 	       | Real 		|	Standard deviation for the probability density function of bubble velocity (only for ``qbmm = 'T'``) |
+| `qbmm` 	       | Logical 		|	Quadrature by  method of moments |
+| `dist_type` 	       | Integer 		|	Joint probability density function for bubble radius and velocity (only for ``qbmm = 'T'``) |
+| `sigR` 	       | Real 		|	Standard deviation for the probability density function of bubble radius (only for ``qbmm = 'T'``) (EE)|
+| `sigV` 	       | Real 		|	Standard deviation for the probability density function of bubble velocity (only for ``qbmm = 'T'``) (EE)|
 | `rhoRV`	       | Real 		|	Correlation coefficient for the joint probability density function of bubble radius and velocity (only for ``qbmm = 'T'``) |
 
-These options work only for gas-liquid two-component flows.
-Component indexes are required to be 1 for liquid and 2 for gas.
-
-- \* These parameters should be prepended with patch index $1$ that is filled with liquid: `fluid_pp(1)%`.
-- †  These parameters should be prepended with patch indexes filled with liquid and gas: `fluid_pp(1)%` and `fluid_pp(2)%`.
-
 This table lists the ensemble-averaged bubble model parameters.
 
-- `bubbles` activates the ensemble-averaged bubble model.
-
-- `bubble_model` specified a model for spherical bubble dynamics by an integer of 1 and 2.
-`bubble_model = 1`, `2`, and `3` correspond to the Gilmore, Keller-Miksis, and Rayleigh-Plesset models.
-
 - `polytropic` activates polytropic gas compression in the bubble.
 When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on ([Preston et al., 2007](references.md#Preston07)).
 
-- `polydisperse` activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius.
-
 - `thermal` specifies a model for heat transfer across the bubble interface by an integer from 1 through 3.
 `thermal = 1`, `2`, and `3` correspond to no heat transfer (adiabatic gas compression), isothermal heat transfer, and heat transfer with a constant heat transfer coefficient based on [Preston et al., 2007](references.md#Preston07), respectively.
 
+- `polydisperse` activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius.
+
 - `R0ref` specifies the reference bubble radius.
 
 - `nb` specifies the number of discrete bins that define the probability density function (PDF) of the equilibrium bubble radius.
@@ -691,11 +708,6 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 
 - `Ca`, `Web`, and `Re_inv` respectively specify the Cavitation number, Weber number, and the inverse Reynolds number that characterize the offset of the gas pressure from the vapor pressure, surface tension, and liquid viscosity when the polytropic gas compression model is used.
 
-- `mu_l0`, `ss`, and `pv`, `gamma_v`, `M_v`, `mu_v`, and `k_v` specify simulation parameters for the non-polytropic gas compression model.
-`mu_l0`, `ss`, and `pv` correspond to the liquid viscosity, surface tension, and vapor pressure, respectively.
-`gamma_v`, `M_v`, `mu_v`, and `k_v` specify the specific heat ratio, molecular weight, viscosity, and thermal conductivity of a chosen component.
-Implementation of the parameters into the model follow [Ando (2010)](references.md#Ando10).
-
 - `qbmm` activates quadrature by method of moments, which assumes a PDF for bubble radius and velocity.
 
 - `dist_type` specifies the initial joint PDF of initial bubble radius and bubble velocity required in qbmm. `dist_type = 1`  and `2` correspond to binormal and lognormal-normal distributions respectively.
@@ -706,6 +718,46 @@ Implementation of the parameters into the model follow [Ando (2010)](references.
 
 - `rhoRV` specifies the correlation coefficient of the joint PDF of bubble radius and bubble velocity required in the QBMM feature.
 
+#### 9.2 Volume-Averaged Bubble Model
+
+| Parameter             | Type    | Description                                               |
+| ---:                  | :---:   | :---                                                      |
+| `bubbles_lagrange`    | Logical | Lagrangian subgrid bubble model switch                    |
+| `nBubs_glb`           | Integer | Global number of bubbles                                  |
+| `solver_approach`     | Integer | 1: One-way coupling, 2: two-way coupling                  |
+| `cluster_type`        | Integer | Method to find p_inf                                      |
+| `pressure_corrector`  | Logical | Cell pressure correction term                             |
+| `smooth_type`         | Integer | Smoothing function. 1: Gaussian, 2:Delta 3x3              |
+| `heatTransfer_model`  | Logical | Activates the interface heat transfer model               |
+| `massTransfer_model`  | Logical | Activates the interface mass transfer model               |
+| `write_bubbles`       | Logical | Write files to track the bubble evolution each time step  |
+| `write_bubbles_stats` | Logical | Write the maximum and minimum radius of each bubble       |
+| `epsilonb`            | Real    | Standard deviation scaling for the gaussian function      |
+| `charwidth`           | Real    | Domain virtual depth (z direction, for 2D simulations)    |
+| `valmaxvoid`          | Real    | Maximum void fraction permitted                           |
+| `c0`                  | Real    | Reference speed                                           |
+| `rho0`                | Real    | Reference density                                         |
+| `T0`                  | Real    | Reference temperature                                     |
+| `x0`                  | Real    | Reference length                                          |
+| `Thost`               | Real    | Temperature of the surrounding liquid (host)              |
+| `diffcoefvap`         | Real    | Vapor diffusivity in the gas                              |
+| `rkck_adap_dt`        | Logical | Activates the adaptive rkck time stepping algorithm       |
+| `rkck_tolerance`      | Real    | Admissible error truncation tolerance in the rkck stepper  |
+
+- `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.
+
+- `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.
+
+- `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 `pressure_corrector` is true.
+
+- `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.
+
+- `heatTransfer_model` Activates the heat transfer model at the bubble's interface based on ([Preston et al., 2007](references.md#Preston07)).
+
+- `massTransfer_model` Activates the mass transfer model at the bubble's interface based on ([Preston et al., 2007](references.md#Preston07)).
+
+- `rkck_adap_dt` Activates the adaptive 4th/5th order Runge—Kutta–Cash–Karp (RKCK) time-stepping algorithm (requires `time_stepper ==4`). 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 (`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.
+
 ### 10. Velocity Field Setup
 
 | Parameter              | Type    | Description |

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/0D_bubblecollapse_adap/case.py

@@ -119,7 +119,7 @@
     # ==========================================================
 
     # Bubbles ==================================================
-    'bubbles'                       : 'T',
+    'bubbles_euler'                 : 'T',
     'bubble_model'                  : 2,
 
     # Nondimensional numbers

examples/1D_bubblescreen/case.py

@@ -176,7 +176,7 @@
     # ==========================================================
 
     # Bubbles ==================================================
-    'bubbles'                      : 'T',
+    'bubbles_euler'                : 'T',
     'bubble_model'                 : 3,
     'polytropic'                   : 'T',
     'polydisperse'                 : 'F',

examples/1D_exp_bubscreen/case.py

@@ -210,7 +210,7 @@
     # ==========================================================
 
     # Bubbles ==================================================
-    'bubbles'               : 'T',
+    'bubbles_euler'               : 'T',
     # in user guide... 1 = gilbert 2 = keller-miksis
     # but gilbert won't work for the equations that you are using... (i think)
     'bubble_model'          : 2,

examples/1D_poly_bubscreen/case.py

@@ -184,7 +184,7 @@
     'rhoref'                : rho0,
     # ==========================================================================
     # Bubbles ==================================================================
-    'bubbles'               : 'T',
+    'bubbles_euler'         : 'T',
     'bubble_model'          : 2,
     'polytropic'            : 'T',
     # 'polydisperse'          : 'T',

examples/1D_qbmm/case.py

@@ -183,7 +183,7 @@
     # ==========================================================
 
     # Bubbles ==================================================
-    'bubbles'               : 'T',
+    'bubbles_euler'         : 'T',
     'bubble_model'          : 2,
     'polytropic'            : 'F',
     'polydisperse'          : 'T',

examples/2D_bubbly_steady_shock/case.py

@@ -204,7 +204,7 @@
     # ==========================================================
 
     # Bubbles ==================================================
-    'bubbles'               : 'T',
+    'bubbles_euler'         : 'T',
     'bubble_model'          : 2,
     'polytropic'            : 'T',
     'polydisperse'          : 'F',

examples/2D_whale_bubble_annulus/case.py

@@ -173,7 +173,7 @@
     # ==========================================================================
 
     # Bubbles ==================================================================
-    'bubbles'               : 'T',
+    'bubbles_euler'         : 'T',
     'bubble_model'          : 3,
     'polytropic'            : 'T',
     'R0_type'               : 1,