Eron T V Dauricio

Carlos Junqueira-Junior

Diego F Abreu
email: mecabreu@yahoo.com.br

João Luiz

F Azevedo

EXTERNAL LAMINAR BOUNDARY LAYER SIMULATIONS USING A HIGH-FIDELITY WALL-MODELING APPROACH

Keywords: Wall-Modeling, Wall-Modeled Large Eddy Simulation, External Laminar Flow, Boundary Layer Flow

Wall-Modeled Large Eddy Simulation (WMLES) is a well-stablished technique for obtaining high-fidelity solutions of turbulent, high Reynolds number flows, with reasonably acceptable computational costs. However, for external flows, the very thin laminar boundary layer developing near the body leading edge imposes quite restrictive mesh resolution requirements, leading to prohibitively high computational costs for practical Reynolds numbers. We propose a wall-modeling approach for the laminar portion of the boundary layer in order to alleviate these costs by reducing the aforementioned mesh resolution requirements. The wall model is based on local self-similar solutions of the boundary layer, and is implemented in the same context of wall-stress models in the WMLES approach. An assessment of the model is done in terms of both pressure and skin friction coefficient distributions along the surface, for an incompressible, fully laminar flow around a NACA 0012 airfoil geometry, with a chord Reynolds number of Re c = 4.5 × 10 3 . The results obtained in the simulations using the proposed model are in good agreement with the reference solution, demonstrating the feasibility of the model for external laminar flows.

INTRODUCTION

Large Eddy Simulations (LES) have now greatly superseded the long-standing Reynolds-Averaged Navier-Stokes (RANS) approach as a research tool for the study of turbulent flow dynamics in academic and research centers, primarily due to the higher accuracy and fidelity of the resulting simulations. The higher accuracy and resolution of the former is achieved by directly resolving the energetic scales of the turbulent motions, whereas in the latter, all the turbulent motions are modeled. As a consequence, the computational cost of the LES approach is much greater than that of RANS, becoming prohibitively high for practical engineering configurations and high Reynolds number flows.

In the LES approach, only the larger, energy-containing scales are resolved, while the effects of the smaller ones are modeled. However, near the wall, where the boundary layer develops and the turbulent scales are very small, proportional to the viscous length, δ ν , there is a peak in turbulent energy production [START_REF] Pope | Turbulent Flows[END_REF][START_REF] Schiavo | Turbulent kinetic energy budgets in wall bounded flows with pressure gradients and separation[END_REF] that must be captured by the LES simulation. Since the smallest scales become smaller as the Reynolds number increases, the computational cost of resolving for the near-wall region becomes prohibitive for high Reynolds number flows. This situation configures the so-called "LES near-wall problem".

