Decomposition of the mean friction drag on a NACA4412 airfoil under uniform
blowing/suction

Yitong Fan1, Marco Atzori2, Ricardo Vinuesa2, Davide Gatti3, Philipp
Schlatter2, and Weipeng Li1***Corresponding author<EMAIL_ADDRESS>
1 School of Aeronautics and Astronautics, Shanghai Jiao Tong University,
Shanghai, China, 200240

2 SimEx/FLOW, Engineering Mechanics, KTH Royal Institute of Technology,
Stockholm, Sweden, SE-100 44

3 Institute of Fluid Mechanics, Karlsruhe Institute of Technology (KIT),
Karlsruhe, Germany, D-76131

###### Abstract

The application of drag-control strategies on canonical wall-bounded
turbulence, such as periodic channel and zero- or adverse-pressure-gradient
boundary layers, raises the question of how to describe control effects
consistently for different reference cases. We employ the RD identity (Renard
& Deck, J. Fluid Mech., 790, 2016, pp. 339–367) to decompose the mean friction
drag and investigate the control effects of uniform blowing and suction
applied to a NACA4412 airfoil at chord Reynolds numbers $Re_{c}$=$200,000$ and
$400,000$. The connection of the drag reduction/increase by using
blowing/suction with the turbulence statistics (including viscous dissipation,
turbulence-kinetic-energy production, and spatial growth of the flow) across
the boundary layer, subjected to adverse or favorable pressure gradients, are
examined. We found that the peaks of the statistics associated with the
friction-drag generation exhibit good scaling in either inner or outer units
throughout the boundary layer. They are also independent of the Reynolds
number, control scheme, and intensity of the blowing/suction. The small- and
large-scale structures are separated with an adaptive scale-decomposition
method, i.e. empirical mode decomposition (EMD), aiming to analyze the scale-
specific contribution of turbulent motions to friction-drag generation.
Results unveil that blowing on the suction side of the airfoil is able to
enhance the contribution of large-scale motions and to suppress that of small-
scales; on the other hand, suction behaves contrarily. The contributions
related to cross-scale interactions remain almost unchanged with different
control strategies.

Key words: turbulent flows, turbulent boundary layers, turbulence control

##  1 Introduction

The friction drag (mostly associated with turbulent boundary layers) accounts
for approximately $50\%$ of the total aerodynamic drag in long-range
commercial aircraft (Gad-el Hak, 1994). It provides great potential for drag
reduction and energy saving. Among diverse control strategies for turbulent
boundary layers, such as addition of long-chain polymers, oscillating walls,
superhydrophobic surfaces, and riblets (White and Mungal, 2008; Touber and
Leschziner, 2012; Rastegari and Akhavan, 2015; Li, 2020; Ran et al., 2021),
mass blowing and suction is a promising method to control the friction drag or
flow transition in wall-bounded turbulence (Kim et al., 2002; Kametani and
Fukagata, 2011).

Experiments have shown that uniform blowing from smooth perforated surfaces
can reduce the turbulent friction drag with net-energy saving (Hwang, 1996,
2004). Given the proper blowing intensity, porosity, and effective roughness,
the net-energy saving holds for a wide range of conditions, including both
zero- and adverse-pressure-gradient (ZPG and APG) turbulent boundary layers
(TBLs) (Welch et al., 2001). On the other hand, uniform suction will increase
the turbulent friction drag but can be employed for separation control. Direct
numerical simulations (DNS) and large-eddy simulations (LES) of uniform
blowing and suction applied to turbulent boundary layers have also been
performed (Park and Choi, 1999; Kim et al., 2002; Kametani and Fukagata, 2011;
Kametani et al., 2015; Bobke et al., 2016). In most cases, these simulations
were carried out in zero-pressure-gradient conditions to investigate the
effects of blowing and suction on dynamics of wall-bounded turbulence. For
instance, Stroh et al. (2016) compared uniform blowing and suction with body-
force damping (as a model of opposition control) and introduced the concept of
virtual origin to describe blowing and suction effects in the downstream of
the control region.

Recently, researchers have taken interests in the blowing and suction control
of TBLs on wing sections. Using mass suction at the leading edge of a Clark-Z
airfoil to provide pressurized air for blowing, Kornilov (2017) studied
uniform blowing on the pressure side of the airfoil at Reynolds number
$Re_{c}$=$U_{\infty}c/\nu$=$840,000$, where $U_{\infty}$ is the incoming flow
velocity, $c$ is the chord length, and $\nu$ is the fluid kinematic viscosity.
Eto et al. (2019) studied the effects of both active and passive blowing on
the suction side of a Clark-Y airfoil at $Re_{c}$=$1,500,000$. Kornilov et al.
(2019) employed blowing on the pressure side and suction on the suction side
of a NACA0012 airfoil, and later they provided an estimation of the control
energy cost under the same conditions (Kornilov, 2021). Mahfoze et al. (2019)
used Bayesian optimization to discuss how to benefit from downstream effects
of blowing when the control region is separated into individual areas. The
first high-fidelity numerical simulation of a wing section with uniform
blowing was reported by Vinuesa and Schlatter (2017), although at a low
Reynolds number ($Re_{c}$=$100,000$). Soon after, Atzori et al. (2020)
presented a dataset of highly-resolved LES of a NACA4412 airfoil at
$Re_{c}$=$200,000$ and $400,000$ with various configurations of uniform
blowing and suction, using as a reference the simulation carried out by
Vinuesa et al. (2018). This dataset has been later employed by Fahland et al.
(2021) to validate Reynolds-Averaged Navier–Stokes (RANS) simulations, and it
is also considered in the present paper.

The key objective of this study is to investigate the control effects on mean
friction drag on a wing section with uniform blowing and suction. Although the
mean friction drag is a wall property, as can be directly calculated from the
normal gradient of the mean tangential velocity at the wall, it is connected
to the statistical turbulence quantities across the wall layer and can be
further decomposed into various physics-informed components according to
different mathematical derivations and physical interpretations (Li et al.,
2019; Fan et al., 2019a, b). So far, there have been three kinds of friction-
drag decomposition methods, derived from the momentum, vorticity, and energy
balance, respectively. The first one is the so-called Fukagata-Iwamoto-Kasagi
(FIK) identity (Fukagata et al., 2002), in which a triple integration is
performed on the mean momentum balance equation and gives a direct
relationship between the skin-friction coefficient and the Reynolds-shear-
stress profile. The FIK identity has been widely used and extended for more
complex situations over the years, e.g. Mehdi and White (2011); Mehdi et al.
(2014); Modesti et al. (2018); Peet and Sagaut (2009); Bannier et al. (2015),
to name a few. Kametani et al. (2015) and Stroh et al. (2015) applied FIK
identity to quantify the variation of skin-friction coefficients caused by
blowing and suction in ZPG-TBLs. Inspired by the mathematical derivation of
FIK identity, Yoon et al. (2016) derived a vorticity-based formula relating
the mean friction-drag generation with the motion of vortical structures, by
performing a triple integration on the mean spanwise vorticity transport
equation. They later used this method to analyze the contribution of outer
large-scale motions to the friction-drag generation in a moderate APG-TBL
(Yoon et al., 2018). Finally, an energy-based decomposition method was
proposed by Renard and Deck (2016) from the perspective of streamwise kinetic
energy balance. Under an absolute reference frame where the wall is moving,
the friction drag develops a non-zero power, which is characterized as the
energy transfered from the wall to the fluid, by means of molecular viscosity
dissipation, turbulence-kinetic-energy production, and spatial growth of the
flow. This method is referred to as RD identity hereafter. The RD identity has
been used to analyze the friction-drag generation in channel flows, ZPG/APG-
TBLs, turbulent square-duct flows, and pipe flows (Fan et al., 2019a, b,
2020a, 2020b; Wei, 2018). Li et al. (2019) and Fan et al. (2019b) generalized
the RD identity to a compressible form to quantify the compressibility effects
on the friction-drag generation. All these three methods (Fukagata et al.,
2002; Yoon et al., 2016; Renard and Deck, 2016) are mathematically correct and
have been widely validated. In the present study, we only adopt the RD
identity, considering that the momentum- and vorticity-based method (by three
successive integrations) are lack of physics-informed interpretations, thus
their decomposed constituents hardly carry the causal relationship for the
friction-drag generation (Deck et al., 2014; Renard and Deck, 2016; Fan et
al., 2019a).

The turbulent boundary layers on the suction/pressure sides of the wing
section are subjected to adverse/favorable pressure gradients. The pressure
gradients have significant impacts on the scales of coherent structures across
the wall layer. For instance, inner-outer scale separation is more evident in
APG-TBLs than in ZPG-TBLs, even at relatively low Reynolds numbers, due to the
enhancement of outer-scale motions (Tanarro et al., 2020). Additionally, a
significant increment of small-scale energy was found in the outer region, as
the vertical motion induced by the APG transports small scales from the near-
wall region to the outer layer (Vinuesa et al., 2018; Tanarro et al., 2020).
In the present study, we also aim to quantify the contribution of structures
with different scales to the generation of friction drag on the wing section
with/without blowing and suction. Studies on such cases are likely to promote
new drag control strategies. To this end, an appropriate approach to separate
the multi-scale coherent structures is in need. Typically, Fourier analysis
might be a tempting tool to decompose the raw signals into modes with given
wavelengths, yet it relies on a-priori definition of cutoff wavelength and
suffers from inflexibility for complex and transient signals (Cheng et al.,
2019). Another frequently-used method is proper orthogonal decomposition (POD)
(Lumley, 1967; Wu and Christensen, 2010), which sorts the contribution of
velocity fluctuations to the turbulence kinetic energy. However, Wang et al.
(2018, 2019) pointed out that, the energy-ranking spatial modes cannot fully
recover the dynamics of turbulent motions in different length scales. In
contrast, empirical mode decomposition (EMD), proposed by Huang et al. (1998),
provides an adaptive, data-driven, and a-posteriori technique to delineate the
transient and local characteristics of signals. It is in principle free from
pre-established basis functions and represents the original signal as a
superposition of several mono-components and a residual, with the
characteristic wavelengths of the signals automatically determined. With EMD,
Huang et al. (2008) studied the scaling properties and intermittency of
homogeneous turbulence, and Ansell and Balajewicz (2017) analyzed the features
of large-scale vortical structures in a turbulent mixing layer. Agostini and
Leschziner (2014, 2016) used bidimensional empirical mode decomposition (BEMD)
to analyze the modulation of large-scale motions on the small-scale eddies in
the near-wall region, and later they discussed the scale-specific
contributions of large- and small-scale structures to the friction-drag
generation by means of FIK and RD identity (Agostini and Leschziner, 2019) in
channel flows. Dogan et al. (2019) have used EMD to characterize the inner-
outer interaction based on the modulation coefficient. Cheng et al. (2019)
adopted BEMD to identify attached eddies in turbulent channel flows and
quantify their relationship with the friction-drag generation. However, no
relevant study has been found in the open literature to analyze the scale-
specific contribution of turbulent motions to friction-drag generation on the
wing section with/without blowing and suction.

This paper is organized as follows. In Sec.2, we introduce the friction-drag
decomposition method and the database of flow over a NACA4412 airfoil. The
decomposition results on the suction and pressure side of the NACA4412 are
discussed in Sec. 3 and 4, respectively. Concluding remarks are given in Sec.
5.

##  2 Friction-drag decomposition method and the database of flow over the
NACA4412 airfoil

With the energy-based RD identity (Renard and Deck, 2016), the skin-friction
coefficient $C_{f}$ of a turbulent boundary layer can be decomposed as:

$\displaystyle
C_{f}=\underbrace{\frac{2}{U_{e}^{3}}\int_{0}^{\infty}\nu\left(\frac{\partial\left<u\right>}{\partial
y}\right)^{2}{\rm
d}y}_{C_{f,V}}+\underbrace{\frac{2}{U_{e}^{3}}\int_{0}^{\infty}-\left<u^{\prime}v^{\prime}\right>\frac{\partial\left<u\right>}{\partial
y}{\rm d}y}_{C_{f,T}}$
$\displaystyle+\underbrace{\frac{2}{U_{e}^{3}}\int_{0}^{\infty}\left(\left<u\right>-U_{e}\right)\frac{\partial}{\partial
y}\left(\nu\frac{\partial\left<u\right>}{\partial
y}-\left<u^{\prime}v^{\prime}\right>\right){\rm d}y}_{C_{f,G}},$ (1)

where $\left<\cdot\right>$ is the Reynolds averaging operator, the prime ′
denotes fluctuations with respect to the Reynolds averages, $U_{e}$ is the
velocity at the boundary-layer edge $\delta_{99}$, $x$ and $y$ represent the
directions tangential and normal to the wall surface respectively, and $u$ and
$v$ are the corresponding velocity components. The derivation of the RD
identity can be retrieved in Renard and Deck (2016).

Three contributive friction constituents are obtained in equation (2).
Thereinto, (i) $C_{f,V}$ represents the direct molecular viscous dissipation,
transforming the mechanical power into heat; (ii) $C_{f,T}$ represents the
power spent for turbulence-kinetic-energy production; (iii) $C_{f,G}$ accounts
for the spatial growth of the flow, which is also interpreted as the rate of
gain of mean streamwise kinetic energy by the fluid in the absolute frame.
Note that the integrand in $C_{f,G}$ has been substituted with local
information which only depends on the well-documented wall-normal profiles
(Renard and Deck, 2016). This is especially applicable for the cases where the
accurate calculation of explicit streamwise derivatives is unfeasible.

For adverse-/favorable-pressure-gradient turbulent boundary layers around an
airfoil, the roles of the wall-normal convection and pressure gradient are of
particular importance and should be individually discussed, thus a further
decomposition of $C_{f,G}$ is carried out, viz.

