Axisymmetric Hypoelasticity

docs/documentation/case.md

@@ -164,6 +164,7 @@ MPI topology is automatically optimized to maximize the parallel efficiency for
 | `alpha_rho(i)` *     | Real    | Supported             | Partial density of fluid $i$.                                |
 | `pres` *             | Real    | Supported             | Pressure.                                                    |
 | `vel(i)` *           | Real    | Supported             | Velocity in direction $i$.                                   |
+| `tau_e(i)` *         | Real    | Supported             | Elastic stresses.                                            |
 | `hcid` *             | Integer | N/A                   | Hard coded patch id                                          |
 | `cf_val` *           | Real    | Supported             | Surface tension color function value                         |
 | `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
 Note that the domain must be fully filled with patche(s).
 The code outputs error messages when an empty region is left in the domain.
 
+- `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)`.
+
 #### Analytical Definition of Primitive Variables
 
 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
 | `qv`   ** | Real   | Stiffened-gas parameter $q$ of fluid.          |
 | `qvp`  ** | Real   | Stiffened-gas parameter $q'$ of fluid.         |
 | `sigma`   | Real   | Surface tension coefficient                    |
+| `G`       | Real   | Shear modulus of solid.                        |
 
 Fluid material's parameters. All parameters except for sigma should be prepended with `fluid_pp(i)` where $i$ is the fluid index.
 
@@ -349,6 +353,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 
 - `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.
 
+- `fluid_pp(i)%%G` is required for `hypoelasticity`.
+
 ### 6. Simulation Algorithm
 
 | Parameter              | Type    | Description                                    |
@@ -391,6 +397,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `t_stop`               | Real    | Simulation stop time |
 | `surface_tension`      | Logical | Activate surface tension |
 | `viscous`              | Logical | Activate viscosity |
+| `hypoelasticity`       | Logical | Activate hypoelasticity* |
 
 - \* Options that work only with `model_eqns = 2`.
 - † 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
 
 - `viscous` activates viscosity when set to ``'T'``. Requires `Re(1)` and `Re(2)` to be set.
 
