Generating an underwater shock wave with well-controlled intensity and curvature radius is important for relevant experimental and numerical studies but still challenging. The present study employs the stiffened-gas equation of state (SG-EOS) for liquid water and proposes a method of theoretically designing a shock tube to transform a planar underwater shock wave into a cylindrical one with the preset intensity and curvature. Two fundamental relations, including the Chester-Chisnell-Whitham relation and the two-dimensional characteristic relation of shock dynamics, are formulated using SG-EOS. The design procedure of the shock tube profile is thus derived. The present study conducts typical numerical tests to examine the effectiveness of the proposed method. To evaluate the errors between numerical shocks and their theoretical cylindrical shapes, the normalized ${\mathcal{L}}^2$ and ${\mathcal{L}}^{\infty }$ distances are defined and found to be controlled within $O\left({10}^{-3}\right)$, or even $O\left({10}^{-4}\right)$ in some cases. And the relative errors of the postshock pressure and velocity magnitude can be controlled within 2%, or even 0.7% in some cases. Pressure profiles are also presented and found to continuously approach those of an ideal planar shock wave as the curvature radius increases.

Data-driven subgrid-scale (SGS) modeling in the large-eddy simulations (LESs) suffers from an inconsistency between the a priori tests and the a posteriori tests, which make training accurate SGS models a difficult task. We study the difference in filtered high-fidelity data and LESs to identify the numerical deviation between the two cases, which is a combined impact of commutation error, numerical errors, and error coupling. The impact of the numerical deviation is examined through two SGS model formulations: the eddy-viscosity and the complex nonlinear models. By incorporating numerical deviations into model training, we enhance consistency, stabilize simulations, and improve predictions of energy spectra in a posteriori tests. Our findings highlight that data-driven methods introduce significant nonlinearity and equation coupling, exacerbating inconsistencies compared to non-data-driven approaches. Finally, while the impact of the numerical deviation can be generalized, achieving accurate model predictions necessitates a physically grounded model form and an optimal filter width.

Gauss's principle of least constraint transforms a dynamics problem into a pure minimization problem, where the total magnitude of the constraint force is the cost function, minimized at each instant. A candidate motion whose evolution minimizes the Gaussian cost at each instant is guaranteed to satisfy Newton's equation. In other words, Newton's equation is the first-order necessary condition for minimizing the Gaussian cost, subject to the given kinematic constraints. The principle of minimum pressure gradient (PMPG) is to incompressible fluid mechanics what Gauss's principle is to particle mechanics. The PMPG asserts that an incompressible flow evolves from one instant to another by minimizing the $L^2$-norm of the pressure gradient force. A candidate flow field whose evolution minimizes the pressure gradient cost at each instant is guaranteed to satisfy the Navier-Stokes equation. Consequently, the PMPG transforms the incompressible fluid mechanics problem into a pure minimization framework, allowing one to determine the evolution of the flow field by solely focusing on minimizing the cost—without directly invoking the Navier-Stokes equation. In this paper, we show that the resulting minimization problem is a convex quadratic programming (QP) problem—one of the most computationally tractable classes in nonlinear optimization, which has a rich literature with many efficient algorithms. Moreover, leveraging tools from analytical mechanics and the Moore-Penrose theory of generalized inverses, we derive an analytical solution for this QP problem. As a result, we present an explicit formula for the projected dynamics of the spatially discretized Navier-Stokes equation on the space of divergence-free fields. The resulting ordinary differential equation (ODE) is ready for direct time integration, eliminating the need for solving the Poisson equation in pressure at each time step. It is typically an explicit nonlinear ODE with constant coefficients. This compact form is expected to be highly valuable for both simulation and theoretical studies, including stability analysis and flow control design. We demonstrate the framework on the benchmark problem of unsteady flow in a lid-driven cavity.

