1% -------------------------------------------------------------------------
2% EXAMPLE: Helmholtz decomposition
3%
4% Perform Helmholtz decomposition of a velocity field and analyze the
5% turbulence spectra of the original, solenoidal, and dilatational fields.
6%
7% Note: The data used is a downsampled velocity field from a Direct Numerical
8% Simulation (DNS) of a Homogeneous Isotropic Turbulence (HIT) case carried
9% out with the HTR solver [1]. Take the data as an example to illustrate
10% the use of the HelmholtzSolver and TurbulenceSpectra classes.
11%
12% References:
13% [1] Di Renzo, M., Fu, L., & Urzay, J. (2020). HTR solver: An
14% open-source exascale-oriented task-based multi-GPU high-order
15% code for hypersonic aerothermodynamics. Computer Physics
16% Communications, 255, 107262.
17%
18%
19% @author: Alberto Cuadra Lara
20% Postdoctoral researcher - Group Fluid Mechanics
21% Universidad Carlos III de Madrid
22%
23% Last update Nov 23 2024
24% -------------------------------------------------------------------------
25
26% Import required packages
27import combustiontoolbox.turbulence.*
28
29% Load field variables
30filename = 'velocityField.h5';
31filePath = which(filename);
32rho = double(h5read(filePath, ['/', 'density']));
33u = double(h5read(filePath, ['/', 'velocity/u']));
34v = double(h5read(filePath, ['/', 'velocity/v']));
35w = double(h5read(filePath, ['/', 'velocity/w']));
36
37% Convert velocity to VelocityField object
38velocity = VelocityField(u, v, w);
39
40% Initialize HelmholtzSolver
41solver = HelmholtzSolver();
42
43% Perform Helmholtz decomposition
44[solenoidal, dilatational, STOP] = solver.solve(velocity, 'density', rho);
45
46% Compute turbulent kinetic energy (TKE)
47K_solenoidal = solenoidal.getTurbulentKineticEnergy(rho);
48K_dilatational = dilatational.getTurbulentKineticEnergy(rho);
49
50% Get dilatational over solenoidal TKE ratio
51eta = K_dilatational / K_solenoidal;
52
53% Analyze turbulence spectra using TurbulenceSpectra
54analyzer = TurbulenceSpectra();
55
56% Compute energy spectra for the original, solenoidal, and dilatational fields
57[EK1, k1] = analyzer.getEnergySpectra(velocity); % Original field
58[EK2, k2] = analyzer.getEnergySpectra(solenoidal); % Solenoidal component
59[EK3, k3] = analyzer.getEnergySpectra(dilatational); % Dilatational component
60
61% Plot results
62analyzer.plot(k1, EK1, EK2, EK3);