During the past few decades, many different solutions have been proposed to tackle the LES near-wall problem, all of which incorporates some way of modeling the effects of the near-wall motions onto the overall turbulent flow, instead of directly resolving it. Among such modeling propositions are the Detached-Eddy Simulation (DES) approach, hybrid RANS/LES approach and wall-stress modeled LES [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large reynolds numbers[END_REF][START_REF] Schumann | Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli[END_REF][START_REF] Piomelli | Wall-layer models for large-eddy simulation[END_REF][START_REF] Spalart | Detached-eddy simulation[END_REF]. These methods are all considered forms of the Wall-Modeled Large Eddy Simulation (WMLES) approach, as opposed to the traditional Wall-Resolved Large Eddy Simulation (WRLES).

Earlier estimates of grid point requirements [START_REF] Choi | Grid-point requirements for large eddy simulation: Chapman's estimates revisited[END_REF] show that, for the WMLES approach, a scaling of N wm ∝ Re L is achieved, compared to N wr ∝ Re 13/7 L and N DN S ∝ Re 37/14 L for WRLES and DNS computations, respectively. In addition, methods within the WMLES approach break the scaling of computational cost with Re τ , the Reynolds number based on the friction velocity, u τ = τ w /ρ, since in these methods the inner layer is completely modeled, and the cost of resolving the outer layer is independent of Re τ [START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF]. Therefore, the WMLES approach have the potential of achieving high fidelity, high accuracy results, with a much lower computational cost, making such simulations feasible for high Reynolds number flows.

As promising as the WMLES approach may be, there are still some challenges to be faced. One problem with the method that has attracted some attention only recently is the cost of accurately resolving the complete boundary layer in external flows. The current common practice is to employ some form of tripping in the body leading edge to force the boundary layer to be fully turbulent along the surface. However, if WMLES is to be used in contexts other than turbulent research, such as in the engineering design process, then, it is paramount that the developing laminar boundary layer is also accurately resolved, including the region where transition to turbulence occurs. [START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF] estimated the cost of resolving the boundary layer flow along the surface of a NACA0012 airfoil. Even if the WMLES is employed, meaning that only the outer boundary layer is resolved, the number of grid points is massive near the leading edge, where the developing laminar boundary layer is very thin, being proportional to 1/δ 2 (x).

In another study, it has been estimated that the cost of accurately solving the very thin laminar boundary layer near the leading edge of an airfoil may be 10 to 100 times higher than the cost of solving the modeled turbulent region [START_REF] Slotnick | CFD Vision 2030 study: A path to revolutionary computational aerosciences[END_REF]. These estimates show that, for high Reynolds number external flows, even the WMLES is still unfeasible due to high computational costs.

In this work, we propose a wall-modeling strategy for the laminar region of the boundary layer in order to alleviate the computational cost of accurately solving the thin boundary layer near the body leading edge in external flows. The laminar modeling procedure described in this work is applied in the same way as wall-stress models in the WMLES approach, since these are intended to be combined as seamlessly as possible into a single, complete flow simulation process.

NUMERICAL METHODOLOGY

Governing Equations

The governing equations to be solved numerically are the compressible Navier-Stokes equations, which may be written in conservation form and compact vector notation as

∂Q ∂t + ∇ • F(Q, ∇Q) = 0 (1)
where Q = [ρ, ρu i , ρE] T is the vector of conserved variables, and F is the flux vector, given by F = F e -F v , with the superscripts e and v indicating advective, or Euler, and viscous fluxes, respectively. These may be further written as

F e i =       ρu i ρu 1 u i + δ 1i p ρu 2 u i + δ 2i p ρu 3 u i + δ 3i p (ρE + p) u i       , F v i =       0 τ 1i τ 2i τ 3i u j τ ij -q i       , i = 1, 2, 3 (2)
In the above equations, ρ is the density, u i are the components of the velocity vector, and ρE is the total energy per unit volume. In order to obtain closure of the system of equations, the total energy per unit volume is written in terms of the thermodynamic pressure using the thermal and caloric equations of state for an ideal gas, that is, p = ρRT and e = c v T , respectively. Then, one can write

ρE = ρe + 1 2 ρu i u i = p (γ -1) + 1 2 ρu i u i
where e is the internal energy per unit mass, and γ and R are the ratio and the difference between specific heats, respectively.

The viscous stress tensor, τ ij , may be written as

τ ij = µ ∂u i ∂x j + ∂u j ∂x i + λ ∂u k ∂x k δ ij (3)
where µ and λ are the dynamic and the second coefficients of viscosity, related by Stokes' hypothesis as λ = -2 3 µ. The heat flux vector, q i , in Eq. ( 2) is given by

q i = -k ∂T ∂x i
with k as the thermal conductivity coefficient.

Numerical Method

In this work, Eq. ( 1) is numerically solved through a nodal implementation of the discontinuous Galerkin (DG) scheme with split formulations for the inviscid fluxes. We start by discretizing the physical domain, Ω, into K non-overlapping elements, Ω k , such that

Ω ≈ Ω h ≡ K k=1 Ω k , K k=1 Ω k = ∅
where Ω h is the computational domain.

The elements in physical domain are then mapped onto a standard element, Ω st , where each shape in physical domain have a corresponding shape in the standard domain. Although the shape of the discretizing elements are arbitrary, in this work we use only hexahedral elements, primarily due to computational efficiency and ease of implementation. For this shape, the standard element is defined by Ω st = [-1, 1] 3 . Applying this mapping to Eq. ( 1) results in

J ∂Q ∂t + ∇ ξ • F ≡ J ∂Q ∂t + ∇ ξ • e ξ • F = 0 (4)
where ∇ ξ is the divergence with respect to the standard element coordinates, ξ = (ξ, η, ζ), J = |∂x/∂ξ| is the Jacobian of the coordinate transformation, e ξ is the contravariant basis vector, and F = e ξ • F the contravariant flux vector.

For the nodal implementation of the discontinuous Galerkin scheme, the solution within each element, in standard coordinates, is approximated by a polynomial interpolation given by

Q ≈ Q h (ξ) = N p,q,r=0 Q h (ξ p , η q , ζ r , t)ϕ pqr (ξ) (5)
where ϕ pqr is the interpolating polynomial. For hexahedral elements, this polynomial is given by a tensor product of 1-D Lagrange polynomials, that is,

ϕ pqr (ξ) = ℓ p (ξ)ℓ q (η)ℓ r (ζ), ℓ p = Np i=0 i̸ =p ξ -ξ i ξ p -ξ i
with equivalent definitions for the other directions.

In the DG formulation, the residual of the approximation, obtained by substitution of Eq. ( 5) into Eq. ( 4), is required to be orthogonal to a test function, ψ, where the space of this test function is the same as the space spanned by the trial, or solution function. This leads to the following formulation

Ωst J ∂Q h ∂t ψdξ + Ωst (∇ ξ • F e ) ψdξ + ∂Ωst (F e,⋆ -F e ) ψdS ξ - Ωst (∇ ξ • F v ) ψdξ = 0, (6)
where ∂Ω st defines the boundary of the element, and F e,⋆ is a numerical flux function defined on this boundary. The formulation in Eq. ( 6) is obtained after integration by parts of the second term of Eq. (4) two times, and is called the strong formulation. The weak formulation is obtained when integration by parts is done only once. The integrals in Eq. ( 6) are carried out through Gaussian quadratures, where the details may be consulted in [START_REF] Hindenlang | Explicit discontinuous Galerkin methods for unsteady problems[END_REF].

As mentioned above, a split formulation of the inviscid fluxes, F e , is implemented in order to enhance stability of the method, especially for higher order interpolations. The details of this formulation are not described here due to brevity, but may be consulted in the comprehensive exposition of [START_REF] Gassner | Split form nodal discontinuous Galerkin schemes with summationby-parts property for the compressible Euler equations[END_REF]. The specific split form of the fluxes employed throughout the simulations in this work is that of [START_REF] Pirozzoli | Numerical methods for high-speed flows[END_REF]. An approximate Riemann solver, namely the Roe solver with an entropy fix [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF], is used to obtain the numerical flux F e,⋆ on the element boundaries.

Since weak variational formulations using discontinuous function spaces are not suitable for directly discretizing higher order derivatives, as those appearing in the divergence term of the viscous fluxes in Eq. ( 4) through the viscous stress tensor, Eq. ( 3), a mixed formulation must be used for these terms. In this work, we make use of the lifting scheme of [START_REF] Bassi | A high-order discontinuous Galerkin method for compressible turbulent flows[END_REF].

The semidiscrete equations obtained after full discretization of Eq. ( 6) according to the aforementioned methods are advanced in time using the five stage, fourth-order, low-storage Runge-Kutta scheme of [START_REF] Carpenter | Fourth-order 2n-storage Runge-Kutta schemes[END_REF]. This work uses the computational implementation of the FLEXI software [START_REF] Krais | FLEXI: A high order discontinuous galerkin framework for hyperbolic-parabolic conservation laws[END_REF] in order to apply the numerical scheme and discretization methods just described.

WALL-STRESS MODELING

Wall-Modeling Framework

This section briefly outlines the wall-modeling procedure employed in wall-stress models for LES computations, since the laminar wall model proposed in this work is implemented within this same context. For a comprehensive exposition and further details, the reader is referred to previous reviews, e.g., [START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF] and [START_REF] Piomelli | Wall-layer models for large-eddy simulation[END_REF].

The main idea of the WMLES approach is to model the inner part of the boundary layer, which scales with the viscous length, δ v = ν/u τ , and directly resolve only the outer layer, which scales with the boundary layer thickness, δ. In doing so, grid requirements near the wall are greatly relaxed, reducing the computational cost by a great amount, as mentioned in previous sections. Moreover, as the dynamics of the inner layer is under-resolved, the approach further reduces the computational cost by reducing the temporal stiffness induced by the viscous sublayer, meaning that larger time-step sizes may be used [START_REF] Bose | Wall-modeled large-eddy simulation for complex turbulent flows[END_REF].

The effects of the inner layer onto the overall boundary layer dynamics are accounted for by prescribing the wall-shear stress, τ w , as a boundary condition. From a numerical and mathematical point of view, the approach basically replaces a Dirichlet boundary condition on velocity, where the no-slip condition states that u w = 0, by a Neumann-type boundary condition on the velocity gradients through the computed wall-shear stress -see Eq. ( 3).

Although the peak turbulence production in the inner layer is not captured by this approach, it is still physically consistent with the fact that, in high Reynolds number flows, the Reynolds stresses are predominantly produced at the same wall distance as they are later dissipated [START_REF] Hoyas | Scaling of the velocity fluctuations in turbulent channels up to Re τ = 2003[END_REF] and, hence, the method is actually resolving the majority of the boundary layer turbulence production and dissipation. The WMLES approach is also consistent with Townsend's similarity hypothesis, which states that the shear stresses and the no-transpiration conditions at the wall are sufficient to define the outer flow of the boundary layer [START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF].

In the context of wall-stress models, in order to calculate a wall-shear stress for boundary condition imposition, the LES solution at a distance h wm away from the wall in the normal direction is used as the input for a wall model, which, in turn, computes a wall-shear stress and provides this information back to the LES simulation. This process is illustrated in Fig. 1. The distance from the wall where the LES solution is taken is termed the wall-model height or, similarly, the exchange location, where information between the wall model and the LES solution is exchanged. Note that the mesh may be unstructured, with arbitrary shape elements near the wall, and a structured-like grid is shown in Fig. 1 The information taken from the LES solution is usually the velocity field, u, but may be any other relevant quantity, such as pressure gradients, temperature, and more. If the energy equation is part of the dynamical system being solved, then the imposed boundary condition is augmented with the heat flux, q w , which is also computed by the model. It should be noted that in this wall-modeling methodology, the velocity at the wall is allowed to slip in the wall-tangent direction, and only the no-transpiration condition is used, where the velocity component normal to the wall vanishes.

Throughout the decades, many wall models have been proposed, from simple algebraic models to more complex, PDE-based ones [START_REF] Bose | Wall-modeled large-eddy simulation for complex turbulent flows[END_REF]. The models included in the latter category make use of the thin boundary layer equations (TBLE), where a PDE must be solved in the region near the wall in order to compute a wall-shear stress. Such models are very general, with only mild underlying assumptions, but are usually difficult to implement. On the other hand, algebraic models have a straightforward implementation, but lacks generality, where it is usually assumed that the boundary layer is in a state of dynamical equilibrium.

Wall Model for Laminar Boundary Layers

As mentioned in a previous section, one of the remaining challenges in WMLES simulations is that of external flows with high Reynolds numbers. In such cases, the laminar boundary layer developing near the body leading edge is very thin, and the grid resolution required to resolve such thin region leads to a computational cost that is prohibitively high.

In this work, we propose a modeling strategy for laminar boundary layers in external flows in order to alleviate the resolution requirements needed to resolve this portion of the boundary layer in high Reynolds number flows. As such, the model presented herein seeks to calculate a wall-shear stress for boundary condition imposition, and it is argued that this condition -with the no-transpiration condition at the wall -is sufficient to recover the skin friction and pressure coefficients along the wall. This section describes the reasoning behind the model and presents its formulation, whereas the next section contains an assessment of the feasibility of the model in the aforementioned context.

The basic underlying assumption of our model is that the dynamics of the boundary layer along the wall may be given according to locally self-similar solutions. More specifically, it is assumed that at each point along the wall-tangent direction, there exists a similarity solution from the general Falkner-Skan family satisfying the boundary layer dynamics. Hence, at each point along this direction, we solve for a local similarity solution given by

d 3 f dε 3 + αf d 2 f dε 2 + β 1 - df dε 2 = 0 (7)
with the boundary conditions

ε = 0 : f = 0, (8a) df dε = 0, (8b) ε → ∞ : df dε = 1, (8c)
where α and β are constant parameters, and ε is the similarity variable, written as

ε = y λ(x s )
with λ(x s ) a function to be defined, and x s the distance along the surface.

In the above ordinary differential equation, the function f (ε) defines the velocity profile within the boundary layer through its first derivative,

df dε = u U e , ( 9
)
where u is the velocity within the boundary layer in the x s direction, and U e is the outer flow velocity, that is, the velocity at the edge of the boundary layer [START_REF] Currie | Fundamental Mechanics of Fluids[END_REF], which is assumed to be a function of x s only. In boundary layer theory, this is the velocity that would be obtained for the solution in the limit of infinite Reynolds number, represented by the Euler equations.

The constant parameters α and β in Eq. ( 7) may be written as

λ ν d dx s (U e λ) = α, (10a)
λ 2 ν dU e dx s = β, ( 10b
)
where ν is the kinematic viscosity. Thus, it can be seen that these two parameters define both the outer flow velocity, U e , and the function λ(x s ). Therefore, they completely define the problem to be solved. If we set α = 1, then, the corresponding similarity variable is

ε = y 1 2 -β U e νx s , (11)
and the outer flow velocity is

U e (x s ) = cx β 2-β s , ( 12
)
which may be identified as the velocity of a potential flow over a wedge of half angle πβ/2 (Currie, 2012), as illustrated in Fig. 2.

U ∞ Wedge x s x n πβ 2
Figure 2: Flow over a wedge of half angle βπ/2, where the inviscid velocity along the surface is given by Eq. ( 12).

In the solution process of Eq. ( 7), we seek a value of d 2 f /dε 2 at ε = 0 that satisfies the boundary condition at ε → ∞, Eq. ( 8c). This value is then used to calculate the wall-shear stress, which, for α = 1, may be written as

τ w = µ 1 2 -β U 3 e νx s d 2 f dε 2 ε=0 , (13)
where the definitions in Eqs. ( 9) and ( 11) are used to arrive at this final form. This wall-shear stress is subsequently used to impose the boundary condition at the wall. Following the wall-modeling procedure described in the previous section, the outer flow velocity, U e , appearing in the equation for the wall-shear stress above, is taken from the solution farther away from the wall, at the wall-model height, h wm . It is then left to determine the parameter β defining the problem to be solved in Eq. ( 7), thus obtaining the value of d 2 f /dε 2 at ε = 0.

In order to determine the constant parameter β, we remind that the basic assumption of the model is that the boundary layer dynamics locally satisfies a solution of Eq. [START_REF] References Bassi | High-order accurate discontinuous finite element solution of the 2D Euler equations[END_REF]. As illustrated in Fig. 2, for α = 1, the problem represented in Eq. ( 7) is that of the flow over a wedge. Thus, we assume that at each location along the body surface, the dynamics of the fluid particles are similar to those in the flow over a wedge with a half angle related to the geometric angle between the surface wall-tangent and the streamwise directions. With these assumptions, we define β as

β = κ -|κ| 1 -χ x c , Re c ( 14
)
where κ is defined such that κπ/2 is the geometric angle between the wall-tangent and streamwise directions, as illustrated in Fig. 3, and, hence, may be computed as

κ = 2 π cos -1 t • ŝ , ( 15
)
with t as the unit wall-tangent vector, and ŝ a unit vector in the streamwise direction. In Eq. ( 14), χ is a correction function used to adjust the angle κ obtained solely from geometric considerations. This correction is justified on physical grounds. If β = κ is used directly, the solution obtained from Eq. ( 7) is that of a flow that accelerates as much as the inviscid flow over a wedge of half angle πκ/2. In this case, the velocity gradient within the boundary layer must be higher than that obtained for a viscous flow over the same wedge and, consequently, the computed wall-shear stress is overpredicted. The boundary layer is actually a dissipation of this gradient through viscosity, thus reducing the wall-shear stress when compared to the inviscid flow. Therefore, since d 2 f /dε 2 is directly proportional to the constant parameter defining the problem [START_REF] Schlichting | Boundary-Layer Theory[END_REF], a reduction in the wall-shear stress as computed from Eq. ( 13) is achieved by decreasing the geometric angle κ in the definition of the parameter β.

Another way to justify this angle correction is to see the practical effect of the boundary layer as an obstruction for the flow field, acting to decrease the acceleration along the wedge. This same decreased acceleration may be obtained for the flow in a wedge with a lower half angle. Since the reasoning above applies regardless of whether κ is positive or negative, the adjustment is written in terms of |κ|.

From the rationale above, which accounts for the effects of the viscous boundary layer in the definition of the parameter β, it is intuitive that the correction function in Eq. ( 14) must be related to a Reynolds number and a characteristic length along the surface. Viscous effects become less pronounced the higher the Reynolds number, while the boundary layer grows thicker the further downstream along the surface. Since this work deals with flows around airfoils, we choose the chord Reynolds number, Re c , and the downstream distance relative to the chord, x/c, when defining the correction function, and write it as

χ x c , Re c = 1 - x/c e -ARec 1 + e -ARec ( 16
)
where A is a coefficient, and the chord Reynolds number is given by

Re c = U ∞ c ν ,
with U ∞ the free-stream velocity magnitude. Note that x in the equation above is the distance in the streamwise direction, as opposed to the distance along the surface, x s , used in the boundary layer equations.

The function in Eq. ( 16) has the following desired property,

Re c → ∞, χ → 1,
that is, in the limit of infinite Reynolds number, no adjustment is made to the geometric angle κ, and the solution obtained by solving Eq. ( 7) refers to that of a flow accelerating as much as the inviscid flow over a wedge of half angle κπ/2, as desired. In addition, we have that Re c → -∞, χ →x/c, and, although the above limit is not physically possible, this property means that, for a given chord Reynolds number, as x/c grows larger, and hence the boundary layer thickens, the correction function χ decreases, such that the adjustment in Eq. ( 14) is more pronounced.

The coefficient A is a free constant parameter to fine tune the angle correction. Preliminary parametric studies indicate that a sufficiently good value for this parameter is A = 2.63 × 10 -5 , and this is used throughout this work. Figure 4 shows the behavior of the correction function χ for increasing chord Reynolds number and multiple x/c locations. After the definition of the parameter β for each location along the surface, Eq. ( 7) is solved using the iterative procedure described in [START_REF] Zhang | An iterative method for solving the Falkner-Skan equation[END_REF]. The method employs a free boundary formulation, and uses a shooting method to transform the initial boundary-value problem into a set of initial-value problems, which are solved using standard Runge-Kutta methods. The solution using this iterative method directly outputs the quantity of interest d 2 f /dε 2 at ε = 0.

RESULTS AND DISCUSSION

This section presents the results obtained for the simulations of laminar flows around a NACA 0012 airfoil at zero angle of attack using the wall model described in Sec. 3. Since the model is only intended to be used in laminar flow regions, the test cases that constitute the validation process must exhibit relatively low Reynolds number in order to prevent transition to turbulence. At such low Reynolds numbers, it is difficult to sustain an entirely attached flow against an adverse pressure gradient, and separation of the boundary layer is likely to take place at some position downstream of the leading edge.

The validation of the wall model is done for a flow around a NACA 0012 geometry, even though separation is only absent for very low Reynolds numbers in this case. In order to keep the separation region as mild and as further away from the leading edge as possible, the current wall modeling strategy is validated primarily in a flow with a chord Reynolds number of Re c = 4.5 × 10 3 , which remains laminar in its entirety. The Mach number for every simulation in this work is M ∞ ≈ 0.2, and thus the flow is essentially incompressible. A fully resolved simulation is performed, under the same flow conditions and a sufficiently fine mesh -capable of correctly capturing the boundary layer dynamics -in a wall-resolved configuration, to provide the benchmark solution for results comparisons.

The performance of the model is assessed in terms of the predicted pressure coefficient, C p , and skin friction coefficient, C f , along the surface. A mesh with 96 × 44 × 1 hexahedral elements is used, where the first number indicates the elements in the direction along the surface, whereas the second indicates the elements in the wall-normal direction. The distribution of the elements along the surface is uniform, meaning that ∆x s is constant throughout. It should be noted that, although the flow is two-dimensional, a full, three-dimensional simulation is performed, where only a single element is generated in the spanwise direction. Figure 5 illustrates this mesh.

Figure 5: Two-dimensional section of the mesh employed in this work, for the region near the airfoil.

In addition, the effects of the wall-model height, h wm , the main parameter of the wall-modeling procedure, are also investigated. We assume a linear growth of the exchange location where the outer flow velocity is extracted, starting from x/c = 0 and increasing along the streamwise direction, x. Three different growth rates are tested, as described in Table 1 and shown in Fig. 6, where it should be noted that the wall model height, h wm , is a distance from the wall in the wall-normal direction. The boundary layer edge, extracted from the fully resolved simulation, is shown for comparison. In the assessment of the model, the pressure and skin friction coefficients are plotted for both the upper (suction) and lower (pressure) surfaces of the airfoil, with filled symbols for the former and open symbols for the latter. As expected for a symmetric airfoil at zero angle of attack, the coefficients along the surface are exactly the same for both upper and lower surface and, hence, they are superimposed in the plots below.

The results are compared to those obtained from a fully resolved simulation, using a mesh sufficiently fine to resolve for all the viscous dynamics of the flow, with no wall modeling. In all the plots presenting the solutions, points along the lines are used only for distinguishing the different results, and are not representative of grid points.

The accuracy of the proposed wall model is first investigated using the wall model height with medium growth rate -see Fig. 6. This exchange location is chosen for the investigation because it is the one that more closely resembles the actual edge of the boundary layer. The results obtained for this simulation are shown in Fig. 7. 3 The pressure and skin friction coefficients predicted by the model are very close to the reference solution, both qualitatively and quantitatively, despite the quite restrictive basic assumption of the model.

´ µ ÓÓÑ Ó C f Ò Ðº Ö ÓÒ -C p x/c ÙÐÐÝ Ê ×ÓÐÚ Ä Ñ Ò Ö ÏÅ x/c ÙÐÐÝ Ê ×ÓÐÚ Ä Ñ Ò Ö ÏÅ C f x/c ÙÐÐÝ Ê ×ÓÐÚ Ä Ñ Ò Ö ÏÅ x/c ÙÐÐÝ Ê ×ÓÐÚ Ä Ñ Ò Ö ÏÅ
Near the leading edge, it is seen that the solution obtained with the wall model accelerates slightly more than the reference viscous solution, resulting in a lower pressure in this region. Moreover, the deceleration in the rear portion of the airfoil is also more pronounced, resulting in a higher pressure for the results obtained with the wall model. This is a direct consequence of imposing a lower wall-shear stress in the acceleration region, and a higher wall-shear stress in the deceleration region, when compared to the fully resolved simulation, as is evident from the results for the skin friction coefficient. This should not be seen as a flaw of the model, since the imposed wall-shear stress may be adjusted by changing the parameters of the model in the correction function χ -cf. Eq. ( 14).

The behavior of the skin friction coefficient very near the leading edge, as predicted by the wall model solution, is far from the one obtained for the fully resolved simulation. It may be seen that the skin friction coefficient is very high for x/c ≈ 0, followed by a narrow, negative skin friction region, from which it catches up with the expected behavior, qualitatively. From Eq. ( 13), defining the wall-shear stress to be imposed as the boundary condition, we see that the leading edge, x/c = 0, is a singularity point. Thus, at this location, the wall-shear stress is ill-defined, and the wall model needs special treatment. The effect of imposing an extremely high wall-shear stress at the leading edge, when no special treatment of this location is employed, is that the flow separates, resulting in a very narrow region of reverse flow, from which it recovers shortly after, given the favorable pressure gradient.

For the solutions obtained with the wall-modeling approach, the cumulative skin friction along the surface -that is, the integrated force -is not sufficient to cause boundary layer separation. In the fully resolved simulation, however, separation occurs when the flow is approaching the trailing edge, at approximately x/c ≈ 0.92. In a slightly higher Reynolds number flow, Re c = 5 × 10 3 , separation is also present in the simulations of [START_REF] Swanson | Comparison of NACA 0012 laminar flow solutions: Structured and unstructured grid methods[END_REF], and it takes place at x/c ≈ 0.8. It should be noted, however, that the simulations in [START_REF] Swanson | Comparison of NACA 0012 laminar flow solutions: Structured and unstructured grid methods[END_REF] concern compressible flows, with a Mach number of M ∞ = 0.5.

The solution for the pressure coefficient, C p , exhibits wiggles near the leading edge. These are spurious oscillations that appear due to the fact that the mesh is too coarse to accurately represent the curvature of the airfoil geometry [START_REF] References Bassi | High-order accurate discontinuous finite element solution of the 2D Euler equations[END_REF]. Such spurious oscillations appear even in the solution for the fully resolved simulation, although with remarkably lower amplitude.

From the results shown in Fig. 7, the feasibility of the model is demonstrated. The discrepancies observed may be adjusted by parameter fine-tuning, and the behavior of the skin friction coefficient in the vicinity of the leading edge may be corrected by an appropriate treatment of the computed wall-shear stress at the singularity point, x/c = 0.

A very important aspect of every wall model employed in the context of wall-stress models is the effect of changes in the exchange location, h wm , onto the predicted results. This parameter is usually defined in terms of the boundary layer thickness, δ, or, similarly, the ratio h wm /δ. Since this quantity is not know prior to the simulation, a model that is insensitive to the exchange location is less restrictive in terms of the initial definition of h wm .

Figure 8 shows the results obtained for the simulations using the three previously described wall model heights, with low, medium and high growth rates along the streamwise direction. These simulations are all performed in the same mesh, depicted in Fig. 5, such that any influence regarding mesh refinement or element distribution is eliminated. 3 Similar to the results shown in Fig. 7 for the medium growth rate of the wall model height, all the results shown in Fig. 8 are in very good agreement with the reference solution, especially from a qualitative point of view. When the solutions for the distinct wall model heights are compared to each other, it may be seen that only the results obtained with the smallest growth rate is distinguishable from the others. Moreover, differences are only seen for the pressure coefficient. We note that the wall model height for the lower growth rate is too close to the wall and, as a consequence, the probed velocity, U e , is not representative of the outer flow velocity. For both medium and high growth rates, on the other hand, the results are essentially the same, since the wall model height in these cases are farther from the wall, and more closely matches the boundary layer edge.

´ µ ÓÓÑ Ó C f Ò Ðº Ö ÓÒ -C p x/c ÙÐÐÝ Ê ×ÓÐÚ ÄÓÛ Ä Ò Ö Å Ä Ò Ö À Ä Ò Ö x/c ÙÐÐÝ Ê ×ÓÐÚ ÄÓÛ Ä Ò Ö Å Ä Ò Ö À Ä Ò Ö C f x/c ÙÐÐÝ Ê ×ÓÐÚ ÄÓÛ Ä Ò Ö Å Ä Ò Ö À Ä Ò Ö x/c ÙÐÐÝ Ê ×ÓÐÚ ÄÓÛ Ä Ò Ö Å Ä Ò Ö À Ä Ò Ö
Since the results are very close to each other, the same behavior discussed in the previous section is observed here. The skin friction coefficient is underpredicted in the region near the leading edge, whilst being overpredicted at the rear portion of the airfoil. Consequently, the pressure is slightly lower near the leading edge due to an excess acceleration for the results using the wall modeling procedure. None of the solutions exhibit boundary layer separation near the trailing edge, however, as is expected from the behavior in the reference simulation.

Spurious oscillations appear in the solutions of the pressure coefficient for every distinct wall model height tested, as briefly discussed for the results shown in Fig. 7. This is another evidence that such wiggles are a consequence of insufficient discrete representation of the airfoil curvature, and is not related to the proposed wall model.

Another aspect of the results is that the wall-model height, h wm , have very little effect on the predicted pressure coefficient. This is to be expected, since the pressure distribution at the wall is imposed from the outer flow solution and, therefore, viscous effects have little influence for fully attached flows. Nonetheless, the same behavior is observed for the predicted skin friction coefficient, where the influence of the wall model height is mild. As mentioned previously, low sensitiveness to the exchange location is an important feature of the wall model, where restrictions on the correct definition of an initial wall model height, h wm , are relaxed.

CLOSING REMARKS

This work proposes a wall-modeling approach for simulations of external laminar boundary layers, based on local selfsimilar solutions. The wall model uses information from the numerical solution farther away from the wall, combined with local self-similar solutions of the Falkner-Skan family in order to provide a wall-shear stress that is subsequently imposed as a boundary condition. The proposed wall model is validated on simulations of an incompressible, low Reynolds number flow, with Re c = 4.5 × 10 3 and M ∞ ≈ 0.2, where the boundary layer flow remains laminar in its entirety.

The results show the feasibility of the proposed model for incompressible, low Reynolds number flows. The solutions obtained with the wall model exhibits excellent behavior for both the pressure and skin friction coefficients, following very closely the solutions obtained in a fully resolved simulation without wall modeling.

The wall model as proposed in Sec. 3., with the constant parameter β as defined by Eqs. ( 14)-( 16), computes a wallshear stress that is too low in the acceleration region, and too high in the deceleration region. Consequently, the predicted pressure distribution along the surface is too low in the acceleration region, and too high in the deceleration region. One way to adjust the computed wall-shear stress to match that obtained in a fully resolved solution is by a corresponding adjustment of the correction function, χ, through parameter fine tuning.

An important feature of the laminar wall model is that the solutions for the pressure and skin friction coefficients are quite insensitive to the wall model height, h wm . Hence, the limitations on the initial definition of the exchange location are relaxed, since locations that are not strictly matching the boundary layer edge provide reasonably accurate results.

Although spurious oscillations in the pressure distribution along the surface may be seen for every result presented in this work, we claim that such wiggles are a consequence of the mesh coarseness and its failure to accurately represent the curved geometry of the airfoil in the acceleration region.

The laminar wall model, as presented in this work, may be readily combined with a turbulent wall model in the context of wall-stress models in WMLES solutions, in order to perform a simulation of the complete developing boundary layer in external flows with reasonable computational costs for all regions and regimes of the boundary layer.

Figure 1 :

 1 Figure 1: Schematic representation of the wall-modeling procedure employed in the context of WMLES, also used in this work.

Figure 3 :

 3 Figure 3: Definition of the geometry-related parameter κ, illustrated for a NACA 0012 airfoil surface.

Figure 4 :

 4 Figure 4: Corretion function χ as a function of the chord Reynolds number, Re c , for multiple dimensionless streamwise locations, x/c.

Figure 6 :

 6 Figure6: Wall model heights with linear growth along the surface, for different growth rates. The boundary layer edge (BLE) obtained from a fully resolved simulation is also shown.

Figure 7 :

 7 Figure 7: Solutions of pressure and skin friction coefficients for the proposed laminar wall model, using the wall model height with medium growth rate. Filled symbols: suction side; open symbols: pressure side. Symbols are superimposed.

Figure 8 :

 8 Figure 8: Solutions of pressure and skin friction coefficients for the proposed laminar wall model, using the wall model height with medium growth rate. Filled symbols: suction side; open symbols: pressure side. Symbols are superimposed.

Table 1 :

 1 Shape and growth rate of the wall model height, h wm , tested in the simulations.

	Growth type Growth rate	h wm (x)/c
		Low	0.04 x/c + 0.01
	Linear	Medium	0.08 x/c + 0.01
		High	0.11 x/c + 0.01

ACKNOWLEDGEMENTS

The authors acknowledge the support for the present research provided by Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq, under the Research Grant No. 309985/2013-7. The work is also supported by the computational resources from the Center for Mathematical Sciences Applied to Industry, CeMEAI, funded by Fundação de Amparo à Pesquisa do Estado de São Paulo, FAPESP, under the Research Grant No. 2013/07375-0. The authors further acknowledge the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer. This work was also granted access to the HPC resources of IDRIS under the allocation 2021-A0112A12067 made by GENCI. The first author acknowledges the financial support in the form of a doctoral scholarship from FAPESP, under the Grant No. 2018/05524-1. The third author acknowledges authorization by his employer, Embraer S.A., which has allowed his participation in the present research effort. Additional support to the fourth author under the FAPESP Research Grant No. 2013/07375-0 is also gratefully acknowledged.


domain: math phys

text:
Raphaël Neymann

Hendrik Meier

Hugo Lhachemi

Christophe Prieur

Antoine Girard

Minimization of Parameter Sensitivity to Pre-Estimation Errors and its Application to the Calibration of Magnetometer Arrays

We consider the problem of parameters estimation in the situation where a subset of parameters of the underlying model has already been estimated. Potential errors in the pre-estimated parameters limit the accuracy in the computation of the still unknown parameters in a way that crucially depends on the symmetries of the data set. In this article, we develop an optimization approach that minimizes the sensitivity to errors in the determination of parameters by determining optimal weights that depend on both the model and the input data in the situation where the data are incomplete. We apply the method to the problem of calibrating the sensor positions in an array of magnetometers whose scale factors and biases have been estimated beforehand.

I. INTRODUCTION

Sensor calibration is an essential step to get acceptable performance from any navigation devices. Consistent measurements are constantly required in navigation algorithms, in particular in practical applications of magnetoinertial dead-reckoning, which consists in a fusion of inertial and magnetic measurements.

Fundamentally, a calibration is an optimization problem, for which the inputs are a dynamic data set and, often, additionally pre-estimated parameters, e.g., physical constants or model parameters known or estimated beforehand. Inaccurate knowledge of these parameters can have a strong impact on the final calibration as the optimization problem can be very sensitive to their errors.

In the context of magnetic sensors, many calibration methods make use of the invariant associated with a homogeneous field, either the Earth's magnetic field [START_REF] Kok | Calibration of a magnetometer in combination with inertial sensors[END_REF], [START_REF] V. Renaudin | Complete Triaxis Magnetometer Calibration in the Magnetic Domain[END_REF], [START_REF] Li | Calibration of tri-axial magnetometer in magnetic compass using vector observations[END_REF], [START_REF] Yang | Complete tri-axis magnetometer calibration with a gyro auxiliary[END_REF] or coil-generated constant field [START_REF] Díaz-Michelena | Calibration of qm-moura three-axis magnetometer and gradiometer[END_REF], [START_REF] Mohamadabadi | An automated indoor scalar calibration method for three-axis vector magnetometers[END_REF], [START_REF] Pang | Calibration of a fluxgate magnetometer array and its application in magnetic object localization[END_REF]. Calibrations in inhomogeneous fields have also been developed [START_REF] Adachi | Calibration for a multichannel magnetic sensor array of a magnetospinography system[END_REF], [START_REF] Chesneau | Calibration of a magnetometer array using motion capture equipment[END_REF], [START_REF] Turner | Calibration of a novel three-axis fluxgate gradiometer for space applications[END_REF], [START_REF] Yoshida | Calibration system for a multichannel squid magnetometer[END_REF]. It has recently [START_REF] Chesneau | Calibration of a magnetometer array using motion capture equipment[END_REF] been proven that in inhomogeneous fields, positions of magnetometers become identifiable using a motion capture reference. The knowledge of this geometry is needed for the computation of the magnetic field gradient, which is the key ingredient within magneto-inertial dead-reckoning [START_REF] Chesneau | Motion estimation of a rigid body with an EKF using magneto-inertial measurements[END_REF], [START_REF] Chesneau | Improving magneto-inertial attitude and position estimation by means of a magnetic heading observer[END_REF], [START_REF] Vissière | Using distributed magnetometers to increase imu-based velocity estimation into perturbed area[END_REF]. More recent developments include the combination of magneto-inertial algorithms with machine learning [START_REF] Zmitri | Magnetic Field Gradient-Based EKF for Velocity Estimation in Indoor Navigation[END_REF], [START_REF] Zmitri | BiLSTM Network-Based Extended Kalman Filter for Magnetic Field Gradient Aided Indoor Navigation[END_REF] and magnetic mapping [START_REF] Huang | A Tightly-Integrated Magnetic-Field aided Inertial Navigation System[END_REF].

We are thus interested in calibrating the positions of the single-axis sensors in an array of magnetometers as accurately as possible. The full calibration model being highly non-linear [START_REF] Chesneau | Calibration of a magnetometer array using motion capture equipment[END_REF], typical calibration methods rely either on dynamic systems [START_REF] Schoukens | Nonlinear System Identification[END_REF] or linearization and iteration, using Least-Squares methods (LS), especially Ordinary Least-Squares (OLS) that assign the same weight to each sample. However, in the latter case, iterating on inaccurate input parameters can lead to amplified errors in the final results. Moreover, trajectories used to generate calibration data often require particular symmetries in their design (needed for reducing the sensitivity in the pre-estimation errors), which may be yet hard to properly define and execute, especially if the sensors are embedded in a larger device.

In this paper, we develop a calibration algorithm that is resilient with respect to errors in the pre-estimated parameters and unexpected perturbations in the experiment by assigning suitable weights to the input samples. Our approach is inspired by the computation of a fictional material in magnetostatics for which the sensitivity in the magnetic permeability and electrical current is minimal [START_REF] Dyck | Automated design of magnetic devices by optimizing material distribution[END_REF], and aims to compute a sample-dependent weight function such that the estimation of a calibration parameter of interest be as insensitive as possible to errors in a pre-estimated set of parameters. A sensitivity minimization strategy, based on data selection, has recently been implemented in the context of Li-ion batteries calibration [START_REF] Fogelquist | Uncertainty-Aware Data Selection Framework for Parameter Estimation with Application to Li-ion Battery[END_REF]. A complementary approach, rebuilding missing points in incomplete data set with polynomials, has also been suggested [START_REF] Shin | Nonadaptive Quasi-Optimal Points Selection for Least Squares Linear Regression[END_REF]. Our method is yet more general, for it consists in a weighting strategy by assigning a quality coefficient to each sample of the data set.

Thus, we introduce a sensitivity measure of the estimated parameters with respect to pre-estimation errors. This sensivity can then be minimized by selecting an optimal weight function of the data set.

To demonstrate the utility, in Sec. II, we first study within a simple toy problem how selecting suitable subsets of the calibration data can significantly reduce the sensitivity in a simple situation. We then establish a general method in Sec. III applicable to situations in which exploitable symmetries are not evident and finally apply this method in Sec. IV to our initial problem of estimating the sensor positions in an array of magnetometers.

II. A TOY PROBLEM STUDY

In this section, we motivate our approach by showing on a simple toy problem that selecting or weighting the calibration can help to reduce the sensitivity to errors in pre-estimated parameters.

A. Model for a 2D magnetometer

Consider a single-axis sensor in a two-dimensional world measuring the projection of an ambiant tangent field, e.g. a magnetic field, onto its axis. The sensor is attached to a (moving) body, to which a (moving) reference frame b has already been assigned. If the sensor response to the ambient field B is linear, a complete model for the measurement M is given by

M = c • B b + b. (1)
Herein, b ∈ R denotes the sensor bias, and the (row) vector c = c x c y ∈ R 1×2 points along the sensitive axis (in body frame b). Its norm is the scale factor. The upper index b in B b indicates that vector components are in body frame.

Our objective is to have at our disposal the best possible estimates for the calibration parameters b and c. Suppose that the bias b has already been estimated to a value b 0 by a separate calibration method, e.g., in a zero-field chamber, with a residual error δb = b 0b * , where b * denotes the true but unknown bias value. In this situation, we wish to estimate the remaining calibration parameter c while being as insensitive as possible to the potential error δb.

B. A simple calibration method

For the estimation of the scale factor vector c, suppose the following experiment has been carried out: the rigid body carrying the sensor is rotated by a known angle θ with respect to the (inertial) navigation frame n. The (uniform) ambient field is oriented along the y-axis of this frame, B n = 0 B . In this situation, the (by hypothesis noise-free) measurement for a given θ amounts to

M (θ) = c • R (θ) B n + b with R (θ) = cos θ -sin θ sin θ cos θ , (2)
cf. the models in [START_REF] Chesneau | Calibration of a magnetometer array using motion capture equipment[END_REF] and [START_REF] V. Renaudin | Complete Triaxis Magnetometer Calibration in the Magnetic Domain[END_REF] and Fig. 1.

If measurements M (θ) are available for homogeneously distributed angles θ ∈ (-π, π), a standard LS optimization of c with b fixed to the possibly erroneous pre-estimate b 0 yields a value that is independent from the actual bias b, cf. Eq. (28) in App. V-A. This value is thus fully insensitive to any residual error δb, i.e. equal to the true value. This independence is a manifestation of the O(2) rotational symmetry in the chosen samples.

If however, for some reason, measurements are available only for θ ∈ -π 4 , 3π 4 , see Fig. 1, the aforementioned symmetry is broken and the standard LS optimization yields an error

δc = δb B 2 2 π 1 1 , (3)
using Eq. (30) of App. V-A. This error reflects how within the standard approach, the scale factor estimation is quite sensitive to residual errors in the pre-estimated bias b. Eqs. ( 4) and ( 5). (For c y , the regions of high and low sensitivities are simply reversed.)

C. Sensitivity in data subsets

For the trajectory θ ∈ -π 4 , 3π 4 , the result obtained by standard LS calibration contains an error of order δb/B . Let us now consider the subset of the initial trajectory given by θ ∈ (-ε, ε) for 0 ≤ ε ≤ π/4. Applying the LS optimization to this interval only, we obtain, using again Eg. (30), the scale factor vector up to an error

δc = δb B 0 ι(ε) , (4)
with ι(ε) = (2 sin ε)/(ε + sin ε cos ε) taking values between 1 and ≈ 1.1. Whereas the error δc y remains of order δb/B , the estimation of c x is remarkably insensitive to the error δb in the pre-estimated bias, see Fig. 1. This insensitivity can be attributed to the reflectional symmetry (w.r.t. θ = 0) of the subset of calibration data.

Analogously, for θ close to π/2, i.e. θ ∈ π 2 -ε, π 2 + ε , we find the error estimate

δc = δb B ι(ε) 0 . (5)
Here, c x , with an error of order δb/B , is very sensitive to δb, see Fig. 1, as the reflectional symmetry (w.r.t. θ = π/2) favors c y , whose error vanishes.

As a result, we find that the LS sensitivity of the estimation of each calibration parameter to the error in the pre-estimated parameters depends crucially on the underlying input data. In particular, symmetries within a well-chosen subset of input data may render a certain parameter insensitive to such errors.

Following the small toy model study of this section, we develop in Sec. III a general approach to minimize the sensitivity of parameter estimation to errors in preestimated parameters in situations in which exploitable symmetries are not apparent.

III. GENERAL APPROACH

In this section, we propose a general approach for an optimization problem estimating (calibration) parameters with minimal sensitivity to errors in predetermined parameters.

A. Setting of the problem

We consider a sensor calibration problem defined by its involved quantities and a calibration model of the form defined below.

Definition 1: Physical quantities. A calibration problem involves

• the calibration parameter vector π π π ∈ R n , n ∈ N * , to be estimated; • a vector of pre-estimated calibration parameters m 0 ∈ R m , m ∈ N * , needed for the estimation of π π π and determined with known precision caracterized by a given covariance matrix;

• a set of data points k ∈ {1, . . . , N } → ξ ξ ξ k ∈ M , where M
is the manifold in which the data points are recorded (i.e. data of the sensor to be calibrated and possibly other sensors used in the calibration experiment).

Definition 2: Calibration model. A calibration model is a real-valued (C 3 -)function ψ : R n × R m × M -→ R (π π π, m,ξ ξ ξ) -→ ψ(π π π; m,ξ ξ ξ). ( 6
)
such that for the true value m * of the pre-estimated parameters, there is exactly one parameter vector π π π * such that ψ(π π π * ; m * ,ξ ξ ξ) = 0 for all ξ ξ ξ ∈ M .

B. Calibration formalism

For fixed pre-estimated parameters m 0 and measurement data points ξ ξ ξ k , k ∈ {1, . . . , N }, we seek to estimate π π π by minimizing loss-function

J m 0 ,ξ ξ ξ, f : R n -→ R + π π π -→ Ψ(π π π, m 0 ,ξ ξ ξ, f ), (7)
with

Ψ(π π π, m 0 ,ξ ξ ξ, f ) = N k=1 f k ψ π π π; m 0 ,ξ ξ ξ k 2 , (8)
Herein, the coefficents f = ( f k ) k=1,...,N denote weights associated to each data point. They are defined by Definition 3:

Weight function. A weight function k ∈ {1, . . . , N } -→ f k ∈ R + (9)
is a nonnegative function that satisfies the L 1 -property,

k f k = 1. 1
In the following, we develop the formalism that allows us to determine an optimal weight function f in order for the minimum of J m 0 ,ξ ξ ξ, f , which we denote π π π| m 0 ,ξ ξ ξ, f , to have minimal sensitivity to errors in m 0 .

We assume that π π π| m 0 ,ξ ξ ξ, f is (at least locally) a unique minimum and that the Hessian matrix of Ψ(π π π, m 0 ,ξ ξ ξ, f ) w.r.t. π π π is positive definite at π π π| m 0 ,ξ ξ ξ, f for all weight functions f involved. 2 1 In particular, Ordinary Least-Squares (OLS) correspond to f = 1/N . 2 In practice, this can usually be accomplished with a sufficiently rich calibration trajectory.

C. Parameter sensitivity

We define the sensitivity of the unique minimum π π π| m 0 ,ξ ξ ξ, f as the linear response to the error m 0m * with respect to the true value m * of the pre-estimated parameter.

Definition 4: Parameter sensitivity. The sensitivity of the estimated parameter π π π| m,ξ ξ ξ, f to errors in the preestimated parameters m, for given data ξ ξ ξ and weight function f , is defined as the derivative (Jacobian)

S m * ,ξ ξ ξ, f ( π π π) = ∇ m π π π| ξ ξ ξ, f (m * ) (10)
such that the error in π π π| m 0 ,ξ ξ ξ, f is given by

S m * ,ξ ξ ξ, f ( π π π) • (m 0 -m * ) + O m 0 -m * 2 .
(11) With full knowledge of the system, the sensitivity can be readily computed from the loss function Ψ:

Proposition 1: The sensitivity S m * ,ξ ξ ξ, f ( π π π) is given by S m * ,ξ ξ ξ, f ( π π π) = -H -1 π π π,π π π • H π π π,m , (12)
with the Hessian matrices

H π π π,π π π = ∇ 2 π π π,π π π Ψ π π π| ξ ξ ξ, f (m * ), m * ,ξ ξ ξ, f ∈ R n×n , H π π π,m = ∇ 2 π π π,m Ψ π π π| ξ ξ ξ, f (m * ), m * ,ξ ξ ξ, f ∈ R n×m . Proof:
The formula follows immediately from the implicit function theorem applied to the mapping

(π π π, m) → ∇ π π π Ψ(π π π, m,ξ ξ ξ, f ) in the neighborhood of π π π| m * ,ξ ξ ξ, f , m * .
The Hessian H π π π,π π π is invertible by assumption, see Sec. III-B.

As the exact value m * is not known in our calibration, we may use the available (pre-)estimation m 0 to obtain the approximation

S m * ,ξ ξ ξ, f ( π π π) S m 0 ,ξ ξ ξ, f ( π π π), (13)
which, if inserted into Eq. ( 11), yields an equal error estimate of the same order. S m 0 ,ξ ξ ξ, f ( π π π) can be calculated from the Hessian matrices at m 0 .

The i -th row of S m 0 ,ξ ξ ξ, f ( π π π) defines the sensititivity of the i -th component of π π π to the error in m 0 . We denote by e i the i -th (row) vector of the standard basis of R n and thus obtain

s (i ) m 0 ,ξ ξ ξ, f ( π π π) = e i • S m 0 ,ξ ξ ξ, f ( π π π) ( 14
)
as the sensitivity of the i -th component of π π π. Definition 5: Scalar sensitivity. The scalar sensitivity of the i -th component of π π π| m,ξ ξ ξ, f to errors in m for data ξ ξ ξ and weight function f is defined as

Λ (i ) f = s (i ) m 0 ,ξ ξ ξ, f (i ) ( π π π) 2 (15)
We assume that all components of m have the same physical dimensions, which can always be achieved if necessary by applying suitable scale factors. We furthermore assume that the covariance matrix, which caracterizes precision of the estimate m 0 , is proportional to the unit matrix. This also is always achievable through a suitable orthogonal transformation and suitable scale factors.

Using Def. 5, we formulate the calibration problem as the following nonlinear variational problem: Find a weight function f (i ) that minimizes Λ (i ) ( f ).

With an optimal weight function f (i ) at hand, the i -th component of π π π can be found by a LS optimization of π π π. Depending on the situation, it may be suitable to treat several components at once by minimizing the sum of the corresponding (possibly weighted) functions Λ (i ) ( f ) with respect to a common f .

D. Steps of the calibration algorithm

We summarize the steps the general method of calibration with minimal sensitivity to pre-estimated parameters.

1) Formulate the calibration problem as in the form of

