comment

# Second Moment Closure Modeling and DNS of Stratified Shear Layers

Naman Jain  Department of Mechanical Engineering, Pennsylvania State
University, Pennsylvania, USA, 16802, USA  Hieu T. Pham  Department of
Mechanical and Aerospace Engineering, University of California, San Diego,
California 92093, USA  Xinyi Huang  Department of Mechanical Engineering,
Pennsylvania State University, Pennsylvania, USA, 16802, USA  Sutanu Sarkar
Department of Mechanical and Aerospace Engineering, University of California,
San Diego, California 92093, USA  Xiang Yang  Department of Mechanical
Engineering, Pennsylvania State University, Pennsylvania, USA, 16802, USA
Robert Kunz  Department of Mechanical Engineering, Pennsylvania State
University, Pennsylvania, USA, 16802, USA

###### Abstract

Buoyant shear layers are encountered in many engineering and environmental
applications, and have been studied by researchers in the context of
experiments and modeling for decades. Often, these flows have high Reynolds
and Richardson numbers, and this leads to significant/intractable space-time
resolution requirements for DNS or LES modeling. On the other hand, many of
the important physical mechanisms in these systems, such as stress anisotropy,
wake stabilization, and regime transition, inherently render eddy viscosity-
based RANS modeling inappropriate. Accordingly, we pursue second-moment
closure (SMC), i.e., full Reynolds stress/flux/variance modeling, for moderate
Reynolds number non-stratified, and stratified shear layers for which DNS is
possible. A range of sub-model complexity is pursued for the diffusion of
stresses, density fluxes and variance, pressure strain and scrambling, and
dissipation. These sub-models are evaluated in terms of how well they are
represented by DNS in comparison to the exact Reynolds averaged terms, and how
well they impact the accuracy of the full RANS closure.

For the non-stratified case, the SMC model predicts the shear layer growth
rate and Reynolds shear stress profiles accurately. Stress anisotropy and
budgets are captured only qualitatively. Comparing DNS of exact and modeled
terms, inconsistencies in model performance and assumptions are observed,
including inaccurate prediction of individual statistics, non-negligible
pressure diffusion, and dissipation anisotropy. For the stratified case, shear
layer and gradient Richardson number growth rates, and stress, flux and
variance decay rates, are captured with less accuracy than corresponding flow
parameters in the non-stratified case. These studies lead to several
recommendations for model improvement.

## Nomenclature

  * $u_{i}$ instantaneous velocity vector

  * $\overline{u_{i}}$ mean velocity vector

  * $u^{\prime}_{i}$ fluctuating velocity vector

  * $x_{i}$ Cartesian coordinate

  * $t$ time

  * $\rho$ instantaneous fluid density

  * $\rho_{0}$ reference fluid density, $1\frac{kg}{m^{3}}$

  * $\overline{\rho}$ mean fluid density

  * $\rho^{\prime}$ fluctuating fluid density

  * $T$ instantaneous fluid temperature

  * $T_{0}$ reference fluid temperature, $1K$

  * $\overline{T}$ mean fluid temperature

  * $T^{\prime}$ fluctuating fluid temperature

  * $P$ instantaneous pressure

  * $\overline{P}$ mean pressure

  * $p^{\prime}$ fluctuating pressure

  * $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ Reynolds stresses

  * $R_{ij}$ Reynolds stress tensor

  * $a_{ij}$ anisotropy tensor, $a_{ij}=\frac{R_{ij}}{k}-\frac{2}{3}\delta_{ij}$

  * $\overline{u^{\prime}_{i}T^{\prime}}$ temperature fluxes

  * $\overline{T^{\prime}T^{\prime}}$ temperature variance

  * $k$ turbulent kinetic energy, $k=\frac{1}{2}\left(\overline{u^{\prime}_{i}u^{\prime}_{i}}\right)$

  * $q^{2}$ turbulence intensity, $q^{2}=\overline{u^{\prime}_{i}u^{\prime}_{i}}$

  * $\varepsilon_{ij}$ turbulence dissipation rate tensor, $\varepsilon_{ij}=2\nu\overline{\frac{\partial u^{\prime}_{i}}{\partial x_{k}}\frac{\partial u^{\prime}_{j}}{\partial x_{k}}}$

  * $\varepsilon$ isotropic turbulence dissipation rate, $\varepsilon=\nu\overline{\frac{\partial u^{\prime}_{i}}{\partial x_{j}}\frac{\partial u^{\prime}_{i}}{\partial x_{j}}}$

  * $\omega$ specific turbulence dissipation rate, $\omega=\frac{\varepsilon}{C_{\mu}k}$

  * $C_{\mu}$ Prandtl-Kolmogorov constant, $C_{\mu}=0.09$

  * $\varepsilon_{iT}$ temperature flux dissipation rate, $\varepsilon_{iT}=\left(\nu+\alpha\right)\overline{\frac{\partial u^{\prime}_{i}}{\partial x_{k}}\frac{\partial T^{\prime}}{\partial x_{k}}}$

  * $\varepsilon_{TT}$ temperature variance dissipation rate, $\varepsilon_{TT}=2\alpha\overline{\frac{\partial T^{\prime}}{\partial x_{k}}\frac{\partial T^{\prime}}{\partial x_{k}}}$

  * $\nu$ fluid kinematic viscosity

  * $\alpha$ fluid thermal diffusivity

  * $g_{i}$ gravitational vector

  * $\beta$ expansion coefficient at constant pressure

  * $U_{1}$ higher streamwise velocity magnitude, $U_{1}=0.5\frac{m}{s}$

  * $U_{2}$ lower streamwise velocity magnitude, $U_{2}=-0.5\frac{m}{s}$

  * $\Delta U$ velocity difference between bottom and top layer, $\Delta U=U_{1}-U_{2}=1\frac{m}{s}$

  * $\Delta T$ temperature difference across $4\delta_{\theta,0}$ at $t=0$

  * $\delta_{\theta}$ momentum thickness at time $t$, $\delta_{\theta}=\int_{-\infty}^{\infty}\left(\frac{1}{4}-\frac{\overline{u_{1}}^{2}}{\Delta U^{2}}\right)dx_{2}$

  * $\delta_{\theta,0}$ momentum thickness at time $t=0$, $0.25m$

  * $L_{i}$ domain length along $x_{i}$ direction

  * $x_{2,0}$ centerline vertical coordinate

  * $N$ Brunt-Vaisala Frequency, $N=\sqrt{-\frac{g_{2}}{\rho_{0}}\frac{\partial\rho}{\partial x_{2}}}$

  * $N_{i}$ grid points in the $i$th Cartesian direction

  * $m$ wavenumber

  * $\delta_{ij}$ Kronecker delta

  * $S$ shear, $S=\frac{\partial\overline{u_{1}}}{\partial x_{2}}$

  * $Re_{0}$ initial Reynolds number, $Re_{0}=\frac{\Delta U\times 4\delta_{\theta,0}}{\nu}=640$

  * $Pr$ Prandtl number, $Pr=\frac{\nu}{\alpha}=1$

  * $Ri_{g,0}$ initial gradient Richardson number, $Ri_{g}=\frac{N^{2}}{S^{2}}$

  * $\mathcal{P}_{ij}$ Reynolds stress production tensor

  * $\mathcal{P}_{iT}$ temperature flux production vector

  * $\mathcal{P}_{TT}$ temperature variance production scalar

  * $\mathcal{D}_{ij}$ Reynolds stress diffusion tensor

  * $\mathcal{D}_{iT}$ temperature flux diffusion vector

  * $\mathcal{D}_{TT}$ temperature variance diffusion scalar

  * $\Pi_{ij}$ pressure-redistribution tensor

  * $\Pi_{iT}$ pressure-scrambling vector

  * $\mathcal{B}_{ij}$ Reynolds stress buoyancy production tensor

  * $\mathcal{B}_{iT}$ temperature flux buoyancy production vector

##  1 Introduction

Wake flows differ markedly in stratified compared to non-stratified
environments. Specifically, numerous complex physical phenomena arise,
including strong Reynolds stress anisotropy, wake flattening/collapse,
counter-gradient fluxes, the coupling between kinetic and potential energies,
internal gravity waves, and regime transition.

Many researchers have studied stratification in homogeneous flows, and these
studies have led to a community focus on salient physics and flow parameters
that characterize these complex turbulent systems.Numerical analysis focusing
on 2D stratified mixing layers have been performed in [1, 2, 3], including
flows with initially turbulent perturbations [4]. However, these turbulence
resolving Direct Numerical Simulations (DNSs) and Large Eddy Simulations
(LESs) require significant computational resources at the higher Reynolds and
Richardson numbers observed in many practical applications. Accordingly,
Reynolds-Averaged Navier-Stokes (RANS) methods are attractive in these flows,
provided they can return sufficient accuracy.

RANS based Eddy Viscosity Models (EVMs) are inherently incapable of
accommodating the physics of stratified flow where the stratification is of
importance due to underlying stress anisotropy and attendant wake
stabilization and regime transition, particularly at high Richardson numbers
[5]. An alternate approach is to use Second Moment Closure (SMC) methods. SMCs
abandon the Boussinesq approximation that underpins eddy viscosity models, and
solve individual transport equations for the Reynolds stresses, fluxes, and
the density variance. SMC based methods have demonstrated superior predictive
capability to EVMs for many flows, see, e.g., Ref. [6]. Historically, these
methods have been less widely applied than EVMs due to increased stiffness,
complexity, and computational cost associated with solving additional
transport equations, and for being numerically non-robust. However,
significant increases in available computational resources have allowed for
adequate mesh density for 3D applications, leading to stable, practical and
robust applications of SMC methods [7, 8].

To date, experimental measurements of important modeling terms (pressure-
redistribution, for example) have been very challenging and unreliable [9],
and detailed DNS/LES studies that yield such terms ”exactly” have been
restricted to simple flow configurations. Consequently, calibration of the
numerous SMC model constants has been limited to flows using reduced forms of
the RANS equations [10], rather than on the accuracy with which the individual
sub-models approximate the exact terms. Although these key simplifications
enable calibration, they do not generalize well to more complex flows, e.g.,
separated flows and simple shear flows [11].

Research on Reynolds stress RANS modeling has focused primarily on improving
model performance in reproducing important first and second moment statistics,
turbulence energetics [12, 13]. Because of the recent availability of DNS and
LES predictions, it is now feasible to conduct detailed comparisons of the
individual term-wise performance of SMC sub-models with the exact terms in the
governing equations, even in complex dynamical systems like stratified flows.
The aim of the current work is to perform such detailed analysis.

Here, the flow dynamics of initially turbulent, temporally evolving stratified
and non-stratified shear layers are studied. Full Reynolds stress RANS and DNS
are applied, with the goal of quantifying model shortcomings and
inconsistencies. The paper is organized as follows. A summary of the governing
equations, numerical schemes, simulation parameters, and the initial and
boundary conditions is provided in section 2. The RANS SMC sub-models
investigated are delineated in section 3. The results from the two different
flow configurations are presented in section 4. Finally, conclusions of the
current work are provided in section 5.

