Skip to content

Commit f3f9a2a

Browse files
sbryngelsonsbryngelson
committed
minor
1 parent e6721b7 commit f3f9a2a

File tree

3 files changed

+628
-0
lines changed

3 files changed

+628
-0
lines changed

‎example_cases/1D_exp_bubscreen/input_shock.py‎

Lines changed: 309 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,309 @@
1+
#!/usr/bin/python
2+
3+
import math
4+
5+
x0 = 61.E-05
6+
p0 = 112900.
7+
rho0 = 960
8+
c0 = math.sqrt( p0/rho0 )
9+
patm = 1.
10+
11+
#water props
12+
## AKA little \gamma (see coralic 2014 eq'n (13))
13+
n_tait = 7.1
14+
## AKA little \pi(see coralic 2014 eq'n (13))
15+
B_tait = 306.E+06 / p0
16+
17+
mul0 = 0.048 #viscosity
18+
# mul0 = 1.E-12
19+
ss = 20.8 #surface tension
20+
# ss = 1.E-12 ## this would turn-off surface tension
21+
pv = 133.322 #vapor pressure
22+
23+
# water
24+
# These _v and _n parameters ONLY correspond to the bubble model of Preston (2010 maybe 2008)
25+
# (this model would replace the usual Rayleigh-plesset or Keller-miksis model (it's more compilcated))
26+
#gamma_v = 1.33
27+
#M_v = 18.02
28+
#mu_v = 0.8816E-05
29+
#k_v = 0.019426
30+
31+
##air props
32+
#gamma_n = 1.4
33+
#M_n = 28.97
34+
#mu_n = 1.8E-05
35+
#k_n = 0.02556
36+
37+
#air props
38+
gamma_gas = 1.09
39+
40+
#reference bubble size
41+
R0ref = 61.E-05
42+
43+
pa = 0.1 * 1.E+06 / 101325.
44+
45+
#Characteristic velocity
46+
uu = math.sqrt( p0/rho0 )
47+
#Cavitation number
48+
# Ca = (p0 - pv)/(rho0*(uu**2.))
49+
Ca = 1
50+
#Weber number
51+
We = rho0*(uu**2.)*R0ref/ss
52+
#Inv. bubble Reynolds number
53+
Re_inv = mul0/(rho0*uu*R0ref)
54+
55+
#IC setup
56+
vf0 = .0024
57+
# vf0 = 1.E-6
58+
n0 = vf0/(math.pi*4.E+00/3.E+00)
59+
60+
cphysical = 981.6
61+
t0 = x0/c0
62+
63+
nbubbles = 1
64+
myr0 = R0ref
65+
66+
67+
# CFL number should be < 1 for numerical stability
68+
# CFL = speed of sound * dt/dx
69+
cfl = 0.2
70+
Nx = 600
71+
Ldomain = 2
72+
L = Ldomain/x0
73+
dx = L/float(Nx)
74+
dt = cfl*dx/(cphysical/c0)
75+
76+
Lpulse = 0.3*Ldomain
77+
Tpulse = Lpulse/cphysical
78+
Tfinal = 0.3*0.25*35.*Tpulse*c0/x0
79+
Nt = int(Tfinal/dt)
80+
81+
Nfiles = 50.
82+
Nout = int(math.ceil(Nt/Nfiles))
83+
Nt = int(Nout*Nfiles)
84+
85+
# Command to navigate between directories
86+
from os import chdir
87+
88+
# Command to acquire directory path
89+
from os.path import dirname
90+
91+
# Command to acquire script name and module search path
92+
from sys import argv, path
93+
94+
# Navigating to script directory
95+
if len(dirname(argv[0])) != 0: chdir(dirname(argv[0]))
96+
97+
# Adding master_scripts directory to module search path
98+
mfc_dir = '../../src'; path[:0] = [mfc_dir + '/master_scripts']
99+
100+
# Command to execute the MFC components
101+
from m_python_proxy import f_execute_mfc_component
102+
103+
# ==============================================================================
104+
105+
# Case Analysis Configuration ==================================================
106+
107+
# Selecting MFC component
108+
comp_name = argv[1].strip()
109+
110+
# Serial or parallel computational engine
111+
engine = 'serial'
112+
if (comp_name=='pre_process'): engine = 'serial'
113+
114+
# Configuring case dictionary
115+
case_dict = \
116+
{ \
117+
# Logistics ================================================
118+
'case_dir' : '\'.\'', \
119+
'run_time_info' : 'T', \
120+
'nodes' : 1, \
121+
# processes per node... > 1 indicates parallel (avoid this for now)
122+
'ppn' : 1, \
123+
'queue' : 'normal', \
124+
'walltime' : '24:00:00', \
125+
'mail_list' : '', \
126+
# ==========================================================
127+
\
128+
# Computational Domain Parameters ==========================
129+
'x_domain%beg' : 0/x0, \
130+
'x_domain%end' : 2/x0, \
131+
'stretch_x' : 'F', \
132+
'cyl_coord' : 'F', \
133+
'm' : Nx, \
134+
'n' : 0, \
135+
'p' : 0, \
136+
'dt' : dt, \
137+
't_step_start' : 0, \
138+
't_step_stop' : Nt, \
139+
't_step_save' : Nout, \
140+
# ==========================================================
141+
\
142+
# Simulation Algorithm Parameters ==========================
143+
'num_patches' : 3, \
144+
'model_eqns' : 2, \
145+
'alt_soundspeed' : 'F', \
146+
'num_fluids' : 1, \
147+
'adv_alphan' : 'T', \
148+
'mpp_lim' : 'F', \
149+
'mixture_err' : 'F', \
150+
'time_stepper' : 3, \
151+
'weno_vars' : 2, \
152+
'weno_order' : 5, \
153+
'weno_eps' : 1.E-16, \
154+
'char_decomp' : 'F', \
155+
'mapped_weno' : 'T', \
156+
'null_weights' : 'F', \
157+
'mp_weno' : 'T', \
158+
'riemann_solver' : 2, \
159+
'wave_speeds' : 1, \
160+
'avg_state' : 2, \
161+
'commute_err' : 'F', \
162+
'split_err' : 'F', \
163+
'bc_x%beg' : -8, \
164+
'bc_x%end' : -8, \
165+
# ==========================================================
166+
\
167+
# Formatted Database Files Structure Parameters ============
168+
'format' : 1, \
169+
'precision' : 2, \
170+
'prim_vars_wrt' :'T', \
171+
'parallel_io' :'F', \
172+
'fd_order' : 1, \
173+
'schlieren_wrt' :'T', \
174+
'probe_wrt' :'T', \
175+
'num_probes' : 1, \
176+
'probe(1)%x' : 1.962, \
177+
# ==========================================================
178+
179+
# Patch 1 _ Background =====================================
180+
# this problem is 1D... so based on the dimension of the problem
181+
# you have different 'geometries' available to you
182+
# e.g. in 3D you might have spherical geometries
183+
# and rectangular ones
184+
# in 1D (like here)... there is only one option {#1}... which is a
185+
# line
186+
'patch_icpp(1)%geometry' : 1, \
187+
'patch_icpp(1)%x_centroid' : 1/x0, \
188+
'patch_icpp(1)%length_x' : 2/x0, \
189+
'patch_icpp(1)%vel(1)' : 0.0, \
190+
'patch_icpp(1)%pres' : patm, \
191+
# \alpha stands for volume fraction of this phase
192+
# so if there are no bubbles, then it is all water (liquid)
193+
# and \alpha_1 = \alpha_liquid \approx 1
194+
'patch_icpp(1)%alpha_rho(1)' : (1.-1.E-12)*(1.E+03/rho0), \
195+
# \alpha_1 here is always (for num_fluids = 1 and bubbles=True)
196+
# \alpha is always the void fraction of bubbles (usually << 1)
197+
'patch_icpp(1)%alpha(1)' : 1.E-12, \
198+
# dimensionless initial bubble radius
199+
'patch_icpp(1)%r0' : 1., \
200+
# dimensionless initial velocity
201+
'patch_icpp(1)%v0' : 0.0E+00, \
202+
# ==========================================================
203+
204+
# Patch 2 Screen ===========================================
205+
'patch_icpp(2)%geometry' : 1, \
206+
#overwrite the part in the middle that was the
207+
#background (no bubble) area
208+
'patch_icpp(2)%alter_patch(1)' : 'T', \
209+
'patch_icpp(2)%x_centroid' : 1/x0, \
210+
'patch_icpp(2)%length_x' : 0.5/x0, \
211+
'patch_icpp(2)%vel(1)' : 0.0, \
212+
'patch_icpp(2)%pres' : patm, \
213+
# \alpha stands for volume fraction of this phase
214+
# so if there are no bubbles, then it is all water (liquid)
215+
# and \alpha_1 = \alpha_liquid \approx 1
216+
# in the screen case, you have \alpha_1 = 1 - \alpha_bubbles = 1 - vf0
217+
'patch_icpp(2)%alpha_rho(1)' : (1.-vf0)*1.E+03/rho0, \
218+
# void fraction of bubbles
219+
'patch_icpp(2)%alpha(1)' : vf0, \
220+
'patch_icpp(2)%r0' : 1., \
221+
'patch_icpp(2)%v0' : 0.0E+00, \
222+
# ==========================================================
223+
224+
# Patch 3 Shock Wave
225+
'patch_icpp(3)%geometry' : 1, \
226+
'patch_icpp(3)%alter_patch(1)' : 'T', \
227+
'patch_icpp(3)%x_centroid' : 0.2/x0, \
228+
'patch_icpp(3)%length_x' : 0.4/x0, \
229+
'patch_icpp(3)%vel(1)' : 1.4*cphysical, \
230+
'patch_icpp(3)%pres' : 2.157*patm, \
231+
'patch_icpp(3)%alpha_rho(1)' : (1.-1.E-12)*(1.E+03/rho0), \
232+
'patch_icpp(3)%alpha(1)' : 1.E-12, \
233+
'patch_icpp(3)%r0' : 1., \
234+
'patch_icpp(3)%v0' : 0.0E+00, \
235+
# ==========================================================
236+
237+
238+
# Fluids Physical Parameters ===============================
239+
# Surrounding liquid
240+
'fluid_pp(1)%gamma' : 1.E+00/(n_tait-1.E+00), \
241+
'fluid_pp(1)%pi_inf' : n_tait*B_tait/(n_tait-1.), \
242+
# 'fluid_pp(1)%mul0' : mul0, \
243+
# 'fluid_pp(1)%ss' : ss, \
244+
# 'fluid_pp(1)%pv' : pv, \
245+
# 'fluid_pp(1)%gamma_v' : gamma_v, \
246+
# 'fluid_pp(1)%M_v' : M_v, \
247+
# 'fluid_pp(1)%mu_v' : mu_v, \
248+
# 'fluid_pp(1)%k_v' : k_v, \
249+
250+
# Last fluid_pp is always reserved for bubble gas state ===
251+
# if applicable ==========================================
252+
'fluid_pp(2)%gamma' : 1./(gamma_gas-1.), \
253+
'fluid_pp(2)%pi_inf' : 0.0E+00, \
254+
# 'fluid_pp(2)%gamma_v' : gamma_n, \
255+
# 'fluid_pp(2)%M_v' : M_n, \
256+
# 'fluid_pp(2)%mu_v' : mu_n, \
257+
# 'fluid_pp(2)%k_v' : k_n, \
258+
# ==========================================================
259+
260+
261+
# Non-polytropic gas compression model AND/OR Tait EOS =====
262+
'pref' : p0, \
263+
'rhoref' : rho0, \
264+
# ==========================================================
265+
266+
# Bubbles ==================================================
267+
'bubbles' : 'T', \
268+
# in user guide... 1 = gilbert 2 = keller-miksis
269+
# but gilbert won't work for the equations that you are using... (i think)
270+
'bubble_model' : 2, \
271+
# polytropic: this is where the different between Rayleigh--Plesset and
272+
# Preston's model shows up. polytropic = False means complicated Preston model
273+
# = True means simpler Rayleigh--Plesset model
274+
# if polytropic == False then you will end up calling s_initialize_nonpoly in
275+
# m_global_parameters.f90 in both the pre_process and simulation_code
276+
'polytropic' : 'T', \
277+
'polydisperse' : 'F', \
278+
#'poly_sigma' : 0.3, \
279+
# only matters if polytropic = False (complicated model)
280+
'thermal' : 3, \
281+
# only matters if polytropic = False (complicated model)
282+
'R0ref' : myr0, \
283+
'nb' : 1, \
284+
# cavitation number (has something to do with the ratio of gas to vapour in the bubble)
285+
# this is usually near 1
286+
# can set = 1 for testing purposes
287+
'Ca' : Ca, \
288+
# weber number (corresponds to surface tension)
289+
'Web' : We, \
290+
# inverse reynolds number (coresponds to viscosity)
291+
'Re_inv' : Re_inv, \
292+
# ==========================================================
293+
294+
# Acoustic source ==========================================
295+
#'Monopole' : 'T', \
296+
'num_mono' : 1, \
297+
'Mono(1)%loc(1)' : -0.3/x0, \
298+
'Mono(1)%npulse' : 1, \
299+
'Mono(1)%dir' : 1., \
300+
'Mono(1)%pulse' : 1, \
301+
'Mono(1)%mag' : 0.001, \
302+
'Mono(1)%length' : (1./(10000.))*cphysical/x0, \
303+
# ==========================================================
304+
}
305+
306+
# Executing MFC component
307+
f_execute_mfc_component(comp_name, case_dict, mfc_dir, engine)
308+
309+
# ==============================================================================

0 commit comments

Comments
 (0)