Turbulence module

The Turbulence Module provides a suite of tools for analyzing and decomposing turbulent flow fields. It includes utilities for computing energy spectra, Helmholtz decomposition, and handling velocity field data, along with specialized methods for processing turbulence-related datasets.

HelmholtzSolver

class HelmholtzSolver(varargin)

Bases: handle

The HelmholtzSolver() class computes the Helmholtz-Hodge decomposition of a 3D velocity field into its solenoidal and dilatational parts using fast Fourier transform [1].

The decomposition is performed with the spectral method, which is only suitable for relatively smooth fields, i.e., with little power on small scales. The code assumes that the grid is uniform with dx = dy = dz.

This code is based on Ref. [2] and has been rewritten in MATLAB with some modifications.

Notes

For even NX, NY, and NZ, decomposed fields can be complex, with the imaginary part coming from the real part of the kmode at Nyquist frequency. In principle, the Nyquist frequency kmode should be dropped when doing the first derivatives to maintain symmetry. See footnote on page 4 of [2]. However, when the field is smooth enough, the imaginary part caused by the Nyquist frequency kmode should be negligible.

References

[1] Johnson, S. G. (2011). Notes on FFT-based differentiation.

MIT Applied Mathematics, Tech. Rep. Available: http://math.mit.edu/~stevenj/fft-deriv.pdf

[2] Xun Shi, Helmholtz-Hodge decomposition using fft (Python),

Available: https://github.com/shixun22/helmholtz

FLAG_CHECKS = 'true'

Flag to perform checks

FLAG_REPORT = 'false'

Flag to print predefined plots

FLAG_TIME = 'true'

Flag to print elapsed time

FLAG_WEIGHTED = 'true'

Flag to compute weighted velocity field as u * sqrt(rho)

HelmholtzSolver(varargin)

Constructor

check(velocity, solenoidal, dilatational)

Check the results of the Helmholtz decomposition

Parameters:

• obj (HelmholtzSolver) – HelmholtzSolver object

• velocity (VelocityField) – VelocityField instance with fields (u, v, w) containing the velocity components (fluctuations)

• solenoidal (VelocityField) – VelocityField instance with fields (u, v, w) containing the solenoidal velocity components (fluctuations)

• dilatational (VelocityField) – VelocityField instance with fields (u, v, w) containing the dilatational velocity components (fluctuations)

Returns:

STOP (float) – Relative error doing the decomposition status (bool): Flag indicating if the checks passed

decomposition(velocity)

Compute Helmholtz decomposition of the velocity field (fluctuations)

Parameters:

velocity (VelocityField) – VelocityField instance with fields (u, v, w) containing the velocity components (fluctuations)

Returns:

solenoidal (VelocityField) – VelocityField instance with fields (u, v, w) containing the solenoidal velocity components (fluctuations) dilatational (VelocityField): VelocityField instance with fields (u, v, w) containing the dilatational velocity components (fluctuations)

static getChi(solenoidal, rho, pressure, sound_mean)

Compute the correlation between the entropic density and the soleonidal part of the velocity field

Parameters:

• solenoidal (VelocityField) – VelocityField instance with fields (u, v, w) containing the solenoidal velocity components (fluctuations)

• rho (float) – Density field

• pressure (float) – Pressure field

• sound_mean (float) – Mean speed of sound

Returns:

chi (float) – Correlation between the entropic density and the soleonidal part of the velocity field chiVariance (float): Variance of the correlation

Note: The components of solenoidal refer to the fluctuations of the density weighted velocity field, i.e.,

u_weighted = u * sqrt(rho), v_weighted = v * sqrt(rho), w_weighted = w * sqrt(rho).

Example

[chi, chiVariance] = getChi(obj, solenoidal, rho, pressure, sound_mean)

static getWaveNumbers(sz)

Compute wave number grids for FFT

Parameters:

sz (float) – Size of the 3D array

Returns:

KX (float) – 3D array with the wave number in the x-direction KY (float): 3D array with the wave number in the y-direction KZ (float): 3D array with the wave number in the z-direction

