Skip to content

Commit 9f830a4

Browse files
ChrisZYJChrisZYJ
committed
rebase upstream
1 parent 4298db0 commit 9f830a4

File tree

19 files changed

+1186
-166
lines changed

19 files changed

+1186
-166
lines changed

‎docs/documentation/case.md‎

Lines changed: 10 additions & 0 deletions
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 |

‎examples/2D_axisym_hypoelasticity/case.py‎

Lines changed: 100 additions & 0 deletions
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+
)

‎src/common/m_boundary_conditions.fpp‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -354,6 +354,13 @@ contains
354354
q_prim_vf(i)%sf(j - 1, k, l)
355355
end do
356356

357+
if (elasticity) then
358+
do i = 1, shear_BC_flip_num
359+
q_prim_vf(shear_BC_flip_indices(1, i))%sf(-j, k, l) = &
360+
-q_prim_vf(shear_BC_flip_indices(1, i))%sf(j - 1, k, l)
361+
end do
362+
end if
363+
357364
if (hyperelasticity) then
358365
q_prim_vf(xibeg)%sf(-j, k, l) = &
359366
-q_prim_vf(xibeg)%sf(j - 1, k, l)
@@ -404,6 +411,13 @@ contains
404411
q_prim_vf(i)%sf(m - (j - 1), k, l)
405412
end do
406413

414+
if (elasticity) then
415+
do i = 1, shear_BC_flip_num
416+
q_prim_vf(shear_BC_flip_indices(1, i))%sf(m + j, k, l) = &
417+
-q_prim_vf(shear_BC_flip_indices(1, i))%sf(m - (j - 1), k, l)
418+
end do
419+
end if
420+
407421
if (hyperelasticity) then
408422
q_prim_vf(xibeg)%sf(m + j, k, l) = &
409423
-q_prim_vf(xibeg)%sf(m - (j - 1), k, l)
@@ -457,6 +471,13 @@ contains
457471
q_prim_vf(i)%sf(l, j - 1, k)
458472
end do
459473

474+
if (elasticity) then
475+
do i = 1, shear_BC_flip_num
476+
q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, -j, k) = &
477+
-q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, j - 1, k)
478+
end do
479+
end if
480+
460481
if (hyperelasticity) then
461482
q_prim_vf(xibeg + 1)%sf(l, -j, k) = &
462483
-q_prim_vf(xibeg + 1)%sf(l, j - 1, k)
@@ -504,6 +525,13 @@ contains
504525
q_prim_vf(i)%sf(l, n - (j - 1), k)
505526
end do
506527

528+
if (elasticity) then
529+
do i = 1, shear_BC_flip_num
530+
q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, n + j, k) = &
531+
-q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, n - (j - 1), k)
532+
end do
533+
end if
534+
507535
if (hyperelasticity) then
508536
q_prim_vf(xibeg + 1)%sf(l, n + j, k) = &
509537
-q_prim_vf(xibeg + 1)%sf(l, n - (j - 1), k)
@@ -556,6 +584,13 @@ contains
556584
q_prim_vf(i)%sf(k, l, j - 1)
557585
end do
558586

587+
if (elasticity) then
588+
do i = 1, shear_BC_flip_num
589+
q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, -j) = &
590+
-q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, j - 1)
591+
end do
592+
end if
593+
559594
if (hyperelasticity) then
560595
q_prim_vf(xiend)%sf(k, l, -j) = &
561596
-q_prim_vf(xiend)%sf(k, l, j - 1)
@@ -603,6 +638,13 @@ contains
603638
q_prim_vf(i)%sf(k, l, p - (j - 1))
604639
end do
605640

641+
if (elasticity) then
642+
do i = 1, shear_BC_flip_num
643+
q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, p + j) = &
644+
-q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, p - (j - 1))
645+
end do
646+
end if
647+
606648
if (hyperelasticity) then
607649
q_prim_vf(xiend)%sf(k, l, p + j) = &
608650
-q_prim_vf(xiend)%sf(k, l, p - (j - 1))

‎src/common/m_variables_conversion.fpp‎

Lines changed: 7 additions & 13 deletions
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
@@ -1589,7 +1583,7 @@ contains
15891583
end subroutine s_finalize_variables_conversion_module
15901584

15911585
#ifndef MFC_PRE_PROCESS
1592-
pure subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c)
1586+
subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c)
15931587
#ifdef _CRAYFTN
15941588
!DIR$ INLINEALWAYS s_compute_speed_of_sound
15951589
#else

0 commit comments

Comments
 (0)