Orlandi et al. (J. Fluid Mech., vol. 770, 2015, pp. 424–441) present direct numerical simulations over a very wide Reynolds number range for plane Couette and Poiseuille flows. The results reveal new information on the abrupt nature of transition in these flows, and the comparisons between Couette and Poiseuille flows help to provide a clearer picture of Reynolds number trends, especially with regard to inner/outer layer interactions. The stress distributions give strong support to Townsend’s attached eddy hypothesis, particularly for the wall-parallel component where there has been little experimental data available. The results pose some intriguing questions regarding the reconciliation of the present results with data at higher Reynolds numbers in different canonical flows.

The near and very near wake of a flat plate with a circular trailing edge, with vigorous vortex shedding, is investigated with data from direct numerical simulations (DNS). Computations were performed for four different combinations of the Reynolds numbers based on plate thickness ($D$) and momentum thickness near the trailing edge (${\it\theta}$). Unlike the case of the cylinder, these Reynolds numbers are independent parameters for the flat plate. The objectives of the study are twofold, to investigate the entrainment process when the separating boundary layers are turbulent and to better understand the instability of the detached shear layers (DSLs). A visualization of the entrainment process, the effect of changing the ratio ${\it\theta}/D$ on entrainment and wake-velocity statistics, and a way of understanding entrainment in a phase-averaged sense via distributions of the turbulent transport rate are provided here. The discussion on shear-layer instability focuses on the role of log-layer eddies in the destabilization process, the effect of high-speed streaks in the turbulent boundary layer in the vicinity of the trailing edge on shear-layer vortex generation rates, and a relationship between the prevalence of shear-layer vortex generation and shedding phase that is a result of an interaction between the shedding process and the shear-layer instability mechanism. A power-law relationship between the ratio of shear-layer and shedding frequencies and the Reynolds numbers mentioned above is obtained. A discussion of the relative magnitudes of the exponents is provided. A second power-law relationship between shed-vortex strength and these two Reynolds numbers is also proposed.

The choriocapillaris is a capillary bed located in a thin layer adjacent to the outer retina and is part of the oxygen delivery system to the photoreceptors of the eye. The blood flow is approximately planar and is serviced by microvessels, which join the choriocapillaris through inlets perpendicular to its plane. Capillaries are densely organised and separated by avascular septal posts, which direct the blood flow. The capillary bed is composed of a juxtaposition of tessellating vascular units called lobules, which are filled and drained independently from each other. A theoretical analysis of the blood flow in an idealised model of a lobule of the choriocapillaris is developed and studied. Lobules are modelled as tessellating polygonal prisms, where the upper and lower surfaces correspond to planar parallel membranes. The septae are modelled as cylinders randomly distributed between the two membranes. Feeding arterioles and draining venules are modelled as inlets and outlets connecting at the lower surface of the prism perpendicularly to the plane of the lobule. An inlet is placed inside the lobule, while an outlet is placed at each of the vertices. The polygonal prism can be formally subdivided into a set of triangular prisms with one inlet and two outlets, each of them located at one of the vertices. The triangular prisms are taken to be isosceles, and are therefore characterised by a vertex angle ${\it\omega}$ at the inlet and a span $L$. The flow is viscously dominated, and is investigated in the lubrication limit, in which the characteristic thickness of the prism is much smaller than the diameter of the cylinders. As a result of the geometry, a stagnation point is located midway between the outlets. A separation streamline joins the inlet and the stagnation point. The pressure drop ${\rm\Delta}\tilde{p}$ and the average fluid particle residence time $\langle \tilde{T}\rangle$ are analysed as a function of the angle at the inlet ${\it\omega}$ and the septae volume fraction ${\it\Phi}$. When no cylinders are present (${\it\Phi}=0$), an analytical expression for the pressure field is calculated by conformal mapping. Close to the triangle walls, the flow is quasi-parallel and characterised by a shorter fluid particle residence time. In the vicinity of the stagnation point, the velocity decreases and the residence time diverges logarithmically with the distance to the stagnation streamline. The minimum in pressure drop corresponds to a maximum in residence time, and is obtained for ${\it\omega}={\rm\pi}/2$. Asymptotic expressions for the pressure drop and average residence time are formulated in both the limits $\Vert {\it\omega}\Vert \ll 1$ and $\Vert {\rm\pi}-{\it\omega}\Vert \ll 1$. The impact of ${\it\Phi}$ on the flow is characterised by solving the equations for the flow numerically and using the Darwin drift framework. We show that the pressure drop is approximately proportional to $1+2{\it\Phi}$ for relatively small ${\it\Phi}$, and that $\langle \tilde{T}\rangle$ is proportional to $1-{\it\Phi}$ regardless of the void fraction or shape of the septae. In the case ${\it\Phi}=0$, the average residence time equals the volume of the domain divided by the volumetric flux. This analysis provides a new perspective on the blood flow dynamics within the choriocapillaris. Lobules form systems, where perfusion and corpuscle transport are a function of the angle that any two venular openings make with an arteriolar opening, the surface area perfused and the void volume fraction. The blood flow velocity and residence time are significantly heterogeneous, which may be responsible for the high degree of selective localisation observed in the pathogenesis of some inflammatory and degenerative diseases of the eye.

