Update initialization for turbulent mixing layer and add validation results

.typos.toml

@@ -17,6 +17,7 @@ gam = "gam"
 Gam = "Gam"
 strang = "strang"
 Strang = "Strang"
+TKE = "TKE"
 
 [files]
 extend-exclude = ["docs/documentation/references*", "tests/", "toolchain/cce_simulation_workgroup_256.sh"]

docs/documentation/case.md

@@ -825,7 +825,6 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 | `mixlayer_vel_profile` | Logical | Set the mean streamwise velocity to hyperbolic tangent profile |
 | `mixlayer_vel_coef`    | Real    | Coefficient for the hyperbolic tangent profile of a mixing layer |
 | `mixlayer_perturb`     | Logical | Perturb the initial velocity field by instability waves |
-| `mixlayer_domain`      | Real    | Domain size of a mixing layer for the linear stability analysis |
 
 The table lists velocity field parameters.
 The parameters are optionally used to define initial velocity profiles and perturbations.
@@ -840,15 +839,13 @@ The parameters are optionally used to define initial velocity profiles and pertu
 
 - `perturb_sph_fluid` specifies the fluid component whose partial density is to be perturbed.
 
-- `mixlayer_vel_profile` activates setting the mean streamwise velocity to a hyperbolic tangent profile. This option works only for 2D and 3D cases.
+- `mixlayer_vel_profile` activates setting the mean streamwise velocity to a hyperbolic tangent profile. This option works only for `n > 0`.
 
 - `mixlayer_vel_coef` is a parameter for the hyperbolic tangent profile of a mixing layer when `mixlayer_vel_profile = 'T'`. The mean streamwise velocity profile is given as:
 
-$$ u = patch\_icpp(1)\%vel(1) * tanh(y\_cc * mixlayer\_vel\_profile) $$
+$$ u = \mbox{patch\_icpp(1)\%vel(1)} * \tanh( y_{cc} * \mbox{mixlayer\_vel\_coef}) $$
 
-- `mixlayer_perturb` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for `n > 0`, `bc_y%[beg,end] = -6`, `num_fluids = 1`, `model_eqns = 2` and `mixlayer_vel_profile = 'T'`.
-
-- `mixlayer_domain` defines the domain size to compute spatial eigenvalues of the linear instability analysis when `mixlayer_perturb = 'T'`. For example, the spatial eigenvalue in `x` direction in 2D problem will be $2 \pi \alpha / (mixlayer\_domain*patch\_icpp(1)\%length\_y)$ for $\alpha = 1$, $2$ and $4$.
+- `mixlayer_perturb` activates the velocity perturbation for a temporal mixing layer with hyperbolic tangent mean streamwise velocity profile, using an inverter version of the spectrum-based synthetic turbulence generation method proposed by Guo et al. (2023, JFM). This option only works for `p > 0` and `mixlayer_vel_profile = 'T'`.
 
 ### 11. Phase Change Model
 | Parameter              | Type    | Description                                    |

docs/documentation/references.md

@@ -20,6 +20,8 @@
 
 - <a id="Gottlieb98">Gottlieb, S. and Shu, C.-W. (1998). Total variation diminishing runge-kutta schemes. Mathematics of computation of the American Mathematical Society, 67(221):73–85.</a>
 
+- <a id="Guo23">Guo, H., Jiang, P., Ye, L., & Zhu, Y. (2023). An efficient and low-divergence method for generating inhomogeneous and anisotropic turbulence with arbitrary spectra. Journal of Fluid Mechanics, 970, A2.</a>
+
 - <a id="Henrick05">Henrick, A. K., Aslam, T. D., and Powers, J. M. (2005). Mapped weighted essentially nonoscillatory schemes: achieving optimal order near critical points. Journal of Computational Physics, 207(2):542–567.</a>
 
 - <a id="Borges08">Borges, R., Carmona, M., Costa, B., and Don, W. S. (2008). An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of computational physics, 227(6):3191–3211.</a>
@@ -34,6 +36,8 @@
 
 - <a id="Meng16">Meng, J. C. C. (2016). Numerical simulations of droplet aerobreakup. PhD thesis, California Institute of Technology.</a>
 
+- <a id="Michalke64"> Michalke, A. (1964). On the inviscid instability of the hyperbolictangent velocity profile. Journal of Fluid Mechanics, 19(4), 543-556.</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>

