2.1. Steady flow

The vortex-force theory developed by Prandtl (Reference Prandtl1918) and refined by Taylor (Reference Taylor1926) and Wu, Ma & Zhou (Reference Wu, Ma and Zhou2006) represents a significant advancement in steady aerodynamics. In inviscid flows, Prandtl (Reference Prandtl1918) derived the aerodynamic force as

(2.1) \begin{align} \boldsymbol{F} = \int _{V_{f\infty }} \rho \boldsymbol{u}\times \boldsymbol{\omega } \,\textrm{d} V. \end{align}

Here, $\rho$ is the constant fluid density, $\boldsymbol{u}$ is the velocity vector, $\boldsymbol{\omega }:=\boldsymbol{\nabla }\times \boldsymbol{u}$ is the vorticity vector and $V_{f\infty }$ is the entire fluid domain. The volume integrand $\rho \boldsymbol{u}\times \boldsymbol{\omega }$ , which is perpendicular both to the local velocity and to the vorticity, is termed vortex force (von Kármán & Burgers Reference von Kármán and Burgers1935; Saffman Reference Saffman1992) and (2.1) is referred to as the vortex-force theory (Alekseenko et al. Reference Alekseenko, Kuibin and Okulov2007). In 2-D flows, (2.1) yields the Kutta–Joukowski theorem (Kutta Reference Kutta1902; Joukowski Reference Joukowski1906), without drag component. In 3-D flows, the streamwise component of (2.1) is generally non-zero, and it is identified as the induced drag, a drag component caused by the wake-vortex-induced downwash on the wing (Prandtl Reference Prandtl1918). In viscous flows, Taylor (Reference Taylor1926) and Wu et al. (Reference Wu, Ma and Zhou2006) generalised the vortex-force theory by employing a finite control volume $V_{\kern-1pt f}$ whose downstream surface is perpendicular to the incoming flow, enabling the identification of the induced drag and other force components. Applications of viscous steady vortex-force theory can be found in the work of Marongiu & Tognaccini (Reference Marongiu and Tognaccini2010) and Zou et al. (Reference Zou, Wu, Gao, Liu, Kang and Shi2019).

2.2. Unsteady flow

Unsteady flows prevail in nature and engineering. They exist not only at high-Reynolds-number flows where the flow is turbulent, containing rich vortex dynamics, but also in flapping wings and unsteady incoming flows such as gusts. Force diagnostics in unsteady flows attracts great research interest due to its relevance to animal locomotion and bio-inspired design. It is also the main focus of the present study.

According to Newton’s laws of motion, the interactive force exerted on the fluid from the solid body can be related to the time rate of change of the fluid momentum. However, the volume integral of momentum over the entire fluid domain does not converge in incompressible flows (Saffman Reference Saffman1992), preventing the direct application of Newton’s second law in the fluid system. This difficulty is overcome by Burgers (Reference Burgers1920), Wu (Reference Wu1981) and Lighthill (Reference Lighthill1986) independently through replacing the fluid momentum with the vorticity moment, which is also named vortical impulse. Using the DMT identities shown in Appendix B, the fluid momentum in a finite control volume $V_{\kern-1pt f}$ equals

(2.2) \begin{align} \int _{V_{\kern-1pt f}}\rho \boldsymbol{u} \, \textrm{d} V = \frac {1}{k}\int _{V_{\kern-1pt f}}\rho \boldsymbol{x}\times \boldsymbol{\omega }\,\textrm{d} V - \frac {1}{k}\oint _{\partial {V_{\kern-1pt f}}}\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{u})\,\textrm{d} S. \end{align}

Here, $k+1$ equals the space dimension and $\partial V_{\kern-1pt f}$ is the boundary of the control volume $V_{\kern-1pt f}$ , as indicated in figure 1. If there is no vorticity outside the control volume $V_{\kern-1pt f}$ , it is proved that the force and moment experienced by the solid body are

(2.3) \begin{align} \boldsymbol{F} &= -\frac {\rho }{k}\frac {\textrm{d}}{\textrm{d} t}\int _{V_{\kern-1pt f}}\boldsymbol{x}\times \boldsymbol{\omega }\,\textrm{d} V - \frac {\rho }{k}\frac {\textrm{d}}{\textrm{d} t}\oint _{\partial\kern-1pt B}\boldsymbol{x}\times (\hat {\boldsymbol{n}}\times \boldsymbol{u})\,\textrm{d} S, \\[-12pt] \nonumber \end{align}

(2.4) \begin{align} \boldsymbol{M} &= \frac {\rho }{2}\frac {\textrm{d}}{\textrm{d} t}\int _{V_{\kern-1pt f}}|\boldsymbol{x}|^2\boldsymbol{\omega }\,\textrm{d} V + \frac {\rho }{2}\frac {\textrm{d}}{\textrm{d} t}\oint _{\partial\kern-1pt B}|\boldsymbol{x}|^2(\hat {\boldsymbol{n}}\times \boldsymbol{u})\,\textrm{d} S. \\[9pt] \nonumber \end{align}

Equation (2.3) decomposes the force into the time derivative of the total vorticity moment and a surface integral term that relates to the body’s motion. Equation (2.4) does a similar thing for the moment. They are referred to as impulse theory or vorticity moment theory in the literature. The requirement that no vorticity exists outside the control volume $V_{\kern-1pt f}$ is demanding in real applications, and therefore, variants of the impulse theory that apply to finite control volumes are developed (Noca Reference Noca1996; Li & Lu Reference Li and Lu2012; Wang et al. Reference Wang, Zhang, He and Liu2013). Notably, Li & Lu (Reference Li and Lu2012) took into account the contribution from vorticity outside $V_{\kern-1pt f}$ and developed a finite-domain vorticity moment theory that applies to an arbitrary control volume. Kang et al. (Reference Kang, Liu, Su and Wu2018) found that the control volume in the work of Li & Lu (Reference Li and Lu2012) can be minimised only to contain vorticity regions that connect the body, which greatly simplifies the force expression.

The impulse theory and its variants have achieved great success in applications, especially in elucidating the high-lift mechanism in insect flight. However, as noted by Wu et al. (Reference Wu, Ma and Zhou2006), the time derivative in the impulse theory is outside the volume integral, obscuring a direct correlation between the force and flow structures. To address this problem, Wu, Lu & Zhuang (Reference Wu, Lu and Zhuang2007) derived three integral force formulas with time derivatives embedded inside the integrals: an advection form, a diffusion form and a boundary form. These formulas allow for the visualisation of force contributions by analysing the spatial distributions of their integrands. Notably, the advection-form formula is also referred to as the unsteady vortex-force theory in Wu et al. (Reference Wu, Liu and Liu2018) because it contains the volume integral of the vortex force $\rho \boldsymbol{u}\times \boldsymbol{\omega }$ , which resembles the inviscid vortex-force theory (Prandtl Reference Prandtl1918) given by (2.1). While all three formulas accurately recover the total force in numerical test cases, their integrands exhibit remarkably distinct spatial distributions. Specifically, the vortex force is concentrated in the boundary layer, the free shear layer and the near-wake vortices. In contrast, the volume integrand of the diffusion-form formula is strongly localised near the solid wall.

By projecting the N-S equations onto a divergence- and curl-free vector function, Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983) developed an alternative approach to associate the fluid force and moment with the velocity field in incompressible flows. First, a set of auxiliary harmonic functions $\eta _{\overline {i}}$ with $\overline {i}=1,2,3,4,5,6$ are introduced, which satisfy ${\nabla} ^2\eta _{\overline {i}} = 0$ in the fluid domain and boundary conditions $\hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta _i = -\hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{e}_i$ and $\hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta _{i+3} = -(\boldsymbol{x}\times \hat {\boldsymbol{n}})\boldsymbol{\cdot }\boldsymbol{e}_i$ with $i=1,2,3$ on the solid surface. At the infinitely far field, $\eta _{\overline {i}}$ goes to zero. The pressure force $\boldsymbol{F}_\textit{p}$ and the moment due to pressure $\boldsymbol{M}_{\!p}$ can be expressed as

(2.5) \begin{align} \boldsymbol{F}_\textit{p}\boldsymbol{\cdot }\boldsymbol{e}_i := \oint _{\partial\kern-1pt B}-p\hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{e}_i\,\textrm{d} S = \oint _{\partial\kern-1pt B} \hat {\boldsymbol{n}}\boldsymbol{\cdot }(p\boldsymbol{\nabla }\eta _i) \,\textrm{d} S = -\int _{V_{f\infty }} \boldsymbol{\nabla } p\boldsymbol{\cdot }\boldsymbol{\nabla }\eta _i \,\textrm{d} V, \\[-28pt] \nonumber \end{align}

(2.6) \begin{align} \boldsymbol{M}_{\!p}\boldsymbol{\cdot }\boldsymbol{e}_i := \oint _{\partial\kern-1pt B}-p(\boldsymbol{x}\times \hat {\boldsymbol{n}})\boldsymbol{\cdot }\boldsymbol{e}_i\,\textrm{d} S = \oint _{\partial\kern-1pt B} \hat {\boldsymbol{n}}\boldsymbol{\cdot }(p\boldsymbol{\nabla }\eta _{i+3}) \,\textrm{d} S = -\int _{V_{f\infty }} \boldsymbol{\nabla } p\boldsymbol{\cdot }\boldsymbol{\nabla }\eta _{i+3} \,\textrm{d} V. \\[0pt] \nonumber \end{align}

Here, $\boldsymbol{e}_i, i=1,2,3$ are the unit basis vectors of the coordinate system. The force and moment can be subsequently obtained by the volume integral of the dot product of the N–S equations and $\boldsymbol{\nabla }\eta _{\overline {i}}$ . After applying the divergence theorem, the remaining volume integrand is $\rho (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla }\eta _{\overline {i}}$ . Protas (Reference Protas2007) further theoretically examined this method and found that the vorticity in the body-surface integral cannot be eliminated, which is problematic in experimental applications. It is noted that the gradient of the auxiliary function $\eta _i$ can be explained as a virtual velocity in the virtual power principle developed by Yu (Reference Yu2014).

Following the work of Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983), a variety of projection-based force diagnostics methods have been developed. Recognising that the vorticity is more localised compared with the velocity gradient, Howe (Reference Howe1989) and Chang (Reference Chang1992) reformulated the advection term as $\rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = \rho \boldsymbol{\omega }\times \boldsymbol{u}+\boldsymbol{\nabla }(\rho |\boldsymbol{u}|^2/2)$ and derived the force and moment expressions with $(\rho \boldsymbol{\omega }\times \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla }\eta _{\overline {i}}$ as the volume integrand. Yu (Reference Yu2014) and Gao & Wu (Reference Gao and Wu2019) noted that the Galilean invariance, which is a fundamental requirement for the force in Newtonian mechanics, is not inherently guaranteed in most force diagnostics methods due to their velocity dependence. To address this limitation, Gao et al. (Reference Gao, Zou, Shi and Wu2019) derived a force expression with a Galilean-invariant volume integrand $2\rho Q\eta _i$ , which is obtained from the identity $\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}) = -2\rho Q$ , with $Q$ being the second invariant of the velocity gradient tensor. The moment expression based on $2\rho Q\eta _{i+3}$ is obtained by Menon & Mittal (Reference Menon and Mittal2021). Li & Wu (Reference Li and Wu2018) and Li et al. (Reference Li, Wang, Graham and Zhao2020) utilised the force expression of Howe (Reference Howe1989) and developed a vortex force and moment map method that links the vortex motion to the force and moment. These theoretical methods have found wide applications, to name a few, in numerical simulations, see Martín-Alcántara & Fernandez-Feria (Reference Martín-Alcántara and Fernandez-Feria2019) and Chiu et al. (Reference Chiu, Tseng, Chang and Chou2023), and in experiments, see Díaz-Arriba et al. (Reference Díaz-Arriba, Protas, Acher, Pons and David2022), Gehlert, Andreu-Angulo & Babinsky (Reference Gehlert, Andreu-Angulo and Babinsky2023) and Ōtomo et al. (Reference Ōtomo, Gehlert, Babinsky and Li2025).