$\displaystyle
C_{f,G}=\underbrace{\frac{2}{U_{e}^{3}}\int_{0}^{\infty}\left(\left<u\right>-U_{e}\right)\left(\left<v\right>\frac{\partial\left<u\right>}{\partial
y}\right){\rm
d}y}_{C_{f,C}}+\underbrace{\frac{2}{U_{e}^{3}}\int_{0}^{\infty}\left(\left<u\right>-U_{e}\right)I_{x}{\rm
d}y}_{C_{f,D}}$
$\displaystyle+\underbrace{\frac{2}{U_{e}^{3}}\int_{0}^{\infty}\left(\left<u\right>-U_{e}\right)\left(\frac{{\rm
d}p/\rho}{{\rm d}x}\right){\rm d}y}_{C_{f,P}},$ (2)

where $I_{x}={\partial\left<u^{\prime}u^{\prime}\right>}/{\partial
x}+\left<u\right>{\partial\left<u\right>}/{\partial
x}-\nu{\partial^{2}\left<u\right>}/{\partial x^{2}}$, $p$ is the static
pressure, and $\rho$ is the density. Contributions of the mean wall-normal
convection ($C_{f,C}$), streamwise development ($C_{f,D}$), and the pressure
gradient ($C_{f,P}$) are separated in equation (2).

We consider a set of well-resolved LESs of a NACA4412 airfoil at angle of
attack of $5^{\circ}$ at two chord Reynolds numbers, i.e. $Re_{c}$=$200,000$
and $400,000$. The simulations were performed with the spectral-element code
Nek5000, developed by Fischer et al. (2008). The spatial derivatives in the
incompressible Navier–Stokes equations are discretized employing a Garlerkin
method, following the $P_{N}-P_{N-2}$ formulation by Patera (1984) and the
solution is expressed within each spectral element in terms of a nodal-base of
Legendre polynomials on the Gauss–Lobatto–Legendre (GLL) quadrature points.
The discretization of the time derivatives is explicit for the non-linear
terms and implicit for the viscous term, employing an extrapolation and a
backward differentiation scheme, respectively, both of the third order. In
order to trigger transition to turbulence, we employed tripping through a
volume force, implemented as proposed by Schlatter and Örlü (2012), at
$x/c=0.1$ on both suction and pressure sides.

The cases with/without control, listed in Table 1, include various
configurations of uniform blowing and suction applied on the suction side and
uniform blowing applied on the pressure side.

Case | Control strategies | Intensity ($V_{wall}$) | $\Delta(c_{l})$ | $\Delta(c_{d})$ | $\Delta(L/D)$ | $Re_{\tau}$ | $\beta$
---|---|---|---|---|---|---|---
$Re200k,ss,ref$ | – | – | – | – | – | $[132,224]$ | $[0.16,11.07]$
$Re200k,ss,blw1$ | blowing, suction s. ($0.25<x/c<0.86$) | $0.1\%U_{\infty}$ | $-4\%$ | $+3\%$ | $-7\%$ | $[126,209]$ | $[0.14,19.62]$
$Re200k,ss,blw2$ | blowing, suction s. ($0.25<x/c<0.86$) | $0.2\%U_{\infty}$ | $-8\%$ | $+8\%$ | $-15\%$ | $[118,193]$ | $[0.12,36.54]$
$Re200k,ss,sct1$ | suction, suction s. ($0.25<x/c<0.86$) | $0.1\%U_{\infty}$ | $+4\%$ | $-2\%$ | $+6\%$ | $[139,238]$ | $[0.17,6.97]$
$Re200k,ss,sct2$ | suction, suction s. ($0.25<x/c<0.86$) | $0.2\%U_{\infty}$ | $+7\%$ | $-4\%$ | $+11\%$ | $[145,248]$ | $[0.18,4.71]$
$Re400k,ss,ref$ | – | – | – | – | – | $[183,363]$ | $[0.15,9.16]$
$Re400k,ss,blw1$ | blowing, suction s. ($0.25<x/c<0.86$) | $0.1\%U_{\infty}$ | $-4\%$ | $+5\%$ | $-9\%$ | $[174,335]$ | $[0.14,17.56]$
$Re400k,ss,sct1$ | suction, suction s. ($0.25<x/c<0.86$) | $0.1\%U_{\infty}$ | $+3\%$ | $-1\%$ | $+4\%$ | $[193,390]$ | $[0.15,5.37]$
$Re200k,ps,ref$ | – | – | – | – | – | $[96,219]$ | $[-0.30,-0.01]$
$Re200k,ps,blw1$ | blowing, pressure s. ($0.2<x/c<1.00$) | $0.1\%U_{\infty}$ | $+0\%$ | $-3\%$ | $+4\%$ | $[87,232]$ | $[-0.42,-0.03]$
$Re200k,ps,blw2$ | blowing, pressure s. ($0.2<x/c<1.00$) | $0.2\%U_{\infty}$ | $+1\%$ | $-5\%$ | $+7\%$ | $[80,234]$ | $[-0.55,-0.04]$

Table 1: Cases with/without control considered in the present paper and the
relative changes of the total lift and drag coefficients (denoted by
$\Delta(c_{l})$ and $\Delta(c_{d})$, respectively) and aerodynamic efficiency
(denoted by $\Delta(L/D)$) in respect to the reference case. Note that
“suction s.” and “pressure s.” are suction and pressure sides, respectively,
$Re200k$ and $Re400k$ are the chord Reynolds number ($Re_{c}$) of each case ,
$Re_{\tau}$ denotes the friction Reynolds number, and $\beta$ is the Rotta-
Clauser pressure-gradient parameter.

The relative proportions of pressure drag and skin-friction drag determine the
control effects on the total drag, denoted by $c_{d}$. At these moderate
Reynolds numbers, pressure drag is relatively high and uniform blowing on the
suction side increases it by an amount that is high enough to overcome the
skin-friction reduction, eventually leading to higher $c_{d}$. On the
contrary, uniform suction increases skin-friction drag, but it decreases the
pressure drag enough to result in lower $c_{d}$. At the same time, uniform
blowing and suction on the suction side also decreases and increases lift
($c_{l}$), respectively. Uniform blowing on the pressure side has different
effects on the pressure distribution around the airfoil, decreasing both skin
friction and pressure drag and increasing lift. The friction Reynolds numbers
($Re_{\tau}=u_{\tau}\delta_{99}/\nu$) and the Rotta-Clauser pressure-gradient
parameters (Rotta, 1950; Clauser, 1954, 1956) ($\beta=\delta^{*}/\tau_{w}{\rm
d}Pe/{\rm d}x$) are also listed in Table 1. Note that
$u_{\tau}=\sqrt{(\tau_{w}/\rho)}$ is the friction velocity, $\delta^{*}$ is
the displacement thickness, $\tau_{w}$ is the wall shear stress, and ${\rm
d}Pe/{\rm d}x$ is the streamwise pressure gradient at the boundary-layer edge.
For a more complete description of the numerical setup and the aerodynamic
effects of control, we refer to Vinuesa et al. (2018) and Atzori et al.
(2020).

##  3 Friction-drag decomposition on the suction side

###  3.1 The control effects

Using the database, we first show the variation of skin-friction coefficients
on the suction side of a NACA4412 wing section in figure 1. It can be easily
found that uniform blowing causes friction-drag reduction whereas suction
causes friction-drag increase, regardless of the Reynolds number and
streamwise position on the control surface. Stronger intensity of
blowing/suction leads to larger drag-reduction/increase rate, as expected.
Such phenomena are in consistence with previous studies (Kametani and
Fukagata, 2011; Kametani et al., 2015; Atzori et al., 2020). The mechanisms of
the drag reduction/increase by blowing/suction are associated with the
interactions between the “cross-stream” and the quasi-streamwise vortical
structures in the near-wall region, which probably enhance or damp the
behavior of sweep/ejection events and yield modifications of the mean velocity
profiles (Park and Choi, 1999; Kim et al., 2002; Hwang, 2004; Kornilov, 2015).
The variation of turbulent dynamics in the near-wall region leads to a
redistribution of the turbulence kinetic energy and alteration of the
turbulent momentum transport across the wall layer. The vortical structures in
the outer layer will also be influenced by the near-wall blowing/suction,
especially in the downstream of the control surface. Consequently, the
generation of the skin-friction drag, which is linked to the turbulence
statistics across the wall layer, will be correspondingly changed.

Figure 1: Variation of skin-friction coefficients on the suction side of a
NACA4412 wing section at ($a$) $Re_{c}=200,000$ and ($b$) $Re_{c}=400,000$.

In order to clarify such control effects, we conduct the decomposition of
skin-friction coefficients on the suction side (within $0.2\leq x/c\leq 0.85$)
of NACA4412 by RD identity. Note that the relative errors,
$(C_{f,V}+C_{f,T}+C_{f,G}-C_{f})/C_{f}$, where $C_{f}$ is directly calculated
with the normal gradient of tangential velocity at the wall surface (i.e.
$C_{f}=(\mu\partial\left<u\right>/\partial y)|_{wall}/(0.5\rho U_{e}^{2})$),
are well confined within $\pm 0.12\%$ for all cases considered, which confirms
the reliability of the decomposition method.

Figure 2 shows the variations of $C_{f,V}$, $C_{f,T}$, and $C_{f,G}$ induced
by uniform blowing and suction with regard to the reference case, at
$Re_{c}$=$200,000$ (figures 2–2) and $Re_{c}$=$400,000$ (figures 2–2). With
blowing, the friction contribution of direct viscous dissipation ($C_{f,V}$)
is reduced at both Reynolds numbers. This is closely related to the
suppression of the near-wall sweep events, which transport high-speed fluid
towards the wall. With blowing, an increase of the generation of turbulence-
kinetic-energy production ($C_{f,T}$) is found in figures 2 and 2. The
variations of $C_{f,V}$ and $C_{f,T}$ are essentially associated with the
influences on the wall-normal profiles of mean viscous shear stress and
Reynolds shear stress across the wall layer, which will be discussed in Sec.
3.2.

As for the generation of the spatial growth ($C_{f,G}$), it is decreased by
the blowing and the influence on $C_{f,G}$ is stronger than that on $C_{f,V}$
and $C_{f,T}$. To further clarify the cause of such variation, we trace back
to its sub-constituents in equation (2) and plot the variations of $C_{f,C}$,
$C_{f,D}$, and $C_{f,P}$ in figure 3. In the cases with blowing, the
generation of wall-normal convection ($C_{f,C}$) and adverse pressure gradient
($C_{f,P}$) is decreased, while that of streamwise development ($C_{f,D}$) is
increased. A slight turnup of $C_{f,P}$ is observed near the trailing edge for
the blowing cases (see figure 3), which probably relates to the fact that the
boundary layer is approaching the condition of mean separation (Atzori et al.,
2020). Generally, the positive variation of $C_{f,D}$ is overcome by the
negative influence on $C_{f,C}$ and $C_{f,P}$, which consequently leads to the
overall reduction of $C_{f,G}$ by blowing (Mahfoze et al., 2019), as shown in
figures 2 and 2. In Sec. 3.2, we will further analyze the wall-normal
distributions of these friction constituents, to reveal the connection of the
$C_{f}-$constituents with the turbulent dynamics across the boundary layer.
With mass suction, the control effects on the $C_{f}-$constituents shown in
figure 2 and 3 are opposite to those with blowing. Here we did not add more
discussions on the suction cases for brevity.

Figure 2: Variation of ($a,d$) $C_{f,V}$, ($b,e$) $C_{f,T}$, and ($c,f$)
$C_{f,G}$ with regard to the reference case on the suction side of a NACA4412
wing section at ($a-c$) $Re_{c}=200,000$ and ($d-f$) $Re_{c}=400,000$. (The
superscript of “ref” represents the reference case without control. Legends in
($a-c$) refer to figure 1, while those in ($d-f$) refer to figure 1.)

Figure 3: Variation of ($a,d$) $C_{f,C}$, ($b,e$) $C_{f,D}$, and ($c,f$)
$C_{f,P}$ with regard to the reference case on the suction side of a NACA4412
wing section at ($a-c$) $Re_{c}=200,000$ and ($d-f$) $Re_{c}=400,000$.

Figure 4: Total friction-drag change rate as a function of control intensity
($V_{wall}/U_{\infty}$) on the suction side of a NACA4412 wing section.

In order to assess the control effects on the total skin-friction drag over
the control surface, a parameter $\bar{D}_{f,i}$ is introduced:

$\bar{D}_{f,i}={\int\limits_{\Omega_{\rm
ctr}}{\tau_{w,i}}\left(\vec{t}\cdot\vec{n}\right){\rm d}\xi},$ (3)

where ${\tau_{w,i}}=C_{f,i}\cdot\left(0.5\rho U_{e}^{2}\right)$ is the
decomposed component of wall-shear stress, with the subscript ‘i’ working as a
label to denote each $C_{f}-$constituent, i.e. ‘V’, ‘T’, ‘G’, ‘C’, ‘D’, and
‘P’ as mentioned in equations (2) and (2), $\vec{t}$ and $\vec{n}$ denote the
unit vectors tangential to the airfoil surface and along the free-stream
direction respectively, $\xi$ is the curvilinear coordinate along the airfoil
surface, and $\Omega_{\rm ctr}$ represents the area of the control surface.
The change rate of $\bar{D}_{f,i}$, with respect to the uncontrolled reference
case, is then defined as:

$R_{i}=\frac{\left(\bar{D}_{f,i}-\bar{D}_{f,i}^{\rm
ref}\right)}{\bar{D}_{f}^{\rm ref}}.$ (4)