Sec. III-A. 2) Calculate the sensitivity S m * ,ξ ξ ξ, f ( π π π) using Proposition 1.

3) For the selected component (or components) i of the parameter vector π π π, compute the optimal weight function f by minimizing the sensitivity Λ (i ) ( f ) of Eq. ( 15). 4) Compute π π π| ξ ξ ξ, f (m 0 ) using the derivation by LS optimization using the optimal weight function f .

For numerical optimizations discussed in this Article, we have used the Levenberg-Marquardt algorithm [11][13].

E. The toy problem revisited

To demonstrate the general approach, let us apply it to the toy-model of Sec. II. The quantities involved are π π π = c (scale factor vector) 3 , m = b (pre-estimated bias), and ξ ξ ξ k = (M (θ k ), θ k ) with M k = M (θ k ) denoting the measurement at orientation angle θ k . The calibration model is given by the function

ψ (c, b, (M , θ)) = M -b -B c x sin θ -B c y cos θ, ( 16
)
with B denoting the external field strength, cf. Eq. ( 2). Calculating the Hessian matrices required for Proposition 1, we find for the sensitivity

∂ ∂b c| (M ,R), f (b 0 ) = -H -1 c,c • H b,c (17)
with

H c,c = - B 2 2 N k=1 f k 1 -cos 2θ k sin 2θ k sin 2θ k 1 + cos 2θ k , ( 18
)
H c,b = B N k=1 f k sin θ k cos θ k . ( 19
)
For a given weight function, the scale factor estimate is