3.1. General weighted method

Assuming that $ \unicode{x1D652}=W_{\textit{ij}}\boldsymbol{e}_i\boldsymbol{e}_j$ is a smooth second-order tensor function defined in the fluid domain as well as on the body surface, the weighted momentum equation is derived as

(3.3) \begin{align} \rho \boldsymbol{a}\boldsymbol{\cdot }\unicode{x1D652} = \boldsymbol{\nabla }\boldsymbol{\cdot } \big({-}p \unicode{x1D652}+ \unicode{x1D64F}\boldsymbol{\cdot }\unicode{x1D652} \kern1.5pt \big) + p\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652} - \unicode{x1D64F}:\boldsymbol{\nabla } \unicode{x1D652}, \end{align}

with $ \unicode{x1D64F}:\boldsymbol{\nabla } \unicode{x1D652}=T_{\textit{ij}}(\partial W_{\textit{ik}}/\partial x_j)\boldsymbol{e}_k$ . The weighted force $\boldsymbol{F}_{\!w}$ is defined by

(3.4) \begin{align} \boldsymbol{F}_w = \boldsymbol{F}_{w,\textit{p}} + \boldsymbol{F}_{w{,}\kern-2pt f} := \oint _{\partial\kern-1pt B}{-p\hat {\boldsymbol{n}}\boldsymbol{\cdot }\unicode{x1D652} \,\textrm{d} S} + \oint _{\partial\kern-1pt B}{\hat {\boldsymbol{\tau }}\boldsymbol{\cdot }\unicode{x1D652} \,\textrm{d} S}, \end{align}

in which $\boldsymbol{F}_{w,\textit{p}}$ and $\boldsymbol{F}_{w{,}\kern-2pt f}$ denote the weighted pressure force and the weighted friction force, respectively.

The weight function $ \unicode{x1D652}$ , interpreted as the dyad of virtual velocity in the virtual power principle (Yu Reference Yu2014), also indicates the relative importance of local spatial regions. By increasing $W_{\textit{ij}}$ on a local surface area while reducing its magnitude $|W_{\textit{ij}}|$ in other regions, the local stress distribution is effectively represented through $\boldsymbol{F}_{\!w}$ . The total force and moment can also be recovered. If $ \unicode{x1D652}=\delta _{\textit{ij}}\boldsymbol{e}_i\boldsymbol{e}_j$ on the body surface with $\delta _{\textit{ij}}$ being the Kronecker delta, $\boldsymbol{F}_{\!w}$ equals the total fluid force, $\boldsymbol{F}$ ; whereas if $ \unicode{x1D652}=\epsilon _{\textit{ijk}}x_i\boldsymbol{e}_j\boldsymbol{e}_k$ on $\partial\! B$ with $\epsilon _{\textit{ijk}}$ being the permutation tensor, $\boldsymbol{F}_{\!w}$ recovers the total moment, $\boldsymbol{M}$ . Furthermore, if the body is rigid and $ \unicode{x1D652}$ equals $\hat {\boldsymbol{n}}\hat {\boldsymbol{n}}$ on $\partial\! B$ , $\boldsymbol{F}_{\!w}$ becomes the pressure force, $\boldsymbol{F}_\textit{p}$ , because $\hat {\boldsymbol{\tau }}$ is perpendicular to $\hat {\boldsymbol{n}}$ . On the contrary, if $ \unicode{x1D652}=\delta _{\textit{ij}}\boldsymbol{e}_i\boldsymbol{e}_j-\hat {\boldsymbol{n}}\hat {\boldsymbol{n}}$ on $\partial\! B$ , $\boldsymbol{F}_{\!w}$ equals the friction force, $\boldsymbol{F}_\textit{f}$ .

To establish quantitative relationships between the weighted force $\boldsymbol{F}_{\!w}$ and the flow field, (3.3) is integrated over the control volume $V_{\kern-1pt f}$ . After applying the divergence theorem to the term $\int _{V_{\kern-1pt f}}\boldsymbol{\nabla }\boldsymbol{\cdot }(-p \unicode{x1D652}+ \unicode{x1D64F}\boldsymbol{\cdot }\unicode{x1D652} \kern1.5pt )\textrm{d} V$ , the general weighted method is obtained as follows:

(3.5) \begin{align} \boldsymbol{F}_w = \int _{V_{\kern-1pt f}}{({-}\rho \boldsymbol{a}\boldsymbol{\cdot }\unicode{x1D652} + p\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652} - \unicode{x1D64F}:\boldsymbol{\nabla } \unicode{x1D652} \kern1.5pt) \,\textrm{d} V} + \oint _{\varSigma }{({-}p\boldsymbol{n}\boldsymbol{\cdot }\unicode{x1D652}+\boldsymbol{\tau }\boldsymbol{\cdot }\unicode{x1D652} \kern1.5pt)\, \textrm{d} S}. \end{align}

The last term in (3.5) is the weighted surface-stress integral over $\varSigma$ and it is denoted by $\boldsymbol{F}_{w,\varSigma }$ ,

(3.6) \begin{align} \boldsymbol{F}_{w,\varSigma } := \oint _{\varSigma }{({-}p\boldsymbol{n}\boldsymbol{\cdot }\unicode{x1D652}+\boldsymbol{\tau }\boldsymbol{\cdot }\unicode{x1D652} \kern1pt)\, \textrm{d} S}. \end{align}

In addition, the volume integral of $- \unicode{x1D64F}:\boldsymbol{\nabla } \unicode{x1D652}$ in (3.5), denoted by $\boldsymbol{F}_{w,{v\textit{is}}}$ , arises from the viscous effect and it can be transformed into the following four forms by applying the divergence theorem:

(3.7a) \begin{align} \boldsymbol{F}_{w,{v\textit{is}}} &:= \int _{V_{\kern-1pt f}}{ - \unicode{x1D64F}:\boldsymbol{\nabla } \unicode{x1D652} \, \textrm{d} V} = -2\mu \int _{V_{\kern-1pt f}}{ \unicode{x1D63F}:\boldsymbol{\nabla } \unicode{x1D652}\, \textrm{d} V}, \\[-12pt] \nonumber \end{align}

(3.7b) \begin{align} \boldsymbol{F}_{w,{v\textit{is}}} &=-\mu \int _{V_{\kern-1pt f}}{ (\boldsymbol{\nabla }\times \boldsymbol{\omega })\boldsymbol{\cdot }\unicode{x1D652} \, \textrm{d} V}- \oint _{\partial V_{\kern-1pt f}}{\boldsymbol{\tau }\boldsymbol{\cdot }\unicode{x1D652}\,\textrm{d} S}, \\[-12pt] \nonumber \end{align}

(3.7c) \begin{align} \boldsymbol{F}_{w,{v\textit{is}}} &= -\mu \int _{V_{\kern-1pt f}}{ \boldsymbol{\omega }\boldsymbol{\cdot }(\boldsymbol{\nabla }\times \unicode{x1D652} \kern1.5pt )\, \textrm{d} V}+ 2\mu \oint _{\partial V_{\kern-1pt f}}{\left [ (\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u})(\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652} \kern1.5pt ) - \boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla } \unicode{x1D652} \kern1.5pt ) \right ]\,\textrm{d} S}, \\[-12pt] \nonumber \end{align}

(3.7d) \begin{align} \boldsymbol{F}_{w,{v\textit{is}}} &= \mu \int _{V_{\kern-1pt f}}{ \boldsymbol{u}\boldsymbol{\cdot }\left [2\boldsymbol{\nabla }(\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652} \kern1.5pt ) - \boldsymbol{\nabla }\times (\boldsymbol{\nabla }\times \unicode{x1D652} \kern1.5pt )\right ] \textrm{d} V} - \mu \oint _{\partial V_{\kern-1pt f}}{(\boldsymbol{u}\boldsymbol{n}+\boldsymbol{n}\boldsymbol{u}):\boldsymbol{\nabla } \unicode{x1D652}\,\textrm{d} S}. \\[9pt] \nonumber \end{align}

The mathematical details are provided in Appendix C. Equations (3.7a )–(3.7d ) involve spatial derivatives of different orders for the velocity and the weight function. A proper choice of the form of $\boldsymbol{F}_{w,{v\textit{is}}}$ offers convenience in applications. For instance, if $ \unicode{x1D652}$ has an analytical expression while $\boldsymbol{u}$ is known only at discrete points, (3.7d ) can be utilised to avoid the numerical differentiation of $\boldsymbol{u}$ . It should be noted that the weighted friction force $\boldsymbol{F}_{\!w,\textit{f}}$ cancels out in (3.5) when (3.7b ) is employed.

In the general weighted method, the weight function $ \unicode{x1D652}$ is set as a second-order tensor to maximise the flexibility of the formula and also to ensure terms in (3.5) are vectors. If (3.5) is projected onto a basis vector $\boldsymbol{e}_j$ , this is equivalent to having a vector weight function $\boldsymbol{w}=W_{\textit{ij}}\boldsymbol{e}_i$ . Additionally, when $ \unicode{x1D652}\kern1pt$ is an isotropic tensor, i.e. $W_{\textit{ij}}=w\delta _{\textit{ij}}$ , this set-up is equivalent to a scalar weight function $w$ .

If the weight function $ \unicode{x1D652}\kern1pt$ is taken as an arbitrary function, the general weighted method is equivalent to the weak-form momentum balance equations, which contain all information of momentum transport in the control volume. However, one drawback of the general weighted method is that the pressure is involved in the volume integration. Because the pressure field is hard to measure in experiments and the pressure only serves as a Lagrangian multiplier that ensures the incompressible condition (Chorin, Marsden & Marsden Reference Chorin, Marsden and Marsden1993), it is preferred to eliminate pressure by restricting the weight function to a specific function space.

3.2. Weighted projection methods

The fluid motion can be decomposed into a transverse process, described by the vorticity transport equation, and a longitudinal process, characterised by pressure wave propagation, as established by Chu & Kovásznay (Reference Chu and Kovásznay1958), Liu et al. (Reference Liu, Shi, Zhu, Su, Zou and Wu2014) and Mao et al. (Reference Mao, Kang, Wu, Yu, Gao, Su and Lu2020). In incompressible flows, the longitudinal process is governed by the pressure Poisson equation owing to the infinite speed of sound. In this section we demonstrate that the general weighted method (3.5) captures the transverse process when $ \unicode{x1D652}$ is divergence free and the longitudinal process when $ \unicode{x1D652}$ is curl free.

3.2.1. Divergence-free weight function

Following the methodology of Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983), where the N–S equations are projected onto a divergence-free vector field, we assume that the weight function satisfies $\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652}={\boldsymbol{0}}$ . Consequently, the term $p\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652}$ vanishes in (3.5). The weighted projection method with a divergence-free weight function is given by

(3.8) \begin{align} \boldsymbol{F}_{{\kern-1pt}w} = \int _{V_{f}}{-\rho \boldsymbol{a}\boldsymbol{\cdot }\unicode{x1D652} \,\textrm{d} V} +\boldsymbol{F}_{{\kern-1pt}w,{v\textit{is}}} + \boldsymbol{F}_{{\kern-1pt}w,\varSigma }. \end{align}

As proposed by Howe (Reference Howe1989) and Chang (Reference Chang1992), replacing the acceleration $\boldsymbol{a}$ with $\partial \boldsymbol{u}/\partial t + \boldsymbol{\omega }\times \boldsymbol{u}+\boldsymbol{\nabla }(0.5|\boldsymbol{u}|^2)$ yields a localised integrand for uniform incoming flows, and hence, an alternative form of the weighted projection method with a divergence-free weight function is derived as

(3.9) \begin{align} \boldsymbol{F}_w = \int _{V_{f}}{\left [-\rho \frac {\partial \boldsymbol{u}}{\partial t}\boldsymbol{\cdot }\unicode{x1D652}+\rho (\boldsymbol{u}\times \boldsymbol{\omega })\boldsymbol{\cdot }\unicode{x1D652}\right ] \textrm{d} V} - \oint _{\partial V_{\kern-1pt f}}{0.5\rho |\boldsymbol{u}|^2\boldsymbol{n}\boldsymbol{\cdot }\unicode{x1D652} \,\textrm{d} S} +\boldsymbol{F}_{w,{v\textit{is}}} + \boldsymbol{F}_{w,\varSigma }. \end{align}