Figure 4 shows the result of $R_{i}$ under different control schemes at
$Re_{c}$=$200,000$ and $400,000$. It can be seen that $R_{i}$ appears to be
linearly dependent on the control intensity of the blowing/suction, within
$-0.2\%\leq V_{\rm wall}/U_{\infty}\leq 0.2\%$. However, note that this might
not always be true when $V_{\rm wall}/U_{\infty}$ becomes much larger, and
needs to be validated in the future work. Among the decomposed constituents,
the most significant control effect lies on the friction constituent of
spatial growth of the flow ($C_{f,G}$), with its sub-constituents primarily
correlated with the convection, streamwise growth, and pressure gradient in
the outer region (Fan et al., 2020a). This reveals that the tremendously
influenced outer-layer dynamics play an important role in the drag control
with blowing/suction. Moreover, weak Reynolds-number effects are found,
especially for $R_{V}$ and $R_{G}$, where a stronger control effect is
achieved at higher Reynolds number which to some extent validates the
theoretical estimation by Kametani and Fukagata (2011).

###  3.2 Wall-normal distributions of the $C_{f}-$constituents

To answer the question that how uniform blowing/suction specifically
influences the sources of skin-friction generation, the wall-normal
distributions of the decomposed $C_{f}-$constituents across the boundary layer
are investigated. As $C_{f}$ and its constituents vary along the streamwise
direction on the wing surface, we only discuss the wall-normal contributions
of the $C_{f}-$constituents at $x/c\approx 0.75$, where the friction Reynolds
numbers are $Re_{\tau}\approx 224$, $205$, $180$, $237$, $245$, $362$, $332$,
and $387$ in the case of “$Re200k,ss,ref$”, “$Re200k,ss,blw1$”,
“$Re200k,ss,blw2$”, “$Re200k,ss,sct1$”, “$Re200k,ss,sct2$”, “$Re400k,ss,ref$”,
“$Re400k,ss,blw1$”, and “$Re400k,ss,sct1$” respectively. Similar conclusions
can be drawn at other positions within $0.2\leq x/c\leq 0.85$, and the results
are not shown here for simplicity.

The $C_{f}-$constituents are expressed in intrinsic scales as:

$\displaystyle{C_{f,V}}$ $\displaystyle=$ $\displaystyle
2/U_{e}^{+3}\int_{0}^{\infty}\left(\frac{\partial\left<u\right>^{+}}{\partial
y^{+}}\right)^{2}{\rm d}y^{+},$ (5) $\displaystyle{C_{f,T}}$ $\displaystyle=$
$\displaystyle
2/U_{e}^{+3}\int_{0}^{\infty}\left<-u^{\prime}v^{\prime}\right>^{+}\frac{\partial\left<u\right>^{+}}{\partial
y^{+}}{\rm d}y^{+},$ (6) $\displaystyle{C_{f,G}}$ $\displaystyle=$
$\displaystyle
2/U_{e}^{+3}\int_{0}^{\infty}\left(\left<u\right>^{+}-U_{e}^{+}\right)\frac{\partial}{\partial
y^{+}}\left(\frac{\partial\left<u\right>^{+}}{\partial
y^{+}}-\left<u^{\prime}v^{\prime}\right>^{+}\right){\rm d}y^{+},$ (7)

where the superscript $+$ denotes normalization by viscous units, i.e.
friction velocity $u_{\tau}$ and viscous length scale
$\delta_{\nu}=\nu/u_{\tau}$.

Figure 5: Pre-multiplied integrands of ($a,d$) $C_{f,V}$, ($b,e$) $C_{f,T}$,
and ($c,f$) $C_{f,G}$ at $x/c\approx 0.75$, as a function of $y^{+}$, on the
suction side of a NACA4412 wing section at ($a-c$) $Re_{c}=200,000$ and
($d-f$) $Re_{c}=400,000$.

Figure 6: Profiles of ($a,c$) wall-normal velocity gradient and ($b,d$)
Reynolds shear stress, on the suction side of a NACA4412 wing section at
($a,b$) $Re_{c}=200,000$ and ($c,d$) $Re_{c}=400,000$.

Figure 5 shows the wall-normal distributions of the pre-multiplied integrand
of $C_{f,V}$, $C_{f,T}$, and $C_{f,G}$ in equations (5)–(7), as a function of
$y^{+}$. The semi-logarithmic plots retain the advantage that the areas
beneath the curves directly yield the total generation of the constituents.
Comparisons between figures 5–5 for the cases at $Re_{c}$=$200,000$ and 5–5
for the cases at $Re_{c}$=$400,000$ confirms that the Reynolds-number
variation does not change the conclusions which will be presented in the
following from a qualitative perspective.

For the distribution of $C_{f,V}-$contributions, two peaks are respectively
observed in the near-wall and outer region of APG-TBLs. Most of the
$C_{f,V}-$contributions come from the inner region ($y^{+}<30$), indicating
that the viscous dissipation is mostly concentrated in the near-wall region,
as expected. In the meantime, a secondary peak appears in the outer region,
which is probably due to the energy enhancement by APG (Tanarro et al., 2020;
Sanmiguel Vila et al., 2020). The secondary peak is absent in the ZPG-TBLs
even at higher friction Reynolds number up to $Re_{\tau}=1270$ (Fan et al.,
2019b). When uniform blowing/suction is applied, the locations of inner peaks
are fixed at a wall-normal distance of $y^{+}\approx 5.0-6.0$, regardless of
the control scheme. In the blowing cases, the inner peak of
$C_{f,V}-$contributions is reduced whereas the outer peak is increased,
suggesting a lowered mean shear in the near-wall region while enhanced in the
outer region due to the lifting-up of boundary layer (Kornilov, 2015), which
is also validated by showing the wall-normal gradient of the tangential
velocity in figures 6 and 6. This reveals that blowing has different actions
in different sub-layers, namely inhibiting the contribution of inner-layer
dynamics to skin-friction generation while promoting that of outer-layer
dynamics. When blowing intensity is up to 0.2%, the outer-layer contributions
seem to be comparable to the inner-layer contributions. On the other hand,
suction behaves quite opposite for all constituents, which will not be
repeated hereafter.

Similar inner and outer peaks are also observed in the pre-multiplied
distribution of $C_{f,T}-$contributions, as shown in figures 5 and 5, with the
former well collapsed at the inner-scaled wall-normal distance $y^{+}\approx
16.0-17.0$. It can be found that the outer-layer motions dominate the
contributions of $C_{f,T}$, although the $Re_{\tau}$ is lower than $400$,
which is much different from the features in ZPG-TBLs (Fan et al., 2019b). The
prominent peak in the outer region suggests the energization of large-scale
outer motions by APGs (Harun et al., 2013). When blowing is applied, the inner
peak of $C_{f,T}-$contributions is reduced, while the outer peak is increased.
This phenomenon is linked to the wall-normal distributions of the wall-normal
velocity gradient and Reynolds shear stress, as shown in figure 6. With
blowing, the wall-normal velocity gradients are suppressed significantly in
the near-wall region, probably resulting from the damping of the near-wall
sweep motions. On the other hand, in the outer layer, the Reynolds shear
stress is amplified by the blowing. These two actions are consequently
responsible for the changes of the inner and outer peaks in the distribution
of $C_{f,T}-$contributions.

As for the distribution of $C_{f,G}-$contributions in figures 5 and 5,
negative contributions are observed within $y^{+}\lesssim 80$ for the
uncontrolled reference case at $Re_{c}$=$200,000$ ($y^{+}\lesssim 130$ at
$Re_{c}$=$400,000$), and positive values beyond this region. This differs from
the result in ZPG-TBLs: $C_{f,G}$ always remains positive across the wall
layer (Fan et al., 2020a). Blowing enhances both the negative and positive
distributions, as the strengthened adverse pressure gradient promotes a more
pronounced growth of the boundary layer and a more prominent outer region
(Vinuesa et al., 2018).

In contrast to the wall-normal distributions as a function of $y^{+}$ shown in
figure 5, figure 7 plots their profiles as a function of $y/\delta_{99}$. The
outer-peak locations of $C_{f,V}-$, $C_{f,T}-$, and $C_{f,G}-$contributions
normalized by the outer scale are well-collapsed at $y/\delta_{99}\approx
0.7$, $0.53$, and $0.65$, respectively, as marked with vertical dashed lines
in figure 7, regardless of the control scheme and Reynolds number. These
phenomena are consistent with our previous finding (Fan et al., 2020a), i.e.
the inner-peak locations (in $C_{f,V}-$ and $C_{f,T}-$contribution) exhibit
good scaling in the inner unit ($\delta_{\nu}$), and the outer-peak locations
in the outer unit ($\delta_{99}$), regardless of the friction Reynolds number,
the magnitude of APG and its development history. This finding suggests that
self similarity is exhibited in inner or outer scales for the turbulence
statistics associated with the friction-drag generation.

Figure 7: Pre-multiplied integrands of ($a,d$) $C_{f,V}$, ($b,e$) $C_{f,T}$,
and ($c,f$) $C_{f,G}$, as a function of $y/\delta_{99}$, on the suction side
of a NACA4412 wing section at ($a-c$) $Re_{c}=200,000$ and ($d-f$)
$Re_{c}=400,000$.

Figure 8: Pre-multiplied integrands of ($a,d$) $C_{f,C}$, ($b,e$) $C_{f,D}$,
and ($c,f$) $C_{f,P}$, as a function of $y/\delta_{99}$, on the suction side
of a NACA4412 wing section at ($a-c$) $Re_{c}=200,000$ and ($d-f$)
$Re_{c}=400,000$.

The generation of $C_{f,G}$ results from a counterbalance between the negative
work done by $C_{f,C}$ and $C_{f,P}$ and the positive work by $C_{f,D}$.
Figure 8 quantifies their wall-normal distributions, and only the outer
scaling by $\delta_{99}$ is applied herein. Good collapses of the peak
locations are also observed at $y/\delta_{99}\approx 0.59$, $0.56$, and $0.43$
for $C_{f,C}-$, $C_{f,D}-$, and $C_{f,P}-$contributions, respectively. In the
blowing cases, the APG effects, a fact that promotes the
population/energization of outer-layer structures (Harun et al., 2013), are
strengthened. Meanwhile, the wall-normal convection and the streamwise
boundary-layer growth are intensified in the outer region (Vinuesa et al.,
2018). Therefore, in an absolute sense, the generations of the components
($C_{f,C}$, $C_{f,D}$, and $C_{f,P}$) are all enhanced. Thereinto, as the
convection and pressure gradient do negative work for the friction-drag
generation, blowing acts to reduce the $C_{f,C}-$ and $C_{f,P}-$contributions.
On the other hand, the positive contribution of $C_{f,D}$ is increased, as
shown in figures 8 and 8. These three components counterbalance each other.
The negative $C_{f,C}-$ and $C_{f,P}-$contributions are responsible for the
$C_{f,G}$ reduction in the near-wall region, whereas the positive
$C_{f,D}-$contributions is responsible for the increase of $C_{f,G}$ in the
outer region, as seen in figures 5 and 5.

###  3.3 Contributions of small- and large-scale structures to the friction-
drag generation

Firstly, we use empirical mode decomposition (EMD) (Huang et al., 1998) to
identity the small- and large-scale turbulence structures. EMD is an adaptive
mode-decomposition technique, which extracts characteristic wavelengths of
non-stationary signals automatically without a-priori basis functions. It has
been applied for wall-bounded turbulence (Agostini and Leschziner, 2014, 2019;
Cheng et al., 2019; Dogan et al., 2019), and details of the methodology of EMD
can be found in Huang et al. (1998). Here we just describe EMD very briefly.

With EMD, a raw temporal or spatial signal $f(t)$ is decomposed into a sum of
multiple intrinsic mode functions (IMFs) with a residual $R(t)$:

$f(t)=\sum_{i=1}^{m}{\rm IMF_{i}(t)}+R(t),$ (8)

where $m$ is the number of IMFs. The IMFs are data-driven functions,
representing components with different wavelengths or scales in the full
field. In this study, the velocity fluctuations ($u^{\prime}$ and
$v^{\prime}$) at the streamwise location $x/a\approx 0.75$ are decomposed into
four modes (three IMFs with a final residual). The first two modes represent
the small-scale structures and the others characterize the large scales, which
is justified based on a preliminary analysis (similar to the studies of
Agostini and Leschziner (2014, 2019)).

Figure 9: Spanwise pre-multiplied spectra of the tangential and normal
velocity fluctuations on the suction side of a NACA4412 wing section at
$Re_{c}=400,000$ ($a,c$) without control, with ($b,e$) uniform blowing and
($c,f$) uniform suction. The nephogram represents the spectra of the full
field; the dashed and solid contour lines represent those of small and large
scales, respectively.

The contours of spanwise pre-multiplied spectra of $u^{\prime}$ and
$v^{\prime}$ scaled with friction velocities are plotted in figure 9 for the
$Re_{c}$=$400,000$ cases. Similar features can also be found in the low-
Reynolds-number ($Re_{c}$=$200,000$) cases, which are not shown here for
brevity. In figure 9, the iso-contour levels marked with dashed and solid
lines represent the spectra of small- and large-scale velocity fluctuations,
respectively. These contour lines indicate 0.12, 0.42, and 0.72 of their
maxima, from outside to inside. For the reference case, as shown in figure 9,
the spectra of small-scale $u^{\prime}-$structures peak at the wall-normal
distance $y^{+}\approx 10$ with a spanwise wavelength $\lambda_{z}^{+}\approx
80$, and those of large-scale $u^{\prime}-$structures peak at $y^{+}\approx
100$ ($y/\delta_{99}\approx 0.28$) with $\lambda_{z}^{+}\approx 310$
($\lambda_{z}/\delta_{99}\approx 0.85$). This observation is consistent with
the study of Cheng et al. (2019), that small-scale $u^{\prime}-$structures
identified by EMD are representative of the near-wall coherent motions,
whereas the large-scale $u^{\prime}-$structures characterize large-scale
motions.