c| (M ,R), f (b 0 ) = -H -1 c,c • B N k=1 f k (M k -b 0 ) sin θ k cos θ k , (20)
which we note to be linear in b 0 . The error estimate δc = δb • (∂ c/∂b) is thus equal to the exact error δc.

In the simulated calibration problem under consideration, a set of N = 400 data points is available for angles θ 3 In this section, we consider c to be a column vector. Applying the approach of the previous section, we optimize the weight function f k for c x and c y separately. For c x , 4 we thus compute the weight function f (x)

k that minimizes Λ (x) = 1 0 • H -1 c,c • H b,c 2 .
For the numeric optimization, we employ a Levenberg-Marquard (LM) algorithm, starting with homogeneous weights f k = 1/N for which the scalar sensitivity Λ (x) 1/2 is of order 1 G -1 , see Table I. As result of the LM optimization, we obtain the non-trivial weight function f (x) k depicted in Fig. 2 and for which the sensitivity practically vanishes. Evaluating c x using Eq. ( 20) with the optimal LM weights, we thus find the true value with practically vanishing error.

We further note that there is more than one optimal weight function. For instance, the much simpler weight function given by f (x) k = 2/N for θ k ∈ (-π/4, π/4) and f (x) k = 0 otherwise, see Fig. 2, yields (practically) zero sensitivity and error as well (column "Symmetry" in Table I). It is this weight function that, chosen by symmetry considerations, has led to Eq. ( 4) in Sec. II-C.