In (3.9) the weighted vortex-force term $\rho (\boldsymbol{u}\times \boldsymbol{\omega })\boldsymbol{\cdot }\unicode{x1D652}$ in the volume integrand underscores the contribution of the vortex force to the generation of aerodynamic or hydrodynamic force.

If we consider a subspace of the divergence-free weight functions whose total flux over the body surface is zero,

(3.10) \begin{align} \oint _{\partial\kern-1pt B}{\hat {\boldsymbol{n}}\boldsymbol{\cdot }\unicode{x1D652} \, \textrm{d} S} = {\boldsymbol{0}}, \end{align}

by vector calculus, this subspace is constituted by $\boldsymbol{\nabla }\times \boldsymbol{\varPsi }$ with $\boldsymbol{\varPsi }$ being an arbitrary second-order tensor function. As demonstrated in Appendix D, (3.8) can be cast into the equation

(3.11) \begin{align} \int _{V_{\kern-1pt f}}{\left (\rho \frac {\partial \boldsymbol{\omega }}{\partial t}+\rho \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega }-\rho \boldsymbol{\omega }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} - \mu {\nabla} ^2\boldsymbol{\omega }\right )\boldsymbol{\cdot }\boldsymbol{\varPsi }\,\textrm{d} V} + & \nonumber \\\oint _{\partial V_{\kern-1pt f}}{\left [\mu \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega } -\big ( \mu (\boldsymbol{n}\times \boldsymbol{\nabla })\times \boldsymbol{\omega } +\rho \boldsymbol{n}\times \boldsymbol{a}+\boldsymbol{n}\times \boldsymbol{\nabla } p\big )\right ]\boldsymbol{\cdot }\boldsymbol{\varPsi }\,\textrm{d} S}&=0, \end{align}

in which the volume integral reveals the vorticity transport equation and the surface integral yields the boundary vorticity flux (BVF) that measures the vorticity creation rate from the boundary (Wu & Wu Reference Wu and Wu1993). Therefore, the weighted projection method with a divergence-free weight function, (3.8) or (3.9), captures all the information of vorticity evolution and vorticity generation on the wall, i.e. the transverse process.

3.2.2. Curl-free weight function

If the weight function is curl-free and it has no circulation for multi-connected fluid domains, the weight function can be expressed as the gradient of a vector potential, $ \unicode{x1D652}=\boldsymbol{\nabla }\boldsymbol{\phi }$ . Here, $\boldsymbol{\phi }=\phi _i\boldsymbol{e}_i$ . Using vector identities, $-\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\phi }=-\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{a}\boldsymbol{\phi })+(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{a})\boldsymbol{\phi }$ and $-(\boldsymbol{\nabla }\times \boldsymbol{\omega })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\phi }=\boldsymbol{\nabla }\boldsymbol{\cdot }[({\nabla} ^2\boldsymbol{u})\boldsymbol{\phi } ]$ , and noting that $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{a}=-2Q$ , where $Q$ is the second invariant of the velocity gradient tensor, the weighted pressure force $\boldsymbol{F}_{w,\textit{p}}$ can be obtained from (3.5) by employing (3.7b ) for $\boldsymbol{F}_{w,{v\textit{is}}}$ as

(3.12) \begin{align} \boldsymbol{F}_{w,\textit{p}} = & \int _{V_{f}}{\big ({-}2\rho Q \boldsymbol{\phi } + p{\nabla} ^2\boldsymbol{\phi }\big )\, \textrm{d} V} + \oint _{\partial V_{\kern-1pt f}}{\boldsymbol{n}\boldsymbol{\cdot } \big({-}\rho \boldsymbol{a}+\mu {\nabla} ^2\boldsymbol{u} \big)\boldsymbol{\phi } \,\textrm{d} S} \nonumber \\&+ \oint _{\varSigma }{-p\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\phi }\, \textrm{d} S}. \end{align}

Equation (3.12) is the weighted projection method with a curl-free weight function. Note that the pressure is not eliminated in the volume integral because (3.12) establishes a quantitative relationship between $2\rho Q$ , which is the source term of the pressure Poisson equation, and the pressure throughout the control volume and on the boundary $\partial V_{\kern-1pt f}$ .

As illustrated in Appendix E, (3.12) can be transformed into

(3.13) \begin{align} \int _{V_{f}}{\phi _i \big({\nabla} ^2p-2\rho Q \big)\, \textrm{d} V}+\oint _{\partial V_{\kern-1pt f}}{\phi _i\boldsymbol{n}\boldsymbol{\cdot } \big({-}\boldsymbol{\nabla } p - \rho \boldsymbol{a}+\mu {\nabla} ^2\boldsymbol{u} \big)\,\textrm{d} S} = 0, \end{align}

where the volume integral yields the pressure Poisson equation and the surface integral reveals the Neumann pressure boundary condition. Therefore, (3.12) is equivalent to the pressure Poisson equation together with the pressure boundary condition, i.e. the longitudinal process.

If the weight function is both divergence- and curl-free, (3.8) and (3.12), representing the transverse process and the longitudinal process, respectively, are both applicable. In this case, the pressure source term $2\rho Q$ in (3.12) remains as the sole volume integrand, revealing the intrinsic relationship that the transverse process serves as the source of the longitudinal process.

3.3. Weighted DMT methods

In this section the three DMT methods (Wu et al. Reference Wu, Lu and Zhuang2007) are generalised to weighted formulations to enhance their capability and elucidate their interconnections.

3.3.1. Advection form

First, three vectors $\boldsymbol{f}_{\!w}$ , $\boldsymbol{l}_w$ and $\boldsymbol{b}_w$ are introduced, which are

(3.14) \begin{align} \rho\! \boldsymbol{f}_{\!w} := \rho \boldsymbol{a}\boldsymbol{\cdot }\unicode{x1D652} - p\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652} + \unicode{x1D64F}:\boldsymbol{\nabla } \unicode{x1D652}, \end{align}

the weighted Lamb vector $\boldsymbol{l}_w:=(\boldsymbol{\omega }\times \boldsymbol{u})\boldsymbol{\cdot }\unicode{x1D652}$ and

(3.15) \begin{align} \rho \boldsymbol{b}_w := \rho\! \boldsymbol{f}_{\!w}-\rho \boldsymbol{l}_w = \rho \frac {\partial \boldsymbol{u}}{\partial t}\boldsymbol{\cdot }\unicode{x1D652} +\boldsymbol{\nabla } \big( 0.5\rho |\boldsymbol{u}|^2 \big)\boldsymbol{\cdot }\unicode{x1D652} - p\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652} + \unicode{x1D64F}:\boldsymbol{\nabla } \unicode{x1D652}. \end{align}

The weighted N–S equations (3.3) can be rewritten as

(3.16) \begin{align} \rho\! \boldsymbol{f}_{\!w} = \boldsymbol{\nabla }\boldsymbol{\cdot }(-p \unicode{x1D652}+ \unicode{x1D64F}\boldsymbol{\cdot }\unicode{x1D652} \kern1.5pt ). \end{align}

It is noted that the volume integrand in the general weighted method (3.5) equals $-\rho\! \boldsymbol{f}_{\!w}$ and that $\rho (\boldsymbol{b}_w-\boldsymbol{f}_{\!w}) = \rho (\boldsymbol{u}\times \boldsymbol{\omega })\boldsymbol{\cdot }\unicode{x1D652}$ is the weighted vortex force.

The advection-form weighted DMT method is obtained by applying DMT identity (B1) to (3.5),

(3.17) \begin{align} \begin{split} \boldsymbol{F}_w = & \int _{V_{\kern-1pt f}}{\left [-\frac {\rho }{k}\boldsymbol{x}\times (\boldsymbol{\nabla }\times \boldsymbol{b}_w)+\rho (\boldsymbol{u}\times \boldsymbol{\omega })\boldsymbol{\cdot }\unicode{x1D652}\right ] \textrm{d} V} -\frac {\rho }{k}\oint _{\partial V_{\kern-1pt f}}{\boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}_w)\, \textrm{d} S} \\&+ \boldsymbol{F}_{w,\textit{DMT},B} + \boldsymbol{F}_{w,\textit{DMT},\varSigma }. \end{split} \end{align}

Here, the surface integrals are

(3.18) \begin{align} \boldsymbol{F}_{w,\textit{DMT},B} := \frac {\rho }{k}\oint _{\partial\kern-1pt B}{\boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{f}_{\!w})\, \textrm{d} S}, \\[-28pt] \nonumber \end{align}

(3.19) \begin{align} \boldsymbol{F}_{w,\textit{DMT},\varSigma } := \frac {\rho }{k}\oint _{\varSigma }{\boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{f}_{\!w})\, \textrm{d} S} + \oint _\varSigma {(-p\boldsymbol{n}\boldsymbol{\cdot }\unicode{x1D652}+\boldsymbol{\tau }\boldsymbol{\cdot }\unicode{x1D652} \kern1.5pt )\, \textrm{d} S}. \\[2pt] \nonumber \end{align}

In (3.17) the weighted vortex force contributes to the volume integral. Additionally, the pressure exists in the volume integrand in the form of $p\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D652}$ , which is generally non-zero for arbitrary weight functions.

Since no simplifications are made in the derivation of (3.17) and the weight function $ \unicode{x1D652}$ is arbitrary, the advection-form weighted DMT method contains the same information as the general weighted method (3.5), which represents the weak-form momentum balance equations.

3.3.2. Diffusion and boundary forms

An isotropic weight function is required to derive the diffusion-form and boundary-form DMT methods, i.e. $W_{\textit{ij}}=w\delta _{\textit{ij}}$ .

Following the work of Wu et al. (Reference Wu, Lu and Zhuang2007), the diffusion-form weighted DMT method is derived by substituting the weighted N–S equations (3.16) into the volume and outer surface integrals of the advection-form equation (3.17). The mathematical details can be found in Appendix F. Here, the weighted force is expressed as

(3.20) \begin{align} \boldsymbol{F}_w = -\frac {1}{k}\int _{V_{\kern-1pt f}}{\boldsymbol{x}\times \left [\boldsymbol{\nabla }\times (\boldsymbol{\nabla }\boldsymbol{\cdot }(w \unicode{x1D64F}\,))\right ] \textrm{d} V} + \boldsymbol{F}_{w,\textit{DMT},B} + \tilde {\boldsymbol{F}}_{w,\textit{DMT},\varSigma }, \end{align}

with $\tilde {\boldsymbol{F}}_{w,\textit{DMT},\varSigma }$ being

(3.21) \begin{align} \tilde {\boldsymbol{F}}_{w,\textit{DMT},\varSigma } := \frac {1}{k}\oint _{\varSigma }{\boldsymbol{x}\times \left [\boldsymbol{n}\times \left (\boldsymbol{\nabla }\boldsymbol{\cdot }(w \unicode{x1D64F}\,)\right )\right ] \textrm{d} S} + \oint _\varSigma {w\boldsymbol{\tau } \textrm{d} S}. \end{align}

In the derivation of (3.20), the weighted N–S equations are used for a second time in the control volume and on the outer surface, which leads to the cancellation of information in $V_{\kern-1pt f}$ and on $\varSigma$ . As demonstrated in Appendix G, the diffusion-form method predicts the same force value even when the flow velocity and pressure are arbitrarily disturbed in the interior of the fluid domain, indicating that (3.20) contains only the information on momentum balance on the body surface $\partial\! B$ .

The boundary-form DMT method is obtained by substituting the weighted N-S equations (3.16) into the surface integral of (3.20). The mathematical details can also be found in Appendix F. Here, the weighted force is

(3.22) \begin{align} \boldsymbol{F}_w = \frac {1}{k}\int _{\partial\kern-1pt B}{\boldsymbol{x}\times [\hat {\boldsymbol{n}}\times \boldsymbol{\nabla }(wp)] \,\textrm{d} S} + \int _{\partial\kern-1pt B}{w\hat {\boldsymbol{\tau }} \,\textrm{d} S}. \end{align}

