GRCBC

docs/documentation/case.md

@@ -815,6 +815,22 @@ The boundary condition supported by the MFC are listed in table [Boundary Condit
 Their number (`#`) corresponds to the input value in `input.py` labeled `bc_[x,y,z]%[beg,end]` (see table [Simulation Algorithm Parameters](#5-simulation-algorithm)).
 The entries labeled "Characteristic." are characteristic boundary conditions based on [Thompson (1987)](references.md#Thompson87) and [Thompson (1990)](references.md#Thompson90).
 
+### Generalized Characteristic Boundary conditions
+
+| Parameter                     | Type    | Description |
+| ---:                          | :----:  | :--- |
+| `bc_[x,y,z]%grcbc_in`         | Logical | Enable grcbc for subsonic inflow |
+| `bc_[x,y,z]%grcbc_out`        | Logical | Enable grcbc for subsonic outflow (pressure)|
+| `bc_[x,y,z]%grcbc_vel_out`    | Logical | Enable grcbc for subsonic outflow (pressure + normal velocity) |
+| `bc_[x,y,z]%vel_in`           | Real Array | Inflow velocities in x, y and z directions |
+| `bc_[x,y,z]%vel_out`          | Real Array | Outflow velocities in x, y and z directions |
+| `bc_[x,y,z]%pres_in`          | Real    | Inflow pressure |
+| `bc_[x,y,z]%pres_out`         | Real    | Outflow pressure |
+| `bc_[x,y,z]%alpha_rho_in`     | Real Array | Inflow density |
+| `bc_[x,y,z]%alpha_in`         | Real Array | Inflow void fraction |
+
+This boundary condition can be used for subsonic inflow (`bc_[x,y,z]%[beg,end]` = -7) and subsonic outflow (`bc_[x,y,z]%[beg,end]` = -8) characteristic boundary conditions. These are based on [Pirozzoli (2013)](references.md#Pirozzoli13). This enables to provide inflow and outflow conditions outside the computational domain.
+
 ### Patch types
 
 | #    | Name               | Dim.  | Smooth | Description |

docs/documentation/references.md

@@ -32,6 +32,8 @@
 
 - <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>
+
 - <a id="Preston07">Preston, A., Colonius, T., and Brennen, C. (2007). A reduced-order model of diffusive effects on the dynamics of bubbles. Physics of Fluids, 19(12):123302.</a>
 
 - <a id="Saurel09">Saurel, R., Petitpas, F., and Berry, R. A. (2009). Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. journal of Computational Physics, 228(5):1678–1712</a>

examples/2D_acoustic_pulse/case.py

@@ -0,0 +1,115 @@
+import math
+import json
+
+# Numerical setup
+Nx = 99
+Ny = 99
+dx = 8./(1.*(Nx+1))
+
+alf_st = 0.4
+
+p_inf = 101325
+rho_inf = 1
+gam = 1.4
+
+c = math.sqrt(gam*(p_inf) / rho_inf)
+cfl = 0.3
+mydt = cfl * dx / c 
+Tfinal = 80*1 / c
+Nt = int(Tfinal/mydt)
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : -4,
+    'x_domain%end'                 : 4,
+    'y_domain%beg'                 : -4,
+    'y_domain%end'                 : 4,
+    'm'                            : Nx,
+    'n'                            : Ny,
+    'p'                            : 0,
+    'dt'                           : mydt,
+    't_step_start'                 : 0,
+    't_step_stop'                  : Nt,
+    't_step_save'                  : int(Nt/100),
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'num_patches'                  : 2,
+    'model_eqns'                   : 2,
+    'alt_soundspeed'               : 'F',
+    'num_fluids'                   : 1,
+    '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'                      : 'F',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     : -8,
+    'bc_x%end'                     : -8,
+    'bc_y%beg'                     : -8,
+    'bc_y%end'                     : -8,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'parallel_io'                  :'T',
+    'omega_wrt(3)'                 :'T',
+    'fd_order'                     : 2,
+    # ==========================================================================
+
+    # Patch 1 ==================================================================
+    'patch_icpp(1)%geometry'       : 3,
+    'patch_icpp(1)%x_centroid'     : 0,
+    'patch_icpp(1)%y_centroid'     : 0,
+    'patch_icpp(1)%length_x'       : 8.,
+    'patch_icpp(1)%length_y'       : 8.,
+    'patch_icpp(1)%vel(1)'         : 0,
+    'patch_icpp(1)%vel(2)'         : 0,
+    'patch_icpp(1)%pres'           : p_inf,
+    'patch_icpp(1)%alpha_rho(1)'   : rho_inf,
+    'patch_icpp(1)%alpha(1)'       : 1.,
+    # ==========================================================================
+
+    # Patch 2 ==================================================================
+    'patch_icpp(2)%geometry'       : 2,
+    'patch_icpp(2)%x_centroid'     : 0,
+    'patch_icpp(2)%y_centroid'     : 0,
+    'patch_icpp(2)%radius'         : 1.,
+    'patch_icpp(2)%vel(1)'         : 0,
+    'patch_icpp(2)%vel(2)'         : 0,
+    'patch_icpp(2)%pres'           : f"{p_inf}*(1 - 0.5*({gam} - 1)*({alf_st})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**({gam} / ({gam} - 1))", 
+    'patch_icpp(2)%alpha_rho(1)'   : f"{rho_inf}*(1 - 0.5*({gam} - 1)*({alf_st})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**(1 / ({gam} - 1))",
+    'patch_icpp(2)%alpha(1)'       : 1.,
+    'patch_icpp(2)%alter_patch(1)' : 'T',
+    # ==========================================================================
+
+    # CBC Inflow / Outflow ========================================
+    'bc_x%grcbc_in'                 : 'F', 
+    'bc_x%grcbc_out'                : 'T', 
+    'bc_x%grcbc_vel_out'            : 'F', 
+    'bc_x%pres_out'                 : p_inf,
+    'bc_y%grcbc_in'                 : 'F', 
+    'bc_y%grcbc_out'                : 'T', 
+    'bc_y%grcbc_vel_out'            : 'F', 
+    'bc_y%pres_out'                 : p_inf,
+    # # ========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.E+00/(gam-1.E+00),
+    'fluid_pp(1)%pi_inf'           : 0.0,
+    # ==========================================================================
+}))
+
+# ==============================================================================

examples/2D_ibm_steady_shock/case.py

@@ -0,0 +1,176 @@
+import json
+import math
+
+Mu = 1.84E-05
+gam_a = 1.4
+gam_b = 1.1
+
+
+x0      = 10E-06
+p0      = 101325
+rho0    = 1.0
+c0      = math.sqrt( p0/rho0 )
+patm    = 1.
+rhoatm    = 1.
+
+#air props
+n_tait  = 1.4
+B_tait  = 0.0
+
+
+vf0 = 0.0 
+
+cact    =  math.sqrt(n_tait*(p0 + p0*B_tait) / ((1 - vf0) * rho0) )
+cfl     = 0.3
+Nx = 400
+Ny = 200
+dx = 6.E-03 / (x0*float(Nx))
+dt      = cfl*dx*c0/cact
+
+
+vel = 1.5
+vel_ac = vel * c0 
+
+Min = (n_tait + 1)*vel_ac + math.sqrt((n_tait + 1)**2 * vel_ac**2 + 16 * cact**2)
+Min = Min / (4 * cact)
+beta = (n_tait + 1) * Min**2 / (Min**2 * (n_tait - 1 + 2*vf0) + 2 * (1 - vf0))
+delta = (1 - vf0) + n_tait * Min**2 * (beta - 1) * (1 + B_tait) / beta
+
+vel = Min * cact * (beta - 1) / (beta * c0)
+
+# Configuring case dictionary
+print(json.dumps({
+    # Logistics ================================================================
+    'run_time_info'                : 'F',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : 0.0E+00,
+    'x_domain%end'                 : 6.0E-03 / x0,
+    'y_domain%beg'                 : 0.0E+00,
+    'y_domain%end'                 : 3.0E-03 / x0,
+    'cyl_coord'                    : 'F',
+    'm'                            : Nx,
+    'n'                            : Ny,
+    'p'                            : 0,
+    'dt'                           : dt,
+    't_step_start'                 : 0,
+    't_step_stop'                  : 1000,  #3000
+    't_step_save'                  : 10,  #10
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'num_patches'                  : 2,
+    # Use the 5 equation model
+    'model_eqns'                   : 2,                  
+    'alt_soundspeed'               : 'F',
+    'num_fluids'                   : 1,
+    # No need to ensure the volume fractions sum to unity at the end of each
+    # time step
+    'mpp_lim'                      : 'F',
+    # Correct errors when computing speed of sound
+    'mixture_err'                  : 'T',
+    # Use TVD RK3 for time marching
+    'time_stepper'                 : 3,
+    # Reconstruct the primitive variables to minimize spurious
+    # Use WENO5
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1.E-16,
+    'weno_Re_flux'                 : 'F',
+    'weno_avg'                     : 'T',
+    'avg_state'                    : 2,
+    # Use the mapped WENO weights to maintain monotinicity
+    'mapped_weno'                  : 'T',
+    'null_weights'                 : 'F',
+    'mp_weno'                      : 'F',
+    # Use the HLLC  Riemann solver
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 1,
+    'bc_x%beg'                     : -7,
+    'bc_x%end'                     : -8,
+    'bc_y%beg'                     : -3,
+    'bc_y%end'                     : -3,
+    # Set IB to True and add 1 patch
+    'ib'                           : 'T',
+    'num_ibs'                      : 1,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    # Export primitive variables in double precision with parallel
+    # I/O to minimize I/O computational time during large simulations
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                :'T',
+    'fd_order'                     : 1,
+    'omega_wrt(3)'                 :'T',
+    'parallel_io'                  :'T',
+    # ==========================================================================
+
+    #Ambient State =====================================
+    'patch_icpp(1)%geometry'       : 3,
+    'patch_icpp(1)%x_centroid'     : 3.0E-03 / x0,
+    'patch_icpp(1)%y_centroid'     : 1.50E-03 / x0,
+    'patch_icpp(1)%length_x'       : 6.0E-03 / x0,
+    'patch_icpp(1)%length_y'       : 3.0E-03 / x0,  
+    'patch_icpp(1)%alpha_rho(1)'   : rhoatm,
+    'patch_icpp(1)%alpha(1)'       : 1,
+    'patch_icpp(1)%vel(1)'         : 0.0,
+    'patch_icpp(1)%vel(2)'         : 0.0E+00,
+    'patch_icpp(1)%pres'           : patm,
+    'patch_icpp(1)%r0'             : 1.,
+    'patch_icpp(1)%v0'             : 0.0E+00,
+    # # ========================================================================
+
+    #Shocked State =====================================   
+    'patch_icpp(2)%geometry'       : 3,
+    'patch_icpp(2)%x_centroid'     : 0.5E-03 / x0,
+    'patch_icpp(2)%y_centroid'     : 1.50E-03 / x0,
+    'patch_icpp(2)%length_x'       : 1.0E-03 / x0,
+    'patch_icpp(2)%length_y'       : 3.0E-03 / x0,  
+    'patch_icpp(2)%alpha_rho(1)'   : beta*rhoatm,
+    'patch_icpp(2)%alpha(1)'       : 1.,
+    'patch_icpp(2)%vel(1)'         : vel,
+    'patch_icpp(2)%vel(2)'         : 0.0E+00,
+    'patch_icpp(2)%pres'           : delta*patm,
+    'patch_icpp(2)%r0'             : 1.,
+    'patch_icpp(2)%v0'             : 0.0E+00,
+    'patch_icpp(2)%alter_patch(1)' : 'T',  
+    # # ========================================================================
+
+
+    # CBC Inflow / Outflow ========================================
+    'bc_x%grcbc_in'                 : 'T', 
+    'bc_x%grcbc_out'                : 'F', 
+    'bc_x%grcbc_vel_out'            : 'F', 
+    'bc_x%vel_in(1)'                : vel, 
+    'bc_x%vel_in(2)'                : 0, 
+    'bc_x%vel_in(3)'                : 0, 
+    'bc_x%pres_in'                  : delta*patm,
+    'bc_x%alpha_rho_in(1)'          : beta*rhoatm,
+    'bc_x%alpha_in(1)'              : 1,
+    'bc_x%vel_out(1)'               : vel, 
+    'bc_x%vel_out(2)'               : 0, 
+    'bc_x%vel_out(3)'               : 0, 
+    'bc_x%pres_out'                 : 1.,
+    # # ========================================================================
+
+    # Patch: Cylinder Immersed Boundary ========================================
+    'patch_ib(1)%geometry'       : 4,
+    'patch_ib(1)%x_centroid'     : 1.5E-03 / x0,
+    'patch_ib(1)%y_centroid'     : 1.5E-03 / x0,
+    'patch_ib(1)%c'     : 1.0E-03 / x0,
+    'patch_ib(1)%t'     : 0.15,
+    'patch_ib(1)%p'     : 0.4,
+    'patch_ib(1)%m'     : 0.02,
+    'patch_ib(1)%slip'     : 'F',
+    'patch_ib(1)%theta'     : 15,
+
+    # 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/(n_tait-1.),
+    #'fluid_pp(1)%Re(1)' 	    : 67567,
+
+
+    # ==========================================================================
+}))