When blowing/suction is applied on the airfoil surface, such energy spectra
are affected, as seen in figures 9 and 9. In the case of blowing, the small-
scale structures are enhanced, and penetrate deeper into the outer region. The
peak of spectra locates at $(\lambda_{z}^{+},y^{+})\approx(80,8)$. Whereas
suction has the opposite influence on the amplitude of small scales, with the
peak location faintly influenced at $(\lambda_{z}^{+},y^{+})\approx(90,10)$.
As for the large-scale structures, they are enhanced by blowing, due to the
energized large-scale motions, which agrees well with the conclusion drawn by
Kametani et al. (2015). The peak of large-scale spectra locates at
$(\lambda_{z}/\delta_{99},y/\delta_{99})\approx(0.96,0.33)$. As shown in
figure 9, suction trends to diminish the secondary peak in the outer region
with the peak at $(\lambda_{z}/\delta_{99},y/\delta_{99})\approx(0.93,0.35)$.

As for the $v^{\prime}-$structures shown in figures 9–9, the small- and large-
scale structures have approximately the same spanwise wavelengths as the
$u^{\prime}-$structures, which is consistent with the EMD results of channel
flows (Cheng et al., 2019), suggesting that the decomposed two scales of
$u^{\prime}-$ and $v^{\prime}-$structures can be characterized with the same
spanwise wavelength. Both the small- and large-scale $v^{\prime}-$structures
are greatly enhanced by blowing while reduced by suction. On the other hand,
the wall-normal locations of the $v^{\prime}-$structures are much different
from those of $u^{\prime}-$structures, because that the presence of wall
prevents the normal velocity fluctuations from extending close to the near-
wall region, in contrast to the wall-parallel component $u^{\prime}$. As shown
in figures 9–9, the small-scale spectra of $v^{\prime}$ are more intense
beyond the buffer layer at $y^{+}\approx 50$, with the spanwise wavelength
scale $\lambda_{z}^{+}$ peaking at around $80$–$100$ for all of the three
cases, and the large-scale $v^{\prime}-$spectra peak around
$(\lambda_{z}/\delta_{99},y/\delta_{99})\approx(0.9,0.35)$.

Hereafter, we denote the small- and large-scale tangential and wall-normal
velocity fluctuations as $u_{s}^{\prime}$, $u_{l}^{\prime}$, $v_{s}^{\prime}$,
and $v_{l}^{\prime}$, respectively. Then the Reynolds stress is decomposed as:

$-\left<u^{\prime}v^{\prime}\right>=-\left<u_{s}^{\prime}v_{s}^{\prime}\right>-\left<u_{s}^{\prime}v_{l}^{\prime}\right>-\left<u_{l}^{\prime}v_{s}^{\prime}\right>-\left<u_{l}^{\prime}v_{l}^{\prime}\right>,$
(9)

where $-\left<u_{s}^{\prime}v_{s}^{\prime}\right>$ represents the Reynolds
stress carried by small-scale structures,
$-\left<u_{l}^{\prime}v_{l}^{\prime}\right>$ represents the Reynolds stress
associated with large-scale structures, and
$-\left<u_{s}^{\prime}v_{l}^{\prime}\right>$ and
$-\left<u_{l}^{\prime}v_{s}^{\prime}\right>$ denote the scale interactions
from small- to large- and from large- to small-scale structures, respectively.
Substituting these Reynolds-stress components into equation (6), the term of
TKE production ($C_{f,T}$) can be further divided into four parts, viz.

$\displaystyle{C_{f,T,ss}}$ $\displaystyle=$ $\displaystyle
2/U_{e}^{+3}\int_{0}^{\infty}\left<-u_{s}^{\prime}v_{s}^{\prime}\right>^{+}\frac{\partial\left<u\right>^{+}}{\partial
y^{+}}{\rm d}y^{+},$ (10) $\displaystyle{C_{f,T,sl}}$ $\displaystyle=$
$\displaystyle
2/U_{e}^{+3}\int_{0}^{\infty}\left<-u_{s}^{\prime}v_{l}^{\prime}\right>^{+}\frac{\partial\left<u\right>^{+}}{\partial
y^{+}}{\rm d}y^{+},$ (11) $\displaystyle{C_{f,T,ls}}$ $\displaystyle=$
$\displaystyle
2/U_{e}^{+3}\int_{0}^{\infty}\left<-u_{l}^{\prime}v_{s}^{\prime}\right>^{+}\frac{\partial\left<u\right>^{+}}{\partial
y^{+}}{\rm d}y^{+},$ (12) $\displaystyle{C_{f,T,ll}}$ $\displaystyle=$
$\displaystyle
2/U_{e}^{+3}\int_{0}^{\infty}\left<-u_{l}^{\prime}v_{l}^{\prime}\right>^{+}\frac{\partial\left<u\right>^{+}}{\partial
y^{+}}{\rm d}y^{+}.$ (13)

Figure 10: Pre-multiplied integrands of $C_{f,T,i}$ as a function of $y^{+}$
on the suction side of a NACA4412 wing section at $Re_{c}=400,000$ ($a$)
without control, with ($b$) uniform blowing and ($c$) uniform suction. Note
that the premultiplication factor is $y^{+}$.

Figure 11: Scale-specific contributions ($C_{f,T,ss}$, $C_{f,T,sl}$,
$C_{f,T,ls}$, and $C_{f,T,ll}$ normalized by $C_{f,T}$ itself) on the suction
side of the NACA4412 wing section at $Re_{c}=400,000$ ($a$) without control,
with ($b$) uniform blowing and ($c$) uniform suction.

Similar to the analysis in Sec. 3.2, we show the wall-normal distributions of
the $y^{+}$-pre-multiplied integrands in equations (10)–(13) in figure 10. For
the reference case without control, the small-scale motions play important
roles both in the inner and outer region, as two comparable peaks are
respectively observed there, as seen in figure 10. This is different from the
phenomenon in channel flows (Agostini and Leschziner, 2019) and ZPG-TBLs,
since the APG strengthens the energy of small-scale structures in the outer
region (Tanarro et al., 2020). Whereas the large-scale structures dominate the
production of $C_{f,T}$ in the outer region, which is associated with the
enhanced generation of large-scale motions in APG-TBLs (Lee and Sung, 2008;
Harun et al., 2013; Vinuesa et al., 2018). When blowing is applied, the
distributions of $C_{f,T,ss}-$contribution and $C_{f,T,ll}-$contribution in
the outer region are shifted upwards, suggesting that the wall-normal mass
flux enhances the outer fluctuations of both small- and large-scale
structures. Comparison between $C_{f,T,ss}$ and $C_{f,T,ll}$ indicates that
the blowing raises the relative importance of large scales with respect to the
small ones. In addition, the scale interactions ($C_{f,T,sl}$ and
$C_{f,T,ls}$) mainly exhibited in the outer region are not very sensitive to
the control schemes.

In order to further check the role of small- and large-scale motions in the
generation of turbulence-kinetic-energy production, figure 11 quantifies their
integrations normalized by $C_{f,T}$ itself. We can find that the blowing and
suction have opposite influences on the contributions of small- and large-
scale motions, that is, blowing is able to enhance the contribution of large-
scales and to suppress the contribution of small-scales, whereas suction
behaves contrarily. As for the scale-interactions, i.e. $C_{f,T,sl}$ and
$C_{f,T,ls}$, they account for approximately 20% of the total $C_{f,T}$ and
remain almost unchanged with different control strategies.

##  4 Friction-drag decomposition on the pressure side

In this section, we pay attention to the pressure side of a NACA4412 wing
section, where the TBLs are subjected to favorable pressure gradients. Figure
12 shows the skin-friction coefficients on the pressure side of the airfoil.
In spite of the favorable pressure gradients, blowing is still able to reduce
the friction drag. Whereas, the skin-friction coefficients are no longer
decreased monotonously with regard to $x/c$ if compared with those on the
suction side of the airfoil, which may be caused by the coupling influences of
$Re_{\tau}$ and FPGs (Atzori et al., 2020). The relative errors
$([C_{f,V}+C_{f,T}+C_{f,G}-C_{f}]/C_{f})$ of the friction-drag decomposition
are limited within $\pm 0.03\%$.

Key points of the decomposition results on the pressure side are listed below:

Figure 12: Variation of skin-friction coefficients on the pressure side of the
NACA4412 wing section at $Re_{c}=200,000$.

Figure 13: Variation of ($a$) $C_{f,V}$, ($b$) $C_{f,T}$, ($c$) $C_{f,G}$,
($d$) $C_{f,C}$, ($e$) $C_{f,D}$, and ($f$) $C_{f,P}$ with regard to the
reference case on the pressure side of NACA4412 at $Re_{c}=200,000$.

$\bullet$ The variations of $C_{f,V}$, $C_{f,T}$, and $C_{f,G}$ (see figures
13–13) are similar to the results on the suction side, that is, by blowing,
$C_{f,V}$ and $C_{f,T}$ are decreased, and $C_{f,G}$ is reduced. As for the
sub-constituents of $C_{f,G}$ in figures 13–13, blowing reduces $C_{f,C}$,
similar to the observations on the suction side; however, differently, blowing
just increases $C_{f,D}$ until $x/c\approx 0.85$ and then decreases; in
contrary to the control effects on the suction side, $C_{f,P}$ is increased,
suggesting that blowing strengthens the FPG on the pressure side of NACA4412.
Stronger intensity of blowing leads to larger increases or decreases of the
constituents.

Figure 14: Total friction-drag change rate as a function of control intensity
($V_{wall}/U_{\infty}$) on the pressure side of NACA4412.

$\bullet$ Again, near-linear dependence of $R_{i}$ on the intensity of blowing
is still observed within the range of $V_{wall}/U_{\infty}$ under scrutiny, as
seen in figure 14.

Figure 15: Pre-multiplied integrands of ($a$) $C_{f,V}$, ($b$) $C_{f,T}$, and
($c$) $C_{f,G}$, as a function of $y^{+}$, on the pressure side of NACA4412 at
$Re_{c}=200,000$.

Figure 16: Pre-multiplied integrands of ($a$) $C_{f,C}$, ($b$) $C_{f,D}$, and
($c$) $C_{f,P}$, as a function of $y/\delta_{99}$, on the pressure side of
NACA4412 at $Re_{c}=200,000$.

$\bullet$ At $x/c\approx 0.75$ (similar features are observed at other
positions), the wall-normal distributions of the $C_{f}-$constituents are
plotted in figure 15. With the presence of FPG, the outer peak of $C_{f,V}$ is
barely legible, and the outer peak of $C_{f,T}$ is comparable to the inner
one, which are different from the results on the suction side. With blowing,
the outer peak of $C_{f,T}$ is also increased however by a lesser amount,
probably because that the FPG attenuates the outer-layer Reynolds shear stress
and the production of turbulence kinetic energy (Harun et al., 2013). In
figure 16, different from the observation in figure 8, the generation of
$C_{f,C}$ mainly arises from the inner-layer dynamics on the pressure side of
the airfoil.

$\bullet$ The peak locations in figures 15 and 16 are well collapsed when
normalized by inner or outer units. It suggests that the self-similar features
of turbulence statistics will not be affected by the presence of FPG. Note
that the distributions in figures 15-15 are not plotted as a function of
$y/\delta_{99}$ since the friction Reynolds numbers of these three cases at
$x/c\approx 0.75$ are very similar, being $Re_{\tau}\approx 177$, $185$, and
$187$ in the case of “$Re200k,ps,ref$”, “$Re200k,ps,blw1$”, and
“$Re200k,ps,blw2$”, respectively.

$\bullet$ At last, EMD was not conducted, since the inner-outer scale
separation is less evident in the FPG-TBLs at such low $Re_{\tau}$ than those
on the suction side of the airfoil.

##  5 Conclusion

We employed the RD identity in conjunction with empirical mode decomposition
(EMD) to study the control effects of uniform blowing and suction on the
generation of mean friction drag on a NACA4412 airfoil at chord Reynolds
numbers $Re_{c}$=$200,000$ and $400,000$. In general, blowing reduces the mean
friction drag, and suction increases the friction drag, both on the suction
and pressure sides of the airfoil. With the RD identity, the mean friction
drag on the airfoil is decomposed into three components, associated with
viscous dissipation ($C_{f,V}$), turbulence-kinetic-energy production
($C_{f,T}$), and spatial growth of the flow ($C_{f,G}$), respectively. The
$C_{f,G}$ component is further decomposed into three terms related to the mean
wall-normal convection ($C_{f,C}$), streamwise development ($C_{f,D}$), and
the pressure gradient ($C_{f,P}$). The effects of suction on these
$C_{f}-$constituents are quite opposite to those of blowing, thus we just
summarize the key conclusions in the blowing cases.

For the adverse-pressure-gradient turbulent boundary layers on the suction
side of the airfoil, concluding remarks are listed as below:

$\bullet$ Blowing reduces the generation of $C_{f,V}$ and $C_{f,G}$, while
increases that of $C_{f,T}$. The drag reduction with blowing is mostly
attributed to $C_{f,G}$, which results from the decreased $C_{f,C}$ and
$C_{f,P}$, which overwhelm the increase in $C_{f,D}$. The integrated
$C_{f}-$constituents over the control surface are observed to be linearly
dependent on the intensity of blowing, and weakly influenced by the chord
Reynolds number (at least within the parameters we considered).