+- `hypoelasticity` activates elastic stress calculations for fluid-solid interactions. Requires `G` to be set in the fluid material's parameters.
+
 #### Constant Time-Stepping
 
 - `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.
 | `omega_wrt(i)`       | Logical | Add the $i$-direction vorticity to the database	 |
 | `schlieren_wrt`      | Logical | Add the numerical schlieren to the database|
 | `qm_wrt`             | Logical | Add the Q-criterion to the database|
+| `tau_wrt`            | Logical | Add the elastic stress components to the database|
 | `fd_order`           | Integer | Order of finite differences for computing the vorticity and the numerical Schlieren function [1,2,4] |
 | `schlieren_alpha(i)` | Real    | Intensity of the numerical Schlieren computed via `alpha(i)` |
 | `probe_wrt`          | Logical | Write the flow chosen probes data files for each time step	|

examples/2D_axisym_hypoelasticity/case.py

@@ -0,0 +1,100 @@
+#!/usr/bin/env python3
+import json
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": 0.001,
+            "y_domain%beg": 0.0,
+            "y_domain%end": 0.0005,
+            "cyl_coord": "T",
+            "m": 100,
+            "n": 50,
+            "p": 0,
+            "dt": 4e-12,
+            "t_step_start": 0,
+            "t_step_stop": 40000,
+            "t_step_save": 2000,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "T",
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 1,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -6,
+            "bc_x%end": -6,
+            "bc_y%beg": -2,
+            "bc_y%end": -6,
+            # Hypoelasticity
+            "hypoelasticity": "T",
+            "fd_order": 4,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Patch 1 Liquid
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.0005,
+            "patch_icpp(1)%y_centroid": 0.00025,
+            "patch_icpp(1)%length_x": 0.001,
+            "patch_icpp(1)%length_y": 0.0005,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 1e05,
+            "patch_icpp(1)%alpha_rho(1)": 1100 * (1.0 - 1e-6),
+            "patch_icpp(1)%alpha(1)": 1.0 - 1e-6,
+            "patch_icpp(1)%alpha_rho(2)": 1100 * 1e-6,
+            "patch_icpp(1)%alpha(2)": 1e-6,
+            # Patch 2 Solid
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%x_centroid": 0.0005,
+            "patch_icpp(2)%y_centroid": 0.000125,
+            "patch_icpp(2)%length_x": 0.0005,
+            "patch_icpp(2)%length_y": 0.00025,
+            "patch_icpp(2)%vel(1)": 0.0,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": 1e05,
+            "patch_icpp(2)%alpha_rho(1)": 1100 * 1e-6,
+            "patch_icpp(2)%alpha(1)": 1e-6,
+            "patch_icpp(2)%alpha_rho(2)": 1100 * (1.0 - 1e-6),
+            "patch_icpp(2)%alpha(2)": 1.0 - 1e-6,
+            # Acoustic source
+            "acoustic_source": "T",
+            "num_source": 1,
+            "acoustic(1)%support": 6,
+            "acoustic(1)%loc(1)": 0.00006,
+            "acoustic(1)%loc(2)": 0.0,
+            "acoustic(1)%pulse": 2,
+            "acoustic(1)%npulse": 1,
+            "acoustic(1)%mag": 1.0,
+            "acoustic(1)%gauss_sigma_time": 2e-8,
+            "acoustic(1)%foc_length": 0.00054,
+            "acoustic(1)%aperture": 0.0008,
+            "acoustic(1)%delay": 1e-7,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 4.4e00 * 5.57e08 / (4.4e00 - 1.0e00),
+            "fluid_pp(1)%G": 0.0,
+            "fluid_pp(2)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
+            "fluid_pp(2)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
+            "fluid_pp(2)%G": 1.0e09,
+        }
+    )
+)

src/common/m_variables_conversion.fpp

@@ -162,10 +162,8 @@ contains
                 do s = stress_idx%beg, stress_idx%end
                     if (G > 0) then
                         E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
-                        ! Additional terms in 2D and 3D
-                        if ((s == stress_idx%beg + 1) .or. &
-                            (s == stress_idx%beg + 3) .or. &
-                            (s == stress_idx%beg + 4)) then
+                        ! Double for shear stresses
+                        if (any(s == shear_indices)) then
                             E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
                         end if
                     end if
@@ -1126,10 +1124,8 @@ contains
                             if (G_K > verysmall) then
                                 qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
                                                                 ((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
-                                ! extra terms in 2 and 3D
-                                if ((i == strxb + 1) .or. &
-                                    (i == strxb + 3) .or. &
-                                    (i == strxb + 4)) then
+                                ! Double for shear stresses
+                                if (any(i == shear_indices)) then
                                     qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
                                                                     ((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
                                 end if
@@ -1396,10 +1392,8 @@ contains
                             if (G > verysmall) then
                                 q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
                                                                (q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
-                                ! extra terms in 2 and 3D
-                                if ((i == stress_idx%beg + 1) .or. &
-                                    (i == stress_idx%beg + 3) .or. &
-                                    (i == stress_idx%beg + 4)) then
+                                ! Double for shear stresses
+                                if (any(i == shear_indices)) then
                                     q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
                                                                    (q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
                                 end if

src/post_process/m_global_parameters.fpp

@@ -674,7 +674,8 @@ contains
                 elasticity = .true.
                 stress_idx%beg = sys_size + 1
                 stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
-                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
+                if (cyl_coord) stress_idx%end = stress_idx%end + 1
+                ! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
                 sys_size = stress_idx%end
 
                 ! shear stress index is 2 for 2D and 2,4,5 for 3D

src/pre_process/m_global_parameters.fpp

@@ -767,7 +767,8 @@ contains
                 elasticity = .true.
                 stress_idx%beg = sys_size + 1
                 stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
-                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
+                if (cyl_coord) stress_idx%end = stress_idx%end + 1
+                ! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
                 sys_size = stress_idx%end
 
                 ! shear stress index is 2 for 2D and 2,4,5 for 3D

src/simulation/m_checker.fpp

@@ -164,9 +164,9 @@
                 "Only acoustic("//trim(jStr)//")%support = 1 is allowed for 1D simulations")
             @:PROHIBIT(dim == 1 .and. acoustic(j)%support == 1 .and. f_is_default(acoustic(j)%loc(1)), &
                 "acoustic("//trim(jStr)//")%loc(1) must be specified for acoustic("//trim(jStr)//")%support = 1")
-            @:PROHIBIT(dim == 2 .and. (.not. any(acoustic(j)%support == (/2, 5, 6, 9, 10/))), &
+            @:PROHIBIT((dim == 2  .and. .not. cyl_coord) .and. (.not. any(acoustic(j)%support == (/2, 5, 9/))), &
                 "Only acoustic("//trim(jStr)//")%support = 2, 5, 6, 9, or 10 is allowed for 2D simulations")
-            @:PROHIBIT(dim == 2 .and. (.not. any(acoustic(j)%support == (/6, 10/))) .and. cyl_coord, &
+            @:PROHIBIT((dim == 2  .and. cyl_coord) .and. (.not. any(acoustic(j)%support == (/2, 6, 10/))), &
                 "Only acoustic("//trim(jStr)//")%support = 6 or 10 is allowed for 2D axisymmetric simulations")
             @:PROHIBIT(dim == 2 .and. any(acoustic(j)%support == (/2, 5, 6, 9, 10/)) .and. &
                 (f_is_default(acoustic(j)%loc(1)) .or. f_is_default(acoustic(j)%loc(2))), &

src/simulation/m_global_parameters.fpp

@@ -1039,7 +1039,8 @@ contains
                 elasticity = .true.
                 stress_idx%beg = sys_size + 1
                 stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
-                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
+                if (cyl_coord) stress_idx%end = stress_idx%end + 1
+                ! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
                 sys_size = stress_idx%end
 
                 ! shear stress index is 2 for 2D and 2,4,5 for 3D

src/simulation/m_hypoelastic.fpp

@@ -331,6 +331,39 @@
             end do
         end if
 
+        if (cyl_coord .and. idir == 2) then
+
+            !$acc parallel loop collapse(3) gang vector default(present)
+            do q = 0, p
+                do l = 0, n
+                    do k = 0, m
+                        ! S_xx -= rho * v/r * (tau_xx + 2/3*G)
+                        rhs_vf(strxb)%sf(k, l, q) = rhs_vf(strxb)%sf(k, l, q) - &
+                                                    rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
+                                                    (q_prim_vf(strxb)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q)) ! tau_xx + 2/3*G
+
+                        ! S_xr -= rho * v/r * tau_xr
+                        rhs_vf(strxb + 1)%sf(k, l, q) = rhs_vf(strxb + 1)%sf(k, l, q) - &
+                                                        rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
+                                                        q_prim_vf(strxb + 1)%sf(k, l, q) ! tau_xx
+
+                        ! S_rr -= rho * v/r * (tau_rr + 2/3*G)
+                        rhs_vf(strxb + 2)%sf(k, l, q) = rhs_vf(strxb + 2)%sf(k, l, q) - &
+                                                        rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
+                                                        (q_prim_vf(strxb + 2)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q)) ! tau_rr + 2/3*G
+
+                        ! S_thetatheta += rho * ( -(tau_thetatheta + 2/3*G)*(du/dx + dv/dr + v/r) + 2*(tau_thetatheta + G)*v/r )
+                        rhs_vf(strxb + 3)%sf(k, l, q) = rhs_vf(strxb + 3)%sf(k, l, q) + &
+                                                        rho_K_field(k, l, q)*( &
+                                                        -(q_prim_vf(strxb + 3)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q))* &
+                                                        (du_dx(k, l, q) + dv_dy(k, l, q) + q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)) &
+                                                        + 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))
+                    end do
+                end do
+            end do
+
+        end if
+
     end subroutine s_compute_hypoelastic_rhs
 
     subroutine s_finalize_hypoelastic_module()

src/simulation/m_riemann_solvers.fpp

@@ -612,8 +612,8 @@
                                     if ((G_L > 1000) .and. (G_R > 1000)) then
                                         E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
                                         E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
-                                        ! Additional terms in 2D and 3D
-                                        if ((i == 2) .or. (i == 4) .or. (i == 5)) then
+                                        ! Double for shear stresses
+                                        if (any(strxb - 1 + i == shear_indices)) then
                                             E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
                                             E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
                                         end if
@@ -1060,17 +1060,26 @@
                                         flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
                                     end do
                                     ! Recalculating the radial momentum geometric source flux
-                                    flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1)) = &
-                                        (s_M*(rho_R*vel_R(dir_idx(1)) &
-                                              *vel_R(dir_idx(1))) &
-                                         - s_P*(rho_L*vel_L(dir_idx(1)) &
-                                                *vel_L(dir_idx(1))) &
-                                         + s_M*s_P*(rho_L*vel_L(dir_idx(1)) &
-                                                    - rho_R*vel_R(dir_idx(1)))) &
-                                        /(s_M - s_P)
+                                    flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) = &
+                                        flux_rs${XYZ}$_vf(j, k, l, contxe + 2) &
+                                        - (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
                                     ! Geometrical source of the void fraction(s) is zero
                                     !$acc loop seq
                                     do i = advxb, advxe
+                                        ! flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
+                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
+                                    end do
+                                end if
+
+                                if (cyl_coord .and. hypoelasticity) then
+                                    ! += tau_sigmasigma using HLL
+                                    flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) = &
+                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) + &
+                                        (s_M*tau_e_R(4) - s_P*tau_e_L(4)) &
+                                        /(s_M - s_P)
+
+                                    !$acc loop seq
+                                    do i = strxb, strxe
                                         flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
                                     end do
                                 end if