examples/2D_mixing_artificial_Ma/case.py

@@ -108,8 +108,6 @@
             # Mixing layer
             "mixlayer_vel_profile": "T",
             "mixlayer_vel_coef": 1.0,
-            "mixlayer_domain": 1.0,
-            "mixlayer_perturb": "T",
             # Artificial Mach number
             "pi_fac": pi_fac,
             # Fluids Physical Parameters

examples/3D_turb_mixing/case.py

@@ -4,39 +4,38 @@
 
 # SURROUNDING FLOW
 # Nondimensional parameters
-Re0 = 50.0  # Reynolds number
+Re0 = 320.0  # Reynolds number
 M0 = 0.2  # Mach number
 
 # Fluid properties
 gamma = 1.4
-pi_inf = 0.0
 
 # Free stream velocity & pressure
 u0 = 1.0
 pres0 = 1.0 / (gamma * M0**2)
 
 # Domain size
-Lx = 59.0
-Ly = 59.0
-Lz = 59.0
+Lx = 160.0
+Ly = 160.0
+Lz = 80.0
 
 # Number of grid cells
-Nx = 191
-Ny = 191
-Nz = 191
+Nx = 1023
+Ny = 1023
+Nz = 511
 
 # Grid spacing
 dx = Lx / float(Nx)
 dy = Ly / float(Ny)
 dz = Lz / float(Nz)
 
 # Time advancement
-cfl = 0.1
-T = 100.0
-dt = cfl * dx / u0
+cfl = 0.5
+T = 150.0
+dt = cfl * dx / (u0 / M0 + 1)
 Ntfinal = int(T / dt)
 Ntstart = 0
-Nfiles = 100
+Nfiles = 30
 t_save = int(math.ceil((Ntfinal - 0) / float(Nfiles)))
 Nt = t_save * Nfiles
 t_step_start = Ntstart
@@ -53,6 +52,11 @@
             "x_domain%end": Lx,
             "y_domain%beg": -Ly / 2.0,
             "y_domain%end": Ly / 2.0,
+            "stretch_y": "T",
+            "a_y": 2,
+            "y_a": -0.3 * Ly,
+            "y_b": 0.3 * Ly,
+            "loops_y": 2,
             "z_domain%beg": 0.0,
             "z_domain%end": Lz,
             "m": Nx,
@@ -68,10 +72,9 @@
             "num_fluids": 1,
             "time_stepper": 3,
             "weno_order": 5,
-            "weno_eps": 1.0e-16,
-            "weno_Re_flux": "T",
-            "weno_avg": "T",
-            "mapped_weno": "T",
+            "weno_eps": 1.0e-40,
+            "weno_Re_flux": "F",
+            "wenoz": "T",
             "riemann_solver": 2,
             "wave_speeds": 1,
             "avg_state": 2,
@@ -85,7 +88,7 @@
             # Formatted Database Files Structure Parameters
             "format": 1,
             "precision": 2,
-            "cons_vars_wrt": "T",
+            "cons_vars_wrt": "F",
             "prim_vars_wrt": "T",
             "parallel_io": "T",
             "fd_order": 1,
@@ -99,7 +102,7 @@
             "patch_icpp(1)%y_centroid": 0.0,
             "patch_icpp(1)%z_centroid": Lz / 2.0,
             "patch_icpp(1)%length_x": Lx,
-            "patch_icpp(1)%length_y": Ly,
+            "patch_icpp(1)%length_y": 1e6,
             "patch_icpp(1)%length_z": Lz,
             "patch_icpp(1)%alpha_rho(1)": 1.0,
             "patch_icpp(1)%alpha(1)": 1.0,
@@ -109,8 +112,6 @@
             "patch_icpp(1)%pres": pres0,
             # Mixing layer
             "mixlayer_vel_profile": "T",
-            "mixlayer_vel_coef": 1.0,
-            "mixlayer_domain": 1.0,
             "mixlayer_perturb": "T",
             # Fluids Physical Parameters
             # Surrounding liquid

examples/3D_turb_mixing/turbulence_stat/README.md