$\bullet$ Wall-normal distributions of the $C_{f}-$constituents at $x/c\approx
0.75$ are checked to clarify how the control schemes impact the sources of
friction-drag generation. With blowing, the generation of $C_{f,V}$, which is
mainly related to the inner-layer dynamics, is suppressed in the inner region
($y^{+}<30$), while it is enhanced in the outer region. Contrarily, the
generation of $C_{f,T}$, mainly originats from the outer-layer motions, is
amplified by the blowing. These phenomena are linked to the variations of the
wall-normal velocity gradients and Reynolds shear stresses in the wall-normal
direction. The generation of $C_{f,G}$ results from a counterbalance between
the negative work done by $C_{f,C}$ and $C_{f,P}$ and the positive work by
$C_{f,D}$. Blowing is able to enhance the generation of all sub-components, as
the strengthened adverse pressure gradient promotes a more pronounced growth
of the boundary layer and a more prominent outer region (Vinuesa et al.,
2018).

$\bullet$ We observed that, in the wall-normal direction, the inner-peak
locations of $C_{f,V}-$ and $C_{f,T}-$contributions scale well in the inner
unit ($\delta_{\nu}$), and the outer-peak locations of $C_{f,V}-$, $C_{f,T}-$,
$C_{f,G}-$ as well as its sub-contributions, collapse well in the outer unit
($\delta_{99}$), regardless of the friction Reynolds number, control scheme,
and the intensity of blowing/suction. This reveals that self similarity is
exhibited in inner or outer scales for the turbulence statistics associated
with the friction-drag generation.

$\bullet$ The small- and large-scale structures are separated with empirical
mode decomposition (EMD), aiming to analyze the scale-specific contribution of
turbulent motions to friction-drag generation. Results unveil that, normalized
by $C_{f,T}$ itself, blowing is able to enhance the contribution of large-
scale motions and to suppress that of small scales; note that suction behaves
contrarily. The contributions related to cross-scale interactions remain
almost unchanged with different control strategies.

For the favorable-pressure-gradient turbulent boundary layers on the pressure
side of the airfoil, the most significant observation is that the outer-layer
motions are of less importance for the generation of $C_{f}-$constituents. In
the case of blowing, the generation of $C_{f,P}$ is increased, which is the
opposite behavior to that on the suction side of the airfoil.

## Acknowledgments

The funding support of the National Natural Science Foundation of China (under
the grant No. 91952302 and 92052101) is acknowledged. Davide Gatti
acknowledges support by the state of Baden-Württemberg through bwHPC. Marco
Atzori, Ricardo Vinuesa and Philipp Schlatter also acknowledge support from
the Swedish Foundation for Strategic Research, project “In-Situ Big Data
Analysis for Flow and Climate Simulations” (ref. number BD15-0082), from the
Knut and Alice Wallenberg Foundation and from the Swedish Research Council
(VR). The simulations were performed on resources provided by the Swedish
National Infrastructure for Computing (SNIC) and within the project CWING on
the national supercomputer Cray XC40 Hazel Hen at the High Performance
Computing Center Stuttgart (HLRS).

## References

  * Agostini and Leschziner (2014) L. Agostini and M. A. Leschziner.  On the influence of outer large-scale structures on near-wall turbulence in channel flow.  _Phys. Fluids_ , 26(7):075107, 2014.
  * Agostini and Leschziner (2016) L. Agostini and M. A. Leschziner.  Predicting the response of small-scale near-wall turbulence to large-scale outer motions.  _Phys. Fluids_ , 28(1):015107, 2016.
  * Agostini and Leschziner (2019) L. Agostini and M. A. Leschziner.  The connection between the spectrum of turbulent scales and the skin-friction statistics in channel flow at ${R}e_{\tau}\approx 1000$.  _J. Fluid Mech._ , 871:22–51, 2019.
  * Ansell and Balajewicz (2017) P. J. Ansell and M. J. Balajewicz.  Separation of unsteady scales in a mixing layer using empirical mode decomposition.  _AIAA J._ , 55(2):419–434, 2017.
  * Atzori et al. (2020) M. Atzori, R. Vinuesa, G. Fahland, A. Stroh, D. Gatti, B. Frohnapfel, and P. Schlatter.  Aerodynamic effects of uniform blowing and suction on a NACA4412 airfoil.  _Flow, Turbul. Combust._ , 105:735–759, 2020.
  * Bannier et al. (2015) A. Bannier, É. Garnier, and P. Sagaut.  Riblet flow model based on an extended FIK identity.  _Flow, Turbul. Combust._ , 95(2-3):351–376, 2015\.
  * Bobke et al. (2016) A. Bobke, R. Örlü, and P. Schlatter.  Simulations of turbulent asymptotic suction boundary layers.  _J. Turbul._ , 17(2):157–180, 2016.
  * Cheng et al. (2019) C. Cheng, W.-P. Li, A. Lozano-Durán, and H. Liu.  Identity of attached eddies in turbulent channel flows with bidimensional empirical mode decomposition.  _J. Fluid Mech._ , 870:1037–1071, 2019.
  * Clauser (1954) F. H. Clauser.  The turbulent boundary layer in adverse pressure gradients.  _J. Aero. Sci._ , 21:91–108, 1954.
  * Clauser (1956) F. H. Clauser.  The turbulent boundary layer.  _Adv. Appl. Mech._ , 4:1–51, 1956.
  * Deck et al. (2014) S. Deck, N. Renard, R. Laraufie, and P. É. Weiss.  Large-scale contribution to mean wall shear stress in high-Reynolds-number flat-plate boundary layers up to ${Re}_{{\theta}}{=}$ 13650\.  _J. Fluid Mech._ , 743:202–248, 2014.
  * Dogan et al. (2019) E. Dogan, R. Örlü, D. Gatti, R. Vinuesa, and P. Schlatter.  Quantification of amplitude modulation in wall-bounded turbulence.  _Fluid Dyn. Res._ , 51:011408, 2019.
  * Eto et al. (2019) K. Eto, Y. Kondo, K. Fukagata, and N. Tokugawa.  Assessment of friction drag reduction on a Clark-Y airfoil by uniform blowing.  _AIAA J._ , 57(7):2774–2782, 2019.
  * Fahland et al. (2021) G. Fahland, A. Stroh, B. Frohnapfel, M. Atzori, R. Vinuesa, P. Schlatter, and D. Gatti.  Investigation of blowing and suction for turbulent flow control on airfoils.  _AIAA J._ , (To Appear), 03 2021.
  * Fan et al. (2019a) Y.-T. Fan, C. Cheng, and W.-P. Li.  Effects of the Reynolds number on the mean skin friction decomposition in turbulent channel flows.  _Appl. Math. Mech. (English Ed.)_ , 40(3):331–342, 2019a.
  * Fan et al. (2019b) Y.-T. Fan, W.-P. Li, and S. Pirozzoli.  Decomposition of the mean friction drag in zero-pressure-gradient turbulent boundary layers.  _Phys. Fluids_ , 31(8):086105, 2019b.
  * Fan et al. (2020a) Y.-T. Fan, W.-P. Li, M. Atzori, R. Pozuelo, P. Schlatter, and R. Vinuesa.  Decomposition of the mean friction drag in adverse-pressure-gradient turbulent boundary layers.  _Phys. Rev. Fluids_ , 5:114608, 2020a.
  * Fan et al. (2020b) Y.-T. Fan, W.-P. Li, and S. Pirozzoli.  Energy-based decomposition of friction drag in turbulent square-duct flows.  _Int. J. Heat Fluid Flow_ , 86:108731, 2020b.  ISSN 0142-727X.
  * Fischer et al. (2008) P. F. Fischer, J. W. Lottes, and S. G. Kerkemeier.  Nek5000: Open source spectral element CFD solver.  Available at: http://nek5000.mcs.anl.gov, 2008.
  * Fukagata et al. (2002) K. Fukagata, K. Iwamoto, and N. Kasagi.  Contribution of reynolds stress distribution to the skin friction in wall-bounded flows.  _Phys. Fluids_ , 14(11):L73–L76, 2002.
  * Gad-el Hak (1994) M. Gad-el Hak.  Interactive control of turbulent boundary layers - A futuristic overview.  _AIAA J._ , 32(9):1753–1765, 1994.
  * Harun et al. (2013) Z. Harun, J. P. Monty, R. Mathis, and I. Marusic.  Pressure gradient effects on the large-scale structure of turbulent boundary layers.  _J. Fluid Mech._ , 715:477–498, 2013.
  * Huang et al. (1998) N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu.  The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis.  _Proc. R. Soc. Lond. A_ , 454:903–995, 1998.
  * Huang et al. (2008) Y. X. Huang, F. G. Schmitt, Z. M. Lu, and Y. L. Liu.  An amplitude-frequency study of turbulent scaling intermittency using empirical mode decomposition and Hilbert spectral analysis.  _Europhys. Lett._ , 84(4):40010, nov 2008.
  * Hwang (1996) D. Hwang.  A proof of concept experiment for reducing skin friction by using a micro-blowing technique.  _NASA_ , pages TM–107315, 1996.
  * Hwang (2004) D. Hwang.  Review of research into the concept of the microblowing technique for turbulent skin friction reduction.  _Prog. Aerosp. Sci._ , 40:559–575, 2004.
  * Kametani and Fukagata (2011) Y. Kametani and K. Fukagata.  Direct numerical simulation of spatially developing turbulent boundary layers with uniform blowing or suction.  _J. Fluid Mech._ , 681:154–172, 2011.
  * Kametani et al. (2015) Y. Kametani, K. Fukagata, R. Örlü, and P. Schlatter.  Effect of uniform blowing/suction in a turbulent boundary layer at moderate Reynolds number.  _Int. J. Heat Fluid Flow_ , 55:132–142, 2015.
  * Kim et al. (2002) K. Kim, H-.J. Sung, and M.-K. Chung.  Assessment of local blowing and suction in a turbulent boundary layer.  _AIAA J._ , 40(1):175–177, 2002.
  * Kornilov (2015) V. I. Kornilov.  Current state and prospects of researches on the control of turbulent boundary layer by air blowing.  _Prog. Aerosp. Sci._ , 76:1–23, 2015.
  * Kornilov (2017) V. I. Kornilov.  Implementation of air injection into the turbulent boundary layer of aircraft wing using external pressurized flow.  _Thermophys. Aeromech._ , 24:175–185, 2017.
  * Kornilov (2021) V. I. Kornilov.  Combined blowing/suction flow control on low-speed airfoils.  _Flow Turbul. Combust._ , 106:81–108, 2021.
  * Kornilov et al. (2019) V. I. Kornilov, I. N. Kavun, and A. N. Popkov.  Modification of turbulent airfoil section flow using a combined control action.  _Thermophys. Aeromech._ , 26:165–178, 2019.
  * Lee and Sung (2008) J.-H. Lee and H.-J. Sung.  Effects of an adverse pressure gradient on a turbulent boundary layer.  _Int. J. Heat Fluid Flow_ , 29(3):568–578, 2008\.
  * Li (2020) W.-P. Li.  Turbulence statistics of flow over a drag-reducing and a drag-increasing riblet-mounted surface.  _Aerosp. Sci. Technol._ , 104:106003, 2020.  ISSN 1270-9638.
  * Li et al. (2019) W.-P. Li, Y.-T. Fan, D. Modesti, and C. Cheng.  Decomposition of the mean skin-friction drag in compressible turbulent channel flows.  _J. Fluid Mech._ , 875:101–123, 2019.
  * Lumley (1967) J. L. Lumley.  The structure of inhomogeneous turbulent flows.  In A. M. Yaglom and V. I. Tartarsky, editors, _Atmospheric Turbulence and Radio Wave Propagation_ , pages 166–177. 1967.
  * Mahfoze et al. (2019) O. A. Mahfoze, A. Moody, A. Wynn, R. D. Whalley, and S. Laizet.  Reducing the skin-friction drag of a turbulent boundary-layer flow with low-amplitude wall-normal blowing within a Bayesian optimization framework.  _Phys. Rev. Fluids_ , 4:094601, 2019.
  * Mehdi and White (2011) F. Mehdi and C. M. White.  Integral form of the skin friction coefficient suitable for experimental data.  _Exp. Fluids_ , 50(1):43–51, 2011.
  * Mehdi et al. (2014) F. Mehdi, T. G. Johansson, C. M. White, and J. W. Naughton.  On determining wall shear stress in spatially developing two-dimensional wall-bounded flows.  _Exp. Fluids_ , 55(1):1656, 2014.
  * Modesti et al. (2018) D. Modesti, S. Pirozzoli, P. Orlandi, and F. Grasso.  On the role of secondary motions in turbulent square duct flow.  _J. Fluid Mech._ , 847:R1, 2018.  doi: 10.1017/jfm.2018.391.
  * Park and Choi (1999) J. Park and H. Choi.  Effects of uniform blowing or suction from a spanwise slot on a turbulent boundary layer flow.  _Phys. Fluids_ , 11(10):3095–3105, 1999.
  * Patera (1984) A. T. Patera.  A spectral element method for fluid dynamics: laminar flow in a channel expansion.  _J. Comput. Phys._ , 54:468–488, 1984.
  * Peet and Sagaut (2009) Y. Peet and P. Sagaut.  Theoretical prediction of turbulent skin friction on geometrically complex surfaces.  _Phys. Fluids_ , 21(10):105105, 2009.
  * Ran et al. (2021) W. Ran, A. Zare, and M. R. Jovanović.  Model-based design of riblets for turbulent drag reduction.  _J. Fluid Mech._ , 906:A7, 2021.  doi: 10.1017/jfm.2020.722.
  * Rastegari and Akhavan (2015) A. Rastegari and R. Akhavan.  On the mechanism of turbulent drag reduction with super-hydrophobic surfaces.  _J. Fluid Mech._ , 773:R4, 2015.  doi: 10.1017/jfm.2015.266.
  * Renard and Deck (2016) N. Renard and S. Deck.  A theoretical decomposition of mean skin friction generation into physical phenomena across the boundary layer.  _J. Fluid Mech._ , 790:339–367, 2016.
  * Rotta (1950) J. C. Rotta.  Über die theorie der turbulenten grenzschichten.  _Mitt. Max Planck Inst. Strömungsforsch., Göttingen_ , 1950\.
  * Sanmiguel Vila et al. (2020) C. Sanmiguel Vila, R. Vinuesa, S. Discetti, A. Ianiro, P. Schlatter, and R. Örlü.  Separating adverse-pressure-gradient and Reynolds-number effects in turbulent boundary layers.  _Phys. Rev. Fluids_ , 5:064609, 2020.
  * Schlatter and Örlü (2012) P. Schlatter and R. Örlü.  Turbulent boundary layers at moderate Reynolds numbers: inflow length and tripping effects.  _J. Fluid Mech._ , 710:5–34, 2012.
  * Stroh et al. (2015) A. Stroh, B. Frohnapfel, P. Schlatter, and Y. Hasegawa.  A comparison of opposition control in turbulent boundary layer and turbulent channel flow.  _Phys. Fluids_ , 27:075101, 2015.
  * Stroh et al. (2016) A. Stroh, Y. Hasegawa, P. Schlatter, and B. Frohnapfel.  Global effect of local skin friction drag reduction in spatially developing turbulent boundary layer.  _J. Fluid Mech._ , 805:303–321, 2016.  doi: 10.1017/jfm.2016.545.
  * Tanarro et al. (2020) Á. Tanarro, R. Vinuesa, and P. Schlatter.  Effect of adverse pressure gradients on turbulent wing boundary layers.  _J. Fluid Mech._ , 883:A8, 2020.
  * Touber and Leschziner (2012) E. Touber and M. A. Leschziner.  Near-wall streak modification by spanwise oscillatory wall motion and drag-reduction mechanisms.  _J. Fluid Mech._ , 693:150–200, 2012.  doi: 10.1017/jfm.2011.507.
  * Vinuesa and Schlatter (2017) R. Vinuesa and P. Schlatter.  Skin-friction control of the flow around a wing section through uniform blowing.  In _Proceedings of European Drag Reduction and Flow Control Meeting (EDRFCM)_ , 2017.
  * Vinuesa et al. (2018) R. Vinuesa, P. S. Negi, M. Atzori, A. Hanifi, D. S. Henningson, and P. Schlatter.  Turbulent boundary layers around wing sections up to ${R}e_{c}$=1,000,000.  _Int. J. Heat Fluid Flow_ , 72:86 – 99, 2018.
  * Wang et al. (2018) W.-K. Wang, C. Pan, and J.-J. Wang.  Quasi-bivariate variational mode decomposition as a tool of scale analysis in wall-bounded turbulence.  _Exp. Fluids_ , 59:1, 2018.
  * Wang et al. (2019) W.-K. Wang, C. Pan, and J.-J. Wang.  Multi-component variational mode decomposition and its application on wall-bounded turbulence.  _Exp. Fluids_ , 60:95, 2019.
  * Wei (2018) T. Wei.  Integral properties of turbulent-kinetic-energy production and dissipation in turbulent wall-bounded flows.  _J. Fluid Mech._ , 854:449–473, 2018.
  * Welch et al. (2001) G. Welch, L. Larosiliere, D. Hwang, and J. Wood.  _Effectiveness of micro-blowing technique in adverse pressure gradients_.  2001\.
  * White and Mungal (2008) C. M. White and M. G. Mungal.  Mechanics and prediction of turbulent drag reduction with polymer additives.  _Annu. Rev. Fluid Mech._ , 40(1):235–256, 2008\.
  * Wu and Christensen (2010) Y. Wu and K. T. Christensen.  Spatial structure of a turbulent boundary layer with irregular surface roughness.  _J. Fluid Mech._ , 655:380–418, 2010.
  * Yoon et al. (2016) M. Yoon, J. Ahn, J. Hwang, and H. J. Sung.  Contribution of velocity-vorticity correlations to the frictional drag in wall-bounded turbulent flows.  _Phys. Fluids_ , 28(8):081702, 2016.
  * Yoon et al. (2018) M. Yoon, J. Hwang, and H. J. Sung.  Contribution of large-scale motions to the skin friction in a moderate adverse pressure gradient turbulent boundary layer.  _J. Fluid Mech._ , 848:288–311, 2018.