The local flow topology is studied using the invariants of the velocity gradient tensor in compressible turbulent mixing layer via direct numerical simulation (DNS) data. The topological and dissipating behaviours of the flow are analysed in two different regions: in proximity of the turbulent/non-turbulent interface (TNTI), and inside the turbulent region. It is found that the distribution of various flow topologies in regions close to the TNTI differs from inside the turbulent region, and in these regions the most probable topologies are non-focal. In order to better understand the behaviour of different flow topologies, the probability distributions of vorticity norm, dissipation and rate of stretching are analysed in incompressible, compressed and expanded regions. It is found that the structures undergoing compression–expansion in axial–radial directions have the highest contraction rate in locally compressed regions, and in locally expanded regions the structures undergoing expansion–compression in axial–radial directions have the highest stretching rate. The occurrence probability of different flow topologies conditioned by the dilatation level is presented and it is shown that the structures in the locally compressed regions tend to have stable topologies while in locally expanded regions the unstable topologies are prevalent.

The self-preserving region of a free round turbulent air jet at high Reynolds number is investigated experimentally (at $x/D=30$, $\mathit{Re}_{D}=1.4\times 10^{5}$ and $\mathit{Re}_{{\it\lambda}}=548$). Air is slightly heated ($20\,^{\circ }\text{C}$ above ambient) in order to use temperature as a passive scalar. Laser doppler velocimetry and simultaneous laser doppler velocimetry–cold-wire thermometry measurements are used to evaluate turbulent kinetic energy and temperature variance budgets in identical flow conditions. Special attention is paid to the control of initial conditions and the statistical convergence of the data acquired. Measurements of the variance, third-order moments and mixed correlations of velocity and temperature are provided (including $\overline{vw^{2}}$, $\overline{u{\it\theta}^{2}}$, $\overline{v{\it\theta}^{2}}$, $\overline{u^{2}{\it\theta}}$, $\overline{v^{2}{\it\theta}}$ and $\overline{uv{\it\theta}}$). The agreement of the present results with the analytical expressions given by the continuity, mean momentum and mean enthalpy equations supports their consistency. The turbulent kinetic energy transport budget is established using Lumley’s model for the pressure diffusion term. Dissipation is inferred as the closing balance. The transport budgets of the $\overline{u_{i}u_{j}}$ components are also determined, which enables analysis of the turbulent kinetic energy redistribution mechanisms. The impact of the surrogacy $\overline{vw^{2}}=\overline{v^{3}}$ is then analysed in detail. In addition, the present data offer an opportunity to evaluate every single term of the passive scalar transport budget, except for the dissipation, which is also inferred as the closing balance. Hence, estimates of the dissipation rates of turbulent kinetic energy and temperature fluctuations (${\it\epsilon}_{k}$ and ${\it\epsilon}_{{\it\theta}}$) are proposed here for use in future studies of the passive scalar in a turbulent round jet. Finally, the budgets of turbulent heat fluxes ($\overline{u_{i}{\it\theta}}$) are presented.