Through DMT identity (B2), there is

(3.23) \begin{align} \boldsymbol{F}_{w,\textit{p}}=\int _{\partial\kern-1pt B}{-w p\hat {\boldsymbol{n}} \,\textrm{d} S} = \frac {1}{k}\int _{\partial\kern-1pt B}{\boldsymbol{x}\times [\hat {\boldsymbol{n}}\times \boldsymbol{\nabla }(w p)] \,\textrm{d} S}. \end{align}

Therefore, the boundary-form DMT method is equivalent to the weighted surface-stress integral (3.4), which reflects the fluid stress distribution on the $\partial\! B$ .

If $B$ is a rigid body, the generalised Caswell formula (Wu et al. Reference Wu, Ma and Zhou2006) shows $2 \unicode{x1D63F}=\hat {\boldsymbol{n}}(\boldsymbol{\omega }_r\times \hat {\boldsymbol{n}})+(\boldsymbol{\omega }_r\times \hat {\boldsymbol{n}})\hat {\boldsymbol{n}}$ , in which $\boldsymbol{\omega }_r:=\boldsymbol{\omega }-2\boldsymbol{\varOmega }$ is the relative vorticity and $\boldsymbol{\varOmega }$ is the angular velocity of the rigid body. Hence, the viscous stress equals $\hat {\boldsymbol{\tau }}=\mu \boldsymbol{\omega }_r\times \hat {\boldsymbol{n}}$ . Using the DMT identities (B3) and (B2), $\boldsymbol{F}_{\!w,\textit{f}}$ can be rewritten as

(3.24a) \begin{align} \boldsymbol{F}_{\!w,\textit{f}} &= \mu \int _{\partial\kern-1pt B}{\boldsymbol{x}\times [(\hat {\boldsymbol{n}}\times \boldsymbol{\nabla })\times (w\boldsymbol{\omega }_r)] \,\textrm{d} S} \text{ in 3-D space}, \\[-12pt] \nonumber \end{align}

(3.24b) \begin{align} \boldsymbol{F}_{\!w,\textit{f}} &= -\mu \int _{\partial\kern-1pt B}{\boldsymbol{x}\left [(\boldsymbol{e}_z\times \hat {\boldsymbol{n}})\boldsymbol{\cdot }\boldsymbol{\nabla }(w\omega _{rz}) \right ]\,\textrm{d} S} \text{ in 2-D space.} \\[9pt] \nonumber \end{align}

In 2-D space, $\boldsymbol{\omega }_r=\omega _{rz}\boldsymbol{e}_z$ .

3.4. Discussions

A summary of the weighted integral methods developed in this work is presented in table 1, along with the information captured by each method. First, the weighted surface-stress integral (3.4) generalises the standard surface-stress integral (1.2) with additional degrees of freedom in the weight function, enabling evaluation of the fluid stress distribution over the body surface $\partial\! B$ . Second, by restricting the weight function to be divergence-free or curl-free, the projection methods (3.8) and (3.12) capture the transverse process and the longitudinal process, governed by the vorticity transport equation and the pressure Poisson equation, respectively. Specifically, the divergence-free weight function eliminates the pressure term in the volume integral. Third, the three forms of weighted DMT methods are not equivalent. The advection form (3.17) contains the momentum balance equations, i.e. the N-S equations, in the control volume $V_{\kern-1pt f}$ and on the boundary $\partial V_{\kern-1pt f}$ , whereas the diffusion form (3.20) captures only the momentum balance on the body surface $\partial\! B$ . The boundary form (3.22) is equivalent to the weighted surface-stress integral (3.4) that reveals the stress distribution on the body surface. Finally, time derivatives in the methods presented in this work are all embedded inside the integral. A weighted generalisation of the impulse theory, which has the time derivative outside the integral, remains unexplored here.

Table 1.

Summary of the weighted integral methods.

The weighted integral methods recover the original control-volume integral methods when special weight functions are adopted. If $W_{\textit{ij}}=-\partial \eta _j/\partial x_i$ , the weighted projection method with divergence-free weight function (3.9) reduces to the force expression developed by Chang (Reference Chang1992), and the weighted projection method with curl-free weight function (3.12) becomes the force expression developed by Gao et al. (Reference Gao, Zou, Shi and Wu2019). The case is also valid for the moment when $\eta _{j+3}$ is employed. Here, $\eta _j$ and $\eta _{j+3}$ that appear in (2.5) and (2.6) are first introduced by Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983). For the weighted DMT methods, if $ \unicode{x1D652}$ is a second-order unit tensor, i.e. $W_{\textit{ij}}=\delta _{\textit{ij}}$ , extra terms involving spatial derivatives of $ \unicode{x1D652}$ vanish and, as a result, the original DMT methods are recovered.

It is widely acknowledged that the force in Newtonian mechanics is independent of the inertial frame of reference, known as Galilean invariance. The Galilean invariance in fluid force decomposition has been addressed by Yu (Reference Yu2014), who claimed that each integral term in the force expression needs to be Galilean invariant, and by Gao & Wu (Reference Gao and Wu2019), who emphasised the Galilean invariance of integrals over local fluid regions. In aerodynamic and hydrodynamic force diagnostics, the Galilean invariance can be satisfied at three levels. First, for a prescribed weight function, the total weighted force $\boldsymbol{F}_{\!w}$ is Galilean invariant, which is already guaranteed owing to the rigorous derivation of all methods. Second, when evaluating force contributions from local flow structures, the local integrals must be independent of the inertial frame of reference. Third, when analysing the spatial distributions of volume or surface integrands, the relevant integrands must be Galilean invariant. For example, integrands in (3.4), (3.5), (3.8) and (3.12) are all Galilean invariant, provided that $ \unicode{x1D652}$ is Galilean invariant. The second or third level is not necessarily satisfied in a specific application when only the total force is concerned. A practical criterion is that the focused force or force contribution should be Galilean invariant.

4.1. 2-D flow past a circular cylinder

The uniform flow past a stationary circular cylinder is employed to assess the performance of various weighted integral methods. The incoming flow velocity $U_\infty$ , the fluid density $\rho$ and the cylinder diameter $d$ are designated as reference quantities. The Reynolds number is defined as $Re=\rho U_\infty d/\mu =100$ . The fluid velocity, expressed as $\boldsymbol{u}=u_x\boldsymbol{e}_x+u_y\boldsymbol{e}_y$ , equals $U_\infty \boldsymbol{e}_x$ at the infinitely far field and zero on the cylinder surface. The pressure at infinity, set to $p_\infty = 0$ , is adopted as the reference pressure.

Figure 4.

Vorticity field in the flow past a circular cylinder at $Re=100$ . (a) Instantaneous flow at the instant of maximum lift. Black lines are contour lines of $Q/(U_\infty ^2d^{-2})=0.1$ . Here V0, V2 and V4 label the most recently formed vortices from the upper surface; V1 and V3 are the most recently formed vortices from the lower surface. (b) Time-averaged flow. Black lines with arrows show the velocity profiles.

The numerical simulation is conducted using a high-order spectral/hp element method, implemented in the open-source code, Nektar++ (Cantwell et al. Reference Cantwell2015). The computational domain is $-100d\le x\le 250d$ , $-100d\le y\le 100d$ , with mesh refinement in the subdomain $-5d\le x \le 160d$ , $-5d\le y\le 5d$ . The minimum element size on the cylinder surface is $0.01d\times 0.02d$ and the maximum element size in the wake is $0.2d\times 0.2d$ . The centre of the cylinder is located at $x=0, y=0$ . The computational domain contains 57 009 elements, with a complete polynomial expansion up to the sixth degree applied within each element. The second-order velocity-correction scheme (Karniadakis, Israeli & Orszag Reference Karniadakis, Israeli and Orszag1991) with a time step of $0.00024 d/U_\infty$ is utilised for the time integration. According to previous studies (Gao et al. Reference Gao, Zou, Shi and Wu2019, Reference Gao, Cantwell, Son and Sherwin2023), this spatial and temporal resolution suffices to fully resolve the cylinder wake in the range of $x\le 160d$ . The flow is simulated for a total duration of $500 d/U_\infty$ to achieve a periodic flow state. With the flow period denoted by $T$ , the dimensionless vortex-shedding frequency is given by $d/(T U_\infty )=0.1636$ . Validation of the numerical results is provided in Appendix H.

The time-dependent drag $F_x$ and lift $F_y$ , non-dimensionalised by $\rho U_\infty ^2d$ , are presented in figure 3. Here, the dimensionless force, rather than the force coefficient, is adopted to facilitate direct comparison with the control-volume integral results. The dimensionless drag fluctuates within a narrow range of $[0.6535, 0.6623]$ , with a mean value of 0.6579, whereas the dimensionless lift varies from −0.1587 to 0.1587, with a mean value of 0. The pressure force, indicated by the red solid line, and the friction force, shown by the black dashed line, are also provided in figure 3. Whilst the pressure force dominates, the friction force remains non-negligible, accounting for 26 % of the total drag and 12 % of the maximum lift.

The instantaneous vorticity field $\omega _z$ , at the moment of maximum lift, is shown in figure 4(a). At this instant, a vortex, identified by the $Q$ criterion and labelled V0, has just formed from the upper surface of the cylinder. The wake exhibits a typical Bérnad-von Kármán vortex street, with concentrated vortices significantly weakened for $x\gt 60d$ under the influence of viscosity. The time-averaged vorticity field is provided in figure 4(b), alongside mean velocity profiles, which are symmetric about the $x$ axis. A recirculation region extends to $x=1.925d$ . The near-wake region ( $x\lt 10d$ ) displays relatively strong nonlinearities, whereas the farther downstream region shows a gradual decay of vorticity. The velocity defect on the centreline persists: at $x=100d$ , the centreline velocity is $0.7780 U_\infty$ .

Unless otherwise specified, the control volume $V_{\kern-1pt f}$ in the following analyses is taken as an Eulerian domain $-5 d\le x \le x_{w\textit{ake}}$ , $-5d\le y\le 5d$ and $x^2+y^2\ge 0.25d^2$ in the body frame, with $x_{w\textit{ake}}$ varying from $0.8d$ to $150d$ . Because the disturbance velocity and pressure decay rapidly in the upstream direction and in the upward and downward directions, the surface integral over $\varSigma$ can be reduced to a wake-plane integral along the vertical line $-5d \le y\le 5d$ at $x=x_{w\textit{ake}}$ .

4.2. Total drag from kinetic energy balance

Analysis of the kinetic energy budget of the disturbance flow provides insight into the power input, energy flux and dissipation in the flow system. By setting the projected weight function $ \unicode{x1D652}\boldsymbol{\cdot }\boldsymbol{e}_x$ as the disturbance velocity $\boldsymbol{u}^{\prime}:= \boldsymbol{u}-U_\infty \boldsymbol{e}_x$ , which becomes zero at infinity and equals $-U_\infty \boldsymbol{e}_x$ on the cylinder surface, the kinetic energy balance relation is obtained from (3.8) as

(4.1) \begin{align} \begin{split} \underbrace {U_\infty F_x}_{\text{Input power}} &= \underbrace {\frac {\textrm{d}}{\textrm{d} t}\int _{V_{f}}{\frac {1}{2}\rho |\boldsymbol{u}^{\prime}|^2 \textrm{d} V}}_{\text{Kinetic energy increase}} + \underbrace {\int _{w\textit{ake}}{\frac {1}{2}\rho |\boldsymbol{u}^{\prime}|^2 u_x\textrm{d} S}}_{\text{Kinetic energy flux}} + \underbrace {\int _{V_{\kern-1pt f}}{ \unicode{x1D64F}:\boldsymbol{\nabla }\boldsymbol{u}^{\prime}\textrm{d} V}}_{\text{Dissipation rate}} \\& \qquad +\underbrace {\int _{w\textit{ake}}{pu^{\prime}_x \textrm{d} S}}_{\text{Pressure power}} +\underbrace {\int _{w\textit{ake}}{-\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{u}^{\prime} \textrm{d} S}}_{\text{Viscous stress power}} . \end{split} \end{align}