The present study investigates turbulence/flame/wall interaction in turbulent boundary layer combustion with varying wall surface reactivity using direct numerical simulation (DNS). The effects of combustion on the turbulent boundary layer were investigated. It was found that the hairpin vortices of the boundary layer turbulence are lifted upward, upstream of the flame, due to the combustion-induced adverse pressure gradient. Variations in wall surface reactivity rarely affect the characteristics of the boundary layer turbulence. Flame/wall interaction with varying wall surface reactivity was analyzed. The configuration of one-dimensional laminar head-on quenching was also studied for reference. Results showed that wall heat flux reduces as the wall surface reactivity enhances in both laminar and turbulent cases, with a more pronounced reduction in turbulent cases. Near the inert wall, heat release overshoots due to radical recombination, while near the reactive wall, the heat release weakens because the recombination reactions are suppressed due to the adsorption and consumption of radicals by the wall surface reactions. The influential zones of wall heat loss and surface reactions reach the logarithmic layer and buffer layer of the turbulent boundary layer, respectively. Near-wall flame quenching was also examined. In the turbulent case with high wall surface reactivity, the quenching distance was found to be larger than with inert wall, whereas in the laminar case, wall surface reactivity has no effect on the quenching distance. The analysis of local quenching characteristics showed that head-on quenching dominates in turbulent cases, and its occurrence increases as wall surface reactivity intensifies.

We analyze the effect of a scale-dependent horizontal hyperdiffusivity in rotating Rayleigh-Bénard convection. The hyperdiffusion (HD) is parametrized by a cutoff wave number $k_0$ and a growth rate $q$ to limit its effects to the smaller scales, enabling simulations to be conducted at more rapid rotation rates and buoyancy forcing (as measured by the Ekman number Ek and Rayleigh number Ra, respectively) than in conventional direct numerical simulations (DNSs). We have performed 107 simulations spanning the ranges $E={10}^{-4}\text{--}{10}^{-7}$ and $\text{Ra}={10}^6\text{--}7\times {10}^{10}$ and directly compared HD and DNS heat transfer (measured by the Nusselt number Nu), flow speeds (quantified by the Péclet number Pe), force balances, and boundary layer thicknesses. We identify two different effects of HD that are particularly clear when $k_0$ is below the dominant scale of the flow $k_{\perp }$: At low supercriticality $\widetilde{\mathrm{Ra}}$, Nu, and Pe are increased compared to DNS due to a weakening of the rotational constraint; at high $\widetilde{\mathrm{Ra}}$, Nu and Pe are decreased compared to DNS owing to suppression of energy at scales greater than $k_{\perp }$. The thermal boundary layer thickness changes in proportion to Nu, while the mechanical boundary layer thickness is marginally affected because it is determined by a vertical balance, whereas HD is applied horizontally. Heat transfer diagnostics are more strongly affected by HD than flow diagnostics, while changing $k_0$ has a greater impact on solution behavior than changing $q$ over the range of values studied. For all HD simulations with $k_0\ge 3k_{\perp }$, all diagnostics remain within the standard deviation of the DNS.

This study investigates the retraction dynamics of viscous droplets on hydrophobic surfaces in the capillary-viscous regime. Using glycerol-water solutions on various hydrophobic substrates, we show that the observed droplet retraction rates can be described using the dynamic receding contact angle rather than the static one. When incorporating the dynamic receding contact angle into existing models, experimental data from different wettability substrates collapse onto a unified curve proportional to the inverse of the viscous timescale. Thus, we reveal that the retraction dynamics in this regime is governed by a balance between capillary forces and viscous dissipation, with the dynamic receding contact angle serving as the key parameter characterizing surface wettability effects. These findings establish a framework for predicting droplet behavior based on nonequilibrium wetting parameters. They extend previous observations from the inertial-capillary regime and provide insights for applications involving viscous droplet interactions with hydrophobic surfaces, including self-cleaning technologies and coating processes.

We study the one-dimensional capillary-driven retraction of a finite planar liquid sheet in the asymptotic regime where both the initial length-to-thickness ratio $l_0/h_0$ and a characteristic Ohnesorge number, Oh, which is inversely related to the Reynolds number of the flow, are large. In this regime, the fluid domain decomposes into two regions: a thin-film region governed by one-dimensional mass and momentum equations and a small tip region near the free edge described by a self-similar Stokes flow. Asymptotic matching between these regions yields an effective boundary condition for the thin-film region, representing a balance between viscous and capillary forces at the free edge. We find that surface tension drives the thin-film flow only through this boundary condition, while the local momentum balance is dominated by viscous and inertial stresses. We show that the thin-film flow possesses a conserved quantity, reducing the governing equations to a one-dimensional free-boundary problem in which the sheet thickness satisfies the heat equation with time-dependent boundary conditions. The reduced problem depends on a single dimensionless parameter $\mathcal{L}=l_0/\left(4h_0\text{Oh}\right)$. Numerical solutions of the reduced model agree well with previous studies based on different formulations, including full Navier-Stokes simulations, and reveal that the sheet undergoes distinct retraction regimes depending on $\mathcal{L}$ and a dimensionless time after rupture $T$. We derive asymptotic approximations for the thickness profile, velocity profile, and retraction speed during the early and late stages of retraction. At early times, the retraction speed grows as $T^{1/2}$, characteristic of diffusive processes, while at late times it decays as $1/T^2$ due to the finiteness of the sheet. Both results are consistent with previous findings. We also analyze an intermediate regime, not previously studied in detail, which arises for very long sheets ($\mathcal{L}\gg 1$). During this phase, the retraction speed approaches the classical Taylor-Culick value, despite the sheet not exhibiting a large circular rim as in classical studies of film rupture. When $T\approx \mathcal{L}$, the speed undergoes a sudden deceleration, which we characterize numerically, going from the Taylor-Culick speed to that predicted by the late-time asymptotics.