Wind forcing of the ocean generates a spectrum of inertia–gravity waves that is sharply peaked near the local inertial (or Coriolis) frequency. The corresponding near-inertial waves (NIWs) are highly energetic and play a significant role in the slow, large-scale dynamics of the ocean. To analyse this role, we develop a new model of the non-dissipative interactions between NIWs and balanced motion. The model is derived using the generalised-Lagrangian-mean (GLM) framework (specifically, the ‘glm’ variant of Soward & Roberts, J. Fluid Mech., vol. 661, 2010, pp. 45–72), taking advantage of the time-scale separation between the two types of motion to average over the short NIW period. We combine Salmon’s (J. Fluid Mech., vol. 719, 2013, pp. 165–182) variational formulation of GLM with Whitham averaging to obtain a system of equations governing the joint evolution of NIWs and mean flow. Assuming that the mean flow is geostrophically balanced reduces this system to a simple model coupling Young & Ben Jelloul’s (J. Mar. Res., vol. 55, 1997, pp. 735–766) equation for NIWs with a modified quasi-geostrophic (QG) equation. In this coupled model, the mean flow affects the NIWs through advection and refraction; conversely, the NIWs affect the mean flow by modifying the potential-vorticity (PV) inversion – the relation between advected PV and advecting mean velocity – through a quadratic wave term, consistent with the GLM results of Bühler & McIntyre (J. Fluid Mech., vol. 354, 1998, pp. 301–343). The coupled model is Hamiltonian and its conservation laws, for wave action and energy in particular, prove illuminating: on their basis, we identify a new interaction mechanism whereby NIWs forced at large scales extract energy from the balanced flow as their horizontal scale is reduced by differential advection and refraction so that their potential energy increases. A rough estimate suggests that this mechanism could provide a significant sink of energy for mesoscale motion and play a part in the global energetics of the ocean. Idealised two-dimensional models are derived and simulated numerically to gain insight into NIW–mean-flow interaction processes. A simulation of a one-dimensional barotropic jet demonstrates how NIWs forced by wind slow down the jet as they propagate into the ocean interior. A simulation assuming plane travelling NIWs in the vertical shows how a vortex dipole is deflected by NIWs, illustrating the irreversible nature of the interactions. In both simulations energy is transferred from the mean flow to the NIWs.

The micro-convection caused by the ponderomotive forces of the self-magnetic field in a magnetic fluid is studied here both numerically and experimentally. The theoretical approach based on the general Brinkman model substantially improves the description with respect to the previously proposed Darcy model. The predictions of both models are here compared to finely controlled experiments. The Brinkman model, in contrast to the Darcy model, allows us to describe the formation of mushrooms on the plumes of the micro-convective flow and the width of the fingers. In the Brinkman approach, excellent quantitative agreement is also obtained for the finger velocity dynamics and the velocity maximal values as a function of the magnetic Rayleigh number. The diffusion coefficient of particles of the water-based magnetic colloid determined by the threshold field strength value of the micro-convection is significantly larger than the diffusion coefficient of individual particles. This result is confirmed by independent measurements of the diffusion coefficient at the smearing of the diffusion front.