On the wake plane, the volume flux is $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{e}_x=u_x$ and the streamwise disturbance velocity is $\boldsymbol{u}^{\prime}\boldsymbol{\cdot }\boldsymbol{e}_x=u^{\prime}_x$ . As the lift force is perpendicular to the cylinder’s velocity and, thus, consumes no power, only the drag force is yielded. Equation (4.1) indicates that the power supplied to the disturbance flow to overcome cylinder drag is partially converted into kinetic energy and dissipated into heat within the control volume $V_{\kern-1pt f}$ . The remainder leaves the control volume via the kinetic energy flux through the wake plane, as well as the power exerted by the stress on the wake plane.

Figure 5 shows the distribution of the kinetic energy, $0.5\rho |\boldsymbol{u}^{\prime}|^2$ , and the dissipation rate, $2\mu \unicode{x1D63F}: \unicode{x1D63F}$ , of the disturbance flow. The dissipation rate is rescaled to $ \unicode{x1D63F}: \unicode{x1D63F}=D_{\textit{ij}}D_{\textit{ij}}$ due to the small viscosity coefficient $\mu$ . Both the kinetic energy and the dissipation rate are mainly distributed in the near wake, and their strengths decrease gradually with the downstream location. The kinetic energy attains its local maximum value around the centreline, i.e. the $x$ axis, whereas the dissipation rate reaches local maximum values in the surroundings of concentrated vortices. In the vortex centre, where the motion of fluid elements is dominated by rigid rotation, the dissipation rate is relatively weak.

Figure 5.

Instantaneous kinetic energy of the disturbance flow and rescaled dissipation rate at the instant of maximum lift. (a) Kinetic energy, $0.5\rho |\boldsymbol{u}^{\prime}|^2/(\rho U_\infty ^2)$ . (b) Rescaled dissipation rate, $ \unicode{x1D63F}: \unicode{x1D63F}/(U_\infty ^2d^{-2})$ . Black contour lines represent dimensionless vorticity, with negative values (−0.1 and −0.03) as dashed lines and positive values (0.03 and 0.1) as solid lines.

Figure 6(a) shows the variation of terms in (4.1) with respect to the wake-plane location, $x_{w\textit{ake}}$ . The total drag calculated by the control-volume integral (4.1), as indicated by the black dash-dotted line, agrees well with the surface-stress integral shown by the star symbol ‘*’, with a maximum relative difference of $0.2\,\,\%$ for $0.8d \le x_{w\textit{ake}} \lt 47d$ and $0.03\,\,\%$ for $47d\le x_{w\textit{ake}}\le 150d$ . The slightly larger error in the near wake is due to the insufficiency of the wake-plane size. The dissipation-rate term (black line with inverted triangles) dominates the right-hand side of (4.1) and it grows monotonically with $x_{w\textit{ake}}$ because its integrand $ \unicode{x1D64F}:\boldsymbol{\nabla }\boldsymbol{u}^{\prime}=2\mu \unicode{x1D63F}: \unicode{x1D63F}$ is always non-negative. The kinetic energy flux is also persistent due to the slow decay of the disturbance velocity in the wake, and it has a value of $0.08F_x U_\infty$ at the downstream location of $x_{w\textit{ake}}=100d$ . Additionally, both the kinetic energy flux (red line with circles) and the kinetic energy increase (blue line with triangles) exhibit significant oscillations when the wake plane cuts vortices in the range of $x_{w\textit{ake}}\lt 60d$ . The pressure power on the wake plane drops quickly and becomes less than $0.05F_x U_\infty$ for $x_{w\textit{ake}}\gt 7d$ . The viscous stress power on the wake plane remains negligible (less than $0.01F_x U_\infty$ ) throughout the wake range examined.

Figure 6.

Drag decomposition based on the kinetic energy balance. (a) Instantaneous drag at the instant of maximum lift. (b) Time-averaged drag.

To eliminate the unsteady effect, the time-averaged form of (4.1) is derived as

(4.2) \begin{align} U_\infty \overline {F_x} = \int _{w\textit{ake}}{\frac {1}{2}\rho \overline {|\boldsymbol{u}^{\prime}|^2 u_x}\,\textrm{d} S} + \int _{V_{\kern-1pt f}}{\overline { \unicode{x1D64F}:\boldsymbol{\nabla }\boldsymbol{u}^{\prime}}\,\textrm{d} V} + \int _{w\textit{ake}}{\overline {pu^{\prime}_x} \,\textrm{d} S} + \int _{w\textit{ake}}{-\overline {\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{u}^{\prime}} \,\textrm{d} S}. \end{align}

Since the flow is periodic, the time-averaged kinetic-energy-increase term vanishes. As shown in figure 6(b), terms in (4.2) exhibit similar trends to the instantaneous terms. However, oscillations disappear for the kinetic energy flux, which first grows due to the increase in the advection velocity, $u_x$ , and then decreases because of the decay of kinetic energy, $0.5\rho |\boldsymbol{u}^{\prime}|^2$ .

In conclusion, the power input that overcomes the drag force of the moving cylinder is mainly dissipated into heat. A small but not negligible portion of power is advected far downstream via the kinetic energy flux. In the near wake of $x_{w\textit{ake}}\lt 7d$ , the power transferred downstream by the pressure stress is also significant.

4.3. Total force from the advection-form weighted DMT method

The advection-form DMT method is characterised by the explicit decomposition of the vortex force. Since the concept of vortex force arises from steady aerodynamics (Prandtl Reference Prandtl1918) where the velocity employed is the relative velocity between the incoming flow and the body, it is natural to fix the frame of reference on the stationary cylinder, which yields the fluid velocity $\boldsymbol{u}$ . By setting $W_{\textit{ij}}=\delta _{\textit{ij}}$ in (3.17), the original advection-form DMT method is obtained as

(4.3) \begin{align} \begin{split} \boldsymbol{F} &= \underbrace {\int _{V_{\kern-1pt f}}{\rho \boldsymbol{u}\times \boldsymbol{\omega } \,\textrm{d} V}}_{\text{Vortex force}} +\underbrace {\int _{V_{\kern-1pt f}}{-\rho \boldsymbol{x}\times \frac {\partial \boldsymbol{\omega }}{\partial t} \,\textrm{d} V}}_{\partial _t\boldsymbol{\omega } \text{ term}}+\underbrace {\oint _{\partial\kern-1pt B}{-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}) \,\textrm{d} S}+ \boldsymbol{F}_{\textit{DMT},B}}_{\text{Body-surface term}}\\& \qquad +\underbrace {\int _{w\textit{ake}}{-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}) \,\textrm{d} S}}_{\text{Wake-plane integral of $-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l})$}} + \boldsymbol{F}_{\textit{DMT},\varSigma }. \end{split} \end{align}

The force expression consists of the volume integration of the vortex force, the negative Eulerian time derivative of vorticity momentum, the body-surface term, the wake-plane integral of $-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l})$ and $\boldsymbol{F}_{\textit{DMT},\varSigma }$ . Here, $\boldsymbol{F}_{\textit{DMT},B}$ and $\boldsymbol{F}_{\textit{DMT},\varSigma }$ denote $\boldsymbol{F}_{w,\textit{DMT},B}$ and $\boldsymbol{F}_{w,\textit{DMT},\varSigma }$ with $W_{\textit{ij}}=\delta _{\textit{ij}}$ , respectively. The Eulerian time derivative is abbreviated to  $\partial _t$ . Because $\boldsymbol{u}$ , as well as $\boldsymbol{l}$ and $\boldsymbol{a}$ , is zero on the cylinder surface, the body-surface term equals zero. Time averaging can also be conducted on (4.3), which eliminates the volume integral of $-\rho \boldsymbol{x}\times \partial \boldsymbol{\omega }/\partial t$ .

Drag decomposition based on (4.3) is shown in figure 7(a,b). The maximum relative difference between the control-volume integral (black dash-dotted line) and the surface-stress integral (black star symbol) is $0.1\,\%$ for $x_{w\textit{ake}}\le 26d$ , $1\,\%$ for $26d\lt x_{w\textit{ake}}\le 111d$ and $4\,\%$ for $111d\lt x_{w\textit{ake}}\le 150d$ . Interestingly, neither the vortex force nor the volume integral of $-\rho \boldsymbol{x}\times \partial \boldsymbol{\omega }/\partial t$ dominates the drag. The vortex force even contributes a small negative drag after experiencing dramatic oscillations in the range of $x_{w\textit{ake}}\lt 60d$ . Instead, the wake-plane integral of $-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l})$ is much closer to the total drag, although it also exhibits oscillations for $x_{w\textit{ake}}\lt 60d$ . Throughout the $x_{w\textit{ake}}$ range examined, $\boldsymbol{F}_{\textit{DMT},\varSigma }$ is negligibly small. The time-averaged curves (figure 7 b) show that the vortex force first has a large negative value and increases rapidly to a slightly positive value, before slowly decreasing to a small negative plateau value. The time-averaged wake-plane integral of $-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l})$ exhibits a contrary trend to the time-averaged vortex force.

Figure 7.

Force decomposition based on the advection-form DMT method. (a) Instantaneous drag at the instant of maximum lift. (b) Time-averaged drag. (c,d) Instantaneous lift at the instant of maximum lift.

Lift decomposition based on (4.3) is shown in figure 7(c,d). Due to the large value of $x$ , the volume integral of $-\rho (\boldsymbol{x}\times \partial \boldsymbol{\omega }/\partial t)\boldsymbol{\cdot }\boldsymbol{e}_y=\rho x\partial \omega _z/\partial t$ and the wake-plane integral of $-\rho [\boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}) ]\boldsymbol{\cdot }\boldsymbol{e}_y=\rho x l_y$ exhibit violent oscillations, with an increasing amplitude for $x_{w\textit{ake}}\lt 30d$ and a decreasing amplitude for $x_{w\textit{ake}}\gt 30d$ . Significant oscillations also exist in the vortex-force term and the $\boldsymbol{F}_{\textit{DMT},\varSigma }$ term. The relative error of the control-volume integral (4.3), compared with the surface-stress integral, is large at far downstream. In the range of $x_{w\textit{ake}}\le 3.8d$ , this relative error is less than $1\,\%$ ; it quickly grows to $5.3\,\%$ at $x_{w\textit{ake}}=9.5d$ and to $10\,\%$ at $x_{w\textit{ake}}= 35d$ . At $x_{w\textit{ake}}=98d$ , this relative error can be as high as $60\,\%$ , indicating deterioration in accuracy due to the large control volume.

To overcome the issue caused by the large $x$ value in (4.3), the divergence- and curl-free auxiliary function introduced in the work of Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983) is adopted but with a negative sign. Here, we set $ \unicode{x1D652}=\boldsymbol{\nabla }\boldsymbol{\phi } = \boldsymbol{\nabla }\phi _x\boldsymbol{e}_x+\boldsymbol{\nabla }\phi _y\boldsymbol{e}_y$ . In the fluid domain, $\boldsymbol{\phi }$ is governed by the Laplacian equation,

(4.4) \begin{align} {\nabla} ^2\boldsymbol{\phi } = {\boldsymbol{0}}. \end{align}

On the cylinder surface, $\boldsymbol{\phi }$ satisfies the Neumann boundary condition

(4.5) \begin{align} \hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\phi } =\hat {\boldsymbol{n}}, \end{align}

and $\boldsymbol{\phi }$ goes to zero at the infinitely far field. Here $\boldsymbol{\nabla }\phi _x$ or $\boldsymbol{\nabla }\phi _y$ can be seen as the potential flow generated by the motion of the body in the $x$ or $y$ direction with a unit speed. It is noted that $\boldsymbol{\nabla }\boldsymbol{\phi }$ is dimensionless and that $\boldsymbol{\phi }$ has a dimension of length. At the far field, $\boldsymbol{\nabla }\boldsymbol{\phi }$ decreases to zero at a rate of $|\boldsymbol{x}|^{-2}$ , which effectively reduces the magnitudes of integrands. The advection-form weighted DMT method is written as

