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%
21% Last update Nov 23 2024
22% -------------------------------------------------------------------------
23
24% Import required packages
25import combustiontoolbox.turbulence.*
26
27% Load field variables
28filename = 'velocityField.h5';
29filePath = which(filename);
30rho = double(h5read(filePath, ['/', 'density']));
31u = double(h5read(filePath, ['/', 'velocity/u']));
32v = double(h5read(filePath, ['/', 'velocity/v']));
33w = double(h5read(filePath, ['/', 'velocity/w']));
34
35% Convert velocity to VelocityField object
36velocity = VelocityField(u, v, w);
37
38% Initialize HelmholtzSolver
39solver = HelmholtzSolver();
40
41% Perform Helmholtz decomposition
42[solenoidal, dilatational, STOP] = solver.solve(velocity, 'density', rho);
43
44% Compute turbulent kinetic energy (TKE)
45K_solenoidal = solenoidal.getTurbulentKineticEnergy(rho);
46K_dilatational = dilatational.getTurbulentKineticEnergy(rho);
47
48% Get dilatational over solenoidal TKE ratio
49eta = K_dilatational / K_solenoidal;
50
51% Analyze turbulence spectra using TurbulenceSpectra
52analyzer = TurbulenceSpectra();
53
54% Compute energy spectra for the original, solenoidal, and dilatational fields
55[EK1, k1] = analyzer.getEnergySpectra(velocity); % Original field
56[EK2, k2] = analyzer.getEnergySpectra(solenoidal); % Solenoidal component
57[EK3, k3] = analyzer.getEnergySpectra(dilatational); % Dilatational component
58
59% Plot results
60analyzer.plot(k1, EK1, EK2, EK3);