##  2 THEORETICAL FORMULATION

The flow configuration considered is a temporally evolving shear layer that
develops when two miscible fluids with velocities equal in magnitude but
opposite signs are brought together. The mixing layer generated is developed
in a stably stratified fluid with a linear density gradient as shown in Fig 1.

Figure 1: Schematic of the flow configuration.

###  2.1 Governing Equations

For the temporal mixing layer, results from two different DNS solvers (labeled
UCSD-DNS [4] and AFiD-DNS [14]) and a RANS SMC solver (labeled NPhase [15])
are presented. The flow is assumed to be incompressible, unsteady and three-
dimensional. UCSD-DNS solves the conservation of mass, momentum and scalar
(density) equations with the Boussinesq approximation invoked:

$\displaystyle\frac{\partial u_{i}}{\partial x_{i}}=$ $\displaystyle 0$ (1)
$\displaystyle\frac{\partial u_{i}}{\partial
t}+\frac{\partial(u_{i}u_{j})}{\partial x_{j}}=$
$\displaystyle-\frac{1}{\rho_{0}}\frac{\partial P}{\partial
x_{i}}+\nu\frac{\partial^{2}u_{i}}{\partial x_{j}\partial
x_{j}}+\frac{\rho}{\rho_{0}}g_{i}$ (2)
$\displaystyle\centering\frac{\partial\rho}{\partial t}+\frac{\partial(\rho
u_{j})}{\partial x_{j}}=\@add@centering$
$\displaystyle\alpha\frac{\partial^{2}\rho}{\partial x_{j}\partial x_{j}}$ (3)

AFiD-DNS along with the instantaneous Eqns. (1) and (2), solves the transport
equation for temperature:

$\centering\frac{\partial T}{\partial t}+\frac{\partial(Tu_{j})}{\partial
x_{j}}=\alpha\frac{\partial^{2}T}{\partial x_{j}\partial
x_{j}}\@add@centering$ (4)

and the density field is estimated from a linear equation of state:

$\centering\frac{\rho-\rho_{0}}{\rho_{0}}=-\beta(T-T_{0});\qquad\beta=-\frac{1}{\rho_{0}}\left(\frac{\partial\rho}{\partial
T}\right)_{P}\@add@centering$ (5)

NPhase solves the Unsteady RANS (URANS) equations under the Boussinesq
approximation and Eqn. (5), along with the exact transport equations for
important second-order statistics
$\left(\overline{u^{\prime}_{i}u^{\prime}_{j}},\overline{u^{\prime}_{i}T^{\prime}},\overline{T^{\prime}T^{\prime}}\right)$:

$\displaystyle\centering\frac{\partial\overline{u_{i}}}{\partial
x_{i}}=\@add@centering$ $\displaystyle 0$ (6)
$\displaystyle\centering\frac{\partial\overline{u_{i}}}{\partial
t}+\overline{u_{j}}\frac{\partial\overline{u_{i}}}{\partial
x_{j}}=\@add@centering$
$\displaystyle-\frac{1}{\rho_{0}}\frac{\partial\overline{P}}{\partial
x_{i}}+\frac{\partial}{\partial
x_{j}}\left(\nu\frac{\partial\overline{u_{i}}}{\partial
x_{j}}-\overline{u^{\prime}_{i}u^{\prime}_{j}}\right)+\frac{\overline{\rho}}{\rho_{0}}g_{i}$
(7) $\displaystyle\centering\frac{\partial\overline{T}}{\partial
t}+\overline{u_{j}}\frac{\partial\overline{T}}{\partial
x_{j}}=\@add@centering$ $\displaystyle\frac{\partial}{\partial
x_{j}}\left(\alpha\frac{\partial\overline{T}}{\partial
x_{j}}-\overline{u^{\prime}_{j}T^{\prime}}\right)$ (8)
$\frac{\partial(\overline{u^{\prime}_{i}u^{\prime}_{j}})}{\partial
t}+\overline{u_{k}}\left(\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial
x_{k}}\right)=\underbrace{\left(-\overline{u^{\prime}_{i}u^{\prime}_{k}}\frac{\partial\overline{u_{j}}}{\partial
x_{k}}-\overline{u^{\prime}_{j}u^{\prime}_{k}}\frac{\partial\overline{u_{i}}}{\partial
x_{k}}\right)}_{\mathcal{P}_{ij}}-\varepsilon_{ij}+\underbrace{\left(-\beta\left(g_{i}\overline{u^{\prime}_{j}T^{\prime}}+g_{j}\overline{u^{\prime}_{i}T^{\prime}}\right)\right)}_{\mathcal{B}_{ij}}\\\
+\underbrace{\left(\frac{\partial}{\partial
x_{k}}\left[-\overline{u^{\prime}_{i}u^{\prime}_{j}u^{\prime}_{k}}-\frac{1}{\rho}\overline{p^{\prime}\left(u^{\prime}_{i}\delta_{jk}+u^{\prime}_{j}\delta_{ik}\right)}+\nu\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial
x_{k}}\right]\right)}_{\mathcal{D}_{ij}}+\underbrace{\left(\overline{\frac{p^{\prime}}{\rho}\left[\frac{\partial
u^{\prime}_{i}}{\partial x_{j}}+\frac{\partial u^{\prime}_{j}}{\partial
x_{i}}\right]}\right)}_{\Pi_{ij}}$ (9)
$\frac{\partial(\overline{u^{\prime}_{i}T^{\prime})}}{\partial
t}+\overline{u_{k}}\left(\frac{\partial\overline{u^{\prime}_{i}T^{\prime}}}{\partial
x_{k}}\right)=\underbrace{\left(-\overline{u^{\prime}_{k}u^{\prime}_{i}}\frac{\partial\overline{T}}{\partial
x_{k}}-\overline{u^{\prime}_{k}T^{\prime}}\frac{\partial\overline{u_{i}}}{\partial
x_{k}}\right)}_{\mathcal{P}_{iT}}-\varepsilon_{iT}+\underbrace{\left(\overline{\frac{p^{\prime}}{\rho}\frac{\partial
T^{\prime}}{\partial
x_{i}}}\right)}_{\Pi_{iT}}+\underbrace{\left(-\beta\overline{T^{\prime}T^{\prime}}g_{i}\right)}_{\mathcal{B}_{iT}}\\\
+\underbrace{\left(\frac{\partial}{\partial
x_{k}}\left[-\overline{u^{\prime}_{i}u^{\prime}_{k}T^{\prime}}-\frac{1}{\rho}\overline{p^{\prime}T^{\prime}\delta_{ik}}+\nu\overline{T^{\prime}\frac{\partial
u^{\prime}_{i}}{\partial
x_{k}}}+\frac{\nu}{Pr}\overline{u^{\prime}_{i}\frac{\partial
T^{\prime}}{\partial x_{k}}}\right]\right)}_{\mathcal{D}_{iT}}$ (10)
$\frac{\partial(\overline{T^{\prime}T^{\prime}})}{\partial
t}+\overline{u_{k}}\left(\frac{\partial\overline{T^{\prime}T^{\prime}}}{\partial
x_{k}}\right)=\underbrace{\left(-2\overline{u^{\prime}_{k}T^{\prime}}\frac{\partial\overline{T}}{\partial
x_{k}}\right)}_{\mathcal{P}_{TT}}-\varepsilon_{TT}+\underbrace{\left(\frac{\partial}{\partial
x_{k}}\left[\frac{\nu}{Pr}\frac{\partial\overline{T^{\prime}T^{\prime}}}{\partial
x_{k}}-\overline{u^{\prime}_{k}T^{\prime}T^{\prime}}\right]\right)}_{\mathcal{D}_{TT}}$
(11)

where in Eqns. (9), (10) and (11), production terms,
$\mathcal{P}_{ij},\mathcal{P}_{iT},\mathcal{P}_{TT}$, and buoyancy terms,
$\mathcal{B}_{ij},\mathcal{B}_{iT}$, are retained in their exact form at the
level of SMC. The remaining terms require modeling as summarized in section 3.

###  2.2 Numerical Schemes

The NPhase solver employs a segregated pressure based methodology with a
collocated variable arrangement and lagged coefficient linearization method
(see [16], for example). A diagonal dominance preserving, finite volume
spatial discretization scheme is selected for the momentum and turbulence
transport equations. Continuity is introduced through a pressure correction
equation, based on the SIMPLE-C algorithm [17]. The cell face fluxes are
generated through a momentum interpolation scheme [18], which introduces
damping in continuity equation. At each iteration, the discrete momentum
equations are solved approximately, followed by a more exact solution of the
pressure correction equation. The turbulence scalar and enthalpy equations are
then solved in succession. For further details, refer to Ref. [15].

The AFiD-DNS solver is an open-source code that uses a second-order finite
difference scheme and a staggered grid to solve for fluid velocities, and a
second-order Adams-Bashforth method for time discretization. Further details
can be found in Ref. [14]. For a better convergence of higher-order
statistics, an ensemble average of 100 AFiD-DNS realizations was found
sufficient and used here for results.

The UCSD-DNS solver employs a staggered grid with normal velocities stored at
cell faces, and density and pressure at cell centers. A second-order central
difference scheme for spatial discretization and low storage third-order
Runge-Kutta method for temporal integration is employed. A sponge region is
employed at the top and bottom boundaries to control spurious reflections. A
detailed description can be found in Refs. [4, 3].

###  2.3 Simulation Parameters

Details of the different simulation parameters including the computational
domain size $(L_{1},L_{2},L_{3})$, node count $(N_{1},N_{2},N_{3})$, Prandtl
number $(Pr)$, initial gradient Richardson number $(Ri_{g,0})$ and Reynolds
number $(Re_{0})$ are listed in the Table 1.

Table 1: Flow parameters for the UCSD-DNS, AFiD-DNS and RANS-SMC (NPhase) for the two test cases, S0 and S1 Case | S0 | S1
---|---|---
UCSD-DNS | AFiD-DNS | NPhase | UCSD-DNS | AFiD-DNS | NPhase
$N_{1}$ | 1024 | 384 | 4 | 1024 | 384 | 4
$N_{2}$ | 512 | 192 | 401 | 512 | 450 | 401
$N_{3}$ | 256 | 128 | 2 | 256 | 128 | 2
$Ri_{g,0}$ | 0 | 0 | 0 | 0.04 | 0.04 | 0.04
$Re_{0}$ | 640 | 640 | 640 | 640 | 640 | 640
$Pr$ | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0
$L_{1}$ | 87.04 | 64.5 | 60 | 87.04 | 64.5 | 60
$L_{2}$ | 48.53 | 32.25 | 65 | 48.53 | 32.25 | 65
$L_{3}$ | 21.76 | 21.5 | 1 | 21.76 | 21.5 | 1

###  2.4 Initial Conditions