We investigate the transport of relatively heavy, small particles through a downward displacement-ventilated space. A flux of particles is supplied to the space from a localised source at a high level and forms a turbulent particle-laden plume which descends through the space. A constant flow of ambient fluid which does not contain particles is supplied to the space at a high level, while an equal amount of fluid is vented from the space at a low level. As a result of the entrainment of ambient fluid into the particle plume, a return flow is produced in the ambient fluid surrounding the plume in the lower part of the space. At steady state, particles are suspended by this return flow. An interface is formed which separates the ambient fluid in the lower part of the space, which contains particles, from the particle-free ambient fluid in the upper part of the space. New laboratory experiments show that the concentration of particles in the ambient fluid below the interface is larger than the average concentration of particles in the plume fluid at the level of the interface. Hence, as the plume fluid crosses the interface and descends through the particle-laden fluid underneath, it becomes relatively buoyant and forms a momentum-driven fountain. If the fountain fluid impinges on the floor, it then spreads radially over the surface until lifting off. We develop a quantitative model which can predict the height of the interface, the concentration of particles in the lower layer, and the partitioning of the particle flux between the fraction which sediments over the floor and that which is ventilated out of the space. We generalise the model to show that when particles and negatively buoyant fluid are supplied at the top of the space, a three-layer stratification develops in the space at steady state: the upper layer contains relatively low-density ambient fluid in which no particles are suspended; the central layer contains a mixture of ambient and plume fluid in which no particles are suspended; and the lower layer contains a suspension of particles in the same mixture of ambient and plume fluid. We quantify the heights of the two interfaces which separate the three layers in the space and the concentration of particles in suspension in the ambient fluid in the lower layer. We then discuss the relevance of the results for the control of airborne infections in buildings. Our experiments show that the three-layer stratification is subject to intermittent large-scale instabilities when the concentration of particles in the plume at the source is sufficiently small, or the rate of ventilation of the space is sufficiently large: we describe the transient concentration of particles in the space during one of these instabilities.

We study the bifurcation structure of zonal flows on a rotating sphere. The setting of our problem is similar to the Kolmogorov problem on a flat torus, where the vorticity forcing is given by a single eigenfunction of the Laplacian. First we prove the global stability of two-jet zonal flow for arbitrary Reynolds number and the rotation rate of the sphere. Then we study the bifurcation structure of steady solutions arising from three-jet zonal flow. In the non-rotating case, we find that two steady travelling-wave solutions bifurcate from a three-jet zonal flow via Hopf bifurcation. As the Reynolds number increases, steady-travelling solutions arise via pitchfork bifurcation from the steady-travelling solutions. On the other hand, in the rotating case, we find saddle-node bifurcations and closed-loop branches. We carry out time integration to study the properties of unsteady solutions at high Reynolds numbers. In the non-rotating case, the unsteady solution is chaotic and it wanders around the steady-travelling solutions bifurcating from three-jet zonal flow. We show that a linear combination of the steady and steady-travelling solutions gives a good approximation of the zonal-mean zonal flow of the unsteady solution, suggesting that the chaotic solution at high Reynolds numbers exists mostly within a relatively low-dimensional space spanned by the steady and steady-travelling solutions, which become unstable at low Reynolds numbers.

The transient deformation of a weakly conducting (‘leaky dielectric’) drop under a uniform DC electric field is computed via an axisymmetric boundary integral method, which accounts for surface charge convection and a finite relaxation time scale over which the drop interface charges. We focus on drops that attain an ultimate oblate (major axis normal to the applied field) steady-state configuration. The computations predict that as the time scale for interfacial charging increases, a shape transition from prolate deformation (major axis parallel to the applied field) to oblate deformation occurs at intermediate times due to the slow buildup of charge at the surface of the drop. Convection of surface charge towards the equator of the drop is shown to weaken the steady-state oblate deformation. Additionally, convection results in sharp shock-like variations in surface charge density near the equator of the drop. Our numerical results are then compared with an experimental system consisting of a millimetre-sized silicone oil drop suspended in castor oil. Agreement in the transient deformation is observed between our numerical results and experimental measurements for moderate electric field strengths. This suggests that both charge relaxation and charge convection are required, in general, to quantify the time-dependent deformation of leaky dielectric drops. Importantly, accurate prediction of the observed modest deformation requires a nonlinear model. Discrepancies between our numerical calculations and experimental results arise as the field strength is increased. We believe that this is due to the observed onset of rotation and three-dimensional flow at such high electric fields in the experiments, which an axisymmetric boundary integral formulation naturally cannot capture.