For a good choice of a proper weight function in real situations, noise can be a factor that may favor some weight functions over others that are equally optimal in the noise-free setting.

IV. GEOMETRY OF AN ARRAY OF MAGNETIC SENSORS

The estimation of the magnetic gradient is an important ingredient of magneto-inertial navigation [START_REF] Vissière | Using distributed magnetometers to increase imu-based velocity estimation into perturbed area[END_REF] and, if realized in an array of magnetometers, requires accurate knowledge of the array geometry, in particular the positions of the single sensors.

As recently shown by Chesneau et al. [START_REF] Chesneau | Calibration of a magnetometer array using motion capture equipment[END_REF], it is possible to calibrate simulationeously for one single-axis magnetometer its scale factor, bias, and position relative to a body frame origin for which motion capture data are available. Being able to calibrate a single sensor in this way, we are a fortiori in a position to calibrate an array of such sensors. The calibration requires a trajectory providing as many independent observations than the number of parameters to estimate, and encompassing a region of space with an inhomogeneous magnetic field. In practice, it may become a difficult task to produce a trajectory that excites all the degrees of freedom necessary to identify all involved calibration parameters at once.

Here, we suppose that sensor scale factors and biases are known beforehand. This may be the case if a separate calibration, e.g. in a homogeneous magnetic field [START_REF] V. Renaudin | Complete Triaxis Magnetometer Calibration in the Magnetic Domain[END_REF], has already been carried out or if an iterative calibration algorithm is used alternating scale factor and bias optimization on one hand and position optimization on the other hand. Furthermore, we assume that a local fit of the ambient magnetic field to a suitable model is available.