# Parameter Estimation for multi-mixed Fractional Ornstein–Uhlenbeck Processes
by Generalized Method of Moments

Hamidreza Maleki Almani  School of Technology and Innovations, University of
Vaasa, P.O. Box 700, FIN-65101 Vaasa, FINLAND<EMAIL_ADDRESS>and  Tommi
Sottinen  School of Technology and Innovations, University of Vaasa, P.O. Box
700, FIN-65101 Vaasa, FINLAND<EMAIL_ADDRESS>
(Date: September 3, 2024)

###### Abstract.

We develope the generalized method of moments (GMM) estimation for the
parameters of the finitely mixed multi-mixed fractional Ornstein–Uhlenbeck
(mmfOU) processes, and analyze the consistency and assymptotic normality of
this estimator. We also include some simulations and provide numerical
observations considering different statistical errors.

###### Key words and phrases:

generalized method of moments, fractional Langevin equation, multi-mixed
fractional Ornstein–Uhlenbeck process, parameter estimation

###### 2020 Mathematics Subject Classification:

60G10, 60G15, 60G22, 62M09

##  1\. Introduction and Preliminaries

Langevin [13] postulated the model

(1) $\frac{\mathrm{d}U_{t}}{\mathrm{d}t}=-\frac{f}{m}U_{t}+\frac{F_{t}}{m},$

to the velocity $U$ of a free particle with mass $m>0$, moving in a liquid by
the alternating force $F$ of surrounding molecules, and a constant friction
coefficient $f$. By imposing probabilistic hypotheses to the force $F$,
Ornstein and Uhlenbeck [18] derived that the solution of (1) has a Gaussian
distribution with exponential mean. Later on, considering the initial variable
$U_{0}$ to be central Gaussian and independent from $F$, Doob [5] showed that
the solution to (1) is stationary and for some constant $c>0$, the time-
changed scaled process $t^{1/2}\,U(c\log t)$ is in fact the well-known
Einstein’s [6] Brownian motion for $t\geq 0$. As the Brownian motion in
nowhere differentiable, Doob understood (1) as

(2) $\mathrm{d}U_{t}=-\lambda U_{t}\mathrm{d}t+\mathrm{d}X_{t},$

where $U_{0}=\xi$ is the initial variable and $X$ is a Lévy process. He
derived the solution

(3) $U_{t}=\mathrm{e}^{-\lambda
t}\xi+\int_{0}^{t}\mathrm{e}^{-\lambda(t-s)}\,\mathrm{d}X_{s},$

which is called the Ornstein–Uhlenbeck (OU) process. This solution was later
extent for all cases that $X$ is a semimartingale (see e.g. [17]). Lamperti
[12] proved the processes $V_{t},t\geq 0$ is $H$-selfsimilar if and only if
$U_{t}=e^{-\lambda t}V(c\exp(\lambda t/H))$ is stationary for all
$\lambda,c>0$. Cheridito and et al. [4] applied Lamperti transform and
verified the solution (3) is valid for $X=B^{H}$, where $B^{H}$ is the
fractional Brownian mtoion (fBm) introduced by Mandelbrot and Van Ness [14].

In [2] the Langevin equation (2) with
$X=\sum_{i=1}^{\infty}\sigma_{i}B^{H_{i}}$, the multi-mixed fractional
Brownian motion (mmfBm), is the dynamic of particle’s velocity when the liquid
is not purely homogeneous, i.e. the local surrounding molecules can have
imposing forces of different roughnesses. That is

(4) $\mathrm{d}U_{t}=-\lambda
U_{t}\mathrm{d}t+\mathrm{d}\left(\sum_{i=1}^{\infty}\sigma_{i}B^{H_{i}}_{t}\right).$

They confirmed the solution (3) for this model is possible if for all $i\geq
1$, $H_{i}\in[H_{\inf},H_{\sup}]\subset(0,1)$ and
$\sum_{i=1}^{\infty}\sigma_{i}^{2}<\infty$. The solution is

(5) $U_{t}=\mathrm{e}^{-\lambda
t}\xi+\int_{0}^{t}\mathrm{e}^{-\lambda(t-s)}\,\mathrm{d}\left(\sum_{i=1}^{\infty}\sigma_{i}B^{H_{i}}_{t}\right).$

When $U_{0}=\xi=\sum_{i}\sigma_{i}\xi_{i}$ is the initial condition, this
solution can also be written as

(6) $U_{t}=\sum_{i=1}^{\infty}\sigma_{i}U^{\lambda,H_{i}}_{t},$

where

(7) $\mathrm{d}U^{\lambda,H_{i}}_{t}=-\lambda
U^{\lambda,H_{i}}_{t}\mathrm{d}t+\mathrm{d}B^{H_{i}}_{t},$

and $U^{\lambda,H_{i}}_{0}=\xi_{i}$. To approximate $U$, one can naturally
consider

(8) $U^{n}_{t}=\sum_{i=1}^{n}\sigma_{i}U^{\lambda,H_{i}}_{t},$

that is convergent to $U$ in $L^{2}(\Omega\times[0,T])$, for any $T\geq 0$
(see [2]). To apply the process (8) for modelling the real data, we need
$2n+1$ real-valued parameters

$\lambda;\,H_{1},\ldots,H_{n};\,\sigma_{1},\ldots,\sigma_{n}$

to be statistically estimated. To do this, we use the generalized method of
moments (GMM) introduced in Barbiza and Viens [3].

The outline of this paper is as follows. In Section 2 we develop the GMM
estimator for the process (8) analytically. Next, in Section 3, we investigate
the consistency of the introduced estimator. In Section 4 we prove the
asymptotic normality of our estimator. Finally, in Section 5, we provide some
simulations.

##  2\. Generalized Methdod of Moments Estimation

The parameters we aim to estimate are

$\lambda>0,\qquad
H_{k}\in[H_{\inf},H_{\sup}]\subset(0,1),\qquad\sigma_{k}>0;\qquad k=1,2,\ldots
n,$

i.e. the following vector will be estimated

$\displaystyle\theta=(\lambda;\mathbf{H};\bm{\sigma})^{\prime}\in\mathbb{R}^{2n+1}$
where $\displaystyle\mathbf{H}=(H_{1},\ldots,H_{n}),$
$\displaystyle\bm{\sigma}=(\sigma_{1},\ldots,\sigma_{n}).$

To develope GMM estimator, the autocovariance function of $U^{n}$ is required.
It is given in [2] by

$\displaystyle\rho_{\theta}(t)=\rho_{\lambda,n}(t)=\mathbb{E}[U^{n}_{s}U^{n}_{s+t}]$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{n}\sigma_{k}^{2}\frac{\Gamma(1+2H_{k})}{4}\frac{\mathrm{e}^{-\lambda
t}}{\lambda^{2H_{k}}}\bigg{\\{}1+{\gamma}_{2H_{k}-1}(\lambda
t)+\mathrm{e}^{2\lambda t}{\Gamma}_{2H_{k}-1}(\lambda t)\bigg{\\}}$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{n}\sigma_{k}^{2}\frac{\Gamma(2H_{k}+1)}{2\lambda^{2H_{k}}}\,\bigg{\\{}\cosh(\lambda
t)-\frac{1}{\Gamma(2H_{k})}\int_{0}^{\lambda t}s^{2H_{k}-1}\cosh(\lambda
t-s)\,ds\bigg{\\}},$

where for $\alpha\in(-1,0)\cup(0,1)$

$\displaystyle{\gamma}_{\alpha}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha)}\int_{0}^{x}s^{\alpha-1}\mathrm{e}^{s}\,\mathrm{d}s,$
$\displaystyle{\Gamma}_{\alpha}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha)}\int_{x}^{\infty}s^{\alpha-1}\mathrm{e}^{-s}\,\mathrm{d}s,$

and ${\gamma}_{0}(x)=1$, ${\Gamma}_{0}(x)=0$. Here ${\gamma}_{\alpha}$,
${\Gamma}_{\alpha}$ are the normalized incomplete gamma functions (see [1]).
The formula (2) is continuous in $t$, and so $U^{n}$ is ergodic (see [8], page
257). By Lemma 5.2 in [11] for sufficiently small $t\searrow 0$ we have

(10)
$\rho_{\theta}(t)\approx\sum_{k=1}^{n}\sigma_{k}^{2}\left(\frac{H\Gamma(2H)}{\lambda^{2H}}-\frac{t^{2H}}{2}\right).$

Further, for our estimation, we need some filters associated to the parameter
$\theta$ and a positive integers $L\geq 2n+1$. We define them as following.

######  Definition 1.

A filter of length $L+1$ and order $l$ is a sequence of real numbers
$\mathbf{a}=(a_{k})_{k=0}^{L}=(a_{0},\ldots,a_{L})$, that for $l>0$

$\displaystyle\sum_{q=0}^{L}a_{q}q^{r}=0\,;\quad 0\leq r\leq l-1,\
r\in\mathbb{Z},$ $\displaystyle\sum_{q=0}^{L}a_{q}q^{l}\neq 0.$

When $l=0$, we put $a_{0}=1$ and $a_{q}=0$ for $0<q\leq L$.

The GMM estimation requires a family of filters by the length $L$ and orders
$l=0,1,\ldots,L$. For each filter $\mathbf{a}$ with order $l$ define

$\mathbf{b}=(b_{0},\ldots,b_{L})^{\prime},$

where

$\displaystyle b_{0}:=\sum_{q=0}^{L}a_{q}^{2},$ $\displaystyle
b_{k}:=2\sum_{q=0}^{L-k}a_{q+k}a_{q}\,;\quad k=1,\ldots,L.$