This paper is focused on the investigation of vortex rings evolving in a tube. A new theoretical model for a confined axisymmetric vortex ring is developed. The predictions of this model are shown to be in agreement with available experimental data and numerical simulations. The model combines the viscous vortex ring model, developed by Kaplanski & Rudi (Phys. Fluids, vol. 17, 2005, 087101), with Brasseur’s (PhD thesis, Stanford University) approach to deriving a wall-induced streamfunction correction. Using the power-law assumption for the time variation of the viscous length of the vortex ring, the time variations of the main integral characteristics, circulation, kinetic energy and translational velocity are obtained. Direct numerical simulation (DNS) is used to test the range of applicability of the model and to investigate new physical features of confined vortex rings recently reported in the experimental study by Stewart et al. (Exp. Fluids, vol. 53, 2012, pp. 163–171). The model is shown to lead to a very good approximation of the spatial distribution of the Stokes streamfunction, obtained by DNS. The vortex signature and the time evolution of the energy of the vortex are also accurately predicted by the model. A procedure for fitting the model with realistic vortex rings, obtained by DNS, is suggested. This opens the way to using the model for practical engineering applications.

The pairing of helical tip vortices in the wake of a two-bladed rotor is investigated experimentally. Time-resolved particle image velocimetry measurements provide a clear temporal and spatial evolution of the vortical structures, highlighting the transition to instability and the effect of tip speed ratio and helical spacing. The temporal growth rate of the vortex system instabilities were measured and are shown to be dependent on helical spacing. The evolution of filaments and their growth rates support the argument that the mutual inductance mode is the driving mechanism behind the transition to an unstable wake. The measurements are in agreement with maximum growth rates predicted by linear stability analyses of single- and double-helix arrangements. In addition, the wake topology due to varying rotor load through tip speed ratio variation is shown to play an important role in the initial symmetry breaking that drives the wake transition.

This paper is concerned with the boundary-layer separation in subsonic and transonic flows caused by a two-dimensional isolated wall roughness. The process of the separation is analysed by means of two approaches: the direct numerical simulation (DNS) of the flow using the Navier–Stokes equations, and the numerical solution of the triple-deck equations. Since the triple-deck theory relies on the assumption that the Reynolds number ($\mathit{Re}$) is large, we performed the Navier–Stokes calculations at $\mathit{Re}=4\times 10^{5}$ based on the distance of the roughness element from the leading edge of the flat plate. This $\mathit{Re}$ is also relevant for aeronautical applications. Two sets of calculation were conducted with the free-stream Mach number $\mathit{Ma}_{\infty }=0.5$ and $\mathit{Ma}_{\infty }=0.87$. We used different roughness element heights, some of which were large enough to cause a well-developed separation region behind the roughness. We found that the two approaches generally compare well with one another in terms of wall shear stress, longitudinal pressure gradient and detachment/reattachment points of the separation bubbles (when present). The main differences were found in proximity to the centre of the roughness element, where the wall shear stress and longitudinal pressure gradient predicted by the triple-deck theory are noticeably different from those predicted by DNS. In addition, DNS predicts slightly longer separation regions.