The mean velocity and density profiles for all three simulations (UCSD-DNS,
AFiD-DNS, NPhase) are initialized as:

$\displaystyle\overline{u_{1}}=$ $\displaystyle-\frac{\Delta
U}{2}\tanh\left(\frac{2\left(x_{2}-x_{2,0}\right)}{4\delta_{\theta,0}}\right)$
(12) $\displaystyle\overline{u_{2}}=$
$\displaystyle\overline{u_{3}}=\overline{P}=0$ (13)
$\displaystyle\overline{\rho}=$
$\displaystyle\overline{\rho_{0}}-\left(\frac{\overline{\rho_{0}}N^{2}_{0}}{\left|g_{2}\right|}\right)\left(x_{2}-x_{2,0}\right)$
(14)

where $N_{0}=N(t=0)$, is set from the value of $Ri_{g,0}$ (Table 1). The
values of $\left(\delta_{\theta,0}=0.25\right)$ and $\left(\Delta
U=1.0\right)$ are set for both the test cases. For the UCSD-DNS and AFiD-DNS
simulations, the flow is initialized with low amplitude velocity perturbations
using a broadband spectrum:

$E(m)\propto m^{4}\exp\left[-2\left(m/m_{0}\right)^{2}\right]$ (15)

where $m_{0}=3.7$, and the peak value of $q^{2}$ is set as $0.03\Delta U^{2}$.
For UCSD-DNS, the density field was initialized with zero density
fluctuations, $\rho^{\prime}=0$. Corresponding to this, AFiD-DNS and NPhase
are initialized to match the UCSD-DNS $\rho^{\prime}=0$ specification exactly
using zero initial temperature fluctuations, $T^{\prime}=0$. Additionally, for
NPhase, the specific dissipation rate, $\omega$, and shear stress,
$(\overline{u^{\prime}_{1}u^{\prime}_{2}})$, are initialized as:

$\displaystyle\omega(t=0)=$ $\displaystyle
2.7125\left(\frac{k^{\frac{1}{2}}}{C_{\mu}^{\frac{1}{4}}\left(4\delta_{\theta,0}\right)}\right)$
(16) $\displaystyle\overline{u^{\prime}_{1}u^{\prime}_{2}}(t=0)=$
$\displaystyle-0.0185\left(C_{\mu}^{\frac{1}{4}}k^{\frac{1}{2}}\left(4\delta_{\theta,0}\right)\frac{\partial\overline{u_{1}}}{\partial
x_{2}}\right)$ (17)

###  2.5 Boundary Conditions

The UCSD-DNS and AFiD-DNS models use periodic boundary conditions in the
streamwise ($x_{1}$) and spanwise ($x_{3}$) directions for all the variables.
In the vertical ($x_{2}$) direction, the following conditions are applied:

$\overline{u_{1}}(0)=\frac{1}{2},\quad\overline{u_{1}}(L_{2})=-\frac{1}{2},\quad\overline{u_{2}}(0)=\overline{u_{2}}(L_{2})=0,\\\
\overline{u_{3}}(0)=\overline{u_{3}}(L_{2})=0,\quad\overline{\frac{\partial\rho}{\partial
x_{2}}}(0)=\overline{\frac{\partial\rho}{\partial
x_{2}}}(L_{2})=-\frac{\overline{\rho_{0}}N^{2}_{0}}{\left|g_{2}\right|}$ (18)

For NPhase, the periodic cyclic and symmetric boundary conditions are
specified along the streamwise ($x_{1}$) and spanwise ($x_{3}$) directions
respectively. In the vertical ($x_{2}$) direction, in addition to Eqn. (18), a
symmetry boundary is specified for Reynolds stresses, temperature fluxes,
variance, and specific dissipation rate.

##  3 RANS SMC MODELS

A short description of the models considered to close the system of Eqns. (9),
(10) and (11) is presented below, and all model coefficients are listed in
Tables 2 and 3.

###  3.1 Pressure-Redistribution Term $(\Pi_{ij})$ Model

The standard approach of modeling the pressure-redistribution term for buoyant
flows comprises decomposition into; slow term $\Pi^{S}_{ij}$, rapid term
$\Pi^{R}_{ij}$, and buoyant-contribution $\Pi^{B}_{ij}$, such that
$\Pi_{ij}=\Pi^{S}_{ij}+\Pi^{R}_{ij}+\Pi^{B}_{ij}$. The non-linear SSG model
[19] for $\Pi^{S}_{ij}+\Pi^{R}_{ij}$, (first five terms below), and the
isotropization of production (IP) strategy [20] for $\Pi^{B}_{ij}$ (last
term), are investigated:

$\Pi_{ij}=-\left(C_{1}\varepsilon+\frac{1}{2}C^{*}_{1}\mathcal{P}_{kk}\right)a_{ij}+C_{2}\varepsilon\left(a_{ik}a_{kj}-\frac{1}{3}a_{kl}a_{kl}\delta_{ij}\right)+C_{5}k\left(a_{ik}W_{jk}+a_{jk}W_{ik}\right)\\\
+\left(C_{3}-C^{*}_{3}\sqrt{a_{kl}a_{kl}}\right)k\left(S_{ij}-\frac{1}{3}S_{kk}\delta_{ij}\right)+C_{4}k\left(a_{ik}S_{jk}+a_{jk}S_{ik}-\frac{2}{3}a_{kl}S_{kl}\delta_{ij}\right)-C_{bf}\left(\mathcal{B}_{ij}-\frac{1}{3}\mathcal{B}_{kk}\delta_{ij}\right)$
(19)

where,

$a_{ij}=\frac{R_{ij}}{k}-\frac{2}{3}\delta_{ij};\qquad
S_{ij}=\frac{1}{2}\left(\frac{\partial\overline{u_{i}}}{\partial
x_{j}}+\frac{\partial\overline{u_{j}}}{\partial x_{i}}\right);\qquad
W_{ij}=\frac{1}{2}\left(\frac{\partial\overline{u_{i}}}{\partial
x_{j}}-\frac{\partial\overline{u_{j}}}{\partial x_{i}}\right)$ (20)

###  3.2 Pressure-Scrambling Term $(\Pi_{iT})$ Model

The pressure-scrambling term $\Pi_{iT}$, is consistently comprised of three
components, $\Pi_{iT,1}$, $\Pi_{iT,2}$, and $\Pi_{iT,3}$. The simple return-
to-isotropy model for $\Pi_{iT,1}$, basic IP model for $\Pi_{iT,2}$ and quasi-
isotropic model for $\Pi_{iT,3}$, detailed in [21], are implemented:

$\Pi_{iT}=\underbrace{\left(-C_{T1}\frac{\varepsilon}{k}\overline{u^{\prime}_{i}T^{\prime}}\right)}_{\Pi_{iT,1}}+\underbrace{\left(C_{T2}\overline{u^{\prime}_{k}T^{\prime}}\frac{\partial\overline{u_{i}}}{\partial
x_{k}}\right)}_{\Pi_{iT,2}}+\underbrace{\left(C_{T3}\beta
g_{i}\overline{T^{\prime}T^{\prime}}\right)}_{\Pi_{iT,3}}$ (21)

###  3.3 Diffusion Term
$\left(\mathcal{D}_{ij},\mathcal{D}_{iT},\mathcal{D}_{TT}\right)$ Models

The standard approach of modeling the diffusion term involves neglecting the
pressure-diffusion contribution and modeling the triple-velocity correlation
terms. For the current case, the simple gradient diffusion (SD) model [22] and
the Daly-Harlow (DH) generalized gradient diffusion model [23] are considered:

$\displaystyle\mathcal{D}^{SD}_{ij}=$ $\displaystyle\frac{\partial}{\partial
x_{k}}\left[\mathcal{A}_{SD}\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial
x_{k}}\right];\mathcal{D}^{DH}_{ij}=\frac{\partial}{\partial
x_{k}}\left[\mathcal{A}_{DH}\frac{\partial\overline{u^{\prime}_{i}u_{j}^{\prime}}}{\partial
x_{l}}\right]$ (22) $\displaystyle\mathcal{D}^{SD}_{iT}=$
$\displaystyle\frac{\partial}{\partial
x_{k}}\left[\mathcal{A}_{SD}\frac{\partial\overline{u^{\prime}_{i}T^{\prime}}}{\partial
x_{k}}\right];\mathcal{D}^{DH}_{iT}=\frac{\partial}{\partial
x_{k}}\left[\mathcal{A}_{DH}\frac{\partial\overline{u^{\prime}_{i}T^{\prime}}}{\partial
x_{l}}\right]$ (23) $\displaystyle\mathcal{D}^{SD}_{TT}=$
$\displaystyle\frac{\partial}{\partial
x_{k}}\left[\mathcal{A}_{SD}\frac{\partial\overline{T^{\prime}T^{\prime}}}{\partial
x_{k}}\right];\mathcal{D}^{DH}_{TT}=\frac{\partial}{\partial
x_{k}}\left[\mathcal{A}_{DH}\frac{\partial\overline{T^{\prime}T^{\prime}}}{\partial
x_{l}}\right]$ (24)

where,

$\mathcal{A}_{SD}=\left(\nu+C_{SD}\frac{k^{2}}{\varepsilon}\right);\hfill\mathcal{A}_{DH}=\left(\nu\delta_{kl}+C_{DH}\frac{k\left(\overline{u^{\prime}_{k}u^{\prime}_{l}}\right)}{\varepsilon}\right)$
(25)

###  3.4 Dissipation Rate Term
$\left(\varepsilon_{ij},\varepsilon_{iT},\varepsilon_{TT}\right)$ Models

The dissipation rate terms, $\varepsilon_{ij}$ and $\varepsilon_{iT}$, are
modeled based on the commonly employed local isotropy assumption [24], and an
algebraic relation is employed for $\varepsilon_{TT}$:

$\varepsilon_{ij}=\frac{2}{3}\varepsilon\delta_{ij};\qquad\varepsilon_{iT}=0;\qquad\varepsilon_{TT}=R\frac{\varepsilon}{k}\overline{T^{\prime}T^{\prime}}$
(26)

where, the isotropic dissipation rate, $\varepsilon$, is obtained by solving a
transport equation for the specific dissipation rate, $\omega$:

$\frac{\partial\omega}{\partial
t}+\overline{u_{k}}\left(\frac{\partial\omega}{\partial
x_{k}}\right)=\frac{\alpha_{\omega}\omega}{2k}\left(\mathcal{P}_{kk}+\mathcal{B}_{kk}\right)-\beta_{\omega}\omega^{2}+\frac{\partial}{\partial
x_{k}}\left(\left[\nu+\sigma_{\omega}\frac{k}{\omega}\right]\frac{\partial\omega}{\partial
x_{k}}\right)+\frac{\sigma_{d}}{\omega}max\left(\frac{\partial k}{\partial
x_{k}}\frac{\partial\omega}{\partial x_{k}},0\right)$ (27)