Example

[KX, KY, KZ] = getWaveNumbers(sz)

plotConfig = None

PlotConfig object

printTime()

Print execution time

Parameters:

obj (EquilibriumSolver) – Object of the class EquilibriumSolver

solve(velocity, varargin)

Compute Helmholtz decomposition of the velocity field

Parameters:

• obj (HelmholtzSolver) – HelmholtzSolver object

• velocity (VelocityField) – Velocity field as a VelocityField object, struct, or 4D matrix

Optional Args:

• density (float): Density field

Returns:

solenoidal (VelocityField) – VelocityField object with fields (u, v, w) containing the solenoidal velocity components (fluctuations) dilatational (VelocityField): VelocityField with fields (u, v, w) containing the dilatational velocity components (fluctuations) velocity (VelocityField): VelocityField with fields (u, v, w) containing the velocity components (fluctuations) STOP (float): Relative error doing the decomposition

Shortcuts:

• [solenoidal, dilatational, velocity, STOP] = solve(obj, velocity, rho)

Examples

• [solenoidal, dilatational, velocity, STOP] = solve(obj, velocity)

• [solenoidal, dilatational, velocity, STOP] = solve(obj, velocity, rho)

• [solenoidal, dilatational, velocity, STOP] = solve(obj, velocity, ‘density’, rho)

time = None

Elapsed time

tol0 = '1e-3'

Tolerance for checks

TurbulenceSpectra

class TurbulenceSpectra(varargin)

Bases: handle

The TurbulenceSpectra() class provides methods for computing turbulence spectra. Supports spherical and cross-plane averaging.

FLAG_REPORT = 'false'

Flag to print predefined plots

FLAG_TIME = 'true'

Flag to print elapsed time

TurbulenceSpectra(varargin)

Constructor for TurbulenceSpectra

Optional Args:

averaging (char): Type of averaging (‘spherical’, ‘spherical2D’, ‘spherical2dPerp’, ‘crossplane’, ‘crossplane2D’) axis (char): Axis for cross-plane averaging (‘x’, ‘y’, ‘z’)

Example

analyzer = TurbulenceSpectra(‘averaging’, ‘crossplane’, ‘axis’, ‘z’);

averaging = "'spherical'"

Type of averaging (‘spherical’, ‘spherical2D’, ‘spherical2dPerp’, ‘crossplane’, ‘crossplane2D’)

axis = "'x'"

Axis for cross-plane averaging (‘x’, ‘y’, ‘z’)

fps = '60'

Frame rate for movie (temporal will be added into PlotConfig)

static getCrossplaneAveragedSpectra(EK)

Compute the cross-plane averaged energy spectra of a 3D field. The homogeneous plane are the first two dimensions of the 3D field.

Parameters:

EK (float) – 3D energy spectra (combined from all velocity components)

Returns:

EK_avg (float) – Radially averaged energy spectra for the slices k (float): Wavenumber along the homogeneous direction

Example

[EK_avg, k] = getCrossplaneAveragedSpectra(EK);

static getCrossplaneAveragedSpectra2D(EK)

Compute the cross-plane averaged 2D energy spectra of a 3D field. The homogeneous plane corresponds to the first two dimensions.

Parameters:

EK (float) – 3D energy spectra (combined from all velocity components)

Returns:

EK_avg2D (float) – 2D radially averaged energy spectra k (float): Wavenumber along homogeneous directions

Example

[EK_avg2D, k] = getCrossplaneAveragedSpectra2D(EK);

static getCrossplaneFFT(fluctuation, axisType)

Compute 2D FFT on homogeneous planes based on axisType

Parameters:

• fluctuation (float) – 3D fluctuation field

• axisType (char) – Axis for cross-plane averaging (‘x’, ‘y’, ‘z’)

Returns:

EK (float) –

3D array representing the 2D energy spectra computed on the

