diff --git a/README.md b/README.md index d7d30fa9a5d1fa456cacd8d83656068e62c5c43e..37036258543400cc15d20a4985c1de7ed96e62cb 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ TODO: ## Features -- Generate solenoidal velocity fields using the method of Rogallo. +- Generate solenoidal velocity fields using Kraichnan's method, improved for mass conservation by Saad. - Compute incompressible pressure fluctuations using the method of Ristorcelli & Blaisdell. - Compute thermodynamic fields (density, pressure, temperature) from the generated velocity and pressure fields. - Plot the energy spectrum of the generated turbulence. @@ -41,8 +41,8 @@ This project is licensed under the MIT License. See the LICENSE file for details ## Acknowledgments -- Ristorcelli, J. R., & Blaisdell, G. A. (1997). Consistent initial conditions for the DNS of compressible turbulence. Journal of Computational Physics, 140(2), 427-466. -- Sagaut, P. (2018). Homogeneous Turbulence Dynamics. Springer. -- Tony Saad, Derek Cline, Rob Stoll, and James C. Sutherland. “Scalable Tools for Generating Synthetic Isotropic Turbulence with Arbitrary Spectraâ€. http://dx.doi.org/10.2514/1.J055230. (available online: Aug 25, 2016). -- Saad, T., & Sutherland, J. C. (2016). Comment on “Diffusion by a random velocity field†[Phys. Fluids 13, 22 (1970)]. Physics of Fluids (1994-Present), 28(11), 119101. https://doi.org/10.1063/1.4968528. -- Austin Richards, Tony Saad, and James C. Sutherland. “A Fast Turbulence Generator using Graphics Processing Unitsâ€, 2018 Fluid Dynamics Conference, AIAA AVIATION Forum, (AIAA 2018-3559). https://doi.org/10.2514/6.2018-3559. \ No newline at end of file +[1] Kraichnan R. H., “Diffusion by a Random Velocity Field,†Physics of Fluids, Vol. 13, No. 1, 1970, Paper 22. doi: https://doi.org/10.1063/1.1692799 + +[2] Tony Saad, Derek Cline, Rob Stoll, and James C. Sutherland. “Scalable Tools for Generating Synthetic Isotropic Turbulence with Arbitrary Spectraâ€. http://dx.doi.org/10.2514/1.J055230. (available online: Aug 25, 2016). + +[3] Ristorcelli, J. R., & Blaisdell, G. A. (1997). Consistent initial conditions for the DNS of compressible turbulence. Journal of Computational Physics, 140(2), 427-466. \ No newline at end of file diff --git a/pyTurbulence/solver.py b/pyTurbulence/solver.py index 08a99c5e914abadd6458f258f3ffc0242a4d6826..ee82f27f15b6f8c74884a18d5a2e6e801a4f80f2 100644 --- a/pyTurbulence/solver.py +++ b/pyTurbulence/solver.py @@ -183,7 +183,7 @@ def _validate_params(params) -> dict: # Define allowed values for specific parameters allowed_values = { "Spectrum": {"PassotPouquet", "Gaussian", "Exponential"}, # Allowed spectrum types - "case": {0, 1, 2}, # Allowed integer cases + "case": {1, 2}, # Allowed integer cases } # Define numeric ranges for specific parameters diff --git a/pyTurbulence/spectrum.py b/pyTurbulence/spectrum.py index a4ea495c49d48cc80ee3de468ab86a089f365822..a6d2ac36cf8cbe09a296f46f5cc8744797276209 100644 --- a/pyTurbulence/spectrum.py +++ b/pyTurbulence/spectrum.py @@ -1,9 +1,7 @@ import numpy as np import matplotlib.pyplot as plt -from numba import njit from scipy import stats -@njit def energy_spectrum(spectrum: str, k: np.ndarray, urms: float, k0: float) -> np.ndarray: """ Compute the energy spectrum E(k) for a given vector of wavenumber k. @@ -32,7 +30,7 @@ def energy_spectrum(spectrum: str, k: np.ndarray, urms: float, k0: float) -> np. E_k = urms*k**4*np.exp(-2.*k**2/k0**2) return E_k -def compute_tke_spectrum(u,v,w,lx,ly,lz): +def compute_tke_spectrum(u,v,w,lx,ly,lz,smooth=True): """ Compute the turbulent kinetic energy (TKE) spectrum from the velocity fields. @@ -87,10 +85,17 @@ def compute_tke_spectrum(u,v,w,lx,ly,lz): # --- Spectrum Abins, _, _ = stats.binned_statistic(knrm, fourier_amplitudes, statistic = "sum", bins = kbins) + tke_spectrum = Abins/(kbins[1:] - kbins[:-1]) knyquist = knrm * min(n0, n1, n2) / 2 + + if smooth: + window = np.ones(int(5)) / float(5) + tkespecsmooth = np.convolve(tke_spectrum, window, 'same') + tkespecsmooth[0:4] = tke_spectrum[0:4] # get the first 4 values from the original data + tke_spectrum = tkespecsmooth - return knyquist, kvals, Abins/(kbins[1:] - kbins[:-1]) + return knyquist, kvals, tke_spectrum def plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=None, save_path=None): plt.figure() diff --git a/pyTurbulence/syntheticTurbulence.py b/pyTurbulence/syntheticTurbulence.py index 89889236d9ed7049dd38e811eb0d1a7a02982a3a..64f0db35a5dba00eec9bb005a90c6fcfb42148cd 100644 --- a/pyTurbulence/syntheticTurbulence.py +++ b/pyTurbulence/syntheticTurbulence.py @@ -21,14 +21,14 @@ def _compute_velocity_field(u_, v_, w_, nx, ny, nz, kx, ky, kz, xc, yc, zc, um, # MIT License # Original function from turboGenPy by Tony Saad -# Modified by Amaury Bilocq in 2024 +# Modified by Amaury Bilocq in 2025 def compute_solenoidal_velocities(spectrum: str, nbModes: int, urms: float, k0: float, domain_size: Tuple[float, float, float], domain_resolution: Tuple[int, int, int] ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float]: """ - Generate a synthetic turbulence velocity field using the Tony Saad methodology. + Generate a synthetic turbulence velocity field using the Tony Saad's methodology. Parameters ---------- @@ -64,12 +64,12 @@ def compute_solenoidal_velocities(spectrum: str, nbModes: int, urms: float, dz = Lz/nz # Proper wavenumber scaling - wmin = 2.*np.pi/Lx - wmax = np.pi/dx + # unique_kk, wave_instances, wave_width = wavespace(nbModes) + wmin = 2.*np.pi/np.min([Lx,Ly,Lz]) + wmax = np.pi/np.min([dx,dy,dz]) dk = (wmax - wmin)/nbModes wn = wmin + np.arange(0,nbModes)*dk - dkn = np.ones(nbModes)*dk # Generate random angles - phases @@ -109,7 +109,7 @@ def compute_solenoidal_velocities(spectrum: str, nbModes: int, urms: float, Ek = energy_spectrum(spectrum, km, urms, k0) Ek = Ek.clip(0.0) - # generate turbulence at cell centers + # generate turbulence um = np.sqrt(Ek * dkn) u_ = np.zeros((nx, ny, nz)) v_ = np.zeros((nx, ny, nz)) @@ -149,6 +149,7 @@ def compute_solenoidal_pressure(u: np.ndarray,v: np.ndarray,w: np.ndarray, pressure : np.ndarray Solenoidal pressure field. """ + nx = domain_resolution[0] ny = domain_resolution[1] nz = domain_resolution[2]