where for $R$, a constant value of 1.5 recommended by [25] was used here.

Table 2: RANS-SMC model coefficients SSG | $\Pi_{iT}$ | SD | DH
---|---|---|---
$C_{1}$ | $C^{*}_{1}$ | $C_{2}$ | $C_{3}$ | $C^{*}_{3}$ | $C_{4}$ | $C_{5}$ | $C_{T1}$ | $C_{T2}$ | $C_{T3}$ | $C_{SD}$ | $C_{bf}$ | $C_{DH}$ | $C_{bf}$
1.7 | 0.9 | 1.05 | 0.8 | 0.65 | 0.625 | 0.2 | 3.5 | 0.5 | 0.4 | $\frac{0.44}{3}$ | 0.3 | 0.22 | 0.15
Table 3: RANS-SMC model coefficients for $\omega$ transport equation $\omega$
---
$\alpha_{\omega}$ | $\beta_{\omega}$ | $\sigma_{\omega}$ | $\sigma_{d}$
0.44 | 0.0828 | 0.856 | 1.712

##  4 RESULTS

The RANS and DNS predictions for two different test cases are considered here:
a non-stratified shear layer (S0) and a stratified shear layer (S1). In the S0
case, the mixing layer is evolved with a linear density gradient (with
Boussinesq approximation). For this case, the buoyancy terms are removed from
the equations to make the flow dynamics identical to a non-stratified system.
The resulting system develops with the density/temperature fields behaving as
a passive scalar. In the stratified shear layer (S1) case, the buoyant effects
are fully included resulting in a complex dynamic system. The simulation
details are listed in Table 1.

###  4.1 Non-Stratified Shear Layer (S0)

For the S0 case, the flow eventually approaches a self-similar state as found
in experimental [26, 27] and DNS [28, 29, 4] studies. The mean velocity and
RMS profiles evolve self-similarly, and the Reynolds stresses attain a maximum
at the centerline and remain constant thereafter. Of the many methods of
estimating the shear layer thickness, the momentum thickness,
$\delta_{\theta}$, is used in the current case. A well established empirical
relation for the growth rate is:

$\centering\frac{d\delta_{\theta}}{dt}=2\frac{C_{\delta}}{D_{\omega}}(U_{1}-U_{2})\@add@centering$
(28)

where $C_{\delta}=0.16$ and $D_{\omega}=5$. In Fig. 2, the non-dimensional
momentum thickness is plotted versus non-dimensional time for the RANS model
(only SSG with DH model results are presented for the sake of brevity), Eqn.
(28), and the two DNS simulations.

Figure 2: Momentum thickness growth rate comparison between RANS (NPhase),
UCSD-DNS, AFiD-DNS and the empirical relation.

The SMC model successfully captures the linear momentum thickness growth rate
and thus the self-similar state (for SD diffusion case as well). In Fig. 3
profiles of the three normalized RMS velocity fluctuations and square root of
shear stress are compared between the RANS (SSG with DH model) and DNS
calculations, and two different experimental results, at a non-dimensional
timestep, $t\Delta U/4\delta_{\theta,0}=50$, which is in the self-similar
region. In these plots the velocity difference, $\Delta U$, and length-scale,
$\delta_{\omega}=D_{\omega}\delta_{\theta}$ are used for non-
dimensionalization.

(a) Streamwise RMS velocity

(b) Vertical RMS velocity

(c) Spanwise RMS velocity

(d) Square root of shear stress

Figure 3: RMS velocities and square root of shear stress comparison between
RANS (NPhase), UCSD-DNS, AFiD-DNS, Bell and Mehta [27], Pantano and Sarkar
[29], Rogers and Moser [28], Spencer and Jones [26].

The shear stress profile, plotted in Fig. 3(d), matches the DNS and
experimental predictions very well, attendant to the linear growth rate
prediction. The $\sqrt{R_{ij}}$ profiles are well-predicted by the SMC model
with a slight over-prediction of $\sqrt{R_{11}}$ and under-prediction of
$\sqrt{R_{22}}$. These results also serve to validate the model implementation
and the ability of SMC model to capture anisotropy and predict the self-
similar state. The centerline evolution of the non-dimensional temperature
fluxes and RMS of temperature fluctuations for the three simulations (SSG with
DH model) are shown in Figs. 4(a) and 4(b). Good agreement between the RANS
and DNS results can be seen, and the flux anisotropy is well captured by the
SMC model.

(a) Temperature fluxes

(b) RMS temperature

Figure 4: Comparison of Temperature fluxes and RMS of temperature evolution at
centerline between RANS (NPhase), UCSD-DNS and AFiD-DNS.

Next, the performance of the SMC sub-models is studied by comparing, in Fig.
5, the non-dimensional budget terms for the
$\overline{u^{\prime}_{1}u^{\prime}_{1}}$,
$\overline{u^{\prime}_{2}u^{\prime}_{2}}$ and
$\overline{u^{\prime}_{1}u^{\prime}_{2}}$ Reynolds stresses with the exact
terms obtained by AFiD-DNS, at $t\Delta U/4\delta_{\theta,0}=30$ (again, SSG
with DH for SMC).

Figure 5: Reynolds stress budget terms comparison with exact DNS (AFiD-DNS)
predicted terms. Pressure diffusion is neglected in the RANS diffusion models.

The model performs well except, importantly, for the terms in the vertical and
streamwise budgets, where the peak production and pressure-redistribution
terms are seen to be over-predicted by RANS. To further understand SSG sub-
model performance, the exact pressure-redistribution terms from DNS are
compared in Fig. 6 (labeled DNS), to what DNS returns for the pressure-
redistribution model, Eqn. (19) (labeled DNS,SSG). Qualitative agreement is
observed with an under-prediction of peak values for the non-zero stress
budget components. This gives an assessment of the accuracy of the SSG
pressure-redistribution model independent of the full RANS implementation.

Figure 6: Comparison of AFiD-DNS exact (DNS) $\Pi_{ij}$ term and AFiD-DNS
predicted model (DNS,SSG) $\Pi_{ij}$ term.

The diffusion term in the current flow configuration does not play an
important role close to the centerline but is dominant at the shear layer
edges. The widely used diffusion models applied in this work (SD and DH) have
historically been formulated to model the exact triple-velocity correlation
(TVC) terms (turbulent-diffusion) in the context of assuming that the
pressure-diffusion term is negligible [24, 22]. In order to assess this
approximation for the present flow, in Fig. 7, the exact turbulent-diffusion
and pressure-diffusion terms returned by AFiD-DNS are presented, with time
extruded onto the horizontal axis.

(a) Turbulent diffusion term

(b) Pressure-diffusion term

Figure 7: AFiD-DNS exact dimensional diffusion terms in
$\left(\overline{u^{\prime}_{1}u^{\prime}_{2}}\right)$ transport equation.

It is observed that the pressure-diffusion term is of the same order of
magnitude as the turbulent-diffusion and hence non-negligible. Despite this,
the SMC models (SSG with DH models) predicted the shear layer growth rate and
Reynolds stresses accurately. Thus, it appears that neglecting pressure-
diffusion and the over-prediction of pressure-redistribution and diffusion
terms offset one another in the calibrated SMC model leading to an overall
good prediction of bulk properties.

###  4.2 Stratified Shear Layer (S1)

The introduction of buoyancy in the transport equations has a significant
impact on the flow dynamics. Compared to the S0 case, the shear layer growth
rate is reduced significantly in the stably stratified S1 case due to the
damping of vertical velocity fluctuations. A comparison of the SMC non-
dimensional momentum thickness growth rate and centerline evolution of
gradient Richardson number are compared with the DNS predictions in Fig. 8.
The non-dimensional momentum thickness from both DNS models, exhibits a
decreasing growth rate over time and eventually reaches a near-asymptotic
state. The SMC models (SD and DH results included here) also return a
decreasing growth rate compared to the S0 case, indicating that the effect of
stratification is being accounted for. However, the growth rate returned by
the SMC is much higher than DNS. Also, RANS under-predicts $Ri_{g}$ and its
growth rate, and fails to capture the plateauing behavior. These S1 shear
layer thickness and gradient Richardson number results are quantitatively far
less accurate than the shear layer evolution results presented for the S0
case. These incorrect predictions by RANS SMC models imply the need to examine
the accuracy of second-order statistics.

Figure 8: RANS (NPhase, DH and SD models), UCSD-DNS, AFiD-DNS momentum
thickness and gradient Richardson number comparison.

The centerline evolution of the non-dimensional temperature fluxes and RMS of
velocity and temperature fluctuations for RANS (SD and DH models) and AFiD-DNS
are compared in Fig. 9. The growth rates predicted by SMCs (both models) match
well with the DNS results. For the peak values of RMS velocity fluctuations
and streamwise temperature flux, the error in SD model predictions is over
$10\%$, while the DH model is under $5\%$. The error in vertical temperature
flux and temperature variance predictions by both models is over $20\%$. The
decay rate following this peak is significantly under-predicted by SMCs. The
DNS results exhibit two different decay rates due to the laminarization of the
flow around, $t\Delta U/4\delta_{\theta,0}\approx 80$. In the SMC predictions,
an initial rapid decay transitions to a slower decay rate around, $t\Delta
U/4\delta_{\theta,0}\approx 30$. Additionally, the decay rates predicted by
the two diffusion models are almost identical, indicating RANS predictions
during the decay period is independent of diffusion model choice. The
incorrect decay rate predictions signify the need to improve the SMC models to
capture the flow dynamics accurately.

(a) Centerline RMS velocities

(b) Centerline temperature flux

(c) Centerline RMS temperature

Figure 9: RMS velocities, temperature fluxes and rms of temperature variations
comparison between RANS (NPhase, DH and SD models) and AFiD-DNS.

The SMC models considered here employ the isotropic dissipation rate,
$\varepsilon$, as a key parameter in formulation of all sub-models considered.
The validity of this assumption is checked by comparing the exact components
of the non-dimensional dissipation rate tensor, $\varepsilon_{ij}$, with the
corresponding isotropic value, $\left(\frac{2}{3}\varepsilon\right)$,
estimated by AFiD-DNS, at $t\Delta U/4\delta_{\theta,0}=60$, in Fig. 10.

Figure 10: AFiD-DNS estimated exact dissipation rate tensor components (S1).

At early time instants, the flow evolves like a non-stratified system, and the
departure from isotropy was found negligible (not presented for the sake of
brevity). The buoyancy term increases with time, and a time instant when its
impact on the flow is significant is chosen for comparison in Fig. 10. A
significant departure from isotropy can be observed, highlighting an important
inconsistency in all sub-model formulations where the isotropic value is used
as a key parameter. This conclusion can be extended to future times where the
buoyancy effects are significant. This invalid assumption is hypothesized to
be a reason behind inaccurate SMC predictions. Departure from isotropy in the
S0 case was found negligible, consistent with the computational and numerical
studies at high Reynolds number [30, 31]. The anisotropy and complex dynamics
introduced by stratification render the assumption invalid, resulting in
incorrect SMC decay rate predictions. Thus, this result implies the need to
account for anisotropy in the dissipation rate tensor to improve RANS
predictions.