Broadband pressure fluctuations induced by tip vortex cavitation (TVC) instabilities pose significant challenges to the safe operation of marine propulsion systems, yet the underlying mechanisms remain poorly understood. To explore these mechanisms, this study employs an Euler-Lagrange hybrid model and develops an analytical model to quantify the instability intensity. Bubble-type vortex bursting is identified as the primary mode of instability under both wetted and cavitating conditions. In wetted conditions, instability intensity correlates strongly with the blade load variation, while in cavitating conditions, it decreases due to the reduced load variation. Nevertheless, sharp pressure rise in the bursting vortex core can still trigger intense TVC instability. This instability satisfies the zero-group-velocity condition, giving rise to significant broadband pressure fluctuations near the corresponding characteristic frequency. These findings provide insights for controlling flow instabilities in marine propulsion.

This work (1) introduces a compensation framework to predict attached-eddy signatures from instantaneous fields and (2) explores the logarithmic decay rate $A_1$ of the streamwise turbulence intensity $\overline{u^{\prime 2}}$, where $u^{\prime }$ is the streamwise velocity fluctuation. While the attached eddy hypothesis implies a constant logarithmic decay rate $A_1$ of $\overline{u^{\prime 2}}$ concerning the wall distance, $A_1$ varies significantly in previous studies. Very recently, energy spectra have been decomposed via spectral linear stochastic estimation (SLSE) and linear coherence function (LCF) to extract the attached-eddy contribution and avoid the contamination of detached eddies. Even after the decomposition, $A_1$ shows an evident increasing trend with Reynolds numbers. This work re-examines the original single-input/single-output (SI/SO) model of SLSE and identifies the near-wall universal signal (noise of the SI/SO system) as one possible reason for lower $A_1\mathrm{s}$ at low and moderate Reynolds numbers. Both rigorous mathematical derivations and intuitive physical understandings are provided to explain the impact of noise on the SLSE and LCF. Then, two theoretically equivalent compensation approaches are proposed and compared, both of which transfer the original SI/SO system to a multi-input/single-output (MI/SO) system by introducing an extra input, near-wall wall-normal velocity fluctuation $v_i^{\prime }$, to eliminate the impact of the near-wall universal signal. The refined attached-eddy spectral subcomponent also complies with the classic model. The decay rate $A_1$ becomes, overall, much larger, reducing the Reynolds number variation. It reaches the asymptotic value 0.98 for ${\text{Re}}_{\tau }\ge 4200$. Although lower-Reynolds-number cases still have smaller $A_1\mathrm{s}$, this work cannot exclude the probability of a constant $A_1$ through qualitative analysis of LCFs between the near-wall $u_i^{\prime }$, $v_i^{\prime }$, and $u_o^{\prime }$ in the logarithmic region using the theory of the MI/SO system. The comprehensive refinement framework for analyzing and eliminating the effect of noise is also applicable to other prediction problems, and promising results in the present work could inspire new frameworks for future investigation to improve the precision of extracting attached-eddy signatures for low-Reynolds-number cases.