For different orders $l_{i}\neq l_{j}$, the corresponding vectors
$\mathbf{b}_{i}$ and $\mathbf{b}_{j}$ are linearly independent. So, we can
take $L$ filters $\mathbf{a}_{1},\ldots,\mathbf{a}_{L}$ with respective orders
$l_{1},\ldots,l_{L}$ such that the $L\times(L+1)$ matrix $B=[\mathbf{b}_{1}\
\mathbf{b}_{2}\ \cdots\ \mathbf{b}_{L}]^{\prime}$ satisfies
$2n+1\leq\mbox{rank}(B)\leq L$. Then, for a fixed step size $\alpha>0$, and
$\theta\in\Theta$ (a closed subset of $\mathbb{R}^{2n+1}$), we can define the
function

$V(\theta):=\sum_{k=0}^{L}b_{k}\rho_{\theta}(\alpha
k)=\sum_{k=0}^{L}b_{k}\rho_{\lambda,n}(\alpha k),$

and the filtered process of the step size $\alpha>0$ for $t\geq 0$

$\varphi(t):=\sum_{q=0}^{L}a_{q}U^{n}_{t-q\alpha}.$

Next, let $\theta_{0}$ be the exact value of the parameter. Due to the
stationarity of $U^{n}$, the random variable $\varphi(t)$ is centered Gaussian
with variance $\mathbb{E}[\varphi(t)^{2}]=V(\theta_{0})$. Following the
notation of [15], for any $\theta\in\Theta$ we have the vector of differences

$\mathbf{g}(t,\theta):=\left(g_{1}(t,\theta),\ldots,g_{L}(t,\theta)\right)^{\prime},$

where

$g_{\ell}(t,\theta)=\varphi_{\ell}(t)^{2}-V_{\ell}(\theta),\quad 1\leq\ell\leq
L,$

and the subscript $\ell$ means that we use the filter $\mathbf{a}_{\ell}$ in
the filtered process and variance. Moreover, $\mathbf{g}(t,\theta)$ satisfies
a population moment condition (see [9] and [15])

$\mathbb{E}[\mathbf{g}(t,\theta_{0})]=0,\quad t\geq 0.$

Also, the vector $\mathbf{g}$ is in the second Wiener chaos (see [16]).

Further, assume that we have observed the process $U^{n}$ at discrete
equidistant times $0=t_{0}<t_{1}<\cdots<t_{N-1}<t_{N}=T$ with fixed time step
$\alpha=t_{i}-t_{i-1}$. Let $A$ be an $L\times L$ symmetric positive-definite
matrix (a suitably chosen $A$ ensures that the GMM estimate is efficient, see
[3] Section 2.3). Now, for $\theta\in\Theta$ and time $t\geq 0$ denote

  1. $\mathbf{g}_{0}(\theta):=\mathbb{E}[\mathbf{g}(t,\theta)]$ (vector of expected differences),

  2. $\hat{\mathbf{g}}_{N}(\theta):=\frac{1}{N-L+1}\sum_{i=L}^{N}\mathbf{g}(t_{i},\theta)$ (vector of sampled differences),

  3. $Q_{0}(\theta):=\mathbf{g}_{0}(\theta)^{\prime}A\mathbf{g}_{0}(\theta)$ (squared distance of expected differences),

  4. $\hat{Q}_{N}(\theta):=\hat{\mathbf{g}}_{N}(\theta)^{\prime}A\hat{\mathbf{g}}_{N}(\theta)$ (sampled version of the above distance).

Note that $\mathbf{g}_{0}(\theta)$ does not depend on time $t$, due to the
stationarity of $U^{n}$. Finally, we define the GMM estimator of
$\theta_{0}=(\lambda_{0};\mathbf{H}_{0};\bm{\sigma}_{0}^{2})$ as

(11)
$\hat{\theta}_{N}=(\hat{\lambda}_{N};\hat{\mathbf{H}}_{N};\hat{\bm{\sigma}}_{N})^{\prime}=\arg\min_{\theta\in\Theta}\hat{Q}_{N}(\theta).$

In fact, if the function $Q_{0}(\theta)$ attains a unique zero at
$\theta_{0}$, then we can find a value $\hat{\theta}_{N}$ such that the
approximation of $Q_{0}$, that is $\hat{Q}_{N}$, should be the smallest
possible (see [9, 10]).

Moreover, all the remaining results are still valid if we substitute the fixed
matrix $A$ with a sequence of random matrices (perhaps depending on the data)
with a deterministically bounded eigenstructure. In particular, we can choose
such sequence in order to attain convergence in probability to the efficient
alternative of $A$ (see [3], Section 2.3).

In the following we will show the consistency and asymptotic normality of the
GMM estimator (11).

##  3\. Consistency

As mentioned in [3], for the consistency, it is sufficient to find conditions
under which the mapping $\theta\mapsto\bm{\rho}_{\theta,2n+1}(\alpha)$ is
smooth and invertible at least in a closed subset
$\theta\in\Theta\subseteq\mathbb{R}^{2n+1}$, where

$\bm{\rho}_{\theta,2n+1}(\alpha)=(\rho_{\theta}(0\cdot\alpha),\ldots,\rho_{\theta}((2n)\cdot\alpha))^{\prime},$

and $\alpha>0$ is the step size. This is equivalent for
$\nabla_{\theta}\bm{\rho}_{\theta,2n+1}(\alpha)$ to be non-singular in a
closed subset $\theta\in\Theta\subseteq\mathbb{R}^{2n+1}$. To show this for
the mmfOU processes, with
$\theta=(\lambda;\mathbf{H};\bm{\sigma})^{\prime}\in\mathbb{R}^{2n+1}$, we
have the next lemma and theorem.

######  Lemma 1.

In the closed rectangle

(12)
$\Upsilon_{3}=\Bigg{\\{}\theta=(\lambda,H,\sigma)^{\prime}\,\Bigg{|}\,\begin{matrix}\lambda,\sigma>0,\,H\in(0,1),\\\
\lambda<\exp\big{(}\Psi(2H+1)\big{)}\end{matrix}\Bigg{\\}}\subset\Theta\subseteq\mathbb{R}^{3},$

where $\Psi(\cdot)$ is the so called digamma function

$\Psi(z)=\frac{\mathrm{d}}{\mathrm{d}z}\log\Gamma(z)=\frac{\Gamma^{\prime}(z)}{\Gamma(z)},$

for $\rho_{\theta}(t)=\sigma^{2}\rho_{\lambda,H}(t)$, the matrix

$\nabla_{\theta}\bm{\rho}_{\theta,3}(\alpha)=\left[\begin{array}[]{lll}\frac{\partial\rho_{\theta}}{\partial\lambda}(0)&\frac{\partial\rho_{\theta}}{\partial
H}(0)&\frac{\partial\rho_{\theta}}{\partial\sigma}(0)\\\ &&\\\
\frac{\partial\rho_{\theta}}{\partial\lambda}(\alpha)&\frac{\partial\rho_{\theta}}{\partial
H}(\alpha)&\frac{\partial\rho_{\theta}}{\partial\sigma}(\alpha)\\\ &&\\\
\frac{\partial\rho_{\theta}}{\partial\lambda}(2\alpha)&\frac{\partial\rho_{\theta}}{\partial
H}(2\alpha)&\frac{\partial\rho_{\theta}}{\partial\sigma}(2\alpha)\\\
\end{array}\right]$

is a P-matrix for sufficiently small $\alpha$.

###### Proof.

Proof is similar to Lemma B.2 in [3]. ∎

######  Theorem 1.

In the closed rectangle

(13)
$\Upsilon_{2n+1}=\Bigg{\\{}\theta=(\lambda;\mathbf{H};\bm{\sigma})^{\prime}\,\Bigg{|}\,\begin{matrix}\lambda,\sigma_{k}>0,\,H_{k}\in(0,1),\\\
\lambda<\exp\big{(}\Psi(2H_{\inf}+1)\big{)}\end{matrix}\Bigg{\\}}\subset\Theta\subseteq\mathbb{R}^{2n+1},$

the mapping

(14)
$\bm{\rho}_{\theta,2n+1}(\alpha)=(\rho_{\theta}(0\cdot\alpha),\ldots,\rho_{\theta}((2n)\cdot\alpha))^{\prime},$

is injective for sufficiently small $\alpha$.

###### Proof.

Let $\theta\in\Upsilon_{2n+1}$, and $x=(x_{1},\ldots,x_{2n+1})^{\prime}>0$ is
a positive vector in $\mathbb{R}^{2n+1}$. Then, the $k$th component of the
vector $\nabla_{\theta}\bm{\rho}_{\theta,2n+1}(\alpha)\,x$ is

(15)
$\displaystyle\Big{[}\nabla_{\theta}\bm{\rho}_{\theta,2n+1}(\alpha)\,x\Big{]}_{k}=\sum_{j=1}^{2n+1}\frac{\partial\rho_{\lambda,n}}{\partial\theta_{j}}((k-1)\alpha)x_{j}$
$\displaystyle=\frac{\partial\rho_{\lambda,n}}{\partial\lambda}((k-1)\alpha)\,x_{1}+\sum_{j=1}^{n}\frac{\partial\rho_{\lambda,n}}{\partial
H_{j}}((k-1)\alpha)\,x_{j+1}$
$\displaystyle+\sum_{j=1}^{n}\frac{\partial\rho_{\lambda,n}}{\partial\sigma_{j}}((k-1)\alpha)\,x_{n+j+1}$
$\displaystyle=\sum_{j=1}^{n}\sigma_{j}^{2}\frac{\partial\rho_{\lambda,H_{j}}}{\partial\lambda}((k-1)\alpha)\,x_{1}+\sum_{j=1}^{n}\sigma_{j}^{2}\frac{\partial\rho_{\lambda,H_{j}}}{\partial
H_{j}}((k-1)\alpha)\,x_{j+1}$
$\displaystyle+\sum_{j=1}^{n}2\sigma_{j}\frac{\partial\rho_{\lambda,H_{j}}}{\partial\sigma_{j}}((k-1)\alpha)\,x_{n+j+1}$
$\displaystyle=\sum_{j=1}^{n}\bigg{[}\sigma_{j}^{2}\frac{\partial\rho_{\lambda,H_{j}}}{\partial\lambda}((k-1)\alpha)\,x_{1}+\sigma_{j}^{2}\frac{\partial\rho_{\lambda,H_{j}}}{\partial
H_{j}}((k-1)\alpha)\,x_{j+1}$
$\displaystyle+2\sigma_{j}\frac{\partial\rho_{\lambda,H_{j}}}{\partial\sigma_{j}}((k-1)\alpha)\,x_{n+j+1}\bigg{]}.$

Now, for all $j=1,\ldots,2n+1$, let
$\tilde{\rho}_{j}(t)=\sigma_{j}^{2}\rho_{\lambda,H_{j}}(t)$ and