##  5 Conclusions

The performance of RANS-SMC models were compared against DNS results of the
same case. Despite the reasonable prediction of first and second-order
statistics for the non-stratified case, individual sub-models were shown to
depart from DNS in the Reynolds stress budget term-wise comparison. The
omission of pressure-diffusion term offset by calibration of model
coefficients is hypothesized to be a reason for the inconsistencies. For the
stratified case, SMC models were shown to capture the complex flow dynamics
less accurately. The local isotropy assumption in dissipation modeling was
shown to be invalid for the stratified case and hypothesized to be a reason
for incorrect RANS predictions.

## Acknowledgement

This work was sponsored by the US Office of Naval Research under contract
N00014-19-1-2057, with Dr. Peter Chang as Technical Monitor.

## References

  * [1] C. Caulfield and W. Peltier, “The anatomy of the mixing transition in homogeneous and stratified free shear layers,” Journal of Fluid Mechanics, vol. 413, pp. 1–47, 2000.
  * [2] W. Smyth, J. Moum, and D. Caldwell, “The efficiency of mixing in turbulent patches: Inferences from direct simulations and microstructure observations,” Journal of Physical Oceanography, vol. 31, no. 8, pp. 1969–1992, 2001.
  * [3] H. T. Pham, S. Sarkar, and K. A. Brucker, “Dynamics of a stratified shear layer above a region of uniform stratification,” Journal of Fluid Mechanics, vol. 630, pp. 191–223, 2009.
  * [4] K. A. Brucker and S. Sarkar, “Evolution of an initially turbulent stratified shear layer,” Physics of Fluids, vol. 19, no. 10, p. 105105, 2007.
  * [5] G. L. Mellor and T. Yamada, “Development of a turbulence closure model for geophysical fluid problems,” Reviews of Geophysics, vol. 20, no. 4, pp. 851–875, 1982.
  * [6] D. Wilcox, “Reynolds averaging,” Turbulence Modeling for CFD, Third Edition, DCW Industries, Inc., La Cañada, California, pp. 34–38, 2006.
  * [7] R. F. Kunz and B. Lakshminarayana, “Three-dimensional navier-stokes computation of turbomachinery flows using an explicit numerical procedure and a coupled k-$\varepsilon$ turbulence model,” in ASME 1991 International Gas Turbine and Aeroengine Congress and Exposition, American Society of Mechanical Engineers Digital Collection, 1991.
  * [8] R. F. Kunz and B. Lakshminarayana, “Stability of explicit navier-stokes procedures using k-$\varepsilon$ and k-$\varepsilon$/algebraic reynolds stress turbulence models,” Journal of Computational Physics, vol. 103, no. 1, pp. 141–159, 1992.
  * [9] X. Liu and J. Katz, “Pressure–rate-of-strain, pressure diffusion, and velocity–pressure-gradient tensor measurements in a cavity flow,” AIAA Journal, vol. 56, no. 10, pp. 3897–3914, 2018.
  * [10] K. Hanjalić and B. Launder, Modelling turbulence in engineering and the environment: second-moment routes to closure.  Cambridge University Press, 2011.
  * [11] R. H. Bush, T. S. Chyczewski, K. Duraisamy, B. Eisfeld, C. L. Rumsey, and B. R. Smith, “Recommendations for future efforts in rans modeling and simulation,” in AIAA Scitech 2019 Forum, p. 0317, 2019.
  * [12] T. Sommer, R. So, and J. Zhang, “Modeling nonequilibrium and history effects of homogeneous turbulence in a stably stratified medium,” International Journal of Heat and Fluid Flow, vol. 18, no. 1, pp. 29–37, 1997.
  * [13] B. Eisfeld, “Characteristics of incompressible free shear flows and implications for turbulence modeling,” AIAA Journal, vol. 59, no. 1, pp. 180–195, 2021.
  * [14] E. P. Van Der Poel, R. Ostilla-Mónico, J. Donners, and R. Verzicco, “A pencil distributed finite difference code for strongly turbulent wall-bounded flows,” Computers & Fluids, vol. 116, pp. 10–16, 2015.
  * [15] R. Kunz, W.-S. Yu, S. Antal, and S. Ettorre, “An unstructured two-fluid method based on the coupled phasic exchange algorithm,” in 15th AIAA Computational Fluid Dynamics Conference, p. 2672, 2001.
  * [16] S. S. Clift and P. A. Forsyth, “Linear and non-linear iterative methods for the incompressible navier-stokes equations,” International Journal for Numerical Methods in Fluids, vol. 18, no. 3, pp. 229–256, 1994.
  * [17] J. P. Van Doormaal and G. D. Raithby, “Enhancements of the simple method for predicting incompressible fluid flows,” Numerical Heat Transfer, vol. 7, no. 2, pp. 147–163, 1984.
  * [18] C. Rhie and W. L. Chow, “Numerical study of the turbulent flow past an airfoil with trailing edge separation,” AIAA Journal, vol. 21, no. 11, pp. 1525–1532, 1983.
  * [19] C. G. Speziale, S. Sarkar, and T. B. Gatski, “Modelling the pressure–strain correlation of turbulence: an invariant dynamical systems approach,” Journal of Fluid Mechanics, vol. 227, pp. 245–272, 1991.
  * [20] B. Launder, “On the effects of a gravitational field on the turbulent transport of heat and momentum,” Journal of Fluid Mechanics, vol. 67, no. 3, pp. 569–581, 1975.
  * [21] M. Gibson and B. Launder, “Ground effects on pressure fluctuations in the atmospheric boundary layer,” Journal of Fluid Mechanics, vol. 86, no. 3, pp. 491–511, 1978.
  * [22] C. Shir, “A preliminary numerical study of atmospheric turbulent flows in the idealized planetary boundary layer,” Journal of Atmospheric Sciences, vol. 30, no. 7, pp. 1327–1339, 1973.
  * [23] B. J. Daly and F. H. Harlow, “Transport equations in turbulence,” Physics of Fluids, vol. 13, no. 11, pp. 2634–2649, 1970.
  * [24] B. E. Launder, G. J. Reece, and W. Rodi, “Progress in the development of a reynolds-stress turbulence closure,” Journal of Fluid Mechanics, vol. 68, no. 3, pp. 537–566, 1975.
  * [25] L. v. Haren, Etude théorique et modélisation de la turbulence en présence d’ondes internes.  PhD thesis, PhD thesis, Ecole Centrale de Lyon, 1993.
  * [26] B. Spencer and B. Jones, “Statistical investigation of pressure and velocity fields in the turbulent two-stream mixing layers,” AIAA Paper, vol. N 71-613, 1971.
  * [27] J. H. Bell and R. D. Mehta, “Development of a two-stream mixing layer from tripped and untripped boundary layers,” AIAA Journal, vol. 28, no. 12, pp. 2034–2042, 1990.
  * [28] M. M. Rogers and R. D. Moser, “Direct simulation of a self-similar turbulent mixing layer,” Physics of Fluids, vol. 6, no. 2, pp. 903–923, 1994.
  * [29] C. Pantano and S. Sarkar, “A study of compressibility effects in the high-speed turbulent shear layer using direct simulation,” Journal of Fluid Mechanics, vol. 451, p. 329, 2002.
  * [30] S. B. Pope, Turbulent Flows.  Cambridge University Press, 2000.
  * [31] S. G. Saddoughi and S. V. Veeravalli, “Local isotropy in turbulent boundary layers at high reynolds number,” Journal of Fluid Mechanics, vol. 268, pp. 333–372, 1994.


# Safety-Driven Battery Charging: A Fisher Information-guided Adaptive MPC
with Real-time Parameter Identification