The wake of a model wind turbine was experimentally investigated to better understand the dynamics of tip, hub, and tower vortices and their interactions. The free-spinning two-blade model turbine was positioned in a turbulent freestream flow downstream of an active grid, resulting in a tip speed ratio of approximately 4. Volumetric particle tracking velocimetry was used to capture the full near-field wake, up to three rotor diameters downstream of the turbine, and the vortex structures are described in terms of the phase-averaged vorticity and Lamb vectors. A detailed characterization was made in the tip vortex frame of reference by studying the evolution of the circulation, size, and orientation of the structures with downstream distance. The two tip vortices shed from each blade have different circulation magnitudes and merge around 1.5 rotor diameters downstream of the turbine, resulting in changes to the characteristics of the stronger vortex. After merging, the tip vortex is more elliptical, its local orientation is different, and its circulation begins to decay more rapidly. Finally, the positional scatter of instantaneous tip vortices and their development with downstream distance are quantified to explain the decreasing circulation of the phase-averaged structures downstream.

A flexible body can interact with an incident water wave passively and extract wave energy to generate a thrust force. To understand the hydrodynamic performance and the wave energy extraction mechanism in this passive propulsion, we use direct numerical simulation to investigate the two-dimensional unsteady motion of a passive flexible foil with a fixed leading edge under an incident wave superposed with a current flow at a Reynolds number of $\mathit{O}\left(1000\right)$. The results of kinematics studies show that the passive flapping motion is up–down asymmetric, with the trailing-edge displacement having a larger magnitude above the fixed leading edge. The peak-to-peak trailing-edge amplitude is a strong function of the reduced frequency, and a maximum trailing-edge amplitude is observed, indicating a resonance phenomenon. Passive propulsion is achieved under the wave with a maximum cruising speed of $4.5$ times of the chord length per unit time at a Reynolds number of 1000. The foil reaches a peak propulsive efficiency between the Strouhal number $\left(\text{St}\right)$ values of 0.2 and 0.4. At low $\text{St}$ values, the scaling relationship of the mean thrust coefficient $\left(\overline{C_T}\right)$ and propulsion efficiency $\left({\eta }_e\right)$ with respect to $\text{St}$ are similar between the passive propulsion of flexible foils under waves and the active heaving and pitching of rigid foils in free streams, and $\overline{C_T}$ scales with ${\text{St}}^2$. At high $\text{St}$ values, $\overline{C_T}$ scales with ${\text{St}}^3$. When $\text{St}$ is greater than a threshold value between 0.2 and 0.27, a passing-over leading-edge vortex is generated, which is responsible for the $\overline{C_T}\sim {\text{St}}^3$ scaling and the enhanced thrust force generated at high $\text{St}$.

The field of bubble dynamics and cavitation, motivated initially by applications in underwater acoustics and later in biomedical ultrasound, continues to evolve. This article describes the recent focus and ongoing developments in understanding and utilizing bubble dynamics in complex fluids. An introductory overview is given of the phenomena of spherical bubble dynamics in complex fluids. The emerging interest in using bubble dynamics to probe high-frequency viscoelastic properties of complex fluids is then highlighted. This interest is motivated by the fact that timescales as short as ${10}^{-6}\mathrm{s}$ can be probed, comparable to the relaxation timescales for complex fluids containing suspended particles or macromolecules. The focus of the article then shifts to specific examples from my research group. First, I present our work on linear rheology of soft solids using ultrasound-driven bubbles. I then discuss how the deformation and flow fields generated by acoustically driven bubbles can modify the microstructure and properties of complex fluids, bearing potential for bubble removal from yield-stress fluids. Finally, I describe our observations of self-assembly of locally ordered microstructures in complex fluids driven by bubble dynamics. In closing, I discuss the growing synergy between the communities of cavitation and rheology, which will be instrumental to address new challenges in the characterization and manipulation of complex fluids.

The modeling and simulation of turbulent combustion must account for the contribution of many chemical species and the effect of turbulence on their transport and reactions. Such reactions must span a wide range of timescales and often present bottlenecks in accelerating simulations. Data from experiments or simulations enables tools to accelerate simulations and develop accurate predictions of important turbulence-chemistry interactions. Several methods designed to exploit this data are presented and discussed. They are motivated by and rooted in traditional paradigms in turbulent combustion that rely heavily on the existence of a low-dimensional manifold for the composition space and its coupling with turbulent transport. The methods also rely on machine learning techniques for model order reduction, the extraction of closure from observations and models, and learning predictions of reaction rates and chemical states from low-dimensional descriptions of these states. These methods include principal component transport for combustion direct numerical simulations, methods to develop closure models and construct closure terms from multiscalar measurements, and deep operator networks for chemistry integration and acceleration. They provide pathways for efficient simulations of turbulent combustion and overcome the inherent limitations of predicting these complex flows.