A. Formulation of the calibration model

We adopt the calibration model of [START_REF] Chesneau | Calibration of a magnetometer array using motion capture equipment[END_REF] for the single-axis magnetometer

M -b -c • R Bn X + Rp = 0. (21)
Herein, p ∈ R 3 is the sensor position in body frame relative to the frame origin, which is the parameter to be estimated: π π π = p. The sensor bias b ∈ R and scale factor (row) vector c ∈ R 1×3 are the pre-estimated parameter: m = b c . Finally, we dispose for each data point of the sensor measurement M ∈ R as well as the position X ∈ R 3 of the body origin in the navigation frame n and the attitude R ∈ SO(3) of the rigid body relative to n. These quantities together form the variable ξ ξ ξ = (M , X, R) ∈ M = R×R 3 ×SO(3). Equation ( 21) with these definitions sets up a calibration problem in the sense of Sec. III-A.

As mentioned before, we suppose the magnetic field to be known in form of a map X → Bn (X). Assuming positions p close to the body origin on the scale defined the (non-linear) variations in the magnetic field, we obtain the calibration model ψ :

R 3 ×R 4 × R × R 3 × SO(3) → R with ψ p, (b, c) , (M , X, R) = M -b -c • R Bn (X) + ∇ Bn (X) Rp (22
) by linearization of Eq. [START_REF] Yoshida | Calibration system for a multichannel squid magnetometer[END_REF].