Jorge Espin  Yuichi Kajiura  Dong Zhang  School of Aerospace and Mechanical
Engineering, The University of Oklahoma, Norman, OK 73019, USA (e-mail:
[jorge.espin, yuichi<EMAIL_ADDRESS>
###### Abstract

Lithium-ion (Li-ion) batteries are ubiquitous in modern energy storage
systems, highlighting the critical need to comprehend and optimize their
performance. Yet, battery models often exhibit poor parameter identifiability
which hinders the development of effective battery management strategies and
impacts their overall performance, longevity, and safety. This manuscript
explores the integration of Fisher Information (FI) theory with Model
Predictive Control (MPC) for battery charging. The study addresses the
inherent hurdles in accurately estimating battery model parameters due to
nonlinear dynamics and uncertainty. Our proposed method aims to ensure safe
battery charging and enhance real-time parameter estimation capabilities by
leveraging adaptive control strategies guided by FI metrics. Simulation
results underscore the effectiveness of our approach in mitigating parameter
identifiability issues, offering promising solutions for improving the control
of batteries during safe charging process.

###### keywords:

Fisher Information, Li-ion Battery, Parameter Identifiability, Charging
Control

##  1 Introduction

In recent years, the widespread adoption of battery-powered technologies,
spanning from electric vehicles to renewable energy storage systems, has
underscored the critical importance of safe and efficient battery charging.
Improper management during charging can lead to inefficient battery
utilization, diminished battery safety, and premature battery degradation
(Komsiyska et al. (2021)). Therefore, there exists a pressing need to design
advanced battery management systems to optimize charging processes. Among the
various techniques available in the literature, model predictive control (MPC)
stands out for its ability to leverage the underlying battery dynamics and
operational constraints within an optimization framework to achieve desired
objectives. This advanced control strategy has been validated by several
studies. Klein et al. (2011) proposed a nonlinear model predictive control
(NMPC) to rapidly charge batteries while ensuring safety throughout their
lifespan. Zou et al. (2018) introduced a battery charging methodology to
achieve an optimal balance between charging speed and battery State-of-Health
(SoH). This approach involves developing a simplified physics-based battery
model, which enables the design of a model-driven charging strategy.
Additionally, Pozzi et al. (2020) outlined a sensitivity-based MPC with
reduced computational cost suitable for battery packs comprising numerous
cells.

However, a common limitation of these charging control methods is their high
reliance on accurate battery models. Battery models, including empirical
models and physics-based models, often struggle with parametric uncertainties
caused by identifiability issues steaming from complex nonlinearities and
inherent parameter correlations. For instance, Forman et al. (2012) employed
Fisher information (FI) metrics to assess parameter identifiability in
electrochemical models, demonstrating that certain battery parameters cannot
be identified using traditional charge and discharge protocols. In this
context, identifiability refers to the ability to uniquely determine model
parameters from available data, a problem exacerbated by limited and noisy
experimental measurements (Guillaume et al. (2019)).

In response to poor battery parameter identifiability, researchers leveraged
optimal experimental design (OED), a strategy aimed at designing experiments
to maximize information about model parameters by carefully selecting input
profiles (Atkinson et al. (2007)). To name a few, Park et al. (2018) employed
a gradient-based algorithm for efficient parameter sensitivity analysis and
optimal fitting. Optimized current profiles demonstrate accelerated parameter
identification compared to standard methods. In the same way, Sharma and Fathy
(2014) used FI to quantify the identifiability of battery parameters for a
first-order nonlinear equivalent-circuit model subjected to periodic cycling.
Similarly, Lai et al. (2021) introduced an OED strategy for an electrochemical
battery model, focusing on enhancing the identifiability of diffusion
parameters through parameter sensitivity analysis. Likewise, Mendoza et al.
(2017) examined the feasibility of a single current cycling experiment to
identify the parameters of electrochemical and thermal models, and
Rothenberger et al. (2015) conducted OED based on input shaping for maximizing
identifiability using FI.

While these studies collectively underscored the potentials of OED in
advancing battery identifiability, it is crucial to acknowledge the inherent
challenges in its practical applicability. First, OED techniques, despite
their promises, often involve solving nonlinear and non-convex optimization
problems in an offline fashion that are both time-consuming and resource-
intensive, leading to practical implementation issues. Moreover, there exist
notable concerns regarding the resulting unrealistic current profiles that
might pose safety risks, such as overcharging, overheating, or battery failure
(Jaguemont and Bardé (2023)). This is largely due to aggressively pushing
batteries beyond safe limits in search of a maximized FI based on assumed
parameters that may not accurately represent the underlying variability and
complexity of battery dynamics.

These drawbacks reverberate into the MPC performance, where the reliance on
accurate models becomes a critical issue. A minor discrepancy between
predicted and actual battery behavior, resulting from inaccurate parameter
assumptions, directly impacts the efficacy of control strategies during
charging. Additionally, the limited adaptability of fixed-parameter models
hinders MPC’s ability to respond effectively to evolving battery conditions.

Given these critical considerations, the transition from offline OED for
battery identifiability to the online implementation of MPC for safe battery
charging presents complex challenges that need to be addressed holistically.
It is imperative to bridge the gap between offline design and online
implementation to optimize battery charging. To achieve this, enhancing
battery model accuracy through refined identifiability methodologies in OED
can serve as a robust foundation for effective MPC implementation; and by
integrating real-time adaptation mechanisms within MPC, such as adaptive
parameter estimation or dynamic model updating, the system can better
accommodate the inherent variability and complexity of battery dynamics.

In this regard, we propose a novel approach that integrates OED into an
adaptive MPC framework to safely fast charge a battery while enhancing
parameter identifiability in real-time. By synthesizing these two concepts and
addressing online computational challenges, we decompose the optimization
problem into two distinct phases: (i) an offline optimization phase aimed at
generating reference charging protocols with maximized parameter
identifiability; and (ii) an online control phase that leverages pre-computed
information for an adaptive charging control. Specifically, the offline
optimization phase focuses on maximizing FI while also fulfilling charging
objectives with safety constraints. This phase produces an optimized state of
charge (SOC) trajectory to be used in the online control phase. Subsequently,
in the online control phase, the adaptive MPC will track the pre-optimized SOC
trajectory from the offline optimization phase to calculate the real-time
current. This strategy aims to replicate in real-time the optimal input
current produced in the offline phase that maximizes FI. Notably, in this
work, instead of simply applying the offline current for battery charging,
MPC-based online control is essential for dynamically adjusting to changes
related to uncertain model parameters (Lin et al. (2014)). Hence, to enhance
adaptability and robustness against parameter uncertainty, we will incorporate
a real-time parameter updating scheme, which tracks model parameters to ensure
that the predictive model can accurately respond to evolving conditions.

To summarize, this work outlines the following technical contributions:

  * •

Novel integration of Fisher information into an adaptive MPC for battery
charging control, enhancing parameter identifiability and real-time
adaptability.

  * •

Emphasis on battery safe operation during charging and input shape flexibility
in maximizing parameter identifiability.

  * •

Inclusion of online parameter estimation within the MPC framework to improve
model predictive accuracy and system robustness.

To the authors’ best knowledge, this paper marks the first attempt in the
literature to integrate OED into a model predictive control framework for
improving real-time battery safe charging while taking parameter
identifiability into consideration. This work contributes to improving battery
safety, adaptability, and predictive accuracy during battery charging
operations.

##  2 Battery Equivalent Circuit Model

The equivalent circuit model (ECM) is a widely adopted battery representation.
This model simplifies the intricate electrochemical reactions into a network
of passive electrical components, such as resistors, capacitors, and a voltage
source (Hu et al. (2012)). By capturing the essential electrochemical
phenomena in a simplified yet effective manner, ECM facilitates real-time
simulation and control implementation in practical applications (Plett
(2015)). The model incorporates a open circuit voltage (OCV), which varies
with the SOC of the battery, an internal resistance ($R_{0}$) and a parallel
RC pair consisting of a resistance ($R_{1}$) and a capacitance ($C_{1}$). The
RC parallel circuit captures the transient behaviors resulting from charge-
transfer reactions at the electrode interface, with the time constant given by
$\tau=R_{1}C_{1}$ (Hu et al. (2012)). The dynamic equations governing the
model are as follows:

$\displaystyle\dot{z}(t)$ $\displaystyle=\frac{1}{Q}I(t),$ (1)
$\displaystyle\dot{Q_{c}}(t)$
$\displaystyle=-\frac{1}{R_{1}C_{1}}Q_{c}(t)+I(t),$ (2) $\displaystyle V(t)$
$\displaystyle=OCV(z(t))+\frac{1}{C_{1}}Q_{c}(t)+R_{0}I(t),$ (3)

where $z(t)$ represents the SOC of the battery, and $Q$ denotes the charge
capacity. Symbol $Q_{c}(t)$ indicates the amount of charge stored in the
parallel capacitor $C_{1}$. The input current is denoted by $I(t)$ (with
charging associated with positive current values), and the output terminal
voltage is represented by $V(t)$.

##  3 Fisher Information-based Adaptive Battery Charging Controller

In this section, we delve into the design of the Fisher information-guided
adaptive battery charging controller. Fig. 1 illustrates the proposed control
scheme, which comprises two main phases: an offline optimization phase (the
blue dashed-line block) and an online control phase (the green dashed-line
block).

This approach leverages both offline optimization and online control
techniques to maximize parameter identifiability and ensure safe battery
operation during charging phase. In addition, the proposed control strategy
allows for real-time parameter updating, see green dashed-line block in Fig.
1, which is beneficial for adaptively optimizing battery charging process
under evolving battery dynamics.

###  3.1 Fisher Information for Parameter Estimation

Battery parameter identification plays a pivotal role in accurately predicting
battery behavior and in turn, enables the establishment of optimal and precise
control policies that enhance battery performance. To maximize battery
parameter identifiability, FI theory has become a widely explored field which
offers a powerful toolset for analyzing the information content of
experimental data about unknown parameters. This theory aids in battery
parameter estimation by assessing the sensitivity of experimental voltage
measurements regarding model parameters. In this context, a higher FI value
leads to greater sensitivity, yielding that certain battery parameters might
be estimated more precisely based on specific experimental data (Fujita et al.
(2022)). To quantify FI, a key mathematical concept known as Fisher
information matrix (FIM) emerges. This metric provides a comprehensive measure
of parameter sensitivity by evaluating how variations in model parameters can
affect the likelihood of observed data. If the likelihood of data remains
relatively unchanged regardless of variations in the parameter value, then
data provide limited information about the parameter. Conversely, if the
likelihood of data is highly sensitive, then data contain substantial
information about the parameter, thereby enabling a more accurate parameter
identification (Chao et al. (2016)). Central to our study, we describe the
battery model (1)-(3) in a standard state-space representation:

$\displaystyle\dot{x}(t)$ $\displaystyle=f(x(t),u(t),\theta),$ (4)
$\displaystyle y(t)$ $\displaystyle=g(x(t),u(t),\theta)+v(t).$ (5)

Here, $x(t)$ represents the internal states of the battery, $u(t)$ denotes the
input current over time, $\theta$ describes the unknown battery parameters,
and $y(t)$ is defined as the measured terminal voltage corrupted by a noise
signal $v(t)$. Technically, the noise signal is inherently unknown; otherwise,
it would be theoretically feasible to draw a perfect voltage signal by
subtracting the noise from the measured voltage. Nevertheless, a crucial
assumption in FI analysis relies on having a prior knowledge of the noise
signal behavior (Sharma and Fathy (2014)). Accordingly, the experimentally
measured terminal voltage typically differs from the simulated values due to
inherent output voltage noise. This enables the calculation of a likelihood
function $p(y(t)|\theta)$, which quantifies the likelihood that the
discrepancy between the battery’s simulated output and the actual measured
voltage is a realization of the assumed noise process (Rothenberger et al.
(2015), and Sharma and Fathy (2014)).

Figure 1: Proposed Fisher-information-guided adaptive battery controller
scheme.

Building upon this likelihood function, the mathematical expression for FIM is
defined as the expected value of the Hessian matrix of the log-likelihood
function with respect to $\theta$, expressed as:

$\textbf{\text{FIM}}=\mathbb{E}\left\\{\left(\frac{\partial}{\partial\theta}\ln{p(y(t)|\theta)}\right)\left(\frac{\partial}{\partial\theta}\ln{p(y(t)|\theta)}\right)^{\intercal}\right\\}$
(6)

where, $\frac{\partial}{\partial\theta}\ln p(y(t)|\theta)$ represents the
score function, offering insights into how changes in $\theta$ impact the log-
likelihood function, while $\mathbb{E}$ denotes the expectation over the noise
distribution.

The battery model described by equations (1)–(3) uses a parameter set
$\theta=\left[R_{0},\frac{1}{C_{1}},\frac{1}{Q},\frac{1}{R_{1}C_{1}}\right]^{\intercal}$,
yielding the estimated output voltage $V(t,\theta)$.

The parameterized battery model of equations (1)–(3) is expressed as follows:

$\displaystyle\dot{z}(t,\theta)$ $\displaystyle=\theta_{3}I(t),$ (7)
$\displaystyle\dot{Q_{c}}(t,\theta)$
$\displaystyle=-\theta_{4}Q_{c}(t,\theta)+I(t),$ (8) $\displaystyle
V(t,\theta)$
$\displaystyle=OCV(z(t,\theta))+\theta_{2}Q_{c}(t,\theta)+\theta_{1}I(t).$ (9)

For mathematical simplicity, assuming the noise signal is Gaussian and white,
we derive the following FIM (Sharma and Fathy (2014)):

$\textbf{\text{FIM}}=\frac{1}{\sigma^{2}}\int_{0}^{T}\left[\left(\frac{\partial
V(t,\theta)}{\partial\theta}\right)\left(\frac{\partial
V(t,\theta)}{\partial\theta}\right)^{\intercal}\right]dt.$ (10)

In the given expression, $\sigma$ represents the standard deviation of the
measured terminal voltage. The sensitivity term $\frac{\partial
V(t,\theta)}{\partial\theta_{j}}=S_{j}$, in which $j=1,2,3,4$ indexes
parameters, and $T$ denotes the period over which OED is conducted.

We now detail the derivations of the sensitivity terms for computing the FIM.
The first two terms related to $\theta_{1}$ and $\theta_{2}$ are
straightforward and can be expressed as:

$\displaystyle S_{1}$ $\displaystyle=\frac{\partial
V(t,\theta)}{\partial\theta_{1}}=I(t),$ (11) $\displaystyle S_{2}$
$\displaystyle=\frac{\partial
V(t,\theta)}{\partial\theta_{2}}=Q_{c}(t,\theta).$ (12)

To compute the sensitivity of the third term ($S_{3}$), we apply the chain
rule:

$S_{3}=\frac{\partial V(t,\theta)}{\partial\theta_{3}}=\frac{\partial
V(t,\theta)}{\partial\text{OCV}(z(t,\theta))}\frac{\partial\text{OCV}(z(t,\theta))}{\partial
z(t,\theta)}\frac{\partial z(t,\theta)}{\partial\theta_{3}}.$ (13)

The first partial derivative is equal to 1, and
$\frac{\partial\text{OCV}(z(t,\theta))}{\partial
z(t,\theta)}=\alpha(z(t,\theta))$, where $\alpha(\cdot)$ represents the slope
of the OCV function. Next, to compute $\frac{\partial
z(t,\theta)}{\partial\theta_{3}}$, given the initial condition $z(0)=z_{0}$,
we use the relationship

$z(t,\theta)=z_{0}+\theta_{3}\int_{0}^{t}u(\tau)d\tau.$ (14)

This yields the partial derivative with respect to $\theta_{3}$ as:

$\frac{\partial z(t,\theta)}{\partial\theta_{3}}=\int_{0}^{t}u(\tau)d\tau.$
(15)

Considering the relation:

$\int_{0}^{t}u(\tau)d\tau=\theta_{3}^{-1}(z(t,\theta)-z_{0}).$ (16)

Combining these expressions, we arrive at:

$S_{3}=\alpha(z(t,\theta))\theta_{3}^{-1}(z(t,\theta)-z_{0}).$ (17)

To obtain the last sensitivity term ($S_{4}$), we apply the chain rule,
resulting in:

$S_{4}=\frac{\partial V(t,\theta)}{\partial\theta_{4}}=\frac{\partial
V(t,\theta)}{\partial Q_{c}(t,\theta)}\frac{\partial
Q_{c}(t,\theta)}{\partial\theta_{4}}.$ (18)

It is straightforward to determine that $\frac{\partial V(t,\theta)}{\partial
Q_{c}(t,\theta)}=\theta_{2}$, and we can handle the last partial derivative
expression using the Laplace transform in equation (8) as follows:

$Q_{c}(s,\theta)=\frac{1}{s+\theta_{4}}I(s).$ (19)

Next, we calculate the partial derivative with respect to $\theta_{4}$:

$\frac{\partial
Q_{c}(s,\theta)}{\partial\theta_{4}}=\frac{-1}{s^{2}+2\theta_{4}s+\theta_{4}^{2}}I(s).$
(20)

The resulting transfer function $G(s,\theta)$ can be represented in the time
domain by a linear state-space model using the controllable canonical form:

$G(s,\theta)=\frac{-1}{s^{2}+2\theta_{4}s+\theta_{4}^{2}}.$ (21)

###  3.2 Offline Optimization Problem Formulation

Within the scope of offline optimization problem formuation, the theory
suggests that a larger FIM is preferred as it enhances the accuracy of
parameter identification (Rothenberger et al. (2015)). To pursue this goal, we
can maximize certain scalar metrics such as the matrix trace or its
determinant. In our study, we focus on optimizing the determinant of the FIM.
This approach is widely acknowledged in the literature as “D-optimal”
experimental design (Atkinson et al. (2007) and Forman et al. (2013)).

We formulate our offline optimization problem as follows:

[3] I^off_t_i,…, I^off_t_N-det(FIM) c_1:battery dynamics in (1)-(3),
c_2:z(t_f, θ)=1, c_3: I_min≤I^off(t_i)≤I_max, i=1,…,N. c_4: V(t_i, θ)≤V_max,
i=1,…,N.

The objective of the offline optimization problem is to maximize the
determinant of FIM to improve battery parameter identifiability during safe
battery charging. However, it is crucial to acknowledge that previous studies
(Rothenberger et al. (2014) and Sharma and Fathy (2014)) often overlook
critical safety constraints, such as battery overcharge (constrained by
$c_{2}$ and $c_{4}$), which is essential for avoiding potential battery damage
like explosions or premature aging.

In our approach, we consider controlled rapid charging within a user-defined
time frame $[0,t_{f}]$, under constraint $c_{2}$, aiming to ensure that
battery reaches full charge at a specific time $t_{f}$ while simultaneously
maximizing FI. Additionally, our method also offers increased flexibility in
charge current shaping, enabling a more dynamic charging profile beyond
traditional input-shaping techniques (Rothenberger et al. (2014), Rothenberger
et al. (2015) and Mendoza et al. (2017)). This enhanced flexibility improves
sensitivity in the observed data, contributing to more accurate parameter
estimation.

Mathematically, maximizing the determinant of FIM can be considered as a non-
convex optimization problem (Rothenberger et al. (2015) and De Cock et al.
(2016)) because of the nonlinearity of OCV curve, leading to an unbounded
current profile. To address this challenge, existing literature suggests
incorporating additional objectives in the cost function, such as including
the l2-norm of the current signal, to ensure that the resulting current
profile remains reasonable and bounded (Rothenberger et al. (2014)), which may
compromise the importance of maximizing FI. Conversely, our proposal mitigates
this issue by establishing constrained charging defined by $c_{2}$ and
bounding the optimal current profile using $c_{3}$. To navigate the non-convex
landscape and avoid local minima, we utilize genetic algorithms employing
crossover, mutation, and selection strategies. These methods efficiently
explore the multi-modal cost function surface, enhancing the search for global
optima (Bajpai and Kumar (2010)). Ultimately, this comprehensive offline
method not only focuses on enhancing identifiability but also prioritizes
battery safety and charging efficiency.

The resulting optimal trajectories from this offline optimization step are
utilized in the subsequent online control to ensure safe and efficient battery
operation.

###  3.3 Online Charging Control Using Adaptive MPC

This section explores adaptive MPC to ensure that battery operational
constraints are met while performing real-time parameter identification. By
leveraging the pre-optimized SOC trajectory from the offline optimization,
this approach is specifically designed to enhance battery parameter
identifiability and safety simultaneously during online operation.

Adaptive MPC (see the green dashed-line block in Fig. 1) is an advanced
control strategy that optimizes control actions over a finite time horizon $H$
using a predictive model (Schwenzer et al. (2021)). In battery charging, MPC
ensures SOC tracking while adhering to safety constraints to prevent battery
damage (Pozzi and Toti (2023)). Mathematically, the battery charging problem
can be framed as a constrained optimization, as illustrated below:

[3] I^on_t_k+1,…,I^on_t_k+H ∑_i=k+1^k+H ∥ z_ref(t_i) - ^z(t_i,^θ) ∥^2
d_1:battery dynamics in (1)-(3), d_2: z(t_i,^θ) ≤1,i=k+1,…,k+H. d_3:
I_min≤I^on(t_i)≤I_max, i=k+1,…,k+H. d_4: V(t_i, ^θ)≤V_max, i=k+1,…,k+H.

Notably, $z_{\text{ref}}(t_{i})$ represents the pre-optimized SOC trajectory
produced from the offline optimization phase in Section 3.2. This pre-
optimized SOC trajectory, as previously discussed, was generated by maximizing
FI over the entire time horizon $[0,t_{f}]$. By tracking
$z_{\text{ref}}(t_{i})$ in (3.3), the solution produces a dynamical online
charging current in an effort to mimic the current profile from the offline
optimization so as to maximize FI and improve parameter estimation in real
time. Different from other existing works, MPC in this paper incorporates
real-time parameter updating. This methodology minimizes the discrepancy
between actual battery voltage measurements and the corresponding estimated
voltages by iteratively updating the estimated parameters $\hat{\theta}$
within specified bounds. This adaptive approach improves model accuracy,
allowing for the adaptation to varying battery conditions over time (e.g.,
aging) while still adhering to operational battery constraints.

Remark. Parametric uncertainties in OED studies can lead to sub-optimal
trajectories when incorrect nominal model parameters are used during the
offline design phase (Rothenberger et al. (2014), Mendoza et al. (2017), Rojas
et al. (2007)). Consequently, the online current profile computed by the
adaptive MPC strategy might differ from the pre-optimized offline profile due
to discrepancies between the true battery parameters and the nominal
parameters assumed during the offline phase. However, for this study, online
parameter adaptation within the MPC framework mitigates the influence of
parametric uncertainties by updating parameter values in real-time to
compensate for model inaccuracies, maintaining optimality despite possible
offline design deficiencies.

Future work aims to enhance pre-optimized trajectory reliability by
investigating advanced offline design techniques like Bayesian optimal design
(Rothenberger and Fathy (2015)), which incorporates prior parameter
distribution statistics into the offline optimization design, potentially
reducing the existing offline-online optimal profile discrepancy.

##  4 Simulation Results

As part of this analysis, we present the results of our adaptive MPC approach
utilizing Fisher information for real-time parameter identification during
safe battery charging.

###  4.1 Problem Setup

As has been widely acknowledged in the literature, a critical issue in OED is
that a set of nominal parameters must be presumed known a priori for
maximizing FI (Rothenberger et al. (2014), Mendoza et al. (2017), Rojas et al.
(2007)). In this work, the nominal battery parameters are given by Table 1.
However, this assumption may not hold valid in practice, as actual battery
parameters may deviate from the nominal values initially assumed. Therefore,
we consider the following true parameters for the battery model:
$R_{0,true}=0.05~{}\Omega,~{}C_{1,true}=950~{}F,~{}Q_{true}=1.9~{}\text{Ah}$,
while $R_{1,true}=0.03~{}\Omega$. This adjustment enables us to demonstrate
how the adaptive charging controller converges in real-time toward these true
battery parameters to ensure the accuracy and reliability of the battery model
over time. Moreover, the nonlinear OCV function used in this study is derived
from experimental data provided by the Center for Advanced Life Cycle
Engineering (CALCE (2024)), using 1/20 C-rate on a Samsung INR-18650 20R cell.
The OCV curve is modeled using a 7th-order polynomial function with respective
coefficients ($a_{0},\dots,a_{7}$) given in Table 1.

Table 1: Nominal parameters and coefficients of OCV curve Parameter | Value | Parameter | Value
---|---|---|---
$Q$ | 2.0 Ah | $a_{2}$ | \- 77.31237098
$R_{1}$ | 0.03 $\Omega$ | $a_{3}$ | 327.446181
$C_{1}$ | 1000 $F$ | $a_{4}$ | \- 763.3324119
$R_{0}$ | 0.06 $\Omega$ | $a_{5}$ | 988.408671
$a_{0}$ | 3.039475779 | $a_{6}$ | \- 662.9843922
$a_{1}$ | 9.620312047 | $a_{7}$ | 179.301862

Key considerations to highlight include the following: We are implementing a
constrained fast charging protocol aimed at achieving a full battery charge
within 30 minutes ($t_{f}=1800$ seconds). The dynamic charging current is
constrained to -3C to 3C, with a voltage limit of 4.3 V to prevent
overcharging. We employ a GA with a population size of 50 individuals. The
optimization is conducted on a high-performance computer with an Intel Core i9
processor with 64 GB of RAM and 24 cores, also leveraging parallel processing
to expedite solution convergence within the MATLAB software environment.
Furthermore, the adaptive MPC utilizes fmincon to solve the optimization
problem. The same computational processor and software environment are
utilized.

###  4.2 Discussion and Future Work

Figure 2: Comparison of offline optimized trajectories with online adaptive
MPC responses for battery charging and real-time parameter updating.

Our results in Fig. 2.a) demonstrate the controller’s capability to safely
charge the battery within a 30-minute interval, maximizing the FI metric for
real-time parameter estimation. The high tracking accuracy depicted in Fig.
2.a) significantly influences the observed discrepancy in Fig. 2.c),
attributed to practical differences in battery model parameters between
offline optimization and online control. Throughout the charging period, Fig.
2.c) confirms consistent adherence to safety constraints, including the
maximum voltage limit, ensuring battery operation remains within a safe zone
to prevent potential battery damage.