$\tilde{\alpha}=\bigg{\\{}\begin{array}[]{lll}(k-1)\alpha&;&1\leq k\leq n,\\\
(k-1)\alpha/2&;&n+1\leq k\leq 2n+1.\end{array}$

Then, since $\exp(\Psi(\cdot))$ is an strictly increasing function,

$\lambda<\exp\big{(}\Psi(2H_{\inf}+1)\big{)}=\min_{1\leq i\leq
2n+1}\exp\big{(}\Psi(2H_{i}+1)\big{)}\leq\exp\big{(}\Psi(2H_{j}+1)\big{)},$

and so by Lemma 1 the matrix

(16)
$\nabla_{{}_{\lambda,H_{j},\sigma_{j}}}\tilde{\rho}_{j}(\tilde{\alpha})=\left[\begin{array}[]{*3{>{\displaystyle}c}p{5cm}}\frac{\partial\tilde{\rho}_{j}}{\partial\lambda}(0)&\frac{\partial\tilde{\rho}_{j}}{\partial
H_{j}}(0)&\frac{\partial\tilde{\rho}_{j}}{\partial\sigma_{j}}(0)\\\
\hfil&\hfil&\hfil\\\
\frac{\partial\tilde{\rho}_{j}}{\partial\lambda}(\tilde{\alpha})&\frac{\partial\tilde{\rho}_{j}}{\partial
H_{j}}(\tilde{\alpha})&\frac{\partial\tilde{\rho}_{j}}{\partial\sigma_{j}}(\tilde{\alpha})\\\
\hfil&\hfil&\hfil\\\
\frac{\partial\tilde{\rho}_{j}}{\partial\lambda}(2\tilde{\alpha})&\frac{\partial\tilde{\rho}_{j}}{\partial
H_{j}}(2\tilde{\alpha})&\frac{\partial\tilde{\rho}_{j}}{\partial\sigma_{j}}(2\tilde{\alpha})\\\
\end{array}\right]$

is a P-matrix for all $\tilde{\alpha}\leq\tilde{\alpha}_{kj}$ sufficiently
small associated to $k$ and $\tilde{\rho}_{j}$. So, (16) is also P-matrix for

$\tilde{\alpha}\leq\tilde{\alpha}^{*}_{2n+1}:=\min_{\begin{matrix}{\scriptstyle
1\leq k\leq 2n+1}\\\ {\scriptstyle 1\leq j\leq
2n+1}\end{matrix}}\tilde{\alpha}_{kj},$

for all $j=1,\ldots,2n+1$.
Hence, by the Theorem 2 of [7], for the subvector
$\tilde{x}_{j}=(x_{1},x_{j+1},x_{n+j+1})^{\prime}$ we have

$\nabla_{{}_{\lambda,H_{j},\sigma_{j}}}\tilde{\rho}_{j}(\tilde{\alpha})\cdot\tilde{x}_{j}>0,$

and so for $\theta\in\Upsilon_{2n+1}$

$\displaystyle\sigma_{j}^{2}\frac{\partial\rho_{\lambda,H_{j}}}{\partial\lambda}((k-1)\alpha)\,x_{1}+\sigma_{j}^{2}\frac{\partial\rho_{\lambda,H_{j}}}{\partial
H_{j}}((k-1)\alpha)\,x_{j+1}$
$\displaystyle+2\sigma_{j}\frac{\partial\rho_{\lambda,H_{j}}}{\partial\sigma_{j}}((k-1)\alpha)\,x_{n+j+1}>0.$

Consequently, for

(17) $\alpha\leq\frac{\tilde{\alpha}^{*}_{2n+1}}{2n+1},$

(15) is positive for all $k=1,\ldots,2n+1$, and so the vector
$\nabla_{\theta}\bm{\rho}_{\theta,2n+1}(\alpha)\,x$ is positive. Hence, by the
Theorem 2 of [7], the gradiant
$\nabla_{\theta}\bm{\rho}_{\theta,2n+1}(\alpha)$ is a P-matrix, and so the
mapping (14) is injective on $\Upsilon_{2n+1}$. ∎

Now, we can obtain the strong consistency for the GMM estimators for the mmfOU
processes.

######  Theorem 2.

For the GMM parameter estimator (11) for the mmfOU process $U^{n}$, we have

$\hat{\theta}_{N}=(\hat{\lambda}_{N};\hat{\mathbf{H}}_{N};\hat{\bm{\sigma}}_{N})^{\prime}\overset{a.s.}{\longrightarrow}\theta_{0}=(\lambda_{0};\mathbf{H}_{0};\bm{\sigma}_{0}).$

###### Proof.

As an straightforward result of Theorem 1, Assumption 2.1 of [3] is valid and
we can apply Theorem 2.1 in [3]. ∎

######  Remark 1.

First, we note that for $\lambda\in[0,\mathrm{e}^{\Psi(1)}]$ the consistency
that we studied above is valid. Second, the proof of Theorem 1 returns an
important result as the nature of numerical estimation. By (17), the larger
the number of parameters we aim to estimate, the smaller the sufficient step
$\alpha>0$ we need. However, in any case, there exists some sufficient
$\alpha$ for which the consistency is garanteed.

##  4\. Asymptotic Normality

For $\theta=(\lambda;\mathbf{H};\bm{\sigma})^{\prime}\in\mathbb{R}^{2n+1}$,
the spectral density function (SDF) of $U^{n}$ is

(18)
$f_{\theta}(x)=f(\theta;x)=\sum_{i=1}^{n}\sigma_{i}^{2}f_{\lambda,H_{i}}(x),$

where

(19) $f_{\lambda,H_{i}}(x)=\frac{\sin(\pi
H_{i})\Gamma(1+2H_{i})}{2\pi}\,\frac{|x|^{1-2H_{i}}}{x^{2}+\lambda^{2}},$

is the SDF of the fOU process $U^{\lambda,H_{i}}$ (see [2]). Denote

$\bar{f}_{\theta}(x):=\sum_{p\in\mathbb{Z}}f_{\theta}\left(x+\frac{2p\pi}{\alpha}\right)\,;\quad
x\in\left[-\frac{\pi}{\alpha},\frac{\pi}{\alpha}\right],$

and for $\l_{i},\l_{j}\in\\{1,\ldots,2n+1\\}$

$R_{\theta}(x|i,j):=\min(1,|x|^{l_{i}+l_{j}})\bar{f}_{\theta}(x).$

Then, as explained in [3], the GMM estimator has asymptotic normality if
$R_{\theta_{0}}(\cdot|i,j)\in L^{2}(-\pi,\pi)$. For the process $U^{n}$ we
have

$\displaystyle R_{\theta}(x|i,j)$ $\displaystyle=$
$\displaystyle\min(1,|x|^{l_{i}+l_{j}})\overline{f}_{\theta}(x)$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{n}\sigma_{k}^{2}R_{\lambda,H_{k}}(x|i,j),$

where

$\displaystyle
R_{\lambda,H_{k}}(x|i,j):=\min(1,|x|^{l_{i}+l_{j}})\overline{f}_{\lambda,H_{k}}(x),$
$\displaystyle\bar{f}_{\lambda,H_{k}}(x):=\sum_{p\in\mathbb{Z}}f_{\lambda,H_{k}}\left(x+\frac{2p\pi}{\alpha}\right)\,;\quad
x\in\left[-\frac{\pi}{\alpha},\frac{\pi}{\alpha}\right].$

So, $R_{\theta}(\cdot|i,j)\in L^{2}(-\pi,\pi)$ if and only if
$R_{\lambda,H_{k}}(\cdot|i,j)\in L^{2}(-\pi,\pi)$ for all $k=1,\ldots,n$. The
following Lemma is a generalization of Lemma 3.1 in [3].

######  Lemma 2.

For the $U^{n}$ process $R_{\theta_{0}}(\cdot|i,j)\in L^{2}(-\pi,\pi)$ for

  * (i)

$l_{i}+l_{j}\geq 1$, if all $H_{k}\in(0,1)$,

  * (ii)

$l_{i}+l_{j}=0$, if all $H_{k}\in(0,\frac{3}{4})$, i.e. $H_{\sup}\leq 3/4$,

where $n\geq 1$ ($L\geq 3$).

Now, for $\theta\in\Theta$ we define

$\displaystyle\hat{G}_{N}(\theta):=\nabla_{\theta}\hat{\mathbf{g}}_{N}(\theta),$
$\displaystyle
G(\theta):=\mathbb{E}[\nabla_{\theta}\mathbf{g}(\cdot,\theta)],$

and note that $G(\theta)$ does not depend on time $t$ because of the
stationarity of $U^{n}$.

Finally, we have the following theorem for asymptotic normality of the GMM
estimator of $U^{n}$.

######  Theorem 3.

The GMM estimator (11) for the mmfOU processes $U^{n}$ satisfies

  * (i)

$\mathbb{E}[\|\hat{\theta}_{N}-\theta_{0}\|^{2}]=O(N^{-1})$,

  * (ii)

$N^{c}\|\hat{\theta}_{N}-\theta_{0}\|\overset{a.s.}{\longrightarrow}0$ for
$c<\frac{1}{2}$,

  * (iii)

$\sqrt{N}(\hat{\theta}_{N}-\theta_{0})\overset{d}{\longrightarrow}\mathscr{N}\left(0,C(\theta_{0})\Lambda
C(\theta_{0})^{\prime}\right)$,

where
$C(\theta_{0})=[G(\theta_{0})^{\prime}AG(\theta_{0})]^{-1}G(\theta_{0})^{\prime}A$,
and $\Lambda$ is a $L\times L$ matrix (for $L\geq 2n+1$) such that

$\displaystyle\Lambda_{ij}=2\sum_{p\in\mathbb{Z}}\left[\sum_{k=0}^{L}b_{k}^{i,j}\rho_{\theta_{0}}\left(\alpha\cdot(k+p)\right)\right]^{2},$
$\displaystyle
b_{k}^{i,j}=\sum_{q=0}^{L-k}a_{q}^{i}a_{q+k}^{j}+\sum_{q=0}^{L-k}a_{q}^{j}a_{q+k}^{i}.$

###### Proof.

Theorem 1 and Lemma 2 yield that in fact the Assumptions 1.2-1.3 of [3] are
valid. So, we can apply Lemma 2.2 and Theorem 2.3 of [3]. ∎

##  5\. Simulation

In this section, we perform a simulation study to test our GMM estimator for
the parameters vector $\theta=(\lambda;\mathbf{H};\bm{\sigma})$ for $n=3$. The
statistical error metrics we use are

$\displaystyle\widehat{MSE}$
$\displaystyle:=\frac{1}{m}\sum_{i=1}^{m}\|\hat{\theta}_{N,i}-\theta\|^{2},$
$\displaystyle e(\widehat{Var})$ $\displaystyle:=\text{maximum eigenvalue of
}\widehat{Var}(\theta_{N}),$ $\displaystyle\widehat{Bias}^{2}$
$\displaystyle:=\left\|\frac{1}{m}\sum_{i=1}^{m}\hat{\theta}_{N,i}-\theta\right\|^{2},$

where $\widehat{Var}(\theta_{N})$ is the empirical covariance matrix of the
estimated parameters based on the $m$ replications. For statistically high
accuracy (cooking), we setup $m=500$ replications in our estimations.

To have the consistency, as mentioned by Remark 1, we take
$\lambda=\exp(\Psi(1))$. For (17) to hold we take $N=200,600,$ and $1000$ in
time interval $[0,1]$. Also, Theorem 1 is valid when $\bm{\sigma}$ has a fixed
sign. So, we choose the positive-sign values
$\bm{\sigma}=(0.5,1,1.5)^{\prime}$.

For asymptotic normality, we are restricted by the Lemma 2. So, we take
filters of positive orders $l=1,2,\ldots,2n+1$. On the other hand, we are
interested to test our estimator for varying Hurst values. Specially, for
applications sense, we are interested to include the standard Brownian motion
(Bm), the case $H_{j}=1/2$ for some $j$. Hence, we take
$\mathbf{H}=(0.3,0.5,0.7)^{\prime}$.

In each estimation, we first simulate replications for the $U^{n}$ process,
and then estimate its parameters $\theta=(\lambda;\mathbf{H};\bm{\sigma})$. To
produce $U^{n}$, we use iterations of the Langevin equation

(20) $\mathrm{d}U^{n}_{t}=-\lambda
U^{n}_{t}\mathrm{d}t+\mathrm{d}\left(\sum_{i=1}^{n}\sigma_{i}B^{H_{i}}_{t}\right),$

at $N$ points of time interval $[0,1]$ with uniform step-sizes. In each
replication, we use the autocovariance function $\rho_{\lambda,n}$ given in
(10) to estimate the parameters
$\hat{\theta}_{N}=(\hat{\lambda}_{N};\hat{\mathbf{H}}_{N};\hat{\bm{\sigma}}_{N})$.
To do this, we employ finite difference filters of orders $L=2n+1$ (number of
parameters) with $A=I_{2n+1}$ the identity matrix to minimize the calculation
cost (see [3]).

The resulting values of the statistical error metrics for our simulations are
given in Table 1. As one can check, the method produces reasonably accurate
estimations and the errors vanishe while inreasing the time step accuracy by
increasing $N$. This fact is nicely observable in Figure 1, confirming our
analytical formula (17) in action.

$N$ | 200 | 600 | 1000
---|---|---|---
$\widehat{MSE}$ |  0.006867 |  0.001149 |  0.000558
$e(\widehat{Var})$ | 0.004257 | 0.000564 | 0.000242
$\widehat{Bias}^{2}$ | 0.00082 | 0.000294 | 0.000184
Table 1. Statistical errors for estimations by $m=500$ replications. Figure 1.
Estimations of $\theta=(\lambda;\mathbf{H};\bm{\sigma})$ for $m=500$
replications.

## References

  * [1] M. Abramowitz, I. A. Stegun, and R. H. Romer, Handbook of mathematical functions with formulas, graphs, and mathematical tables, 1988.
  * [2] H. M. Almani and T. Sottinen, Multi-mixed fractional brownian motions and ornstein–uhlenbeck processes, Modern Stochastics: Theory and Applications, 10 (2023), pp. 343–366.
  * [3] L. A. Barboza and F. G. Viens, Parameter estimation of Gaussian stationary processes using the generalized method of moments, Electron. J. Stat., 11 (2017), pp. 401–439.
  * [4] P. Cheridito, H. Kawaguchi, and M. Maejima, Fractional Ornstein-Uhlenbeck processes, Electron. J. Probab., 8 (2003), pp. no. 3, 14 pp. (electronic).
  * [5] J. L. Doob, The Brownian movement and stochastic equations, Ann. of Math. (2), 43 (1942), pp. 351–369.
  * [6] A. Einstein, Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen, Annalen der Physik, 322 (1905), pp. 549–560.
  * [7] D. Gale and H. Nikaido, The jacobian matrix and global univalence of mappings, Mathematische Annalen, 159 (1965), pp. 81–93.
  * [8] U. Grenander, Stochastic processes and statistical inference, Arkiv för matematik, 1 (1950), pp. 195–277.
  * [9] A. R. Hall, Generalized method of moments, Oxford university press, 2005\.
  * [10] L. P. Hansen, Large sample properties of generalized method of moments estimators, Econometrica: Journal of the Econometric Society, (1982), pp. 1029–1054.
  * [11] Y. Hu and D. Nualart, Parameter estimation for fractional ornstein–uhlenbeck processes, Statistics & probability letters, 80 (2010), pp. 1030–1038.
  * [12] J. Lamperti, Semi-stable stochastic processes, Trans. Amer. Math. Soc., 104 (1962), pp. 62–78.
  * [13] P. Langevin, Sur la théorie du mouvement brownien, C. R. Acad. Sci. Paris, 146 (1908), pp. 530–533.
  * [14] B. B. Mandelbrot and J. W. Van Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev., 10 (1968), pp. 422–437.
  * [15] K. Newey and D. McFadden, Large sample estimation and hypothesis, Handbook of Econometrics, IV, Edited by RF Engle and DL McFadden, (1994), pp. 2112–2245.
  * [16] I. Nourdin and G. Peccati, Normal approximations with Malliavin calculus: from Stein’s method to universality, vol. 192, Cambridge University Press, 2012.
  * [17] P. Protter, Stochastic Integration and Differential Equation, Springer-Verlag, Berlin, Heidelberg, second ed., 1992.
  * [18] G. E. Uhlenbeck and L. S. Ornstein, On the theory of the Brownian motion, Phys. Rev., 36 (1930), pp. 823–841.