Mixing at the interface between a convection zone and an adjacent, stably-stratified layer plays a crucial role in shaping the structure and evolution of stars and planets. In this work, we present a suite of 2D and 3D Boussinesq simulations that explore how bottom-driven convection penetrates a compositionally stratified region. Our results reveal two distinct regimes: a penetrative regime, where the convection zone steadily grows by entraining fluid from above, and a stalled regime, where growth halts and transitions to overshooting convection. We extend classical entrainment theory by incorporating thermal and compositional diffusion and by deriving a modified entrainment law that predicts interface speeds in the weak-diffusion limit. We show that convection stalls when the interface speed becomes comparable to the compositional diffusion speed and validate the transition between behaviors across a wide parameter space of Richardson and Lewis numbers. Such diffusively controlled stalling is unlikely to occur in stellar and planetary interiors, where the Lewis number is typically large and compositional diffusion is extremely slow. In these environments, compositional diffusion will merely slow the growth of the convection zone, and convective boundaries are expected to stall only in the presence of other curtailing mechanisms such as strong radiative diffusion or rapid rotation.

The dispersion of Lagrangian particle pairs is a fundamental process in turbulence, with implications for mixing, transport, and the statistical properties of particles in geophysical and environmental flows. While classical theories describe pair dispersion through scaling laws related to energy cascades, extreme events in turbulent flows can significantly alter these dynamics. This is especially important in stratified flows, where intermittency manifests itself also as strong updrafts and downdrafts. In this study, we investigate the influence of extreme events on the relative dispersion of particle pairs in stably stratified turbulence. Using numerical simulations we analyze the statistical properties of pair separation across different regimes, and quantify deviations from classical Richardson scaling. Our results highlight the role of extreme drafts in accelerating dispersion. These findings have important implications for turbulent mixing in natural systems, including atmospheric and oceanic flows, as well as applications in cloud microphysics and pollutant transport.

Surface bubbles in the environment or engineering configurations, such as the ocean-atmosphere interface, sparkling wine, or during volcanic eruptions typically live on contaminated surfaces. A particularly common type of contamination is surface active agents (surfactants). We consider the effect of insoluble surfactant on jet drop formation by bubble bursting. Contrary to the observed trend that surfactants decrease the ejected drop radius for bubbles with precursor capillary waves, we find that surfactants increase the ejected drop radius for bubbles without precursor capillary waves—a regime characteristic of small bubbles. Consequently, the results have fundamental implications for understanding aerosol distributions in contaminated conditions. We find that the trend reversal is due to the effect of Marangoni stresses on the focusing of the collapsing cavity. We demonstrate quantitative agreement on the jet velocity and drop size between laboratory experiments and numerical simulations by using the measured surface tension dependence on surfactant concentration as the equation of state for the simulations.

Understanding the water exit dynamics of solid bodies is of fundamental importance for trans-medium vehicles, projectiles, and other ocean engineering applications. The water exit of a neutrally buoyant sphere is investigated using large-eddy simulation combined with an overset mesh technique. The sphere is launched from a fixed depth with varying initial velocities and launch angles. The results reveal distinct features of velocity attenuation, rotational response, and trajectory deflection. The critical exit velocity exhibits a nonlinear dependence on the launch angle (${\theta }_0$). The exit process consists of two stages: in the fully submerged phase, the relative velocity loss decreases with increasing initial velocity, whereas in the partially submerged phase, the loss approaches a ${\theta }_0$-dependent constant. The sphere undergoes rotation whose magnitude converges to a ${\theta }_0$-dependent limit, with the rotation direction switching from clockwise to counterclockwise as ${\theta }_0$ increases. Lateral trajectory deflection occurs in most cases, except at ${\theta }_0=0^{\circ }$ and ${45}^{\circ }$, and is primarily governed by the interaction between the shed vortices and the free surface together with the image effect. Secondary contributions, including the sphere's rotation and asymmetric wave-induced pressure distribution, amplify the deflection magnitude but do not alter its direction. These findings provide new physical insight into the mechanisms of water exit and underscore the critical role of launch conditions, with implications for predicting and controlling the motion and stability of underwater vehicles, projectiles, and marine structures.

