Skip to content

Commit a2a8800

Browse files
ChrisZYJChrisZYJ
authored
Axisymmetric Hypoelasticity (#815)
1 parent 775d08d commit a2a8800

File tree

14 files changed

+573
-29
lines changed

14 files changed

+573
-29
lines changed

Diff for: ‎docs/documentation/case.md

+10
Original file line numberDiff line numberDiff line change
@@ -164,6 +164,7 @@ MPI topology is automatically optimized to maximize the parallel efficiency for
164164
| `alpha_rho(i)` * | Real | Supported | Partial density of fluid $i$. |
165165
| `pres` * | Real | Supported | Pressure. |
166166
| `vel(i)` * | Real | Supported | Velocity in direction $i$. |
167+
| `tau_e(i)` * | Real | Supported | Elastic stresses. |
167168
| `hcid` * | Integer | N/A | Hard coded patch id |
168169
| `cf_val` * | Real | Supported | Surface tension color function value |
169170
| `model_filepath` | String | Not Supported | Path to an STL or OBJ file (not all OBJs are supported). |
@@ -180,6 +181,8 @@ The parameters define the geometries and physical parameters of fluid components
180181
Note that the domain must be fully filled with patche(s).
181182
The code outputs error messages when an empty region is left in the domain.
182183

184+
- `tau_e(i)` is the `i`-th component of the elastic stress tensor, ordered as `tau_xx`, `tau_xy`, `tau_yy`, `tau_xz`, `tau_yz`, and `tau_zz`. 1D simulation requires `tau(1)`, 2D `tau(1:3)`, and 3D `tau(1:6)`.
185+
183186
#### Analytical Definition of Primitive Variables
184187

185188
Some parameters, as described above, can be defined by analytical functions in the input file. For example, one can define the following patch:
@@ -330,6 +333,7 @@ Additional details on this specification can be found in [The Naca Airfoil Serie
330333
| `qv` ** | Real | Stiffened-gas parameter $q$ of fluid. |
331334
| `qvp` ** | Real | Stiffened-gas parameter $q'$ of fluid. |
332335
| `sigma` | Real | Surface tension coefficient |
336+
| `G` | Real | Shear modulus of solid. |
333337

334338
Fluid material's parameters. All parameters except for sigma should be prepended with `fluid_pp(i)` where $i$ is the fluid index.
335339

@@ -349,6 +353,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
349353

350354
- `fluid_pp(i)%%cv`, `fluid_pp(i)%%qv`, and `fluid_pp(i)%%qvp` define $c_v$, $q$, and $q'$ as parameters of $i$-th fluid that are used in stiffened gas equation of state.
351355

356+
- `fluid_pp(i)%%G` is required for `hypoelasticity`.
357+
352358
### 6. Simulation Algorithm
353359

354360
| Parameter | Type | Description |
@@ -391,6 +397,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
391397
| `t_stop` | Real | Simulation stop time |
392398
| `surface_tension` | Logical | Activate surface tension |
393399
| `viscous` | Logical | Activate viscosity |
400+
| `hypoelasticity` | Logical | Activate hypoelasticity* |
394401

395402
- \* Options that work only with `model_eqns = 2`.
396403
- † Options that work only with ``cyl_coord = 'F'``.
@@ -470,6 +477,8 @@ This option requires `weno_Re_flux` to be true because cell boundary values are
470477

471478
- `viscous` activates viscosity when set to ``'T'``. Requires `Re(1)` and `Re(2)` to be set.
472479

480+
- `hypoelasticity` activates elastic stress calculations for fluid-solid interactions. Requires `G` to be set in the fluid material's parameters.
481+
473482
#### Constant Time-Stepping
474483

475484
- `dt` specifies the constant time step size used in the simulation.
@@ -523,6 +532,7 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running.
523532
| `omega_wrt(i)` | Logical | Add the $i$-direction vorticity to the database |
524533
| `schlieren_wrt` | Logical | Add the numerical schlieren to the database|
525534
| `qm_wrt` | Logical | Add the Q-criterion to the database|
535+
| `tau_wrt` | Logical | Add the elastic stress components to the database|
526536
| `fd_order` | Integer | Order of finite differences for computing the vorticity and the numerical Schlieren function [1,2,4] |
527537
| `schlieren_alpha(i)` | Real | Intensity of the numerical Schlieren computed via `alpha(i)` |
528538
| `probe_wrt` | Logical | Write the flow chosen probes data files for each time step |

Diff for: ‎examples/2D_axisym_hypoelasticity/case.py

+100
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
#!/usr/bin/env python3
2+
import json
3+
4+
# Configuring case dictionary
5+
print(
6+
json.dumps(
7+
{
8+
# Logistics
9+
"run_time_info": "T",
10+
# Computational Domain Parameters
11+
"x_domain%beg": 0.0,
12+
"x_domain%end": 0.001,
13+
"y_domain%beg": 0.0,
14+
"y_domain%end": 0.0005,
15+
"cyl_coord": "T",
16+
"m": 100,
17+
"n": 50,
18+
"p": 0,
19+
"dt": 4e-12,
20+
"t_step_start": 0,
21+
"t_step_stop": 40000,
22+
"t_step_save": 2000,
23+
# Simulation Algorithm Parameters
24+
"num_patches": 2,
25+
"model_eqns": 2,
26+
"alt_soundspeed": "F",
27+
"num_fluids": 2,
28+
"mpp_lim": "F",
29+
"mixture_err": "F",
30+
"time_stepper": 3,
31+
"weno_order": 5,
32+
"weno_eps": 1.0e-16,
33+
"mapped_weno": "T",
34+
"null_weights": "F",
35+
"mp_weno": "F",
36+
"riemann_solver": 1,
37+
"wave_speeds": 1,
38+
"avg_state": 2,
39+
"bc_x%beg": -6,
40+
"bc_x%end": -6,
41+
"bc_y%beg": -2,
42+
"bc_y%end": -6,
43+
# Hypoelasticity
44+
"hypoelasticity": "T",
45+
"fd_order": 4,
46+
# Formatted Database Files Structure Parameters
47+
"format": 1,
48+
"precision": 2,
49+
"prim_vars_wrt": "T",
50+
"parallel_io": "T",
51+
# Patch 1 Liquid
52+
"patch_icpp(1)%geometry": 3,
53+
"patch_icpp(1)%x_centroid": 0.0005,
54+
"patch_icpp(1)%y_centroid": 0.00025,
55+
"patch_icpp(1)%length_x": 0.001,
56+
"patch_icpp(1)%length_y": 0.0005,
57+
"patch_icpp(1)%vel(1)": 0.0,
58+
"patch_icpp(1)%vel(2)": 0.0,
59+
"patch_icpp(1)%pres": 1e05,
60+
"patch_icpp(1)%alpha_rho(1)": 1100 * (1.0 - 1e-6),
61+
"patch_icpp(1)%alpha(1)": 1.0 - 1e-6,
62+
"patch_icpp(1)%alpha_rho(2)": 1100 * 1e-6,
63+
"patch_icpp(1)%alpha(2)": 1e-6,
64+
# Patch 2 Solid
65+
"patch_icpp(2)%alter_patch(1)": "T",
66+
"patch_icpp(2)%geometry": 3,
67+
"patch_icpp(2)%x_centroid": 0.0005,
68+
"patch_icpp(2)%y_centroid": 0.000125,
69+
"patch_icpp(2)%length_x": 0.0005,
70+
"patch_icpp(2)%length_y": 0.00025,
71+
"patch_icpp(2)%vel(1)": 0.0,
72+
"patch_icpp(2)%vel(2)": 0.0,
73+
"patch_icpp(2)%pres": 1e05,
74+
"patch_icpp(2)%alpha_rho(1)": 1100 * 1e-6,
75+
"patch_icpp(2)%alpha(1)": 1e-6,
76+
"patch_icpp(2)%alpha_rho(2)": 1100 * (1.0 - 1e-6),
77+
"patch_icpp(2)%alpha(2)": 1.0 - 1e-6,
78+
# Acoustic source
79+
"acoustic_source": "T",
80+
"num_source": 1,
81+
"acoustic(1)%support": 6,
82+
"acoustic(1)%loc(1)": 0.00006,
83+
"acoustic(1)%loc(2)": 0.0,
84+
"acoustic(1)%pulse": 2,
85+
"acoustic(1)%npulse": 1,
86+
"acoustic(1)%mag": 1.0,
87+
"acoustic(1)%gauss_sigma_time": 2e-8,
88+
"acoustic(1)%foc_length": 0.00054,
89+
"acoustic(1)%aperture": 0.0008,
90+
"acoustic(1)%delay": 1e-7,
91+
# Fluids Physical Parameters
92+
"fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
93+
"fluid_pp(1)%pi_inf": 4.4e00 * 5.57e08 / (4.4e00 - 1.0e00),
94+
"fluid_pp(1)%G": 0.0,
95+
"fluid_pp(2)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
96+
"fluid_pp(2)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
97+
"fluid_pp(2)%G": 1.0e09,
98+
}
99+
)
100+
)

Diff for: ‎src/common/m_variables_conversion.fpp

+6-12
Original file line numberDiff line numberDiff line change
@@ -162,10 +162,8 @@ contains
162162
do s = stress_idx%beg, stress_idx%end
163163
if (G > 0) then
164164
E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
165-
! Additional terms in 2D and 3D
166-
if ((s == stress_idx%beg + 1) .or. &
167-
(s == stress_idx%beg + 3) .or. &
168-
(s == stress_idx%beg + 4)) then
165+
! Double for shear stresses
166+
if (any(s == shear_indices)) then
169167
E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
170168
end if
171169
end if
@@ -1126,10 +1124,8 @@ contains
11261124
if (G_K > verysmall) then
11271125
qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
11281126
((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
1129-
! extra terms in 2 and 3D
1130-
if ((i == strxb + 1) .or. &
1131-
(i == strxb + 3) .or. &
1132-
(i == strxb + 4)) then
1127+
! Double for shear stresses
1128+
if (any(i == shear_indices)) then
11331129
qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
11341130
((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
11351131
end if
@@ -1396,10 +1392,8 @@ contains
13961392
if (G > verysmall) then
13971393
q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
13981394
(q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
1399-
! extra terms in 2 and 3D
1400-
if ((i == stress_idx%beg + 1) .or. &
1401-
(i == stress_idx%beg + 3) .or. &
1402-
(i == stress_idx%beg + 4)) then
1395+
! Double for shear stresses
1396+
if (any(i == shear_indices)) then
14031397
q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
14041398
(q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
14051399
end if

Diff for: ‎src/post_process/m_global_parameters.fpp

+2-1
Original file line numberDiff line numberDiff line change
@@ -674,7 +674,8 @@ contains
674674
elasticity = .true.
675675
stress_idx%beg = sys_size + 1
676676
stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
677-
! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
677+
if (cyl_coord) stress_idx%end = stress_idx%end + 1
678+
! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
678679
sys_size = stress_idx%end
679680

680681
! shear stress index is 2 for 2D and 2,4,5 for 3D

Diff for: ‎src/pre_process/m_global_parameters.fpp

+2-1
Original file line numberDiff line numberDiff line change
@@ -767,7 +767,8 @@ contains
767767
elasticity = .true.
768768
stress_idx%beg = sys_size + 1
769769
stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
770-
! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
770+
if (cyl_coord) stress_idx%end = stress_idx%end + 1
771+
! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
771772
sys_size = stress_idx%end
772773

773774
! shear stress index is 2 for 2D and 2,4,5 for 3D

Diff for: ‎src/simulation/m_checker.fpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -164,9 +164,9 @@ contains
164164
"Only acoustic("//trim(jStr)//")%support = 1 is allowed for 1D simulations")
165165
@:PROHIBIT(dim == 1 .and. acoustic(j)%support == 1 .and. f_is_default(acoustic(j)%loc(1)), &
166166
"acoustic("//trim(jStr)//")%loc(1) must be specified for acoustic("//trim(jStr)//")%support = 1")
167-
@:PROHIBIT(dim == 2 .and. (.not. any(acoustic(j)%support == (/2, 5, 6, 9, 10/))), &
167+
@:PROHIBIT((dim == 2 .and. .not. cyl_coord) .and. (.not. any(acoustic(j)%support == (/2, 5, 9/))), &
168168
"Only acoustic("//trim(jStr)//")%support = 2, 5, 6, 9, or 10 is allowed for 2D simulations")
169-
@:PROHIBIT(dim == 2 .and. (.not. any(acoustic(j)%support == (/6, 10/))) .and. cyl_coord, &
169+
@:PROHIBIT((dim == 2 .and. cyl_coord) .and. (.not. any(acoustic(j)%support == (/2, 6, 10/))), &
170170
"Only acoustic("//trim(jStr)//")%support = 6 or 10 is allowed for 2D axisymmetric simulations")
171171
@:PROHIBIT(dim == 2 .and. any(acoustic(j)%support == (/2, 5, 6, 9, 10/)) .and. &
172172
(f_is_default(acoustic(j)%loc(1)) .or. f_is_default(acoustic(j)%loc(2))), &

Diff for: ‎src/simulation/m_global_parameters.fpp

+2-1
Original file line numberDiff line numberDiff line change
@@ -1039,7 +1039,8 @@ contains
10391039
elasticity = .true.
10401040
stress_idx%beg = sys_size + 1
10411041
stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
1042-
! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
1042+
if (cyl_coord) stress_idx%end = stress_idx%end + 1
1043+
! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
10431044
sys_size = stress_idx%end
10441045
10451046
! shear stress index is 2 for 2D and 2,4,5 for 3D

Diff for: ‎src/simulation/m_hypoelastic.fpp

+33
Original file line numberDiff line numberDiff line change
@@ -331,6 +331,39 @@ contains
331331
end do
332332
end if
333333

334+
if (cyl_coord .and. idir == 2) then
335+
336+
!$acc parallel loop collapse(3) gang vector default(present)
337+
do q = 0, p
338+
do l = 0, n
339+
do k = 0, m
340+
! S_xx -= rho * v/r * (tau_xx + 2/3*G)
341+
rhs_vf(strxb)%sf(k, l, q) = rhs_vf(strxb)%sf(k, l, q) - &
342+
rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
343+
(q_prim_vf(strxb)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q)) ! tau_xx + 2/3*G
344+
345+
! S_xr -= rho * v/r * tau_xr
346+
rhs_vf(strxb + 1)%sf(k, l, q) = rhs_vf(strxb + 1)%sf(k, l, q) - &
347+
rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
348+
q_prim_vf(strxb + 1)%sf(k, l, q) ! tau_xx
349+
350+
! S_rr -= rho * v/r * (tau_rr + 2/3*G)
351+
rhs_vf(strxb + 2)%sf(k, l, q) = rhs_vf(strxb + 2)%sf(k, l, q) - &
352+
rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
353+
(q_prim_vf(strxb + 2)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q)) ! tau_rr + 2/3*G
354+
355+
! S_thetatheta += rho * ( -(tau_thetatheta + 2/3*G)*(du/dx + dv/dr + v/r) + 2*(tau_thetatheta + G)*v/r )
356+
rhs_vf(strxb + 3)%sf(k, l, q) = rhs_vf(strxb + 3)%sf(k, l, q) + &
357+
rho_K_field(k, l, q)*( &
358+
-(q_prim_vf(strxb + 3)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q))* &
359+
(du_dx(k, l, q) + dv_dy(k, l, q) + q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)) &
360+
+ 2._wp*(q_prim_vf(strxb + 3)%sf(k, l, q) + G_K_field(k, l, q))*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l))
361+
end do
362+
end do
363+
end do
364+
365+
end if
366+
334367
end subroutine s_compute_hypoelastic_rhs
335368

336369
subroutine s_finalize_hypoelastic_module()

Diff for: ‎src/simulation/m_riemann_solvers.fpp

+18-10
Original file line numberDiff line numberDiff line change
@@ -612,8 +612,8 @@ contains
612612
if ((G_L > 1000) .and. (G_R > 1000)) then
613613
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
614614
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
615-
! Additional terms in 2D and 3D
616-
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
615+
! Double for shear stresses
616+
if (any(strxb - 1 + i == shear_indices)) then
617617
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
618618
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
619619
end if
@@ -1060,20 +1060,28 @@ contains
10601060
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
10611061
end do
10621062
! Recalculating the radial momentum geometric source flux
1063-
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1)) = &
1064-
(s_M*(rho_R*vel_R(dir_idx(1)) &
1065-
*vel_R(dir_idx(1))) &
1066-
- s_P*(rho_L*vel_L(dir_idx(1)) &
1067-
*vel_L(dir_idx(1))) &
1068-
+ s_M*s_P*(rho_L*vel_L(dir_idx(1)) &
1069-
- rho_R*vel_R(dir_idx(1)))) &
1070-
/(s_M - s_P)
1063+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) = &
1064+
flux_rs${XYZ}$_vf(j, k, l, contxe + 2) &
1065+
- (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
10711066
! Geometrical source of the void fraction(s) is zero
10721067
!$acc loop seq
10731068
do i = advxb, advxe
10741069
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
10751070
end do
10761071
end if
1072+
1073+
if (cyl_coord .and. hypoelasticity) then
1074+
! += tau_sigmasigma using HLL
1075+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) = &
1076+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) + &
1077+
(s_M*tau_e_R(4) - s_P*tau_e_L(4)) &
1078+
/(s_M - s_P)
1079+
1080+
!$acc loop seq
1081+
do i = strxb, strxe
1082+
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
1083+
end do
1084+
end if
10771085
#:endif
10781086
10791087
end do

0 commit comments

Comments
 (0)