Fig. 2.b) shows the difference between the offline optimized current profile
and the online adaptive MPC current profile due to the discrepancy between the
assumed nominal battery parameters in the offline optimization and the true
parameter values used online. This discrepancy is further evident from the
contrast in the determinant of the FIM calculated offline ($1.6275\times
10^{40}$) and online ($3.8596\times 10^{39}$) using the true parameter values.
The smaller online FI value is expected due to practical variability in
battery parameters during real-time operation. However, the adaptability of
online adaptive MPC, which dynamically adjusts to evolving parameter values,
may compensate for the reduced FI value, leading to improved parameter
estimation accuracy over time, demonstrating the method’s practical utility in
refining model estimates and enhancing system performance despite parameter
variability.

The proposed FI-based controller demonstrates its capability to accurately
identify and converge towards the true battery model parameters in real-time,
as illustrated in Figs. 2.d), 2.e), and 2.f). This real-time parameter
updating feature contributes to enhancing the reliability and accuracy of the
battery model over time, ensuring that the controller can effectively mitigate
parametric uncertainties by adapting to evolving conditions.

Future directions of our work highlight the following insights: (i) accounting
for the dependency of battery model parameters on temperature, SOC, or aging;
(ii) enhancing the offline optimization phase by incorporating more advanced
optimal design techniques, rather than relying solely on nominal values, to
better capture the underlying parameter variability and uncertainty; (iii)
transitioning towards more sophisticated battery models, such as physics-based
models, to leverage additional features of this approach, including the
consideration of battery degradation during charging; and (iv) validating our
method with experimental data to demonstrate its robustness in a more
demanding environment.