homogeneous planes, sliced along the inhomogeneous axis. The first two dimensions correspond to the homogeneous plane, while the third dimension represents the slices along the inhomogeneous axis.

Examples

• EK = getCrossplaneFFT(fluctuation, ‘x’);

• EK = getCrossplaneFFT(fluctuation, ‘y’);

• EK = getCrossplaneFFT(fluctuation, ‘z’);

static getCrossplaneFFTVelocity(velocity, axisType)

Compute 2D FFT on homogeneous planes based on axisType

Parameters:

• u (float) – 3D velocity field in the x-direction

• v (float) – 3D velocity field in the y-direction

• w (float) – 3D velocity field in the z-direction

• axisType (char) – Axis for cross-plane averaging (‘x’, ‘y’, ‘z’)

Returns:

EK (float) –

3D array representing the 2D energy spectra computed on the

homogeneous planes, sliced along the inhomogeneous axis. The first two dimensions correspond to the homogeneous plane, while the third dimension represents the slices along the inhomogeneous axis.

Examples

• EK = getCrossplaneFFTVelocity(velocity, ‘x’);

• EK = getCrossplaneFFTVelocity(velocity, ‘y’);

• EK = getCrossplaneFFTVelocity(velocity, ‘z’);

getEnergySpectra(fluctuation)

Compute the energy spectra of a 3D fluctuation field

Parameters:

• obj (TurbulenceSpectra) – TurbulenceSpectra object

• fluctuation (float) – 3D fluctuation field

Returns:

EK_avg (float) – Averaged energy spectra k (float): Wavenumber along the homogeneous direction (‘spherical’, ‘spherical2D’, ‘crossplane’, ‘crossplane2D’) kParalell (float): Wavenumber along the parallel direction (‘spherical2dPerp’) kPerp (float): Wavenumber along the perpendicular direction (‘spherical2dPerp’)

Examples

• [EK_avg, k] = getEnergySpectra(TurbulenceSpectra(), fluctuation);

• [EK_avg, kParallel, kPerp] = getEnergySpectra(TurbulenceSpectra(‘averaging’, ‘spherical2dPerp’), fluctuation);

getEnergySpectraVelocity(velocity)

Compute the energy spectrum of a 3D velocity field

Parameters:

• obj (TurbulenceSpectra) – TurbulenceSpectra object

• velocity (VelocityField) – Velocity field as a VelocityField object, struct, or 4D matrix

Returns:

EK_avg (float) – Averaged energy spectra k (float): Wavenumber along the homogeneous direction (‘spherical’, ‘spherical2D’, ‘crossplane’, ‘crossplane2D’) kParalell (float): Wavenumber along the parallel direction (‘spherical2dPerp’) kPerp (float): Wavenumber along the perpendicular direction (‘spherical2dPerp’)

Examples

• [EK_avg, k] = getEnergySpectraVelocity(TurbulenceSpectra(), velocity);

• [EK_avg, kParallel, kPerp] = getEnergySpectraVelocity(TurbulenceSpectra(‘averaging’, ‘spherical2dPerp’), velocity);

static getIntegralLengthScale(EK, k)

Compute the integral length scale from the energy spectra

Parameters:

• EK (float) – Energy spectra

• k (float) – Wavenumber along the homogeneous direction

Returns:

L (float) – Integral length scale

Example

L = getIntegralLengthScale(EK, k);

static getJointPDF(fluctuationX, fluctuationY)

Compute the joint probability density function (PDF) of two 3D fluctuation fields

Parameters:

• fluctuationX (float) – 3D fluctuation field

• fluctuationY (float) – 3D fluctuation field

Returns:

pdf (float) – Joint probability density function edgesX (float): Bin edges for fluctuationX edgesY (float): Bin edges for fluctuationY

Example

pdf = getJointPDF(TurbulenceSpectra(), fluctuationX, fluctuationY);

static getPDF(fluctuation)

Compute the probability density function (PDF) of a 3D fluctuation field

Parameters:

fluctuation (float) – 3D fluctuation field

Returns:

pdf (float) – Probability density function edges (float): Bin edges