@@ -0,0 +1,23 @@
+# Temporally evolving turbulent mixing layer (3D)
+
+Reference:
+> Bell, J. H., & Mehta, R. D. (1990). Development of a two-stream mixing layer from tripped and untripped boundary layers. AIAA Journal, 28(12), 2034-2042.
+> Rogers, M. M., & Moser, R. D. (1994). Direct simulation of a self‐similar turbulent mixing layer. Physics of Fluids, 6(2), 903-923.
+> Pantano, C., & Sarkar, S. (2002). A study of compressibility effects in the high-speed turbulent shear layer using direct simulation. Journal of Fluid Mechanics, 451, 329-371.
+> Vaghefi, S. N. S. (2014). Simulation and modeling of compressible turbulent mixing layer. State University of New York at Buffalo.
+> Wang, X., Wang, J., & Chen, S. (2022). Compressibility effects on statistics and coherent structures of compressible turbulent mixing layers. Journal of Fluid Mechanics, 947, A38.
+
+## Description
+This directory (`turbulence_stat`) contains sub-directories and Matlab scripts for computing turbulence statistics of 3D temporally evolving turbulent mixing layer as described below:
+
+Files:
+- `set_user_inputs.m`: User input parameters are defined in this file.
+- `run_turbulence.m`: Turbulence statistics (growth rate, Reynolds stress, TKE budget) are computed and plots are generated in this file.
+- `average_tke_over_self_similar.m`: Averaging over self-similar period of TKE budget is performed in this file.
+- `submit_batch_job_*.sh`: Example scripts for batch job submission on DoD Carpenter (PBS system) and NCSA Delta (Slurm system)
+
+Directories:
+- `log`: Log files from batch job are stored in this directory.
+- `reference_data`: Data from reference papers are provided in this directory.
+- `variables`: User inputs from `set_user_inputs.m` are saved as Matlab data (`user_inputs.mat`) in this directory.
+- `results`: Outputs are stored in sub-directories in the directory.

examples/3D_turb_mixing/turbulence_stat/average_tke_over_self_similar.m