(4.6) \begin{align} \begin{split} \boldsymbol{F} &= \underbrace {\int _{V_{\kern-1pt f}}{\rho (\boldsymbol{u}\times \boldsymbol{\omega })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\phi } \,\textrm{d} V}}_{\text{Weighted vortex force}} +\underbrace {\int _{V_{\kern-1pt f}}{-\rho \boldsymbol{x}\times (\boldsymbol{\nabla }\times \boldsymbol{b}_w) \,\textrm{d} V}}_{\boldsymbol{b}_w\text{ term}} \\& \qquad +\underbrace {\oint _{\partial\kern-1pt B}{-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}_w) \,\textrm{d} S} + \boldsymbol{F}_{w,\textit{DMT},B} + \big(\boldsymbol{F}_\textit{f} - \boldsymbol{F}_{\!w,\textit{f}}\big)}_{\text{Body-surface term}} \\& \qquad + \underbrace {\int _{w\textit{ake}}{-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}_w) \,\textrm{d} S}}_{\text{Wake-plane integral of $-\rho \boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}_w)$}} + \boldsymbol{F}_{w,\textit{DMT},\varSigma }. \end{split} \end{align}

The force expression consists of the volume integral of the weighted vortex force, the volume integral of $-\rho \boldsymbol{x}\times (\boldsymbol{\nabla }\times \boldsymbol{b}_w)$ , the body-surface term, the wake-plane integral of $-\boldsymbol{x}\times (\boldsymbol{n}\times \boldsymbol{l}_w)$ and $\boldsymbol{F}_{w,\textit{DMT},\varSigma }$ . It is noted that the weighted vortex force is exactly the volume force element in the work of Chang (Reference Chang1992) and Chiu et al. (Reference Chiu, Tseng, Chang and Chou2023). The $\boldsymbol{b}_w$ term contains not only the unsteadiness effect but also spatial derivatives of the weight function; as a result, the time averaging of $-\rho \boldsymbol{x}\times (\boldsymbol{\nabla }\times \boldsymbol{b}_w)$ is not zero. Here $\boldsymbol{f}_{\!w}$ is not zero on the cylinder surface, either, and hence, $\boldsymbol{F}_{w,\textit{DMT},B}$ also contributes to the force. Additionally, although $\boldsymbol{F}_{w,\textit{p}}$ equals $\boldsymbol{F}_\textit{p}$ , $\boldsymbol{F}_{\!w,\textit{f}}$ does not equal $\boldsymbol{F}_\textit{f}$ . Therefore, a term $\boldsymbol{F}_\textit{f}-\boldsymbol{F}_{\!w,\textit{f}}$ is introduced in the body-surface term.

Figure 8(a,b) shows the force decomposition at the maximum-lift instant based on the advection-form weighted DMT method (4.6). The curves converge to constant values quickly with the increasing $x_{w\textit{ake}}$ . Specifically, the drag components (figure 8 a) remain almost invariant for $x_{w\textit{ake}}\gt 10d$ . The lift components (figure 8 b), which contain $x$ , show damped oscillations and reach almost constant values for $x_{w\textit{ake}}\gt 70d$ . The wake-plane integrals (red line with circle and green line) both decay to zero. The numerical accuracy of the control-volume integral (black dash-dotted line) is also significantly improved through the introduction of the weight function. The maximum relative error in the instantaneous drag is $0.28\,\%$ for $x_{w\textit{ake}}\lt 5d$ and drops below $0.07\,\%$ for $5d\le x_{w\textit{ake}}\le 150d$ . The maximum relative error in the lift is slightly larger, which is $2.6\,\%$ for $x_{w\textit{ake}}\le 28d$ and $0.7\,\%$ for $28d\lt x_{w\textit{ake}}\le 150d$ .

Figure 8.

Force decomposition based on the advection-form weighted DMT method. (a) Instantaneous drag and (b) instantaneous lift at the instant of maximum lift. (c) Time-averaged drag.

In the weighted DMT method (4.6), the volume integral of the weighted vortex force and the body-surface term become the main contributions to both the drag and the lift. The volume integral of $-\rho \boldsymbol{x}\times (\boldsymbol{\nabla }\times \boldsymbol{b}_w)$ , on the other hand, has a strong negative force contribution. When the final constant values are reached, the body-surface term equals $1.53F_x$ and $0.73F_y$ , the volume integral of the weighted vortex force equals $0.49F_x$ and $0.76F_y$ , and the volume integral of $-\rho \boldsymbol{x}\times (\boldsymbol{\nabla }\times \boldsymbol{b}_w)$ equals $-1.02F_x$ and $-0.49F_y$ , at the maximum-lift instant. It is also noted that the instantaneous drag components are almost the same as the time-averaged drag components (figure 8 c), which is consistent with the small oscillation amplitude of the time-dependent drag.

The spatial distribution of the vortex force is shown in figure 9(a,b) to gain insights into its integrated behaviour presented in figure 7. The vorticity field is also shown by black contour lines with positive values indicated by solid lines and negative values in dashed lines. The three most recently formed vortices are labelled as V0, V1 and V2 in figure 9(b). By definition, the vortex force is distributed in high-vorticity regions, such as the boundary layers, free shear layers and wake vortices. Within each concentrated vortex except for V0, the $x$ component of the vortex force $u_y\omega _z$ is negative on the left side and positive on the right side, owing to the opposite signs of the vertical velocity across the vortex centre. In contrast, the $y$ component of the vortex force $-u_x\omega _z$ exhibits a distribution similar to that of the negative vorticity field, because $u_x$ is generally positive outside the recirculation region. Due to the slow decay of concentrated vortices in the wake, the volume integral of vortex force inevitably oscillates when the wake plane passes through vortices, and the oscillation frequency of the $x$ component is twice that of the $y$ component, as indicated by the black lines with inverted triangles in figure 7(a,c).

Figure 9.

Colour shows the distribution of (a,b) vortex force and (c,d) weighted vortex force at the instant of maximum lift. Black contour lines represent dimensionless vorticity, with negative values (–0.1 and –0.03) as dashed lines and positive values (0.03 and 0.1) as solid lines.

The spatial distribution of the weighted vortex force is shown in figure 9(c,d). In contrast to the original vortex force, the weighted vortex force is mainly distributed in the boundary layers and the free shear layers within the range $x\lt 3d$ . Its magnitude is small even in the concentrated vortex V2. Because $\phi _x=-0.25 d^2 x/(x^2+y^2)$ changes sign across the $y$ axis, the $x$ component of the weighted vortex force is negative in a small area in the boundary layer of the left half-cylinder, whereas it becomes positive in larger areas in the shear layers and vortices with $x\gt 0$ ; as a result, its volume integral increases monotonically to a positive value of $0.32\rho U_\infty ^2 d$ in the range $0.8d\le x_{w\textit{ake}}\lt 5d$ (see the black line with inverted triangles in figure 8 a). Similarly, $\phi _y=-0.25 d^2 y/(x^2+y^2)$ has different signs across the $x$ axis; therefore, the $y$ component of the weighted vortex force is positive in the boundary layer of the upper surface, while it becomes negative in the boundary layer of the lower surface. Additionally, the $y$ component of the weighted vortex force is negative in vortex V0, whereas it has a relatively strong positive value in vortex V1; as a result, its volume integral grows to a positive value in the range $0.8d\lt x_{w\textit{ake}}\lt 2.9d$ before exhibiting small-amplitude and damped oscillations for $x_{w\textit{ake}}\gt 2.9d$ ; see the black line with inverted triangles in figure 8(b).

In conclusion, the advection-form weighted DMT method (4.6) recovers the volume force element proposed by Chang (Reference Chang1992) when the negative gradient of the auxiliary potential introduced by Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983) is employed as the weight function. The introducing of a rapidly decaying weight function offers advantages in improving the numerical accuracy and the control-volume convergence of the DMT method.

4.4. Pressure distribution over the cylinder surface

A quantitative relationship between the pressure distribution over the cylinder surface and flow structures is established and evaluated in this section. A series of harmonic functions that satisfy boundary conditions

(4.7) \begin{align} \hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{cm} = \cos (m\theta ) ,\quad \hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{\textit{sm}} = \sin (m\theta ), \quad m=0,1,2,\ldots \end{align}

on the cylinder surface, with $\theta$ being the polar angle, is considered. The corresponding weighted pressure forces are

(4.8a) \begin{align} F_{cm,p}(t) &:= \oint _{\partial\kern-1pt B}{ -p\hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{cm}\,\textrm{d} S} = -0.5 d \int _0^{2\pi }{p\cos (m\theta )\,\textrm{d} \theta }, \\[-12pt] \nonumber \end{align}

(4.8b) \begin{align} F_{sm,p}(t) &:= \oint _{\partial\kern-1pt B}{ -p\hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{\textit{sm}}\,\textrm{d} S} = -0.5 d \int _0^{2\pi }{p\sin (m\theta )\,\textrm{d} \theta }. \\[9pt] \nonumber\end{align}

The $F_{cm,p}$ and $F_{sm,p}$ , which are termed pressure distribution coefficients, are proportional to, but with opposite signs to, the Fourier expansion coefficients of the surface pressure. By definition, $F_{s0,p}=0$ . It is also noted that the total pressure force is $\boldsymbol{F}_\textit{p}=F_{c1,p}\boldsymbol{e}_x+ F_{s1,p}\boldsymbol{e}_y$ . The Fourier expansion of the surface pressure is given by

(4.9) \begin{align} p_N(t, \theta ) = -\frac {1}{\pi d}F_{c0,p}(t)-\frac {2}{\pi d}\sum _{m=1}^{N}\left [F_{cm,p}(t)\cos (m\theta )+F_{sm,p}\sin (m\theta )\right ]\!, \end{align}

which recovers the surface pressure distribution when the expansion order $N$ goes to infinity,

(4.10) \begin{align} p(t, \theta ) = \lim _{N\to +\infty }{p_N(t, \theta )}. \end{align}

Since ${\nabla} ^2\varphi _{cm}=0$ and ${\nabla} ^2\varphi _{\textit{sm}}=0$ , $F_{cm,p}$ and $F_{sm,p}$ can be expressed as control-volume integrals through the weighted projection method with curl-free weight functions (3.12). These are given by

(4.11a) \begin{align} F_{cm,p} & = \underbrace {\int _{V_{f}}{-2\rho Q \varphi _{cm} \,\textrm{d} V}}_{\text{Weighted integral of $Q$}} + \underbrace {\oint _{\partial\kern-1pt B}{\boldsymbol{n}\boldsymbol{\cdot }(-\rho \boldsymbol{a}+\mu {\nabla} ^2\boldsymbol{u})\varphi _{cm} \,\textrm{d} S}}_{\text{Body-surface term}} \nonumber \\& \qquad + \underbrace {\int _{w\textit{ake}}{\boldsymbol{n}\boldsymbol{\cdot }(-\rho \boldsymbol{a}+\mu {\nabla} ^2\boldsymbol{u})\varphi _{cm} \,\textrm{d} S} + \int _{w\textit{ake}}{-p\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{cm} \,\textrm{d} S}}_{\text{Wake-plane integral}}, \\[-12pt] \nonumber \end{align}

(4.11b) \begin{align} F_{sm,p} & = \int _{V_{f}}{-2\rho Q \varphi _{\textit{sm}} \,\textrm{d} V} + \oint _{\partial\kern-1pt B}{\boldsymbol{n}\boldsymbol{\cdot }(-\rho \boldsymbol{a}+\mu {\nabla} ^2\boldsymbol{u})\varphi _{\textit{sm}} \,\textrm{d} S} \nonumber \\ & \qquad + \int _{w\textit{ake}}{\boldsymbol{n}\boldsymbol{\cdot }(-\rho \boldsymbol{a}+\mu {\nabla} ^2\boldsymbol{u})\varphi _{\textit{sm}} \,\textrm{d} S} + \int _{w\textit{ake}}{-p\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{\textit{sm}} \,\textrm{d} S}. \\[6pt] \nonumber \end{align}

Therefore, (4.9) and (4.11) establish a quantitative relationship between the $Q$ field and the surface pressure distribution.