The spectral model of Perry et al. (J. Fluid Mech., vol. 165, 1986, pp. 163–199) predicts that the integral length scale varies very slowly with distance to the wall in the intermediate layer. The only way for the integral length scale’s variation to be more realistic while keeping with the Townsend–Perry attached eddy spectrum is to add a new wavenumber range to the model at wavenumbers smaller than that spectrum. This necessary addition can also account for the high-Reynolds-number outer peak of the turbulent kinetic energy in the intermediate layer. An analytic expression is obtained for this outer peak in agreement with extremely high-Reynolds-number data by Hultmark et al. (Phys. Rev. Lett., vol. 108, 2012, 094501; J. Fluid Mech., vol. 728, 2013, pp. 376–395). Townsend’s (The Structure of Turbulent Shear Flows, 1976, Cambridge University Press) production–dissipation balance and the finding of Dallas et al. (Phys. Rev. E, vol. 80, 2009, 046306) that, in the intermediate layer, the eddy turnover time scales with skin friction velocity and distance to the wall implies that the logarithmic derivative of the mean flow has an outer peak at the same location as the turbulent kinetic energy. This is seen in the data of Hultmark et al. (Phys. Rev. Lett., vol. 108, 2012, 094501; J. Fluid Mech., vol. 728, 2013, pp. 376–395). The same approach also predicts that the logarithmic derivative of the mean flow has a logarithmic decay at distances to the wall larger than the position of the outer peak. This qualitative prediction is also supported by the aforementioned data.

We present azimuthal velocity profiles measured in a Taylor–Couette apparatus, which has been used as a model of stellar and planetary accretion disks. The apparatus has a cylinder radius ratio of ${\it\eta}=0.716$, an aspect ratio of ${\it\Gamma}=11.74$, and the plates closing the cylinders in the axial direction are attached to the outer cylinder. We investigate angular momentum transport and Ekman pumping in the Rayleigh-stable regime. This regime is linearly stable and is characterized by radially increasing specific angular momentum. We present several Rayleigh-stable profiles for shear Reynolds numbers $\mathit{Re}_{S}\sim O(10^{5})$, for both ${\it\Omega}_{i}>{\it\Omega}_{o}>0$ (quasi-Keplerian regime) and ${\it\Omega}_{o}>{\it\Omega}_{i}>0$ (sub-rotating regime), where ${\it\Omega}_{i,o}$ is the inner/outer cylinder rotation rate. None of the velocity profiles match the non-vortical laminar Taylor–Couette profile. The deviation from that profile increases as solid-body rotation is approached at fixed $\mathit{Re}_{S}$. Flow super-rotation, an angular velocity greater than those of both cylinders, is observed in the sub-rotating regime. The velocity profiles give lower bounds for the torques required to rotate the inner cylinder that are larger than the torques for the case of laminar Taylor–Couette flow. The quasi-Keplerian profiles are composed of a well-mixed inner region, having approximately constant angular momentum, connected to an outer region in solid-body rotation with the outer cylinder and attached axial boundaries. These regions suggest that the angular momentum is transported axially to the axial boundaries. Therefore, Taylor–Couette flow with closing plates attached to the outer cylinder is an imperfect model for accretion disk flows, especially with regard to their stability.

A slow uniform flow of a rarefied gas past a sphere with a uniform temperature is considered. The steady behaviour of the gas is investigated on the basis of the Boltzmann equation by a systematic asymptotic analysis for small Mach numbers in the case where the Knudsen number is finite. Introducing a slowly varying solution whose length scale of variation is much larger than the sphere dimension, the fluid-dynamic-type equations describing the overall behaviour of the gas in the far region are derived. Then, the solution in the near region which varies on the scale of the sphere size, described by the linearised Boltzmann equation, and the solution in the far region, described by the fluid-dynamic-type equations, are sought in the form of a Mach number expansion up to the second order, in a way that they are joined in the intermediate overlapping region. As a result, the drag is derived up to the second order of the Mach number, which formally extends the linear drag obtained by Takata et al. (Phys. Fluids A, vol. 5, 1993, pp. 716–737) to a weakly nonlinear case. Numerical results for the drag on the basis of the Bhatnagar–Gross–Krook (BGK) model are also presented.