In order to obtain a uniform physical dimension in the tuple m = b c , we may choose to measure the bias in units of the average magnetic field norm B . 5

B. Computation of sensitivity

Using the method presented in Sec. III, our objective is to estimate the sensor position, obtained by LS minimization as

p (M ,X,R), f (b 0 , c 0 ) = -H -1 p,p • N k=1 f k M k -b 0 -c 0 • R k Bn (X k ) R k ∇ Bn (X k ) R k • c 0 ( 23
)
for an optimal weight function f k . 6Separating the sensitivities errors in bias and errors in the scale factor,

∂ ∂b p (M ,X,R), f (b 0 , c 0 ) = -H -1 p,p • H p,b (24)
∇ c p (M ,X,R), f (b 0 , c 0 ) = -H -1 p,p • H p,c , (25)
we set up three cost functions, for each of the components of p = (p x , p y , p z ) in the form

Λ (i ) ( f ) = σ 2 b ∂ p i ∂b (M ,X,R), f 2 + σ 2 c ∇ c p i (M ,X,R), f 2 (26)
for i = x, y, z. The coefficients σ 2 b and σ 2 c , which designate the (known) variances of the pre-estimated parameters, have been chosen in accordance with the remarks following Def. 5.

C. Settings of the calibration

We apply the calibration method to a single-axis magnetometer for which estimates of bias and scale factor are available to a precision characterized by standard deviations σ b ≈ 1.4 mG and σ c ≈ 1.4 • 10 -3 . We confront the method to both simulated and real data, minimizing the (scalar) sensitivity separately for each of the three components of p.