An alternative interpretation of (4.11a ), as well as (4.11b ), is given in Appendix I. The pressure field can be decomposed into four components, $p=p_Q+p_a+p_{\textit{vis}}+p_\infty$ . Here, $p_Q$ represents the partial pressure field generated by the volume pressure source $2\rho Q$ , $p_a$ is generated by the body-surface normal acceleration, $p_{\textit{vis}}$ is generated by the viscosity effect on $\partial\! B$ and $p_\infty$ is the pressure at infinity. Correspondingly, $F_{cm,p}=F_{cm,p_Q}+F_{cm,p_a}+F_{cm,p_{\textit{vis}}}+F_{cm,p_\infty }$ . Auxiliary functions $\varphi _{cm}$ and $\varphi _{\textit{sm}}$ are the solutions to the adjoint equation of the pressure Poisson equation. Therefore, $F_{cm,p_Q}$ , the body-surface integral of $-p_Q \cos (m\theta )$ , can be calculated by the volume integration of $-2\rho Q\varphi _{cm}$ over the whole fluid domain. Here, the function $\varphi _{cm}$ represents the sensitivity of $F_{cm,p_Q}$ to the volume pressure source. A vortical structure with $Q\gt 0$ produces a low-pressure region through ‘vortex suction’ on the nearby body surface, whereas the strain rate region with $Q\lt 0$ generates a repulsive force and a high-pressure area on the neighbouring solid surface. Similarly, $F_{cm,p_a}$ and $F_{cm, p_{\textit{vis}}}$ can be calculated by the surface integration of $-\rho \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{a}\varphi _{cm}$ and $\mu \boldsymbol{n}\boldsymbol{\cdot }{\nabla} ^2\boldsymbol{u}\varphi _{cm}$ on $\partial\! B$ , respectively. The outer surface integral of $-p\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{c0}$ converges to $\pi dp_\infty$ , when $\varSigma$ approaches infinity.

For the circular cylinder, the auxiliary functions $\varphi _{cm}$ and $\varphi _{\textit{sm}}$ can be solved analytically as

(4.12a) \begin{align} \varphi _{c0} &= 0.5 d\ln (2r/d),\quad \varphi _{s0} = 0, \\[-12pt] \nonumber \end{align}

(4.12b) \begin{align} \varphi _{cm} &= -m^{-1} (d/2)^{m+1} r^{-m}\cos (m\theta ),\quad m=1,2,3,\ldots, \\[-12pt] \nonumber \end{align}

(4.12c) \begin{align} \varphi _{\textit{sm}}& = -m^{-1} (d/2)^{m+1} r^{-m} \sin (m\theta ),\quad m=1,2,3,\ldots . \\[9pt] \nonumber \end{align}

Here, $r$ and $\theta$ are the polar coordinates, with $x=r\cos \theta$ and $y=r\sin \theta$ . For $m=0$ , the net flux of $\boldsymbol{\nabla }\varphi _{c0}$ from the cylinder surface is non-zero, and therefore, $\varphi _{c0}$ , which can be seen as the velocity field induced by a point source of strength $\pi d$ in the 2-D potential flow, grows as $O(\ln r)$ at the far field. Here, $\varphi _{c0}$ is set to zero on the cylinder surface. For $m\ge 1$ , $\varphi _{cm}$ and $\varphi _{\textit{sm}}$ decay as $O(r^{-m})$ when $r$ goes to infinity.

Figure 10 shows the pressure distribution over the cylinder surface at two instants. Here, $\theta =\pm \pi$ is the leading edge and $\theta =0$ is the trailing edge of the cylinder. Meanwhile, $\theta =\pi /2$ is the highest point and $\theta =-\pi /2$ is the lowest point on the cylinder surface. The leading-edge stagnation point has a higher pressure than the leeward side of the cylinder, while two local pressure minima occur slightly upstream of the highest and the lowest points. At the maximum-lift instant, the upper-right quarter ( $0\lt \theta \lt 0.5\pi$ ) has a lower pressure than the lower-right quarter ( $-0.5\pi \lt \theta \lt 0$ ); whereas at the zero-lift instant, the pressure distribution is more symmetric between the upper and the lower surfaces.

Figure 10.

Pressure distribution on the cylinder surface and its Fourier expansion $p_N$ with orders $N=0,1,2,3$ . (a) At the instant of maximum lift. (b) At the instant of zero lift.

The truncated Fourier expansion (4.9) is also shown in figure 10, with the dimensionless $F_{cm,p}$ and $F_{sm,p}$ at the maximum-lift instant presented in table 2. With the increasing expansion order $N$ , the Fourier expansion $p_N$ gradually converges to the surface pressure  $p$ . The zeroth-order Fourier expansion gives the surface-averaged pressure, which is a constant lower than $p_\infty = 0$ as indicated by the green solid line. However, due to the minus sign in (4.8a ), $F_{c0,p}$ has a positive value of $0.76\rho U_\infty ^2 d$ . The first-order expansion, as indicated by the blue-dotted line, captures the higher pressure near the leading edge and the lower pressure on the leeward side because of the contribution from $F_{c1,p}$ , i.e. the pressure drag. At the maximum-lift instant (figure 10 a), the line with $N=1$ exhibits asymmetry about $\theta =0$ owing to the significance of the antisymmetric term $F_{s1,p}=0.1395\rho U_\infty ^2 d$ , which equals the pressure lift. Since $F_{s2,p}$ and $F_{s3,p}$ are negligibly small, the second-order coefficient $F_{c2,p}$ with a negative value of $-0.5035\rho U_\infty ^2 d$ captures the even higher pressure at the leading edge and the local pressure minima near the highest and the lowest points; see the red-dashed line. With $N=3$ as indicated by the black dash-dotted line, the third-order coefficient $F_{c3,p}$ further raises the pressure at the leading edge and reduces $p$ near the trailing edge, resulting in a satisfactory approximation to the surface pressure.

Table 2.

The decomposition of $F_{cm,p}/(\rho U_\infty ^2 d)$ and $F_{sm,p}/(\rho U_\infty ^2 d)$ when the integrals are independent of the control-volume size. Here $x_{w\textit{ake}}=150d$ for $m=0$ ; $x_{w\textit{ake}}=36 d$ for $m=1,2,3$ , 4.

The decomposition of dimensionless $F_{cm,p}$ and $F_{sm,p}$ at the maximum-lift instant, based on (4.11), is provided in table 2. A large enough control volume, which is $x_{w\textit{ake}}=150d$ for $m=0$ and $x_{w\textit{ake}}=36d$ for $m\ge 1$ , is adopted to eliminate the $x_{w\textit{ake}}$  dependence. The results calculated by the control-volume integral (4.11) and the cylinder surface integral (4.8) are almost identical. The weighted integral of $Q$ over the control volume is the major contribution to the pressure distribution coefficients. Although the cylinder is stationary and, thus, $\boldsymbol{a}={\boldsymbol{0}}$ , the body-surface term with $m\ge 1$ is not zero due to the viscosity effect in the Neumann pressure boundary condition on the cylinder surface, $\hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla } p = -\rho \hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{a}+\mu {\nabla} ^2\boldsymbol{u}\boldsymbol{\cdot }\hat {\boldsymbol{n}}$ . For prominent coefficients $F_{c1,p}$ and $F_{c2,p}$ , the contribution from the body-surface term is about half that of the weighted integral of $Q$ , representing a non-negligible viscosity effect on the surface pressure distribution.

The dependence of components in $F_{c0,p}$ on the wake-plane location $x_{w\textit{ake}}$ is shown in figure 11(a). The body-surface term is zero since $\varphi _{c0}$ is set to zero on the cylinder surface. Because the total flux of $\boldsymbol{\nabla }\varphi _{c0}$ through $\varSigma$ remains a constant value of $\pi d$ , the outer surface integral of $-p\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _{c0}$ converges to $-\pi d p_\infty$ , with $p_\infty$ set to zero, when  $\varSigma$ goes to infinity. Since $\varphi _{c0}$ grows with $r$ , the decay of $-2\rho Q\varphi _{c0}$ at the far field solely depends on the decrease of $Q$ . As a result, the volume integral of $-2\rho Q\varphi _{c0}$ (blue line with inverted triangles) converges slowly. For $0.8d\le x_{w\textit{ake}}\lt 10d$ , whilst experiencing strong oscillations, this volume integral exhibits an increasing trend. Further downstream, it exhibits damped oscillations before reaching a constant value of $0.76\rho U_\infty ^2 d$ at $x_{w\textit{ake}}=75d$ .

Figure 11.

(a) Decomposition of the zeroth-order pressure distribution coefficient, $F_{c0,p}$ . (b) Spatial distribution of $Q$ at the instant of maximum lift. Black contour lines represent dimensionless vorticity, with negative values (−0.1 and −0.03) as dashed lines and positive values (0.03 and 0.1) as solid lines.

The spatial distribution of $Q$ is shown in figure 11(b) to gain further insights. Due to the singularity in 2-D space, a constant in $\varphi _{c0}$ is undetermined. Although $\varphi _{c0}$ has been artificially set to zero on $\partial\! B$ , the volume integrand $-2\rho Q\varphi _{c0}$ cannot objectively reflect the contribution from local fluid domains. When the total volume integration of $2\rho Q$ is zero, which is satisfied in the present flow, Appendix I provides another explanation for the cause of the average surface pressure: if the region with $Q\gt 0$ is distributed closer to the solid body than the region with $Q\lt 0$ , the average surface pressure will become lower than $p_\infty$ . Since $Q$ measures the relative strength between rotation and deformation of fluid elements, $Q$ is positive at the centre of concentrated vortices, such as the most recently formed vortices V0, V1 and V2, and in the rolling-up shear layers; whereas it is negative between two successive vortices and in front of the leading edge. Clearly, the rolling-up shear layers near the cylinder and the vortex centre of V0 show strong positive values of $Q$ , and they are closest to the cylinder. The strain rate region below V0 is weaker than the vortex centre of V0 in magnitude; the strain rate region between V0 and V2 is also weaker than V1. Therefore, the stronger shear layers and vortices together lead to the lower average pressure on the cylinder surface.

Figure 12.

Decomposition of the pressure distribution coefficients, $F_{cm,p}$ and $F_{sm,p}$ with $m=1,2,3$ , based on the weighted projection method with curl-free weight function. Panels (a,c,e) are $F_{cm,p}$ . Panels (b,d,f) are $F_{sm,p}$ .

Figure 13.

The spatial distribution of (a,c,e) $-2\rho Q\varphi _{cm}/(\rho U_\infty ^2 d^{-1})$ and (b,d,f) $-2\rho Q\varphi _{\textit{sm}}/(\rho U_\infty ^2 d^{-1})$ at the instant of maximum lift. Results are shown for (a,b) $m=1$ , (c,d) $m=2$ , (e,f) $m=3$ . Black contour lines represent dimensionless vorticity, with negative values (−0.1 and −0.03) as dashed lines and positive values (0.03 and 0.1) as solid lines.

The dependence of $F_{cm,p}$ and $F_{sm,p}$ with $m=1,2,3$ on the $x_{w\textit{ake}}$ is shown in figure 12. For $m\ge 1$ , the total fluxes of $\boldsymbol{\nabla }\varphi _{cm}$ and $\boldsymbol{\nabla }\varphi _{\textit{sm}}$ over $\varSigma$ are both zero; therefore, the wake-plane integral (red line) decays to zero quickly with the increasing $x_{w\textit{ake}}$ . For $F_{c1,p}$ and $F_{s1,p}$ , the weighted integral of $Q$ (blue line) grows mainly in the range $x_{w\textit{ake}}\le 1.5d$ , a region that just includes the vortex V0, and converges to a steady value for $x_{w\textit{ake}}\ge 3d$ . For $F_{c2,p}$ and $F_{c3,p}$ , their curves remain almost constant for $x_{w\textit{ake}}\ge 0.8d$ , and therefore, they receive main contributions from the $Q$ structures upstream of $x=0.8d$ . Here $F_{s2,p}$ and $F_{s3,p}$ are small in magnitude and not important to the surface pressure.