Example

pdf = getPDF(TurbulenceSpectra(), fluctuation);

static getSphericallyAveragedSpectra(EK)

Compute the spherically averaged energy spectra of a 3D field

Parameters:

EK (float) – 3D energy spectra

Returns:

EK_avg (float) – Spherically averaged energy spectra

Example

EK_avg = getSphericallyAveragedSpectra(EK);

static getSphericallyAveragedSpectra2D(EK)

Compute the spherically averaged 2D energy spectra of a 3D field

Parameters:

EK (float) – 3D energy spectra

Returns:

EK_avg2D (float) – Spherically averaged 2D energy spectra k (float): Wavenumber along the homogeneous direction

Example

[EK_avg2D, k] = getSphericallyAveragedSpectra2D(EK);

static getSphericallyAveragedSpectra2DPerp(EK, axisType)

Compute the 2D spherical average of a 3D energy spectrum as a function of kParalell and kPerp, where kPerp = sqrt(k2^2 + k3^2) is calculated from the two dimensions orthogonal to the specified axis.

Parameters:

• EK (float) – 3D energy spectra

• axisType (char) – Axis for cross-plane averaging (‘x’, ‘y’, ‘z’)

Returns:

EK_avg2D (float) – Spherically averaged 2D energy spectra E(kParalell, kPerp) kParalell (float): Wavenumber along the parallel direction kPerp (float): Wavenumber along the perpendicular direction

Example

[EK_avg2D, kParalell, kPerp] = getSphericallyAveragedSpectra2DPerp(EK, ‘x’);

movie(k, EK, varargin)

Create a movie by plotting slices of the spectra

Parameters:

• obj (TurbulenceSpectra) – TurbulenceSpectra object

• k (float) – Wavenumber along the homogeneous direction

• EK (float) – Energy spectra (2D array)

Example

obj.movie(k, EK);

plot(k, EK, varargin)

Plot results

Parameters:

• obj (TurbluenceSpectra) – TurbulenceSpectra object

• k (float) – Wavenumber along the homogeneous direction

• EK (float) – Energy spectra

Optional Args:

• EK_i (float): Array of energy spectra

Returns:

ax (Axes) – Axes object

Examples

• ax = plot(TurbulenceSpectra(), k, EK);

• ax = plot(TurbulenceSpectra(), k, EK1, EK2);

plotConfig = None

PlotConfig object

plotJointPDF(pdf, edgesX, edgesY)

Plot the joint probability density function (PDF)

Parameters:

• obj (TurbulenceSpectra) – TurbulenceSpectra object

• pdf (float) – Joint probability density function

• edgesX (float) – Bin edges for fluctuationX

• edgesY (float) – Bin edges for fluctuationY

Returns:

ax (Axes) – Axes object

Example

ax = plotJointPDF(obj, pdf, edgesX, edgesY);

plotPDF(pdf, edges)

Plot the probability density function (PDF)

Parameters:

• obj (TurbulenceSpectra) – TurbulenceSpectra object

• pdf (float) – Probability density function

• edges (float) – Bin edges

Returns:

ax (Axes) – Axes object

Example

ax = plotPDF(obj, pdf, edges);

printTime()

Print execution time

Parameters:

obj (TurbulenceSpectra) – TurbulenceSpectra object

time = None

Elapsed time

VelocityField

class VelocityField(u, v, w, varargin)

Bases: handle, matlab.mixin.Copyable

This VelocityField() class provides methods for computing properties of a velocity field.

VelocityField(u, v, w, varargin)

Constructor for VelocityField

Parameters:

• u (float) – x-component of velocity

• v (float) – y-component of velocity

• w (float) – z-component of velocity

Returns:

obj (VelocityField) – VelocityField object

getCrossplaneFluctuations(axisType, varargin)

Compute fluctuating velocity components

Parameters:

• obj (VelocityField) – VelocityField instance with fields (u, v, w) containing the velocity components