1) Simulation: We simulate the ambient magnetic field in navigation frame as modelled by its value at a chosen origin, B n 0 = (0.46, 0.0, 0.0) G, and a uniform gradient ∇B n 0 whose (Frobenius) norm is of order 45 mG/m. The model for the magnetic field is supposed to be perfectly known in the calibration algorithm. The calibration trajectory consists of 35 samples corresponding to non-uniformly chosen positions and attitudes of the body frame (to which the magnetometer is attached). Positions and attitudes are also supposed to be known without error in the calibration algorithm.

For bias and scale factor, we choose the true values b * ≈ -3.0 mG and c * = 1.0 0.0 0.0 . Their pre-estimates b 0 and c 0 are simulated by adding errors randomly chosen according to a Gaussian law with standard deviations as specified above.

2) Real data: Real data is obtained from a calibration trajectory confined to an indoor region allowing for position and attitude tracking by an optical motion capture. For this considered volume, the magnetic model has been determined beforehand. Harmonizing the data of motion capture with an on-board inertial measurement unit in the body reference frame, we are in the position to define the body frame origin as the position of the accelerometer and the body frame orientation by the accelerometer axes. The calibration parameters to estimate are then the position coordinates of the magnetometer in this frame. A reference position (at millimeter precision) is available by conventional metric measurements.

Errors in the bias b and the scale factors c considered in this study dominate over the errors due to noise or systematic errors in the motion capture data and the magnetic field model. Focussing on the sensitivity to b and c, a detailled characterization of the errors in position, attitude, and model parameters is not necessary. (Such errors may yet be included in the sensitivity measure in a separate study.)

D. Results

1) Simulations:

Table II shows the results of the simulation study. The OLS method (i.e. a uniform weight function) yields a position p with an error of 15 mm, in accordance with the sensitivity to bias and scale factor errors of also 15 mm (in units of their respective standard deviations). Minimizing this sensitivity to both bias and scale factor errors, it is reduced to below 0.15 mm, which is reflected by a position error of less than 0.05 mm, a precision that is difficult to achieve using conventional means.

2) Real data: The results for real-life data are shown in Table III. The position error using the conventional OLS method amounts to 43 mm, which is in magnitude twice the OLS error obtained in simulations. We may attribute this increased error to noise and systematic errors that have not been included in the simulation model. Applying the weights obtained from our sensitivity minimization algorithms, we obtain a position that is accurate up to 5 mm in comparison with the (nominal) reference position. Minimization of sensititivity has thus been reduced the error by a factor of order 10. The residual position error of 5 mm can again be attributed to the errors in the motion capture data and the magnetic model.

E. Discussion

Our simulation shows that for sufficiently rich calibration trajectories, the sensitivity of the position calibration to errors in the pre-estimated parameters can be reduced to arbitrarily small values. Applying the method to data obtained in real experiments, we have achieved a net error reduction (by a factor of order 10), yet a residual position error of several millimeters remains. At this precision, our approach is reliable in the situation where conventional position measurements (including the simple use of a ruler) are unavailable, e.g. if the magneto-inertial measurement unit is hidden in a "black box". In this case, we provide a non-invasive method to estimate the magnetometer position with an accuracy that is optimal with respect to errors in the pre-estimated calibration parameters.

For a further reduction of the position error, we may on one hand try and improve the performance of the reference data sets used in the calibration algorithm. These are the trajectory data set, which depends on the performance of the motion capture system, and the magnetic field model, which may be developed to a higher order and, in our case, would also benefit from a more performing motion capture system, we have used to map the local ambient field.

In a complementary approach, errors of the motion data and the magnetic field model may, too, be included into our sensitivity measure. In doing so, our approach may be extended beyond static errors (such those in bias and scale factor) to include also noise.

V. CONCLUSIONS

We have developed an optimized least-squares approach that minimizes sensitivities to errors in preestimated parameters, leading sometimes to complete vanishing. Such situations arise in calibration problems in which a certain subset of parameters is already known (with limited accuracy) at a given step of the procedure.

We have tested our sensitivity-based approach in both simulation and experiment in the context of position calibration of magnetometers in an inhomogeneous magnetic field. In simulation, we have proven the validity of the approach. Our experimental result shows a significant error reduction. We have discussed how the residual error may further be reduced in future studies.

Our approach may finally be helpful in the design of a calibration process such as the conception of a calibration trajectory with minimal sensitivity to errors.

4 c x very sensitive to bias error c x insensitive to bias errorFig. 1 :

 41 Fig. 1: In red the angles covered by the calibration experiment. The arrows indicate the regions within this trajectory where the scale factor c x is very sensitive or insensitive to an error in the pre-estimated bias b, cf. Eqs. (4) and (5). (For c y , the regions of high and low sensitivities are simply reversed.)

Fig. 2 :

 2 Fig.2: OLS (red, solid) and optimized weights LS for the calibration of c x in the toymodel. The optimal functions depicted are the one suggested by symmetry considerations (green, dashed) and the result of the LM optimization (blue, dotted).

TABLE I :

 I Numerical results for errors and sensitivities for the calibration of c i , i = x, y, for the homogeneous weights (OLS) and optimized weights LS obtained by symmetry considerations or by Levenberg-Marquardt (LM) optimization.

		OLS	Optimized weights
			Symmetry	LM
	Error δc i (-)	1.8 • 10 -4	< 10 -9	< 2 • 10 -13

TABLE II :

 II Simulation results for the position error and sensitivity to the errors with respect to pre-estimated bias and scale factor values. The indicated scalar sensibility is the square root of Λ = i =x,y,z Λ (i ) f (i ) with Λ (i ) as in Eq. (15) and f (i ) the weight function minimizing Λ (i ) .

		OLS	Optimized weights
	Error δp (mm)	15	0.05
	Scalar sensitivity Λ 1/2 (mm)	15	0.15

TABLE III :

 III Experimental results for the position error and sensitivity to errors in pre-estimated bias and scale factor values. The indicated position errors is obtained by comparison with the conventionally measured magnetometer position (relative to the on-board accelerometer).

		OLS	Optimized weights
	Error δp (mm)	43	4.6
	Scalar sensitivity Λ 1/2 (mm)	20	0.1

Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des Signaux et Systèmes, 91190 Gif-sur-Yvette,

France 2 SYSNAV, 27200 Vernon

, France[START_REF] Chesneau | Improving magneto-inertial attitude and position estimation by means of a magnetic heading observer[END_REF] Université Grenoble Alpes, CNRS, Grenoble-INP, GIPSA-lab, 38000 Grenoble, France. The work of the second author has been partially supported by MIAI@Grenoble Alpes (ANR-19-P3IA-0003).

For c y , a completely analogous treatment yields the same results for sensitivity and error.

Having in mind a typical indoor environment, norm and direction of the magnetic field vary relatively by far less 1 along the calibration trajectory.

For Hessian matrices H p,p , we use the notations of in Sec. III-C.

ACKNOWLEDGMENTS

We are grateful to ESIGELEC School of Engineering in Rouen-Madrillet, F-BP 10024 76801, for their expertise and support in the use of their optical tracking device.

APPENDIX

A. Least-squares optimization for the toy model of Sec. II

For a set measurements that are homogeneously distributed in the interval (θ 1 , θ 2 ), we consider for the calibration model of Eq. ( 2) the cost function

The bias in this cost function is considered fixed and preestimated with an error δb. Minimization of (27) yields the estimate

with

Note that if θ 1 and θ 2 span over a full interval of length 2π, Eq. ( 28) shows that ĉ is independant from b 0 and therefore insensitive to any possible error in it.

For c * the true value of the scale factor vector, the (noise-free) measurement is M (θ) = c * R (θ) B n + b * and the LS estimate (28) contains an error

The remaining integrations and matrix products are elementary and yield the results reported in Sec. II.