The spatial distributions of the volume integrands $-2\rho Q\varphi _{cm}$ and $-2\rho Q\varphi _{\textit{sm}}$ , which represent contributions from local flow structures, are shown in figure 13, where the vorticity field is also shown by black contour lines. Because $\varphi _{c1}$ changes sign across $x=0$ , the strain rate in front of the leading edge contributes to drag, i.e. a positive value to $F_{c1,p}$ , whereas vortical structures in the wake also contribute to drag. In contrast, the boundary layer on the left cylinder surface and the strain region below vortex V0 contribute to thrust. Similarly, since $\varphi _{s1}$ changes sign across $y=0$ , the upstream shear layer and the boundary layers on the left half-cylinder show opposite signs in the upper and lower planes, leading to cancellation when evaluating their net contribution to $F_{s1,p}$ . On the other hand, both the vortex centre of V0 and the strain rate region below it contribute to lift, which is the main reason for the increase of the volume integral of $-2\rho Q\varphi _{s1}$ in $0.8d\lt x_{w\textit{ake}}\lt 1.5d$ , as indicated in figure 12(b). Therefore, the formation of V0 and the strain rate region below it is responsible for the asymmetric pressure distribution on the leeward side of the cylinder surface. The $-2\rho Q\varphi _{c2}$ is mainly distributed with negative values in front of the leading edge and near the highest and lowest points. This is consistent with the fact that the effect of $F_{c2,p}$ , which is negative, raises the pressure at the leading edge and reduces $p$ near the highest and lowest points, as illustrated by the Fourier expansion in figure 10(a). In addition, the $-2\rho Q\varphi _{c3}$ is mainly distributed with positive values in front of the leading edge and in the boundary layers on the left half-cylinder. These regions are responsible for the further increased pressure at the leading edge through the additional term $F_{c3,p}$ in the Fourier expansion.

4.5. Force estimation from a discrete and noisy flow field

Estimating the fluid force from the flow field provides an approach for non-intrusive force measurements. However, its accuracy and applicability are impeded by the requirement for near-wall velocity derivatives or the pressure field (Díaz-Arriba et al. Reference Díaz-Arriba, Protas, Acher, Pons and David2022). In the weighted projection method (3.8), these inconveniences can be circumvented by a specially designed weight function.

Consider a divergence-free vector weight function $\boldsymbol{w}$ that satisfies $\boldsymbol{w}=\boldsymbol{w}_B$ on the body surface $\partial\! B$ and $\boldsymbol{w}={\boldsymbol{0}}$ on the outer surface $\varSigma$ . Here, $\boldsymbol{w}_B$ is a constant vector. By substituting $ \unicode{x1D652}\boldsymbol{\cdot }\boldsymbol{e}_x=\boldsymbol{w}$ into (3.8) and using (3.7d ) for $\boldsymbol{F}_{w,{v\textit{is}}}$ , the total force exerted on the body is obtained as

(4.13) \begin{align} \boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{w}_B &= \int _{V_{f}}{-\rho \boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{w} \textrm{d} V} +\int _{V_{\kern-1pt f}}{\mu \boldsymbol{u}\boldsymbol{\cdot }{\nabla} ^2\boldsymbol{w} \,\textrm{d} V} + \oint _{\partial\kern-1pt B}{-\mu \boldsymbol{u}\boldsymbol{\cdot } \big[\boldsymbol{\nabla }\boldsymbol{w}+(\boldsymbol{\nabla }\boldsymbol{w})^T \big]\boldsymbol{\cdot }\boldsymbol{n}\,\textrm{d} S} \nonumber \\ & \quad + \oint _{\varSigma }{-\mu \boldsymbol{u}\boldsymbol{\cdot } \big[\boldsymbol{\nabla }\boldsymbol{w}+(\boldsymbol{\nabla }\boldsymbol{w})^T \big]\boldsymbol{\cdot }\boldsymbol{n} \textrm{d} S}. \end{align}

Although the acceleration is involved in the volume integration, (4.13) does not contain spatial derivatives of the velocity field, nor does it require the pressure field. As demonstrated in the work of Lynch & Scarano (Reference Lynch and Scarano2014) and Jiang et al. (Reference Jiang, Lee, Smith, Chen and Linden2020), both the fluid velocity and the acceleration fields can be accurately acquired through time-resolved PIV techniques; therefore, (4.13) provides a feasible approach for force estimation.

The discretisation of (4.13) can be performed by splitting the control volume and its boundary into volume elements $V_i$ and surface elements $\partial\kern-1pt B_i$ and $\varSigma _i$ , respectively. These are

(4.14) \begin{align} V_{\kern-1pt f} = \bigcup _i V_i, \quad \partial\kern-1pt B = \bigcup _i \partial\kern-1pt B_i, \quad \varSigma = \bigcup _i \varSigma _i. \end{align}

The following coefficients can be introduced as

(4.15) \begin{align} \boldsymbol{f}_\textit{V,i} &= \int _{V_i}{-\rho \boldsymbol{w}\textrm{d} V}, \quad \boldsymbol{g}_\textit{V,i} = \int _{V_i}{\mu {\nabla} ^2\boldsymbol{w}\textrm{d} V},\quad \boldsymbol{\gamma}_\textit{B,i} = \int _{\partial\kern-1pt B_i}{-\mu \big[\boldsymbol{\nabla }\boldsymbol{w}+(\boldsymbol{\nabla }\boldsymbol{w})^T \big]\boldsymbol{\cdot }\boldsymbol{n}\textrm{d} S}, \nonumber \\ \boldsymbol{\gamma }_{\varSigma ,i} &= \int _{\varSigma _i}{-\mu \big[\boldsymbol{\nabla }\boldsymbol{w}+(\boldsymbol{\nabla }\boldsymbol{w})^T \big]\boldsymbol{\cdot }\boldsymbol{n}\textrm{d} S}. \end{align}

It is noted that the coefficients $\boldsymbol{f}_{\!\textit{V},\textit{i}}$ , $\boldsymbol{g}_{\!\textit{V},\textit{i}}$ , $\boldsymbol{\gamma}_{\!\textit{B},\textit{i}}$ and $\boldsymbol{\gamma}_{\varSigma ,i}$ are fixed once the weight function $\boldsymbol{w}$ and domain partition are provided. Assuming the $V_i$ , $\partial\kern-1pt B_i$ and $\varSigma _i$ are small enough that the integrands in (4.15) do not change sign in the corresponding integral domains, from the generalised mean-value theorem for integrals, the following holds:

(4.16) \begin{align} \boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{w}_B = \sum _i \bar {\boldsymbol{a}}_i\boldsymbol{\cdot }\boldsymbol{f}_{\!\textit{V,i}} + \sum _i \bar {\boldsymbol{u}}_i\boldsymbol{\cdot }\boldsymbol{g}_{\!\textit{V,i}} + \sum _i \tilde {\boldsymbol{u}}_i\boldsymbol{\cdot }\boldsymbol{\gamma}_\textit{B,i} + \sum _i \tilde {\boldsymbol{u}}_i\boldsymbol{\cdot }\boldsymbol{\gamma }_{\!{\varSigma ,i}}. \end{align}

Here, $\overline {(\;)}_i$ denotes a special value in the volume domain $V_i$ and $\widetilde{(\;)}_i$ denotes a special value on the surface domain $\partial B_i$ or $\varSigma _i$ . In applications, these special values can be replaced by values at an arbitrary point in each element, resulting in a first-order numerical quadrature method. Additionally, if the weight function $\boldsymbol{w}$ , as well as its first- and second-order derivatives, is smooth, the sign-preserving assumption can be satisfied by most elements. Violation of this assumption by a small group of elements leads to a first-order numerical error at worst. Therefore, the weighted summation method (4.16) enables force estimation with first-order spatial accuracy.

Equation (4.16) is tested using the flow past a circular cylinder at $Re=100$ . Although the control volume $V_{\kern-1pt f}$ can be set as an arbitrary domain that encloses the cylinder, a medium-sized domain ensures a better spatial resolution in practical applications. Here, $V_{\kern-1pt f}$ , which is also the domain of definition of $\boldsymbol{w}$ , is taken as $-3d\le x\le 3d$ , $-3d\le y\le 3d$ and $x^2+y^2\ge 0.25d^2$ ; see figure 14. The control volume $V_{\kern-1pt f}$ is divided into $N\times N$ uniform elements. The sample point is placed at the centre of each volume or surface element, as indicated by the black dots. For sample points lying inside the cylinder, their velocity and acceleration are taken as zero. The weight function $\boldsymbol{w}$ , governed by the divergence-free condition (3.1) and the N-S equations (3.2), is obtained numerically using the incompressible flow solver in Nektar++ (Cantwell et al. Reference Cantwell2015), with boundary conditions $\boldsymbol{w}={\boldsymbol{0}}$ on $\varSigma$ and $\boldsymbol{w}=\boldsymbol{w}_B$ on $\partial\! B$ , an initial condition $\boldsymbol{w}={\boldsymbol{0}}$ and a total evolution time of $d/|\boldsymbol{w}_B|$ . By setting $\boldsymbol{w}_B=\boldsymbol{e}_x$ or $\boldsymbol{w}_B=\boldsymbol{e}_y$ , the drag and lift can be obtained, respectively. Figure 15 shows the distribution of $\boldsymbol{w}$ with $\boldsymbol{w}_B=\boldsymbol{e}_x$ in the control volume. Since $\boldsymbol{w}$ is generated by a unit-speed blowing on the right half-cylinder and a suction on the left half-cylinder within a time duration of $d/|\boldsymbol{w}_B|$ , its streamlines reveal two recirculation regions and its magnitude is mainly distributed within the range of $\sqrt {x^2+y^2}\le 1.5d$ .

Figure 14.

(a) Computational domain and boundary conditions for calculating the vector weight function $\boldsymbol{w}$ . (b) Grid lines in a grey colour show the domain partition and black dots show the locations of velocity and acceleration sample points. The black circle shows the boundary of the circular cylinder.

Figure 15.

Colour shows the strength of the vector weight function $\boldsymbol{w}$ with $\boldsymbol{w}_B=\boldsymbol{e}_x$ ; (a) $\boldsymbol{w}\boldsymbol{\cdot }\boldsymbol{e}_x$ , (b) $\boldsymbol{w}\boldsymbol{\cdot }\boldsymbol{e}_y$ . Black lines with arrows show the streamlines of $\boldsymbol{w}$ .

Figure 16.

Dimensionless force calculated by the weighted summation method (4.16) with different spatial resolutions and noise levels. Here $\sigma$ is the dimensionless standard deviation of the white noise. (a) Drag. (b) Lift.

To mimic the data measured from experiments, a white noise with a dimensionless standard deviation of $\sigma$ is added to the sampled velocity and acceleration data. Forces calculated by the weighted summation method (4.16) with different grid resolutions and $\sigma$ are shown in figure 16. The maximum relative errors are also provided in table 3. Here, the relative error is defined based on the maximum force magnitude obtained from the surface-stress integral. At a spatial resolution of $128\times 128$ , the estimated force shows good agreement with the surface-stress integral. The maximum relative error is 1.5 % for the drag and 0.8 % for the lift. At a higher resolution of $256\times 256$ , the relative errors in the drag and lift are both halved, in agreement with the first-order spatial accuracy of (4.16). Although the white noise reduces the accuracy of the estimated force, the maximum relative error remains well below 5 % even with a dimensionless standard deviation of $\sigma =5\,\%$ . However, if the spatial resolution is insufficient, such as $64\times 64$ , significant errors of $11.6\,\%$ in the drag and $6.5\,\%$ in the lift are still encountered.

Table 3.

The maximum relative error of the force calculated using the weighted summation method (4.16).

It should be noted that there are still considerable degrees of freedom in the choice of control volume and the weight function $\boldsymbol{w}$ . The optimisation of the integration domain, element partition and the weight function $\boldsymbol{w}$ helps reduce the numerical error. Additionally, a higher-order discretisation of (4.13) may also improve the numerical accuracy. These explorations will be addressed in the future.