A direct numerical simulation of incompressible channel flow at a friction Reynolds number ($\mathit{Re}_{{\it\tau}}$) of 5186 has been performed, and the flow exhibits a number of the characteristics of high-Reynolds-number wall-bounded turbulent flows. For example, a region where the mean velocity has a logarithmic variation is observed, with von Kármán constant ${\it\kappa}=0.384\pm 0.004$. There is also a logarithmic dependence of the variance of the spanwise velocity component, though not the streamwise component. A distinct separation of scales exists between the large outer-layer structures and small inner-layer structures. At intermediate distances from the wall, the one-dimensional spectrum of the streamwise velocity fluctuation in both the streamwise and spanwise directions exhibits $k^{-1}$ dependence over a short range in wavenumber $(k)$. Further, consistent with previous experimental observations, when these spectra are multiplied by $k$ (premultiplied spectra), they have a bimodal structure with local peaks located at wavenumbers on either side of the $k^{-1}$ range.

Detailed velocity measurements in the wake of a body of revolution are reported for pitch angles up to $12^{\circ }$, over an unprecedented range of Reynolds numbers ($2.4\times 10^{6}\leqslant \mathit{Re}_{L}\leqslant 30\times 10^{6}$). The body of revolution, an idealized submarine shape (DARPA SUBOFF), is mounted using a support that mimics a semi-infinite sail. The wake measurements at all pitch angles and Reynolds numbers reveal the presence of a pair of streamwise vortices of unequal strengths which tend to rotate around each other as they evolve downstream. Various attempts to perturb the upstream conditions on the body had no significant impact on the relative strength of the vortices. In addition, two different models, tested in two different wind tunnels, show similar asymmetries, and we propose that wake asymmetry appears to be a robust feature of this flow, a result previously only seen for sharp-nosed bodies at high angles of attack. It is also shown that the wake behaviour for $x/D>5$, in terms of the streamwise mean velocity and turbulence intensity distributions, appears to become invariant with Reynolds number for $\mathit{Re}_{L}>4.8\times 10^{6}$.

A method is developed to adjust a vorticity field to satisfy specified values for a finite number of Casimir invariants. The developed method is tested numerically for a neutral fluid in two dimensions. The adjusted vorticity field is adopted as an initial condition for simulated annealing (SA) of an incompressible, ideal neutral fluid and its magnetohydrodynamics (MHD), where SA enables us to obtain a stationary state of the fluid. Since the Casimir invariants are kept unchanged during the annealing process, the obtained stationary state has the required values of the Casimir invariants specified by our method.

The influence of the airfoil thickness on aerodynamic loads is investigated numerically for harmonically pitching airfoils at low incidence, under the incompressible and inviscid flow approximation. Force coefficients obtained from finite-volume unsteady simulations of symmetrical 4-digit NACA airfoils are found to depart from the linear Theodorsen model of an oscillating flat plate. In particular, the value of the reduced frequency resulting in the inversion – from clockwise to counter-clockwise – of the lift/angle-of-attack hysteresis curve is found to increase with the airfoil thickness. Both the magnitude and direction of the velocity vector due to pitching over the airfoil surface differ from their flat-plate values. During the upstroke, namely nose-up rotation, phase, this results in a decrease (increase) of the normal velocity magnitude over the upper (lower) surface of the airfoil. The opposite occurs during the downstroke phase. This is confirmed by comparing the computed pressure distribution to the flat-plate linear Küssner model. Therefore, beyond the inversion frequency, the lift coefficient of a finite-thickness airfoil is higher during upstroke and lower during downstroke than its flat-plate counterpart. A similar dependence is also found for the quarter-chord moment coefficient. Accordingly, a modification to the classical Theodorsen model is proposed to take into account the effects of the airfoil thickness on unsteady loads. The new model is found to accurately predict the unsteady aerodynamics of a thick symmetric and a slightly cambered airfoil with a maximum thickness in the range 4–24 %. The limits of the present inviscid flow analysis are assessed by means of numerical simulation of high Reynolds number ($\mathit{Re}=10^{6}$) flows.