• axis (char) – Axis for cross-plane averaging (‘x’, ‘y’, or ‘z’)

Optional Name-Pair Args:

• density (float): Density field

• weighted (bool): Flag to indicate if the velocity should be weighted by sqrt(density)

Returns:

velocity (VelocityField) – Struct with fields (u, v, w) containing the fluctuating velocity components (fluctuations)

Notes

• The velocity field is assumed to be periodic in the direction perpendicular to the axis.

• By default, the velocity is weighted by density, if provided.

Examples

• velocity = getCrossplaneFluctuations(obj, ‘x’);

• velocity = getCrossplaneFluctuations(obj, ‘x’, ‘density’, rho);

• velocity = getCrossplaneFluctuations(obj, ‘x’, ‘density’, rho, ‘weighted’, true);

getDissipation(temperature, dynamicViscosity, varargin)

Compute dissipation energy

Parameters:

• obj (VelocityField) – VelocityField object

• temperature (float) – 3D array with the temperature field

• dynamicViscosity (float) – 3D array with the dynamic viscosity field

Optional Args:

Returns:

dissipation (float)

getFluctuations(varargin)

Compute fluctuating velocity components

Parameters:

obj (VelocityField) – VelocityField instance with fields (u, v, w) containing the velocity components

Optional Name-Pair Args:

• density (float): Density field

• weighted (bool): Flag to indicate if the velocity should be weighted by sqrt(density)

Returns:

velocity (VelocityField) – Struct with fields (u, v, w) containing the fluctuating velocity components (fluctuations)

Notes

• The velocity field is assumed to be periodic in all three dimensions.

• By default, the velocity is weighted by density, if provided.

Shortcuts:

• velocity = getFluctuations(obj, rho);

Examples

• velocity = getFluctuations(obj);

• velocity = getFluctuations(obj, rho);

• velocity = getFluctuations(obj, ‘density’, rho);

• velocity = getFluctuations(obj, ‘density’, rho, ‘weighted’, true);

getTurbulentKineticEnergy(varargin)

Compute Turbulent Kinetic Energy (TKE)

Parameters:

obj (VelocityField) – VelocityField object

Optional Args:

• rho (float): 3D array with the density field

Returns:

K (float) – 3D array with the TKE

Example

K = getTurbulentKineticEnergy(u, v, w, rho);

getTurbulentMachNumberAverage(soundVelocityAverage)

Compute mean turbulent Mach number

Parameters:

• obj (VelocityField) – VelocityField object

• soundVelocityAverage (float) – Mean sound velocity

Returns:

Mt (float) – Mean turbulent Mach number

getVelocityRMS()

Compute root mean square (rms) of the velocity field

Parameters:

obj (VelocityField) – VelocityField object

Returns:

velocityRMS (float) – Root mean square of the velocity field

globalMagnitude()

Compute the global magnitude (Frobenius norm) of the velocity field

Parameters:

obj (VelocityField) – VelocityField object

Returns:

globalMagnitude (float) – Scalar magnitude of the velocity field

Example

globalMagnitude = globalMagnitude(obj);

normalize()

Normalize the velocity field by its magnitude

Returns:

obj – Normalized velocity field

plus(obj2)

Overload the plus operator to add the velocity field from two VelocityField objects

Parameters:

• obj (VelocityField) – VelocityField object

• obj2 (VelocityField) – VelocityField object

Returns:

obj (VelocityField) – VelocityField object

pointwiseMagnitude()

Compute the pointwise magnitude of the velocity field

Parameters:

obj (VelocityField) – VelocityField object

Returns:

magnitudeField (float) – 3D matrix of velocity magnitudes

Example

magnitudeField = pointwiseMagnitude(obj);

scale(factor)

Scale the velocity field by a factor

Parameters:

factor – Scaling factor

Returns:

obj – Scaled velocity field

size = None

Size of the velocity field

u = None

x-component of the velocity field

v = None

y-component of the velocity field

w = None

z-component of the velocity field

x = None

x-component field

y = None

y-component field

z = None

z-component field