We use diffusing-wave spectroscopy to perform direct space- and time-resolved measurement of the dissipation rate $\epsilon$ at the boundary of a turbulent flow. We have shown in a previous publication that this technique provides maps of the dissipation rate of Newtonian fluids [E. Francisco, et al., Exp. Fluids 64, 156 (2023)]. Here, we apply the technique at the boundary of a turbulent flow generated in a square box by an impeller stirring the fluid. Although the measurement is made on a small region near the boundary, we show that the dissipation remains proportional to the injected power and follows the turbulent scaling $\epsilon \propto {\text{Re}}^3$, with Re being the Reynolds number ranging from $1.5\times {10}^4$ to $6\times {10}^5$. With this flow, there is no need for logarithmic corrections to reproduce the dissipation near the flat boundary. In addition, our setup allows us to measure the spatiotemporal fluctuations of the dissipation near the boundary. These fluctuations are quite large (the relative fluctuations are about 50%) and are well described by a log-normal distribution, as expected for the dissipation rate in the bulk of homogeneous and isotropic turbulence (HIT) but power density spectra do not correspond to those expected for HIT [A. N. Kolmogorov, et al., J. Fluid Mech. 13, 82 (1962)]. However, we observe the same trends in direct numerical simulation of turbulent channel flow of the John Hopkins tubulent database [Y. Li, et al., J. Turbul. 9, N31 (2008); E. Perlman, et al. (ACM, New York, 2007); J. Graham, et al., J. Turbul. 17, 181 (2015)].

Flows at planetary scales are generally driven by buoyancy and influenced by rotation. Rotating Rayleigh-Bénard convection (RRBC) is a practical and simple model that can be used to describe these systems. In RRBC, thermally induced convection occurs, which is influenced by the constant rotation it experiences. We study RRBC in a cylinder in the transition region between rotation-affected and rotation-dominated (also called geostrophic) convection. Experiments are performed to assess the dependence of the Nusselt number $\mathrm{Nu}$ (efficiency of convective heat transfer) on the Prandtl number $\mathrm{Pr}$ (ratio of kinematic viscosity over thermal diffusivity), a relation that is not explored much for geostrophic convection. By using water at different mean temperatures we can reach $2.8\le \mathrm{Pr}\le 6$. We study the relation between $\mathrm{Pr}$ and $\mathrm{Nu}$ at constant Ekman number $\mathrm{Ek}=3\times {10}^{-7}$ (an inverse measure for strength of rotation) for two different diameter-to-height aspect ratios ($\mathrm{\Gamma }=1/5$ and $1/2$) of the setup. The corresponding constant Rayleigh numbers (strength of thermal forcing) are $\mathrm{Ra}=1.1\times {10}^{12}$ and $1\times {10}^{11}$, respectively. Additionally, we measure the relation between the Rayleigh number $\mathrm{Ra}$ and $\mathrm{Nu}$ for $4\times {10}^{10}\le \mathrm{Ra}\le 7\times {10}^{11}, \mathrm{Ek}=3\times {10}^{-7}$, and $\mathrm{Pr}=3.7$. It is found that $\mathrm{Nu}$ exhibits a significant dependence on $\mathrm{Pr}$, even within this limited range. Increasing $\mathrm{Pr}$ by a factor 2 resulted in a decrease of $\mathrm{Nu}$ of about $25\%$. We hypothesize that the decrease of $\mathrm{Nu}$ is caused by the changing ratio of the thermal and kinetic boundary layer thicknesses as a result of increasing $\mathrm{Pr}$. We also consider the anticipated contributions of the wall mode to the heat transfer using sidewall temperature measurements.

The current state-of-the-art in wall-modeled large-eddy simulation for high-speed flows is relatively accurate for the wall shear stress (errors within 5%) but much less accurate for the wall heat flux (errors up to 20% for some conditions). The present work improves the heat flux prediction of the ODE-based equilibrium wall model using the eddy viscosity by Hasan, et al. [Phys. Rev. Fluids 8, L112601 (2023)] by introducing two new modeling components: a simple model for the diffusion of turbulence kinetic energy, and altered near-wall damping of the thermal eddy diffusivity. The proposed wall model is tested both a priori and a posteriori at Mach numbers up to 10, and is found to reduce the a priori error in the predicted heat flux to within 5%. The commutation error between averaging and solving the wall model is found to be about 2% and 5% for the shear stress and heat flux, respectively, causing the a posteriori errors to be up to 10%.