##  5 Conclusions

This manuscript explores a pioneering solution of integrating FI into an
online control framework for safe battery charging. With the objective of
minimizing computational complexity we consider partitioning the optimization
process into offline and online phases. The findings illustrate a constrained
safe charging protocol capable of maximizing FI in real time. This underscores
the viability of this approach for real-time parameter identification in Li-
ion batteries characterized by inherently poor parameter identifiability.
Future directions for this research involve transitioning towards more
sophisticated electrochemical models to capture detailed battery behavior and
conducting experimental validations to showcase the effectiveness of Fisher
Information-guided adaptive battery controller.

## Acknowledgement

This material is based upon work supported by the National Science Foundation
(NSF) under Award No. 2327327. Any opinions, findings, and conclusions or
recommendations in this paper are those of the author(s) and do not
necessarily reflect the views of the NSF.

## References

  * Atkinson et al. (2007) Atkinson, A., Donev, A., and Tobias, R. (2007).  _Optimum experimental designs, with SAS_ , volume 34.  OUP Oxford.
  * Bajpai and Kumar (2010) Bajpai, P. and Kumar, M. (2010).  Genetic algorithm–an approach to solve global optimization problems.  _Indian Journal of computer science and engineering_ , 1(3), 199–206.
  * CALCE (2024) CALCE (2024).  Battery data.  URL https://calce.umd.edu/battery-data.
  * Chao et al. (2016) Chao, J., Ward, E.S., and Ober, R.J. (2016).  Fisher information theory for parameter estimation in single molecule microscopy: tutorial.  _JOSA A_ , 33(7), B36–B57.
  * De Cock et al. (2016) De Cock, A., Gevers, M., and Schoukens, J. (2016).  D-optimal input design for nonlinear fir-type systems: A dispersion-based approach.  _Automatica_ , 73, 88–100.
  * Forman et al. (2013) Forman, J., Stein, J., and Fathy, H. (2013).  Optimization of dynamic battery paramter characterization experiments via differential evolution.  In _2013 American Control Conference_ , 867–874. IEEE.
  * Forman et al. (2012) Forman, J.C., Moura, S.J., Stein, J.L., and Fathy, H.K. (2012).  Genetic identification and fisher identifiability analysis of the doyle–fuller–newman model from experimental cycling of a lifepo4 cell.  _Journal of Power Sources_ , 210, 263–275.
  * Fujita et al. (2022) Fujita, K., Okada, K., and Katahira, K. (2022).  The fisher information matrix: A tutorial for calculation for decision making models.  _PsyArXiv_.
  * Guillaume et al. (2019) Guillaume, J.H., Jakeman, J.D., Marsili-Libelli, S., Asher, M., Brunner, P., Croke, B., Hill, M.C., Jakeman, A.J., Keesman, K.J., Razavi, S., et al. (2019).  Introductory overview of identifiability analysis: A guide to evaluating whether you have the right type of data for your modeling purpose.  _Environmental Modelling & Software_, 119, 418–432.
  * Hu et al. (2012) Hu, X., Li, S., and Peng, H. (2012).  A comparative study of equivalent circuit models for li-ion batteries.  _Journal of Power Sources_ , 198, 359–367.
  * Jaguemont and Bardé (2023) Jaguemont, J. and Bardé, F. (2023).  A critical review of lithium-ion battery safety testing and standards.  _Applied Thermal Engineering_ , 121014.
  * Klein et al. (2011) Klein, R., Chaturvedi, N.A., Christensen, J., Ahmed, J., Findeisen, R., and Kojic, A. (2011).  Optimal charging strategies in lithium-ion battery.  In _Proceedings of the 2011 american Control Conference_ , 382–387. IEEE.
  * Komsiyska et al. (2021) Komsiyska, L., Buchberger, T., Diehl, S., Ehrensberger, M., Hanzl, C., Hartmann, C., Hölzle, M., Kleiner, J., Lewerenz, M., Liebhart, B., et al. (2021).  Critical review of intelligent battery systems: Challenges, implementation, and potential for electric vehicles.  _Energies_ , 14(18), 5989.
  * Lai et al. (2021) Lai, Q., Ahn, H.J., Kim, Y., Kim, Y.N., and Lin, X. (2021).  New data optimization framework for parameter estimation under uncertainties with application to lithium-ion battery.  _Applied Energy_ , 295, 117034.
  * Lin et al. (2014) Lin, X., Perez, H.E., Mohan, S., Siegel, J.B., Stefanopoulou, A.G., Ding, Y., and Castanier, M.P. (2014).  A lumped-parameter electro-thermal model for cylindrical batteries.  _Journal of Power Sources_ , 257, 1–11.
  * Mendoza et al. (2017) Mendoza, S., Rothenberger, M., Liu, J., and Fathy, H.K. (2017).  Maximizing parameter identifiability of a combined thermal and electrochemical battery model via periodic current input optimization.  _IFAC-PapersOnLine_ , 50(1), 7314–7320.
  * Park et al. (2018) Park, S., Kato, D., Gima, Z., Klein, R., and Moura, S. (2018).  Optimal experimental design for parameterization of an electrochemical lithium-ion battery model.  _Journal of The Electrochemical Society_ , 165(7), A1309.
  * Plett (2015) Plett, G.L. (2015).  _Battery management systems, Volume I: Battery modeling_.  Artech House.
  * Pozzi et al. (2020) Pozzi, A., Torchio, M., Braatz, R.D., and Raimondo, D.M. (2020).  Optimal charging of an electric vehicle battery pack: A real-time sensitivity-based model predictive control approach.  _Journal of Power Sources_ , 461, 228133.
  * Pozzi and Toti (2023) Pozzi, A. and Toti, D. (2023).  Imitation learning for agnostic battery charging: A dagger-based approach.  _IEEE Access_.
  * Rojas et al. (2007) Rojas, C.R., Welsh, J.S., Goodwin, G.C., and Feuer, A. (2007).  Robust optimal experiment design for system identification.  _Automatica_ , 43(6), 993–1008.
  * Rothenberger et al. (2014) Rothenberger, M.J., Anstrom, J., Brennan, S., and Fathy, H.K. (2014).  Maximizing parameter identifiability of an equivalent-circuit battery model using optimal periodic input shaping.  In _Dynamic Systems and Control Conference_ , volume 46186, V001T19A004. American Society of Mechanical Engineers.
  * Rothenberger et al. (2015) Rothenberger, M.J., Docimo, D.J., Ghanaatpishe, M., and Fathy, H.K. (2015).  Genetic optimization and experimental validation of a test cycle that maximizes parameter identifiability for a li-ion equivalent-circuit battery model.  _Journal of Energy Storage_ , 4, 156–166.
  * Rothenberger and Fathy (2015) Rothenberger, M.J. and Fathy, H.K. (2015).  Robust bayesian sequential input shaping for optimal li-ion battery model parameter identifiability.  In _Dynamic Systems and Control Conference_ , volume 57250, V002T23A008. American Society of Mechanical Engineers.
  * Schwenzer et al. (2021) Schwenzer, M., Ay, M., Bergs, T., and Abel, D. (2021).  Review on model predictive control: An engineering perspective.  _The International Journal of Advanced Manufacturing Technology_ , 117(5), 1327–1349.
  * Sharma and Fathy (2014) Sharma, A. and Fathy, H.K. (2014).  Fisher identifiability analysis for a periodically-excited equivalent-circuit lithium-ion battery model.  In _2014 American control conference_ , 274–280. IEEE.
  * Zou et al. (2018) Zou, C., Manzie, C., and Nešić, D. (2018).  Model predictive control for lithium-ion battery optimal charging.  _IEEE/ASME Transactions on Mechatronics_ , 23(2), 947–957.