@@ -0,0 +1,118 @@
+close all; clear all;
+
+
+% Setup
+disp("Start average_tke_over_self_similar ..."); tic;
+set_user_inputs(); load variables/user_inputs.mat;
+
+ybeg = -5; yend = 5; ny = 101;
+y = linspace(ybeg,yend,ny);
+
+% Array
+T0_averaged = zeros(ny,1);
+P_averaged = zeros(ny,1);
+D_averaged = zeros(ny,1);
+
+% Compute averaged TKE budget
+for q = 1:Nfiles
+    load("results/tke_budget_data/tstep_"+string(timesteps(q))+".mat");
+
+    % Normalization
+    T0 = T0 / (8/mth);  % T / (Delta U^3 / mth)
+    P = P / (8/mth);    % P / (Delta U^3 / mth)
+    D = D / (8/mth);    % D / (Delta U^3 / mth)
+
+    % Interpolation
+    i_start = 1;
+    for j = 1:ny
+        for i = i_start:length(y_norm_mth) - 1
+            if (y_norm_mth(i) <= y(j) && y_norm_mth(i+1) > y(j))
+                T0_averaged(j) = T0_averaged(j) + ((T0(i+1) - T0(i))/(y_norm_mth(i+1) - y_norm_mth(i))*(y(j) - y_norm_mth(i)) + T0(i))/Nfiles;
+                P_averaged(j) = P_averaged(j) + ((P(i+1) - P(i))/(y_norm_mth(i+1) - y_norm_mth(i))*(y(j) - y_norm_mth(i)) + P(i))/Nfiles;
+                D_averaged(j) = D_averaged(j) + ((D(i+1) - D(i))/(y_norm_mth(i+1) - y_norm_mth(i))*(y(j) - y_norm_mth(i)) + D(i))/Nfiles;
+                i_start = i;
+                break;
+            end
+        end
+    end
+end 
+
+T_averaged = f_compute_derivative_1d(T0_averaged,y*mth); 
+
+% Plot
+plot_tke_budget(T_averaged, P_averaged, D_averaged, y, mth);
+
+disp("End of program"); toc;
+
+%% FUNCTIONS
+% Compute the wall-normal derivative of a discretized function, fun(y)
+function dfunds = f_compute_derivative_1d(fun,s)
+
+    dfunds = zeros(size(fun));    % initialize discrete derivative vector
+
+    % Compute one-sided derivative at the bottom boundary
+    dfunds(1) = (fun(2) - fun(1)) / (s(2) - s(1));
+
+    % Compute one-sided derivative at the top boundary
+    dfunds(end) = (fun(end) - fun(end-1)) / (s(end) - s(end-1));
+
+    % Compute two-sided derivatives for interior points
+    for i = 2:length(s)-1
+        dfunds(i) = (fun(i+1) - fun(i-1)) / (s(i+1) - s(i-1));
+    end
+end
+
+% Plot TKE budget
+function plot_tke_budget(T, P, D, y_norm_mth, mth)
+
+    load variables/user_inputs.mat;
+    load reference_data/reference.mat;
+
+    % Plot
+    f1 = figure("DefaultAxesFontSize",18);
+    set(f1,"Position",[200 200 1000 700]);
+
+    % Present
+    h1 = plot([-100 -100],[-100 -100],'-k','LineWidth',2); hold on; grid on;
+    plot(y_norm_mth,T,'-b','LineWidth',2); 
+    plot(y_norm_mth,P,'-g','LineWidth',2);
+    plot(y_norm_mth,D,'-r','LineWidth',2);
+    xlim([-5 5]); xticks([-5:1:5]);
+    ylim([-0.002 0.003]);
+    xlabel('$y/\delta_\theta$','interpreter','latex');
+    set(gca,'TickLabelInterpreter','latex');
+
+    % Pantano & Sarkar (2002)
+    h2 = plot([-100 -100],[-100 -100],'ko','LineWidth',2,'MarkerSize',8);
+    plot(p2002_tke_transport(:,1),p2002_tke_transport(:,2),'bo','LineWidth',2,'MarkerSize',8);
+    plot(p2002_tke_production(:,1),p2002_tke_production(:,2),'go','LineWidth',2,'MarkerSize',8);
+    plot(p2002_tke_dissipation(:,1),p2002_tke_dissipation(:,2),'ro','LineWidth',2,'MarkerSize',8);
+
+    % Rogers & Moser (1994)
+    h3 = plot([-100 -100],[-100 -100],'k^','LineWidth',2,'MarkerSize',8);
+    plot(r1994_tke_transport(:,1),r1994_tke_transport(:,2),'b^','LineWidth',2,'MarkerSize',8);
+    plot(r1994_tke_production(:,1),r1994_tke_production(:,2),'g^','LineWidth',2,'MarkerSize',8);
+    plot(r1994_tke_dissipation(:,1),r1994_tke_dissipation(:,2),'r^','LineWidth',2,'MarkerSize',8);
+
+    % Vaghefi (2014)
+    h4 = plot([-100 -100],[-100 -100],'k+','LineWidth',2,'MarkerSize',8);
+    plot(v2014_tke_transport(:,1),v2014_tke_transport(:,2),'b+','LineWidth',2,'MarkerSize',8);
+    plot(v2014_tke_production(:,1),v2014_tke_production(:,2),'g+','LineWidth',2,'MarkerSize',8);
+    plot(v2014_tke_dissipation(:,1),v2014_tke_dissipation(:,2),'r+','LineWidth',2,'MarkerSize',8);
+
+    % Wang et al. (2022)
+    h5 = plot([-100 -100],[-100 -100],'k*','LineWidth',2,'MarkerSize',8);
+    plot(w2022_tke_transport(:,1),w2022_tke_transport(:,2),'b*','LineWidth',2,'MarkerSize',8);
+    plot(w2022_tke_production(:,1),w2022_tke_production(:,2),'g*','LineWidth',2,'MarkerSize',8);
+    plot(w2022_tke_dissipation(:,1),w2022_tke_dissipation(:,2),'r*','LineWidth',2,'MarkerSize',8);
+
+    legend([h1,h2,h3,h4,h5], {"$\mbox{Present}$", ...
+            "$\mbox{Pantano \& Sarkar (2002)}$", ...
+            "$\mbox{Rogers \& Moser (1994)}$", ...
+            "$\mbox{Vaghefi (2014)}$", ...
+            "$\mbox{Wang et al. (2022)}$"}, ...
+            'interpreter','latex','location','northeast');
+
+    saveas(f1, "results/tke_budget/avg_self_similar","png");
+    close(f1);
+end