Skip to content

Commit

Permalink
Operator splitting, adaptive time stepping, and other fixes for mixin…
Browse files Browse the repository at this point in the history
…g layer and bubbles (#285)

Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Anand Radhakrishnan <[email protected]>
Co-authored-by: Anand Radhakrishnan <[email protected]>
Co-authored-by: Anand Radhakrishnan <[email protected]>
Co-authored-by: Anand <[email protected]>
  • Loading branch information
12 people authored Apr 17, 2024
1 parent b3761e3 commit 437e2d1
Show file tree
Hide file tree
Showing 52 changed files with 3,383 additions and 227 deletions.
20 changes: 17 additions & 3 deletions docs/documentation/case.md
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@ There are multiple sets of parameters that must be specified in the python input
8. [(Optional) Ensemble-Averaged Bubble Model Parameters](#8-ensemble-averaged-bubble-model)
9. [(Optional) Velocity Field Setup Parameters](#9-velocity-field-setup)
10. [(Optional) Phase Change Parameters](#10-Phase-Change-Model)
11. [(Optional) Artificial Mach Number Parameters](#11-artificial-Mach-number)

Items 7, 8, 9, and 10 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, initial velocity field setup, and phase change, respectively.
Definition of the parameters is described in the following subsections.
Expand Down Expand Up @@ -348,9 +349,11 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
| `model_eqns` | Integer | Multicomponent model: [1] $\Gamma/\Pi_\infty$; [2] 5-equation; [3] 6-equation\\%;%[4] 4-equation |
| `alt_soundspeed` * | Logical | Alternate sound speed and $K \nabla \cdot u$ for 5-equation model |
| `adv_alphan` | Logical | Equations for all $N$ volume fractions (instead of $N-1$) |
| `adv_n` | Logical | Solving directly for the number density (in the method of classes) and compute void fraction from the number density |
| `mpp_lim` | Logical | Mixture physical parameters limits |
| `mixture_err` | Logical | Mixture properties correction |
| `time_stepper` | Integer | Runge-Kutta order [1-3] |
| `time_stepper` | Integer | Runge--Kutta order [1-3] |
| `adap_dt` | Loginal | Strang splitting scheme with adaptive time stepping |
| `weno_order` | Integer | WENO order [1,3,5] |
| `weno_eps` | Real | WENO perturbation (avoid division by zero) |
| `mapped_weno` | Logical | WENO with mapping of nonlinear weights |
Expand Down Expand Up @@ -396,13 +399,17 @@ $$ \alpha_N=1-\sum^{N-1}_{i=1} \alpha_i $$
where $\alpha_i$ is the void fraction of $i$-th component.
When a single-component flow is simulated, it requires that `adv_alphan = 'T'`.

- `adv_n` activates the direct computation of number density by the Riemann solver instead of computing number density from the void fraction in the method of classes.

- `mpp_lim` activates correction of solutions to avoid a negative void fraction of each component in each grid cell, such that $\alpha_i>\varepsilon$ is satisfied at each time step.

- `mixture_err` activates correction of solutions to avoid imaginary speed of sound at each grid cell.

- `time_stepper` specifies the order of the Runge-Kutta (RK) time integration scheme that is used for temporal integration in simulation, from the 1st to 5th order by corresponding integer.
Note that `time_stepper = 3` specifies the total variation diminishing (TVD), third order RK scheme ([Gottlieb and Shu, 1998](references.md#Gottlieb98)).

- `adap_dt` activates the Strang operator splitting scheme which splits flux and source terms in time marching, and an adaptive time stepping strategy is implemented for the source term. It requires `bubbles = 'T'`, `polytropic = 'T'`, `adv_n = 'T'` and `time_stepper = 3`.

- `weno_order` specifies the order of WENO scheme that is used for spatial reconstruction of variables by an integer of 1, 3, and 5, that correspond to the 1st, 3rd, and 5th order, respectively.

- `weno_eps` specifies the lower bound of the WENO nonlinear weights.
Expand Down Expand Up @@ -636,8 +643,7 @@ The parameters are optionally used to define initial velocity profiles and pertu

- `vel_profile` activates setting the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases.

- `instability_wave` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile.
This option only works for 2D and 3D cases, together with `vel_profile = 'T'`.
- `instability_wave` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for `n > 0`, `bc_y%[beg,end] = -5`, and `vel_profile = 'T'`.

### 11. Phase Change Model
| Parameter | Type | Description |
Expand All @@ -655,6 +661,14 @@ This option only works for 2D and 3D cases, together with `vel_profile = 'T'`.

- `ptgalpha_eps` Specifies the tolerance used for the Newton Solvers used in the pTg-equilibrium model.

### 11. Artificial Mach Number
| Parameter | Type | Description |
| ---: | :----: | :--- |
| `pi_fac` | Real | Ratio of artificial and true `pi_\infty` values|

- `pi_fac` specifies the ratio of artificial and true `pi_\infty` values (`=` artificial `pi_\infty` / true `pi_\infty`). This parameter enables the use of true `pi_\infty` in bubble dynamics models, when the `pi_\infty` given in the `case.py` file is an artificial value.


## Enumerations

### Boundary conditions
Expand Down
164 changes: 164 additions & 0 deletions examples/0D_bubblecollapse_adap/case.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,164 @@
#!/usr/bin/env python2

import math
import json

# FLUID PROPERTIES ============================================================
# Water
n_tait = 7.1
B_tait = 306.E+06
rho0 = 1.E+03
mul0 = 1.002E-03
ss = 0.07275
pv = 2.3388E+03

# Vapor
gamma_v = 1.33
M_v = 18.02
mu_v = 0.8816E-05
k_v = 0.019426

# Air
gamma_n = 1.4
M_n = 28.97
mu_n = 1.8E-05
k_n = 0.02556

# REFERENCE VALUES ============================================================
R0ref = 50.E-06
x0 = R0ref
p0 = 8236. # for Ca = 1 in mixing layer scale
u0 = math.sqrt( p0/rho0 )
patm = 1.
cact = math.sqrt(n_tait*(p0+B_tait)/rho0)

# NONDIMENSIONAL NUMBERS ======================================================
Ca = (p0 - pv)/(rho0*(u0**2.)) # Cavitation number
We = rho0*(u0**2.)*R0ref/ss # Weber number
Re_inv = mul0/(rho0*u0*R0ref) # Inv. bubble Reynolds number

# BUBBLES =====================================================================
vf0 = 1e-5
nb = 1

# DOMAIN ======================================================================
Nx = 30
Ldomain = 20.E-03
L = Ldomain/x0
dx = L/float(Nx+1)

# TIME STEPS ==================================================================
Tfinal = 0.05
Nt = int(5e2+1)
t_save = 1
dt = Tfinal/(Nt-1)

# Configuring case dictionary
print(json.dumps({
# Logistics ================================================
'run_time_info' : 'T',
# ==========================================================

# Computational Domain Parameters ==========================
'x_domain%beg' : -0.5*L,
'x_domain%end' : 0.5*L,
'stretch_x' : 'F',
'cyl_coord' : 'F',
'm' : Nx,
'n' : 0,
'p' : 0,
'dt' : dt,
't_step_start' : 0,
't_step_stop' : Nt,
't_step_save' : t_save,
# ==========================================================

# Simulation Algorithm Parameters ==========================
'num_patches' : 1,
'model_eqns' : 2,
'alt_soundspeed' : 'F',
'num_fluids' : 1,
'adv_alphan' : 'T',
'mpp_lim' : 'F',
'mixture_err' : 'F',
'time_stepper' : 3,
'weno_order' : 5,
'weno_eps' : 1.E-16,
'mapped_weno' : 'T',
'null_weights' : 'F',
'mp_weno' : 'T',
'riemann_solver' : 2,
'wave_speeds' : 1,
'avg_state' : 2,
'bc_x%beg' : -1,
'bc_x%end' : -1,
# ==========================================================

# Formatted Database Files Structure Parameters ============
'format' : 1,
'precision' : 2,
'prim_vars_wrt' :'T',
'parallel_io' :'T',
'fd_order' : 1,
# ==========================================================

# Patch 1 _ Background =====================================
'patch_icpp(1)%geometry' : 1,
'patch_icpp(1)%x_centroid' : 0.,
'patch_icpp(1)%length_x' : L,
'patch_icpp(1)%vel(1)' : 0.,
'patch_icpp(1)%pres' : 1000.,
'patch_icpp(1)%alpha_rho(1)' : (1.-vf0),
'patch_icpp(1)%alpha(1)' : vf0,
'patch_icpp(1)%r0' : 1.,
'patch_icpp(1)%v0' : 0.,
# ==========================================================

# Non-polytropic gas compression model AND/OR Tait EOS =====
'pref' : p0,
'rhoref' : rho0,
# ==========================================================

# Bubbles ==================================================
'bubbles' : 'T',
'bubble_model' : 2,

# Nondimensional numbers
'Ca' : Ca,
'Web' : We,
'Re_inv' : Re_inv,

# adv_n
'adv_n' : 'T',

# adap_dt
'adap_dt' : 'T',

# Gas compression model
'polytropic' : 'T',
'thermal' : 1,

# Polydispersity
'polydisperse' : 'F',
'nb' : nb,

# QBMM
'qbmm' : 'F',
# ==========================================================

# Fluids Physical Parameters ===============================
# Surrounding liquid
'fluid_pp(1)%gamma' : 1.E+00/(n_tait-1.E+00),
'fluid_pp(1)%pi_inf' : n_tait*(B_tait/p0)/(n_tait-1.),
'fluid_pp(1)%ss' : ss,
'fluid_pp(1)%pv' : pv,

# Last fluid_pp is always reserved for bubble gas state ===
# if applicable ==========================================
'fluid_pp(2)%gamma' : 1./(gamma_n-1.),
'fluid_pp(2)%pi_inf' : 0.0E+00,
# ==========================================================

}))

# ==============================================================================
140 changes: 140 additions & 0 deletions examples/2D_mixing_artificial_Ma/case.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
#!/usr/bin/env python3

import math
import json

# FLUID PROPERTIES (WATER, VAPOR & AIR) ========================================
# Water
rho_w = 1.E+03 # [kg/m3] density of water
gamma_w = 7.1 # [1] specific heat ratio
pi_inf_w = 3.06E+08 # [N/m2] water stiffness
# ==============================================================================

# REFERENCE VALUES =============================================================
x_ref = 0.002475 # [m] initial vorticity thickness
u_ref = 3.4343 # [m/s] upper stream velocity
# ==============================================================================

# NON-DIMENSIONAL NUMBERS ======================================================
Re = 50. # [1] Reynolds number based on the upper stream
# velocity and the initial vorticity thickness
# ==============================================================================

# SPEED OF SOUND ===============================================================
pres = 8236. # [N/m2] pressure of water
c_w = math.sqrt(gamma_w*(pres+pi_inf_w)/rho_w)
# ==============================================================================

# MODIFIED PROPERTIES FOR ARTIFICIAL MACH NUMBER ==============================
Ma_t = 0.1 # [1] target Mach number
pi_fac = (rho_w*u_ref**2/(gamma_w*Ma_t**2)-pres)/pi_inf_w
pi_inf_w = pi_inf_w*pi_fac
c_w = math.sqrt(gamma_w*(pres+pi_inf_w)/rho_w)
# =============================================================================

# SIMULATION SETUP ============================================================
# Domain size
Lx = 59.0
Ly = 59.0

# Number of grid cells
Nx = 191
Ny = 191

# Grid spacing
dx = Lx/float(Nx+1)
dy = Ly/float(Ny+1)

# Time advancement
cfl = 5e-1
T = 20.
dt = cfl*dx/(1.+c_w/u_ref)
Ntfinal = int(T/dt)
Ntstart = 0
Nfiles = 20
t_save = int(math.ceil((Ntfinal-0)/float(Nfiles)))
Nt = t_save*Nfiles
t_step_start = Ntstart
t_step_stop = int(Nt)

# ==============================================================================


# Configuring case dictionary
print(json.dumps({
# Logistics ================================================================
'run_time_info' : 'T',
# ==========================================================================

# Computational Domain Parameters ==========================================
'x_domain%beg' : 0.,
'x_domain%end' : Lx,
'y_domain%beg' : -Ly/2.,
'y_domain%end' : Ly/2.,
'm' : Nx,
'n' : Ny,
'p' : 0,
'dt' : dt,
't_step_start' : t_step_start,
't_step_stop' : t_step_stop,
't_step_save' : t_save,
# ==========================================================================

# Simulation Algorithm Parameters ==========================================
'num_patches' : 1,
'model_eqns' : 2,
'num_fluids' : 1,
'time_stepper' : 3,
'weno_order' : 5,
'weno_eps' : 1.E-16,
'weno_Re_flux' : 'T',
'mapped_weno' : 'T',
'riemann_solver' : 2,
'wave_speeds' : 1,
'avg_state' : 2,
'bc_x%beg' : -1,
'bc_x%end' : -1,
'bc_y%beg' : -5,
'bc_y%end' : -5,
# ==========================================================================

# Formatted Database Files Structure Parameters ============================
'format' : 1,
'precision' : 2,
'cons_vars_wrt' :'T',
'prim_vars_wrt' :'T',
'parallel_io' :'T',
'fd_order' : 1,
'omega_wrt(3)' :'T',
# ==========================================================================

# Patch 1 ==================================================================
'patch_icpp(1)%geometry' : 3,
'patch_icpp(1)%x_centroid' : Lx/2.,
'patch_icpp(1)%y_centroid' : 0.,
'patch_icpp(1)%length_x' : Lx,
'patch_icpp(1)%length_y' : Ly,
'patch_icpp(1)%alpha_rho(1)' : 1.,
'patch_icpp(1)%alpha(1)' : 1.,
'patch_icpp(1)%vel(1)' : 1.,
'patch_icpp(1)%vel(2)' : 0.,
'patch_icpp(1)%pres' : pres/(rho_w*u_ref**2),
# ==========================================================================

# Mixing layer === =========================================================
'vel_profile' : 'T',
'instability_wave' : 'T',
# ==========================================================================

# Artificial Mach number ===================================================
'pi_fac' : pi_fac,
# ==========================================================================

# Fluids Physical Parameters ===============================================
'fluid_pp(1)%gamma' : 1.E+00/(gamma_w-1.E+00),
'fluid_pp(1)%pi_inf' : gamma_w*(pi_inf_w/(rho_w*u_ref**2))/(gamma_w-1.E+00),
'fluid_pp(1)%Re(1)' : Re,
# =========================================================================
}))

# ==============================================================================
4 changes: 4 additions & 0 deletions src/common/m_helper.fpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,10 @@ module m_helper

use m_global_parameters !< Definitions of the global parameters

use m_mpi_common !< MPI modules

use ieee_arithmetic !< For checking NaN

! ==========================================================================

implicit none
Expand Down
Loading

0 comments on commit 437e2d1

Please sign in to comment.