Elise Grosjean
email: grosjean@mathematik.uni-kl.de

Yvon Maday
email: yvon.maday@sorbonne-universite.fr

A doubly reduced approximation for the solution to PDE's based on a domain truncation and a reduced basis method: Application to Navier-Stokes equations

This paper focuses on the non-intrusive reduced basis (NIRB) method called the two grids method. It is used for the simulation of parametric partial differential equations to reduce the associated computational costs of a High-Fidelity code when such problems must be solved for a large number of parameter values or provide a solution in "real time".

As other reduced basis approaches, the "offline step" relies on the High-Fidelity method with a fine enough grid. On the contrary, the non-intrusiveness of the original two grids method is based on the use, in the "online step", of the same method with a much coarser grid, which considerably reduces the cost of this step.

We extend here this idea by further reducing the "online step" and further simplify the High-Fidelity method. As an example of application we consider a classical fluid problem, the 2D Backward Facing Step (BFS). We simplify the model by i) truncating the outflow part of the channel at extreme and ii) using a coarse uniform mesh instead of refining it at the re-entrant corner, both choices that contradict what is required to get a high fidelity representation of the flow. To accomplish this, we create two reduced bases and a deterministic process that allows us to pass from one to the other.

Several numerical simulations illustrate the ability of this new approach.

Introduction

Reduced Basis Methods (RBM) are widely used to approximate the solutions of parameter-dependent problems for a large number of parameter values within a parameter space. The latter enables the definition of the solutions manifold, which contains the solutions for all parameters within the parameter space. RBM are based on the assumption that the solutions manifold has a small Kolmogorov N-width [START_REF] Kolmogoroff | Über die beste annaherung von funktionen einer gegebenen funktionenklasse[END_REF]. In other words, they take advantage of the fact that the solutions for a wide range of parameters share similar physical behaviors, so that only a few well-chosen solutions are sufficient to approximate well any element of the solution manifold by linear combination. These carefully selected elements are known as snapshots. They are exploited to construct a reduced space that must be close to the solution manifold. Then for a new parameter, the associated solution may be approached very quickly but accurately by projecting it onto this space.

Several RBM exist [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations: an introduction[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Iliev | A Two-grid Infinite-volume/reduced Basis Scheme for the Approximation of the Solution of Parameter Dependent PDE with Applications the AFFFF Devices[END_REF][START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Stabile | Pod-galerkin reduced order methods for cfd using finite volume discretisation: vortex shedding around a circular cylinder[END_REF][START_REF] Casenave | A nonintrusive reduced basis method applied to aeroacoustic simulations[END_REF]. They usually take advantage of a procedure decomposition into two steps, one "offline" and one "online". The snapshots are computed using a High-Fidelity (HF) code during the offline portion. Any classical method, such as the Finite Element Method (FEM) or Finite Volume schemes (FVM), may be used in this code. Then, during the online part, a reduced problem is solved.

The ability to modify the HF code is frequently required for an efficient implementation of these methods. Yet, if the code has been purchased, as it is frequently the case in industry, this may be impossible. This is why we focus on the two grids method, which is a Non-Intrusive Reduced Basis (NIRB) algorithm [START_REF] Chakir | A two-grid finite-element/reduced basis scheme for the approximation of the solution of parametric dependent p[END_REF][START_REF] Chakir | Contribution à l'analyse numérique de quelques problèmes en chimie quantique et mécanique[END_REF][START_REF] Chakir | A non-intrusive reduced basis approach for parametrized heat transfer problems[END_REF][START_REF] Grosjean | Error estimate of the non-intrusive reduced basis method with finite volume schemes[END_REF] (see also different NIRB methods from the two-grid algorithm [START_REF] Casenave | A nonintrusive reduced basis method applied to aeroacoustic simulations[END_REF][START_REF] Geelen | Operator inference for non-intrusive model reduction with quadratic manifolds[END_REF][START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF]). Indeed, it does not require any change to the HF code and instead employs it solely as a "black-box" solver. We refer to [START_REF] Guo | Bayesian operator inference for data-driven reduced-order modeling[END_REF][START_REF] Thomas | Learning high-dimensional parametric maps via reduced basis adaptive residual networks[END_REF] for another approach of non-intrusiveness based on the learning of reduced operators through neural network framework.

Reminders on the classical NIRB two-grid method

Let Ω be a bounded domain in R d , with d ≤ 3 and a smooth enough boundary ∂Ω, and consider a parametric problem P on Ω. Here, we consider time-independent problems (see [START_REF] Grosjean | Error estimate of the non-intrusive reduced basis (nirb) two-grid method with parabolic equations[END_REF] for parabolic equations). Let µ represent a parameter of P chosen in a given set G and u(µ) the associated solution. The manifold of all solutions is defined as S = {u(µ), µ ∈ G}.

Thanks to numerous HF solutions of P already computed for various parameter values, we want to approximate quickly the solution for a new parameter µ ∈ G.

In what follows, we consider a FEM HF code to approximate any solution in S.

The particularity of the two-grid method is to be based on two FEM meshes:

• one fine mesh T h , with which the HF snapshots during the offline part are generated,

• one coarse mesh T H , where the problem solution is computed online with the new parameter µ ∈ G.

The latter is what makes the method so efficient. Because it is based on a coarse FEM approximation, the NIRB solution computation for the new parameter µ is very inexpensive and thus allows to recover the full efficiency of the classical RBM. The approximation is then improved in a very brief runtime by projection and rectification onto the reduced space using a Reduced Basis (RB) of this space created independently during the offline part.

The FEM HF solutions calculated on T h are denoted by u h (µ), while the FEM solution approximations computed on the coarse mesh T H are denoted by u H (µ).

The offline-online decomposition of this method is as follows:

• "Offline step"

The offline part can be thought as a learning procedure and allows us to construct the reduced space. All expensive parts of the algorithm are performed during this step.

1. We begin by defining a large enough training set of parameters

G train := ∪ i∈{1,...,N train } µ i .
Then, in order to generate the reduced space, we use a greedy technique that selects appropriate parameters (µ i ) i=1,...,N within G train . As it is standard in RBM framework, at each greedy iteration, the new RB parameter in G train is chosen such that its associated solution is the poorest approached by the former basis space. The algorithm 1 describes this greedy procedure. Note that another way for generating the RB is to use the Snapshots POD [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations: an introduction[END_REF][START_REF] Ullmann | Pod-galerkin reduced-order modeling with adaptive finite element snapshots[END_REF][START_REF] Stabile | Pod-galerkin reduced order methods for cfd using finite volume discretisation: vortex shedding around a circular cylinder[END_REF][START_REF] Luo | A reduced stabilized mixed finite element formulation based on proper orthogonal decomposition for the non-stationary navier-stokes equations[END_REF].

In algorithm 1, inexpensive a-posteriori error estimates may replace u h (µ) -P k-1 (u h (µ)) L 2 and reduce its computational time [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF][START_REF] Manzoni | An efficient computational framework for reduced basis approximation and a posteriori error estimation of parametrized navier-stokes flows[END_REF][START_REF] Burda | A posteriori error estimates applied to flow in a channel with corners[END_REF][START_REF] Veroy | Reduced-basis approximation of the viscous Burgers equation: rigorous a posteriori error bounds[END_REF][START_REF] Huynh | A static condensation reduced basis element method : approximation and a posteriori error estimation[END_REF].

In practice, the greedy algorithm is halted with a stopping criterion such as an error threshold or a maximum number of basis functions that can be generated. At the end of the algorithm (as in algorithm 1 above), the space X N h := Span{Φ h 1 , . . . , Φ h N } is thus provided with an orthonormal basis. 2. Then, in order to stabilize the approximation we can propose, as an option, a new basis set: we solve the following eigenvalue problem:

Find (Ψ h , λ) ∈ (X N h , R) such that: ∀v ∈ X N h , Ω ∇Ψ h • ∇v dx = λ Ω Ψ h • v dx, (1)
We get an increasing sequence of eigenvalues λ i , and orthogonal eigenfunctions (Φ h i ) i=1,••• ,N , orthogonalized both in L 2 (Ω) and H 1 (Ω) that we choose to normalize in L 2 (Ω). It is worth noting that the Gram-Schmidt method of the traditional greedy algorithm 1 yields only an L 2 (Ω)-orthonormalized RB. As a result, as detailed in [START_REF] Grosjean | Variations and further developments on the Non-Intrusive Reduced Basis two-grid method[END_REF], it stabilizes the NIRB approximation with respect to N. For the sake of simplicity, we continue to name Φ this improved reduce basis i.e. we set

Φ h k = Ψ h k . 1: Input: tol, {u h (µ 1 ), • • • , u h (µ N train ) with µ i ∈ G train }. 2: Output: Reduced basis {Φ h 1 , • • • , Φ h N }. 3: Choose µ 1 = argmax µ∈G train u h (µ) L 2 (Ω) , 4: Set Φ h 1 = u h (µ 1 ) u h (µ 1 ) L 2 (Ω) . 5: Set G 1 = {µ 1 } and X 1 h = span{Φ h 1 }. 6: for k = 2 to N do 7: µ k = arg max µ∈G train \G k-1 u h (µ) -P k-1 (u h (µ)) L 2 (Ω) , 8: with P k-1 (u h (µ)) := k-1 ∑ i=1 (u h (µ), Φ h i ) Φ h i . 9: Compute Φ h k = u h (µ k ) -P k-1 (u h (µ k )) and set Φ h k = Φ h k Φ h k L 2 (Ω)
.

10:

Set G k = G k-1 ∪ {µ k } and X k h = X k-1 h ⊕ span{Φ h k }.
11:

Stop when u h (µ)

-P k-1 (u h (µ)) L 2 (Ω)
≤ tol, ∀µ ∈ G train .

12: end for Algorithm 1: Greedy algorithm 3. Finally, as a further option, in order to improve the NIRB two-grid approximation, we can now employ a post-treatment step. We remark that the projection of the coarse FEM solution onto X N h differs from the projection of the HF snapshot u h (µ k ) computed in the offline stage for any parameter of the RB space µ k , k = 1, . . . , N. Thus, we resort to the "rectification post-processing", as introduced in [START_REF] Chakir | A two-grid finite-element/reduced basis scheme for the approximation of the solution of parametric dependent p[END_REF], to improve the NIRB accuracy. To this purpose, we create a rectification matrix denoted R. This matrix is constructed from the HF snapshots {u h (µ k )} k=1,...,N and coarse snapshots generated on T H with the same parameters µ k as for the fine ones. For all i = 1, • • • , N, and for all µ k ∈ G train , we define the coarse coefficients as

A k,i = Ω u H (µ k ) • Φ h i dx, (2)
and the fine coefficients as

B k,i = Ω u h (µ k ) • Φ h i dx. (3)
Then, we compute the vectors

R i = (A T A + δI N ) -1 A T B i , i = 1, • • • , N, (4)
where I N refers to the identity matrix and δ is a regularization parameter.

• "Online step"

The online part of the algorithm is much faster than a fine HF FEM evaluation.

1. We solve the problem P on the coarse mesh T H with the new parameter µ ∈ G we are interested in.

2. Finally, the approximation used in the two-grid method (with the rectification post-treatment step) is

R[u N Hh ]( µ) := N ∑ i,j=1 R ij (u H ( µ), Φ h j ) Φ h i , (5)
where (R) ij = R ij is the rectification matrix, given by (4).

Main contributions

The original version of the NIRB method still suffers from the cost of the online step which is certainly reduced by the use of the coarse grid but may still represent a high cost.

To summarize our contributions, we propose and implement two new ideas to reduce the costs of the algorithm's online component.

• The first one relies on the truncation of the coarse computational domain and the use of a second reduced basis set. The method is described in section 2 and illustrated in 2.3 with several numerical results on a model problem which is a classical fluid problem that of the Backward-Facing Step, described in section 2.2.

• Then, we propose a new strategy for improving NIRB approximation accuracy when dealing with domain singularities without increasing online stage computational times. This is the purpose of the section 3.1 where we present new results on the BFS problem.

Domain truncation

The underlying idea is to consider that, for a parameterized PDE, posed on a domain Ω, some main characteristics of its solution depend little on the exact shape of the domain. Some specificities of the domain must of course be respected (e.g. distinctive aspects of the boundary, corners, edges, or aspect ratios) but one can think of reducing the domain (allowing to reduce the number of degrees of freedom of the discretization). If, having thus reduced the size of the domain, the two solutions for the same parameters of the PDE can be associated in a bijective way, one can approach the first one (called virtual solution on the reduced domain) to intuit the second one (called real solution on the initial domain). We propose to learn how to go from the virtual solution to the real solution.

To allow this learning, we will look at the solutions under the contracted form of their writing in reduced space... hence the appearance of a second reduced basis associated to the solutions on the reduced domain.

Description of the Non-Intrusive Reduced Basis tool with domain truncations

Let us consider a "truth" domain Ω ∈ R d with d ≤ 3 employed for the simulation, adapted to the parameterized problem.

We then introduce ω ⊂ Ω, a truncation of Ω, reduced with fictitious boundaries. The NIRB method with the domain truncation involves two partitioned meshes:

• the classical fine mesh T h on Ω,

• and one coarse mesh T H on ω,

where h and H are the respective sizes of the meshes and h << H. The size h (resp. H) is defined as

h = max K∈M h h K (resp. H = max K∈M H H K ), (6)
where the diameter h K (or H K ) of any element K in a mesh is equal to sup

x,y∈K |x -y|.
Let µ ∈ G be the new parameter we are interested in.

• The "offline" part of the algorithm is costly in time but only done once.

1. Point 1 of the classical NIRB two-grid method remains the same (see section 1.1 and algorithm 1). We generate the RB (Φ h i ) i=1,...,N 1 from suitable parameters (µ 1 , . . . , µ N 1 ) ∈ G N 1 train and the reduced space [START_REF] Chakir | A non-intrusive reduced basis method for elastoplasticity problems in geotechnics[END_REF], where two POD RB are used to reduce the condition number of the rectification matrix (4), we create another RB. We generate this new RB denoted (Φ

X N 1 h = span{Φ h 1 , . . . , Φ h N 1 }. Taking inspiration from
H i ) i=1,...,N 2 from several parameters (µ ′ 1 , . . . , µ ′ N 2 ) ∈ G N 2 train .
It is constructed on the coarse mesh T H by repeating point 1 again but with the coarse snapshots instead of the fine ones. We denote by

X N 2 H = span{Φ H 1 , . . . , Φ H N 2 } the
new reduced space approaching the manifold of the coarse solutions on the truncated domain ω ⊂ Ω.

The coarse solutions are represented by the coarse RB well (respectively the fine solutions with the fine RB). The number of necessary basis functions N 1 and N 2 may be chosen thanks to a small a priori given treshold (tol in algorithm 1). Note than point 2 of the classical NIRB two-grid method may be applied for both RB to obtain H 1 orthogonalized functions.

2. We deal with two kinds of coefficients:

the fine ones, which correspond to the L 2 (Ω)-inner product of the fine snapshots with the fine RB on the whole domain Ω, -and the coarse coefficients, which are the L 2 (ω)-inner product of the coarse snapshots with the coarse RB on the subdomain ω.

In analogy with the two-grid rectification post-treatment, we introduce a rectification matrix T which allows us to go from the new coarse coefficients to the fine ones. For all i = 1, • • • , N 2 , and for all µ k ∈ G train , we define these new coarse coefficients as

C k,i = ω u H (µ k ) • Φ H i dx, (7)
and the fine coefficients are defined by ( 3) with N = N 1 .

The aim of this processing stage is to seek

T i minimizing CT i -B i 2 .
The solution of this minimization problem with a regularization parameter [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF] are the rectification vectors:

T i = (C T C + λI N 2 ) -1 C T B i , ∀i = 1, • • • , N 1 , (8)
where λ is the regularization parameter.

• The "online" step is then performed on the coarse mesh T H , which is much less expensive than an HF computation. The online part is even less computationally expensive than the usual NIRB online algorithm because the coarse mesh is computed on ω ⊂ Ω and no interpolation is required to compute the coarse coefficients [START_REF] Carpinteri | Stress-singularity and generalized fracture toughness at the vertex of re-entrant corners[END_REF]. The coarse solution coefficient is then corrected using the matrix T before being projected onto the fine RB.

Thus, it consists in:

1. solving the parameterized problem with the FEM solver on the coarse mesh T H for the new parameter µ ∈ G.

setting the new NIRB approximation as

T[u N Hh ](µ) = N 1 ∑ i=1 N 2 ∑ j=1 T ij (u H ( µ), Φ H j ) L 2 (ω) Φ h i . ( 9
)

A problem model: the 2D Backward Facing Step (BFS)

We consider the 2D Backward Facing Step (BFS) with steady Navier-Stokes. It simulates a fluid flowing through a channel with a descending stair and constitutes an important test for numerical simulations in fluid mechanics. It has been widely studied [START_REF] Erturk | Numerical solutions of 2-d steady incompressible flow over a backward-facing step, part i: High reynolds number solutions[END_REF][START_REF] Hossain | Numerical investigation of fluid flow through a 2d backward facing step channel[END_REF][START_REF] Durst | Flows over two-dimensional backward-facing steps[END_REF][START_REF] Lee | Experimental and numerical investigation of 2-d backward-facing step flow[END_REF][START_REF] Gresho | Is the steady viscous incompressible two-dimensional flow over a backward-facing step at re= stable?[END_REF].

The flow is laminar and incompressible in our case, and the channel domain Ω ∞ is assumed to be infinite, with upper and lower walls. To cope with this infinite domain in simulations, we introduce artificial boundaries and use Neumann boundary conditions on the channel's exit border and Dirichlet conditions at the entrance. Thus let us consider as our spatial domain Ω ⊂ Ω ∞ , the channel used for the simulations. The length of the channel to simulate the flow is chosen such that we observe • a Poiseuille flow and a maximum velocity at the upstream of the channel,

• a vortex at the concave corner behind the step,

• and a Poiseuille flow at the backstream of the channel.

The artificial boundaries must not be too close, so as to not interfere with the modeling, and not too far, to avoid too expensive calculations. The conditions on these artificial boundaries result from the modeling phase.

The channel Ω, which is already a truncation of the infinite channel Ω ∞ , is described by Ω = (0, L 1 ] × (H 1 -H 0 , H 1 ) ∪ (L 1 , L) × (0, H 1 ) as shown in Figure 1.

H 0 Γ in L 1 L H 1 Γ out X Y Z Figure 1:
The flow channel Therefore, we consider the following parameter-dependent problem P as our model problem: Find the velocity of the fluid u ∈ H 1 (Ω) 2 and the pressure p ∈ L 2 (Ω), solutions of the steady Navier-Stokes equation

                       u.∇u -ν∆u + ∇p = 0 on Ω, ∇.u = 0 on Ω, u(x) =   4 H 2 0 H 1 -x 2 x 2 -H 1 + H 0 0   on Γ in , ν∂ n u |Γ -pn = 0 on Γ out , u(x) = 0, on ∂Ω\(Γ out ∪ Γ in ), (10)
where x = (x 1 , x 2 ) ∈ Ω, u = (u 1 , u 2 ), and ν = 1 Re , Re being the Reynolds number. The solution of this problem is uniquely defined by the data [ν, H 0 , L 1 , H 1 , L]. We consider µ := Re ∈ G = [START_REF] Iliev | A Two-grid Infinite-volume/reduced Basis Scheme for the Approximation of the Solution of Parameter Dependent PDE with Applications the AFFFF Devices[END_REF]300] as the problem variable parameter. Note that, for Re > 300, there is no stable steady solution with the geometric data we consider H 0 = 0.5, H 1 = 1, and L 1 = 1.

Numerical results on the model problem

We have carried out numerical simulations for equations [START_REF] Chakir | A non-intrusive reduced basis method for elastoplasticity problems in geotechnics[END_REF] in FreeFem++ (version 4.9) [START_REF] Hecht | Bamg: bidimensional anisotropic mesh generator[END_REF] with the Newton algorithm for the numerical treatment of the nonlinear term, and classical Taylor-Hood finite elements P 2 -P 1 . The fine mesh size is equal to 0.03, and the coarse one to 0.21. We have used a finer mesh with a size equal to 0.016 as a surrogate for the solutions of reference in order to compute the errors accurately.

For the fine mesh T h , we have set the length of the channel (which does not change) to L = 5. Thus, in what follows, L refers to the length of the subdomain ω. We recall that the Reynolds number Re ∈ G = [START_REF] Iliev | A Two-grid Infinite-volume/reduced Basis Scheme for the Approximation of the Solution of Parameter Dependent PDE with Applications the AFFFF Devices[END_REF]300] is the problem parameter, and the domain Ω is truncated according to L ∈ [START_REF] Abdulle | Reduced basis finite element heterogeneous multiscale method for high-order discretizations of elliptic homogenization problems[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF] (remember that L ≤ L 1 = 1). The size of the coarse mesh H < 1 is such that H 2 ≃ h. For instance, Figure 2 represents the two partitioned meshes with L = 3 for the length of the coarse mesh. When the Reynolds number increases, the recirculation after the step is expanded, as shown in Figure 3. Thus, we also consider two subgroups of G: • between the u re f and the coarse and fine HF solutions,

G 1 = [30, 150[ with G 1,train = {30 + 15i, i = 0, . . . , 7}, (11)
• obtained by the classical NIRB algorithm of section 1.1, with the rectification postprocessing step (4), as illustrated in Figure 5. With the classical NIRB (5), the error is given by

u re f ( µ) -u N Hh ( µ) H 1 (Ω) u re f ( µ) H 1 (Ω)
. We find that the error between u re f ( µ) and u h ( µ) can be recovered if the new Reynolds parameter belongs to the parameters group used for the RB generation. It requires 4 modes for an RB generated from the group G, whereas 3 modes are sufficient for the parameter subgroups G 1 and G 2 .

We present several tests with 3 different values for L:

• The first situation is with Ω = ω (no truncation). Figure 6 depicts the new NIRB error. With the group G, only 3 and 4 modes are required to retrieve the fine solutions accuracy for Re = 52 and Re = 233, respectively. With the subgroups G 1 and G 2 3 modes are sufficient in both cases. To make the presentation of these results easier, we have first taken N 1 = N 2 . This test provides more accurate results with the new method (i.e. with ( 9)) than the classical rectified NIRB (i.e. with ( 5)) in the case of parameter extrapolation (for example, Re = 52 with an RB generated from the subgroup G 2 ). In fact, unlike the rectified classical version, the error is almost recovered with Re = 52 (with 9 modes).

• Then, we have tested the new algorithm with L = 3. We get the same precision as with ω = Ω as shown in Figure 7. However, if the new parameter does not belong to the RB subgroup range (G 1 or G 2 ), the error is not recovered with Re = 52 and some instabilities are observed.

• Finally, the last tests have been carried out with L = 1.2. Because the recirculation area behind the step is cut, this case is much more complex. As a consequence, the results show more instabilities. The importance of λ in the new rectification process ( 8) is highlighted in this test. Indeed, we compare results obtained with λ = 0 in Figure 8a to those carried out with λ = 10 -10 in Figure 8b. With this regularization parameter, all the errors remain stable, demonstrating the ability of this new algorithm to retrieve accurate approximations even with highly truncated domains.

To highlight the influence of N 1 and N 2 , we compare the error between u h ( µ) and the new NIRB approximation with L = 1.2, L = 3 and L = 5. The errors between the fine solution and the new NIRB approximation are given by

u h ( µ) -T[u N Hh ]( µ) H 1 (Ω) u h ( µ) H 1 (Ω) . ( 13
)
They are presented using heatmaps. Figure 9 presents the results with ω = Ω. Although there may be some instabilities when the matrix C T C (8) is inverted, as we can see with N 2 = 18 with λ = 0, the error reaches a lower threshold than with the classical rectified NIRB approximation. As shown in Figure 10, the term λ set to 10 -10 allows us to stabilize the errors. We remark in Figure 11 that with L = 3, the errors (13) remain smaller than the obtained from the classical NIRB with rectification, with both Reynolds and still quite low when L = 1.2 (Figure 12 and Figure 13). Table 1 shows its great ability to reduce computational times. We end this section with the new NIRB errors (for u 1 and u 2 ) and the new NIRB approximations with L = 3 and Re = 52, N 1 = 7, N 2 = 8 (Figure 14).

We present the FEM and NIRB runtimes in table 1.

NIRB with domain singularities

Adaptation of the new method to the problem

There are numerous studies on the behaviour of the solutions to elliptic PDE's in domain with re-entrant corners and domain singularities in the literature [START_REF] Carpinteri | Stress-singularity and generalized fracture toughness at the vertex of re-entrant corners[END_REF][START_REF] Brenner | Multigrid methods and stress intensity factors i: Corner singularities[END_REF][START_REF] Wahlbin | Local behavior in finite element methods[END_REF][START_REF] Wahlbin | On the sharpness of certain local estimates for • H 1 projections into finite element spaces: Influence of a reentrant corner[END_REF][START_REF] Kim | A finite element method for computing accurate solutions for poisson equations with corner singularities using the stress intensity factor[END_REF][START_REF] Igarashi | A boundary element method for potential fields with corner singularities[END_REF][START_REF] Dauge | Elliptic boundary value problems on corner domains: smoothness and asymptotics of solutions[END_REF][START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF]. The solution lacks regularity and thus, classical uniform regular triangulations do not allow for an approximation with the optimal rate of convergence. In the case of our numerical analysis in the previous section, because the channel has one re-entrant corner in the BFS problem, it generates a local singularity on the FEM solution. It implies that the latter does not globally lie in H 2 (Ω). In fact, the solution belongs to H 1.54448 (Ω) (see, e.g. [START_REF] Bernardi | Méthodes d'éléments finis mixtes pour les équations de stokes et de navier-stokes dans un polygone non convexe[END_REF]) . Thus, a standard FEM yields reduced convergence rates and the following estimates hold for any µ ∈ G with an uniform mesh of size h :

   u(µ) -u h (µ) H 1 (Ω) ≤ Ch 0.54448 f L 2 (Ω) , u(µ) -u h (µ) L 2 (Ω) ≤ Ch 1.08896 f L 2 (Ω) ,
as observed in Figure 15 with the H 1 semi-norm and the L 2 error between the reference and the coarse meshes with the Reynolds number Re = 52.

To better approach the solution of such problems and compensate the effects of singularities, most methods rely on adaptive refinements in the vicinity of the singularities, which complicates the resolution of the discrete system, which may not be desirable in the online phase of a NIRB method.

We propose here to follow the same paradigm as in the previous section based on the fair representation of the coarse approximation in terms of its coefficients in a reduced basis that, by a learning appraoch, allows to retrieve a fine and accurate solution.

The strategy is thus to build, in the first "offline" stage of the NIRB method, an accurate reduced basis relying on an adapted mesh, to correctly capture the singularity and the recirculation zone near the domain reentrant corner. Other methods can also be used. For instance, we may supplement the classical FE basis with singular functions (xFEM) [START_REF] Foucard | An x-fem-based numerical-asymptotic expansion for simulating a stokes flow near a sharp corner[END_REF][START_REF] Rabinovich | Xfem-based crack detection scheme using a genetic algorithm[END_REF]. Another approach employs boundary element methods [START_REF] Igarashi | A boundary element method for potential fields with corner singularities[END_REF]. We have used an adaptative mesh refinement strategy because they are widely used [START_REF] Abdulle | Reduced basis finite element heterogeneous multiscale method for high-order discretizations of elliptic homogenization problems[END_REF][START_REF] Brenner | Multigrid methods and stress intensity factors i: Corner singularities[END_REF]. They may be based on a-posteriori error estimates [START_REF] Burda | A posteriori error estimates applied to flow in a channel with corners[END_REF][START_REF] Verf | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF][START_REF] De Sr Gago | A posteriori error analysis and adaptive processes in the finite element method: Part ii-adaptive mesh refinement[END_REF].

Then we use a coarse uniform mesh "online" to compute an approximate solution (that suffers from the above lack of optimal convergence) that is transformed in a reduced basis representation.

We summarize the offline/online strategy:

1. During the "offline" stage, the fine mesh is refined around the re-entrant corners in order to obtain accurate results. Then the offline part of the classical NIRB two-grid with rectification (section 1.1) or its new version with two RB and an appropriate domain truncation (section 2) is applied with this mesh.

2.

Then, during the "online" part, the coarse approximation is computed on a uniform coarse mesh, and projected and rectified onto the reduced space X N h .

Numerical results

To summarize, in our tests [START_REF] Chakir | A non-intrusive reduced basis method for elastoplasticity problems in geotechnics[END_REF], we have employed three meshes:

• a uniform coarse mesh for the coarse FEM approximation,

• a fine mesh refined around the re-entrant corner for the fine snapshots generation,

• and a highly refined mesh for the reference solution to represent the exact solution.

These refinements have been done in FreeFem++ with the "adaptmesh" function. It uses a variable metric/Delaunay automatic meshing algorithm [START_REF] Hecht | Bamg: bidimensional anisotropic mesh generator[END_REF][START_REF] Hecht | New development in freefem++[END_REF] (bamg software). With the refinement, the fine mesh size is now defined as

h = min K∈M h h K ,
where the diameter h K of any element K in a mesh is equal to sup

x,y∈K |x -y|, K ∈ T h . With this definition, the reference mesh size is equal to 2.54 × 10 -5 .

We have generated N 1 = N 2 = N = 6 snapshots (Re = 5 + 20i, i = 0, . . . , 5) and have tested the new algorithm detailed in section 2 with the truncation parameter L = 3 (ω ⊂ Ω).

We have experimented four fine meshes with a refinement in the vicinity of the re-entrant corner, as in Figure 17a. We have then used a uniform coarse mesh of size H ≃ 0.32 and then an even coarser uniform coarse mesh of size H ≃ 0.70 as in Figure 16a. We compare in table 2 the classical NIRB approximation (with the Figure 15: Display of the convergence in H 1 (left) and L 2 (right) of the discrete solution error computed on different uniform grids rectification post-treatment) with the new approach (Ω = ω and ω ⊂ Ω). The three NIRB methods recover the fine (refined) mesh accuracy (note that the original NIRB method benefits from the rectification step, indeed, without it, the inaccuracy resulting from the use of a coarse uniform mesh would prevent to recover the fine accuracy). These results are thus promising and illustrate the relevance of this approach.

= 3 (ω ⊂ Ω) NIRB + Rectifi- cation New NIRB ω = Ω New NIRB with L = 3 (ω ⊂ Ω) 1.99E -4 1.93E -3 1.93E -3 1.93E -3 1.93E -3 1.93E -3 1.93E -3 1.95E -3 1.79E -4 1.70E -3 1.71E -3 1.71E -3 1.71E -3 1.71E -3 1.71E -4 1.73E -3

Time execution (min,sec)

We present the FEM and NIRB runtimes with the meshes as in table 16 for the fine mesh and table 17 for the coarse one, in table 3.

Conclusion

We have presented in this paper a new NIRB approach based

• on the construction of two reduced basis: one from some fine and accurate approximation of the solution manifold and one from a coarse approximation,

• and on a truncation of the spatial domain for the coarse approximation.

This new method allows us to use a much coarser mesh during the online stage of the two-grid algorithm, compared to the original method. Indeed, the truncation of the coarse mesh allows to further reduce its complexity. This is what makes this approach so efficient in time. Yet, its size and the domain truncations must be physically acceptable so that there still exist a bijection between the coarse approximations and the associated fine and accurate ones. From the reduced basis coarse components we recover the components of the solution in the fine reduced basis through a deterministic processing stage.

We have applied this new process to the well-known 2D BFS, and obtained the same accuracy than the one given by the fine solutions, while using several truncations.

In the original NIRB two-grid algorithm, the coarse solution needs to be interpolated onto the fine mesh in order to compute the L 2 -projection onto the reduced space X N h . Note that the interpolating matrix can be computed during the offline stage but is nevertheless complex. Here with this new algorithm, no interpolation is needed, and in particular, the L 2 -inner product between the coarse solution and the coarse RB is performed on the coarse mesh only.

We even closely consider the effect of the HF simulation involving a mesh refinement to further deal with the singularity of the BFS domain (section 3.1). We numerically show that both the NIRB approaches (original NIRB with rectification and the new tool) allow us to recover the fine FEM accuracy, despite the fact that the coarse approximation is using a uniform (coarse) mesh as observed in table 3.

We believe that this new method may be successfully applied to a wide variety of problems, with complex domains having large number of degrees of freedom, and therefore the computational times will be significantly reduced with a truncated domain.

  (a) Fine mesh on Ω, h = 0.03 (b) Coarse mesh on ω, H = 0.21

Figure 2 :

 2 Figure 2: Fine and coarse meshes

G 2 =

 2 [150, 300] with G 2,train = {30 + 15i, i = 8, . . . , 18},[START_REF] Chakir | A non-intrusive reduced basis approach for parametrized heat transfer problems[END_REF] and G =[START_REF] Iliev | A Two-grid Infinite-volume/reduced Basis Scheme for the Approximation of the Solution of Parameter Dependent PDE with Applications the AFFFF Devices[END_REF] 300], G train = G 1,train ∪ G 2,train , from which we extract, with three greedy methods, three reduced bases of size N 1 = 1, . . . , 8, N 2 = 1, . . . , 11, and N = 1, . . . ,19.

Figure 3 :

 3 Figure 3: Velocities (u 1 and u 2 ) for Reynolds=52 (left) and Reynolds=233 (right) Then we have generated the FEM solutions of two new parameters in G: Re = 52 and Re = 233. The coarse solutions for Re = 52 on the truncated subdomains with L = 1.2 and L = 3 are presented in Figure 4. We have tried several values for the parameter L ∈ [1, 5]. In the following figures, we present the results of the new NIRB algorithm in order to analyze its limits and observe the effects of domain truncations on the phenomena. Our results diplay the H 1 relative error

Figure 4 :

 4 Figure 4: Coarse velocities u 1 (up), u 2 (bottom), Re = 52 on the truncated domain with L = 3 (left), and with L = 1.2 (right)

(a) Re = 52 (b) Re = 233 Figure 5 :

 522335 Figure 5: classical NIRB errors with the rectification postprocessing step

Figure 6 :

 6 Figure 6: FEM vs NIRB errors: between the reference solution and the NIRB approximations (classical algorithm with rectification vs new NIRB algorithm with ω = Ω) compared to FEM projection errors (dotted lines), tests with the new parameters Re = 52 (left) and Re = 233 (right) (λ = 0).

Figure 7 :

 7 Figure 7: FEM vs NIRB errors: between the reference solution and the NIRB approximations (classical algorithm with rectification vs new NIRB algorithm with ω ⊂ Ω) with L = 3 compared to FEM projection errors (dotted lines), tests with the new parameters Re = 52 (left) and Re = 233 (right) (λ = 0).

(a) λ = 0 (

 0 b) λ = 10 -10

Figure 8 :

 8 Figure 8: FEM vs NIRB errors: between the reference solution and the NIRB approximations (classical algorithm with rectification vs new NIRB algorithm with ω ⊂ Ω) with L = 1.2 compared to FEM projection errors (dotted lines), tests with λ = 0 (left) and λ = 10 -10 (right).

Figure 9 :

 9 Figure 9: Error of the new NIRB algorithm with ω = Ω (no truncation, L = 5)

Figure 10 :

 10 Figure 10: NIRB relative H 1 errors with the new algorithm ω = Ω (L = 5), λ = 10 -10 , Re = 52.

Figure 11 :Figure 12 :

 1112 Figure 11: Error of the new NIRB algorithm with ω ⊂ Ω λ = 0 and G: L = 3

2 Figure 13 :Figure 14 :

 21314 Figure 13: Error of the new NIRB algorithm for the new Reynolds number Re = 233 with G 2 : L = 3 and λ = 0 (left) and with L = 1.2and λ = 10 -10 (right)

Figure 16 :Figure 17 :

 1617 Figure 16: FEM coarse mesh and solution

Table 1 :

 1 FEM runtimes min:sec

		FEM high fidelity solver	FEM coarse solution
		00:43	00:01
		NIRB Offline (N = 18)	classical NIRB online
		14:30	00:10
		new NIRB offline (N 1 = N 2 = 18)	new NIRB online
	ω = Ω	14:52	00:10
	L = 3	14:52	00:09
	L = 1.2	14:50	00:09

Table 2 :

 2 relative errors with the semi H 1 -norm, refined meshes N 1 = N 2 = N = 6 with λ = 0

			NIRB with H = 0.320156	NIRB with H = 0.707107
	h	FEM	NIRB +	New	New
		with	Rectifi-	NIRB	NIRB
		the fine	cation	ω = Ω	with
		refined			L
		mesh

Table 3 :

 3 FEM runtimes min:sec

	FEM HF solver	FEM coarse solution H = 0.32
	01:24	00:01
	NIRB Offline (N = 6)	NIRB online
	08:34	00:16

Acknowledgments

This research was partially funded by the French Fonds Unique Interministeriel (MOR DICUS).


Sara Behforootan
email: s.behforootan@imperial.ac.uk

Madelaine Thorniley

Jean-Gabriel Minonzio

Oliver Boughton

Monil Karia

Rajarshi Bhattacharya

Ulrich Hansen

Justin Cobb

Richard Abel

Can guided wave ultrasound predict bone mechanical properties at the femoral neck in patients undergoing hip arthroplasty?

Keywords: patients undergoing hip arthroplasty? Ultrasound, Guided wave, Osteoarthritis, Bone, cortical bone, mechanical properties

Introduction

A surgeon with reliable information about a patient's bone quality can make informed decisions on the type of implant and the impaction force used during orthopaedic surgery [START_REF] Boughton | Measuring Bone Stiffness Using Spherical Indentation[END_REF]. Current techniques, however, cannot accurately characterise the mechanical properties of the bone. A portable, non-expensive and non-invasive tool is needed to assess bone health. Ultrasound has been investigated extensively for bone assessment, but the technology has not been established in clinical practice yet [START_REF] Raum | Clinical Devices for Bone Assessment[END_REF][START_REF] Tran | Ultrasonic Guided Waves in Bone: A Decade of Advancement in Review[END_REF].

Important negative outcomes from surgery are intraoperative, and post-operative fractures, and implant loosening [START_REF] Boughton | Measuring Bone Stiffness Using Spherical Indentation[END_REF]. Intraoperative and post-operative fractures are relatively common in hip arthroplasty [START_REF] Abdel | Epidemiology of Periprosthetic Fracture of the Femur in 32 644 Primary Total Hip Arthroplasties: A 40-Year Experience[END_REF], and orthopaedic surgeons need to make more informed decisions to reduce the risk [START_REF] Issack | Fat Embolism and Respiratory Distress Associated with Cemented Femoral Arthroplasty[END_REF][START_REF] Rhyu | Does Osteoporosis Increase Early Subsidence of Cementless Double-Tapered Femoral Stem in Hip Arthroplasty?[END_REF]. Intraoperative and post-operative fractures can either be a consequence of high loads applied during the surgery or implant loosening and bone fractures around the implant after surgery [START_REF] Doyle | Impaction Technique Influences Implant Stability in Low-Density Bone Model[END_REF][START_REF] Haidukewych | Intraoperative Fractures during Primary Total Hip Arthroplasty[END_REF][START_REF] Valle | Acetabular Reconstruction in Revision Total Hip Arthroplasty[END_REF][START_REF] Winter | Basic Considerations for Determining the Amount of Press Fit in Acetabular Cup Endoprostheses as a Function of the Elastic Bone Behavior[END_REF].

The common technologies used for assessing bone health are X-ray, dual X-ray absorptiometry bone mineral density (DEXA BMD) measurements, or haptic assessment of bone during surgery [START_REF] Tulner | Haptic Assessment of Bone Quality in Orthopedic Surgery: No Consensus but Perspective for High Training Potential[END_REF]. Usually, radiographs are not calibrated against the bone mineral content and therefore, the information that is provided is mainly representative of the shape of the organ [START_REF] Jonasson | Skeletal Bone Mineral Density in Relation to Thickness, Bone Mass, and Structure of the Mandibular Alveolar Process in Dentate Men and Women[END_REF]. Only 4% of patients who undergo hip arthroplasty attend for BMD measurements using DEXA, which is currently the gold standard for bone assessment [START_REF] Maier | The Importance of Bone Mineral Density in Hip Arthroplasty: Results of a Survey Asking Orthopaedic Surgeons about Their Opinions and Attitudes Concerning Osteoporosis and Hip Arthroplasty[END_REF]. Haptic assessment by orthopaedic surgeons during surgery is a subjective method and does not provide accurate, quantitative J o u r n a l P r e -p r o o f measurements [START_REF] Tulner | Haptic Assessment of Bone Quality in Orthopedic Surgery: No Consensus but Perspective for High Training Potential[END_REF]. Hence there is a surgical need for a portable, inexpensive, non-invasive tool to assess bone health before surgery.

Characterising the properties of the femoral neck with ultrasound is challenging, due to the anatomical positioning within the body and the amount of surrounding soft tissue. Therefore, ultrasound at peripheral sites has been investigated extensively. While many devices have been developed and tested, the method has not been established yet [START_REF] Behforootan | Ultrasound and Bone Disease: A Systematic Review OPEN ACCESS[END_REF]. Previous devices primarily tested calcaneus sites, which are too distal from the femoral neck and consist of predominantly cancellous bone [START_REF] Rufus-Membere | Associations between Bone Material Strength Index, Calcaneal Quantitative Ultrasound, and Bone Mineral Density in Men[END_REF].

Additionally, previous devices estimating cortical bone have only been able to acquire a single scan of either the cortical thickness or the cortical porosity [START_REF] Behforootan | Ultrasound and Bone Disease: A Systematic Review OPEN ACCESS[END_REF].

The Azalee ultrasound device is the first device able to measure both the cortical thickness and porosity (J.-G. [START_REF] Minonzio | Ultrasound -Based Estimates of Cortical Bone Thickness and Porosity Are Associated With Nontraumatic Fractures in Postmenopausal Women : A J o u r n a l P r e -p r o o f Pilot Study[END_REF]). After extensive review, a bi-directional axial transmission ultrasound device was selected for our study to investigate the ability of ultrasound to predict the mechanical properties of bone at the neck of the femur. The ultrasound has been developed and validated in-vitro (J.-G. [START_REF] Minonzio | Bone Cortical Thickness and Porosity Assessment Using Ultrasound Guided Waves: An Ex Vivo Validation Study[END_REF][START_REF] Schneider | Ex Vivo Cortical Porosity and Thickness Predictions at the Tibia Using Full-Spectrum Ultrasonic Guided-Wave Analysis[END_REF]) and in-vivo [START_REF] Ramiandrisoa | In Vivo Estimation of Cortical Thickness and Porosity by Axial Transmission: Comparison with High Resolution Computed Tomography[END_REF][START_REF] Schneider | In Vivo Measurements of Cortical Thickness and Porosity at the Proximal Third of the Tibia Using Guided Waves: Comparison with Site-Matched Peripheral Quantitative Computed Tomography and Distal High-Resolution Peripheral Quantitative Computed Tomography[END_REF].

The reproducibility and validity of Azalee ultrasound in scanning cortical thickness and porosity were assessed on 27 radii and 13 tibiae (J.-G. [START_REF] Minonzio | Bone Cortical Thickness and Porosity Assessment Using Ultrasound Guided Waves: An Ex Vivo Validation Study[END_REF]. The reproducibility of Azalee ultrasound for cortical thickness and porosity, given by the root mean square error of the standard deviation of estimates, was 0.11 mm and 1.9% respectively. The ultrasound has been utilised before to distinguish between fractured and non-fractured patients. However, the ability of ultrasound to predict the mechanical properties of bone at the clinically relevant site has not yet been investigated.

J o u r n a l P r e -p r o o f

Aims and objectives

The aim of this study was to determine the ability of Azalee guided wave ultrasound to estimate the mechanical properties of cortical bone at the neck of the femur in patients undergoing elective hip arthroplasty. We tested the hypothesis that guided wave ultrasound scans from peripheral sites correlated with the mechanical properties of bone at the neck of the femur.

J o u r n a l P r e -p r o o f

Materials and Method

Data collection

Femoral heads and necks were retrieved from 19 patients (12 male, 7 female) who had elective total hip arthroplasty for osteoarthritis. All patients were given an information sheet and a consent form, and data were collected from the patients who gave written consent to participate in the study. Ethical approval for this research was received (Imperial Tissue Bank number R13004a, Wales Research Ethics Committee number 18/LO/0277; IRAS project ID: 214772). Tibia and radius of patients were scanned using the guided wave ultrasound on the day of surgery or within two days after surgery.

Sample preparation

The bone samples collected from patients were stored at -80 °C after surgery. They were thawed in a 0.9% saline solution for one hour before cutting. It has been shown that thawing and refreezing bone several times does not change the bone properties significantly [START_REF] Linde | The Effect of Specimen Geometry on the Mechanical Behavior of Trabecular Bone Specimens[END_REF].

The femoral neck was parted from the femoral head using a water-cooled bandsaw (Exakt, Scotland). The cortical bone of the femoral neck was separated from the trabecular bone using the same bandsaw (Fig 1a). Using a high precision, low speed saw with a diamond wafering blade (Isomet, Buehler, Germany), cortical bone specimens were cut from the medial calcar region of the femoral neck, in the osteon's direction. The cortical specimen was cut into 3*3*6 mm rectangular parallelepiped samples [START_REF] Helgason | Mathematical Relationships between Bone Density and Mechanical Properties: A Literature Review[END_REF][START_REF] Keaveny | Systematic and Random Errors in Compression Testing of Trabecular Bone[END_REF][START_REF] Linde | The Effect of Specimen Geometry on the Mechanical Behavior of Trabecular Bone Specimens[END_REF]. A digital caliper was utilised to measure the dimensions with a precision of ±0.05 mm. Three measurements were taken along the length of the sample J o u r n a l P r e -p r o o f in each dimension, and the mean was used to calculate the volume. Fully hydrated samples were weighed just before mechanical testing.

Mechanical testing

Bone specimens were subjected to bulk tissue compression testing under a mechanical testing device (Instron 5565, Instron Ltd, High Wycombe, UK). The bone samples were set down between two polished, parallel, stainless-steel platens and compressive loads were applied to the specimens. The load was measured with a 5kN load cell, and the displacement was measured based on plate displacement. The samples were unconstrained, and the load was applied to the longitudinal axis of the bone specimen. Samples were moistened before testing.

Destructive testing was performed with a strain rate of 0.005/s, which is assumed to be the quasi-static strain rate for bone [START_REF] Hansen | The Effect of Strain Rate on the Mechanical Properties of Human Cortical Bone[END_REF][START_REF] Keaveny | Systematic and Random Errors in Compression Testing of Trabecular Bone[END_REF][START_REF] Linde | Mechanical Properties of Trabecular Bone by a Non-Destructive Compression Testing Approach[END_REF].

The load-displacement data was transformed to stress-strain data based on the geometries of the bone samples. The elastic modulus was calculated as the maximum slope of the stress strain curve's linear (elastic) region, using 1% strain intervals [START_REF] Keaveny | Systematic and Random Errors in Compression Testing of Trabecular Bone[END_REF][START_REF] Reed | Elastic Modulus and Strength of Emu Cortical Bone[END_REF]. The ultimate strength of bone specimens was identified as the maximum stress applied to the bone specimens before the failure point (Fig 1b). Finally, the area under the stress-strain curves were approximated using the trapezoidal method in MATLAB software, whereby the area under the curve is divided into smaller trapezoids and the area of these are summed.

Ultrasound scans

J o u r n a l P r e -p r o o f

A bi-directional Quantitative ultrasound device (Azalee, France) was utilised for scanning the patients (J.-G. [START_REF] Minonzio | Bone Cortical Thickness and Porosity Assessment Using Ultrasound Guided Waves: An Ex Vivo Validation Study[END_REF][START_REF] J O U R N A L P R E -P R O O F | In Vivo Characterization of Cortical Bone Using Guided Waves Measured by Axial Transmission[END_REF]. The ultrasound device consists of two probes, including a 1-MHz probe adapted for radius and a 0.5-MHz probe adapted for tibia (Vermont, Tours, France). Each probe consists of a linear array of piezo composite elements, divided into an array of 24 receivers, and two surrounding arrays of five transmitters each.

The probes were bidirectional, measuring the propagation of guided waves in two opposite directions. This was to correct for errors introduced by inclination angles between the probe and the bone due to uneven overlying soft tissues (Fig 2).

An electronic device was used to transmit and receive signals (Althais, Tours, France). Each transmitter was successively excited with a wideband pulse (170 V, 1-MHz for radius, and 0.5-MHz for tibia as the central frequency). The received signals were averaged 16 times and then sampled with a 20 MHz frequency.

A Graphic User Interface (HMI; Bleu Solid, Paris, France) was developed to display the spectrum of guided waves in real-time, to guide the operator in placing the probe in the optimal position, relative to the main bone axis (Fig 2). When the probe is not aligned with the bone axis, the guided wave spectrum image appears blurry. On the contrary when correctly aligned, guided modes appear as continuous lines. Two particular regions of the spectrum are visually followed: the upper part, associated with the slowest wave or A 0 mode, and the lower part, associated with the highest phase velocities and the modes traveling within the two cortical bone interfaces (J. [START_REF] Minonzio | Bi-Directional Axial Transmission Measurements Is an Easy To Apply Methodology Allowing Risk Assessment of Fracture in Elderly[END_REF].

Ultrasound scans for both radii and tibiae were obtained for each patient. The ultrasound device was placed at the anteromedial surface of the distal one-third of the radius and J o u r n a l P r e -p r o o f midpoint of the tibia. Three acquisitions were performed in each location and 10 scans were collected in each acquisition.

Signal Processing

The ultrasound signal was post-processed to extract the cortical thickness and porosity from the guided wave signal. The velocity of the first arriving signal (v FAS ) and the velocity of the A0 mode (v A0 ) were also calculated. This was done using software previously researched and developed (J.-G. [START_REF] Minonzio | Ultrasound -Based Estimates of Cortical Bone Thickness and Porosity Are Associated With Nontraumatic Fractures in Postmenopausal Women : A J o u r n a l P r e -p r o o f Pilot Study[END_REF].

Initially, the signal processing was designed to obtain the guided wave spectrum image, or Norm (f, k), using the so called SVD based method, interpreted as an enhanced spatiotemporal Fourier transform (J.-G. [START_REF] Minonzio | Bone Cortical Thickness and Porosity Assessment Using Ultrasound Guided Waves: An Ex Vivo Validation Study[END_REF]. This method was expanded to the inverse problem by projecting waveguide models associated with series of guided mode theorical curves in the f-k plane. Finally, each pixel of the inverse problem image in the Ct.Th-Ct.Po corresponds to a particular model, obtained with 2D transverse isotropic free plate, with porosity parametrised elasticity. The best fitting model corresponds to the maximum of the inverse problem image, expressed in a 0-1 scale due to normalised projection (Fig 3).

The ultrasound results were then filtered, firstly by discarding any individual scans where either velocity parameter (v FAS or v A0 ) was unobtainable, then by discarding any acquisitions where the remaining number of scans was less than 3, and finally by discarding acquisitions whose measurements had a standard deviation in thickness greater than 0.3 mm or a standard deviation in porosity greater than 2%. The final value for each site was calculated based on the average of the remaining measurements values of the same site.

Statistical Analysis

J o u r n a l P r e -p r o o f

Linear regression analysis was performed for each of the ultrasound parameters and each of the mechanical testing parameters, using the least square fit. Pearson correlation coefficients were calculated to show the degree of association between the two parameters. Correlations were considered statistically significant for p < 0.05. 'Match side' ultrasound scans, i.e., measurements taken from the same side the bone specimen was collected from, were used in the correlations. Correlations for each of the 4 scanned sites (radii, tibia, left, right) and for the mean between left and right sites were also calculated.

Multiple linear regression models were used to investigate the association between the mechanical testing properties and the ultrasound measurements at both radial and tibial sites. The association between mechanical properties and multiple ultrasound parameters at a single site were also investigated. 36 multiple linear regression models were created in MATLAB, using combinations of radial and tibial measurements, as well as combinations of different ultrasound parameters at a single site. The coefficient of determination (R 2 ) and statistical significance (p) were calculated for each model, and models were considered statistically significant for p < 0.05. The models were used to predict the mechanical properties, and the predicted values were compared to the values obtained experimentally.

Wilcoxon signed-rank tests were performed on all ultrasound parameters to identify any significant differences between scans obtained from the left and right sides.

J o u r n a l P r e -p r o o f

Results

From 19 patients, Azalee was only able to collect a complete set of ultrasound scans (i.e., all 4 sites) from 15 patients. For the remaining 4 patients, scans could not be obtained for one or more sites due to the high amount of soft tissue, resulting in a total of 7 missing data points (i.e., 4 radii and 3 tibiae).

The correlations for guided waves versus mechanical properties are shown in table 2. Most of the correlations were found to be statistically insignificant (p > 0.05). Three statistically significant correlations, shown in figures 4 a, b, and c (R 2 = 0.34, p < 0.05, R 2 = 0.26, p < 0.05, and R 2 = 0.31, p < 0.05, respectively), demonstrate a potential link between the porosity, as estimated by the ultrasound scans, and the elasticity, density, and stress (as indicated by the area under the stress-strain curve) of the femur, however, these correlations were not consistent across all sites.

Two of the 36 multiple linear regression models were found to be statistically significant, the details of which can be found in table 3. Negative linear associations were found between density at the femoral neck and the porosity and v FAS ultrasound measurements at the tibia, and the elasticity at the femoral neck and the porosity ultrasound measurements at both the radius and the tibia. The comparison of the predicted values with the actual values for these models can be found in figure 5a andb (R 2 = 0.48, p < 0.003, and R 2 = 0.41, respectively). This shows that combining multiple ultrasound measurements may have the potential to better predict properties at the femoral neck.

Left and Right Comparison

The Wilcoxon signed-rank test revealed there was no statistically significant difference between the left and right sides for any ultrasound scans.

J o u r n a l P r e -p r o o f

Discussion

Orthopaedic Surgeons need to characterise the mechanical properties of bone at the femoral neck to make more informed clinical decisions in arthroplasty surgeries and achieve better patient outcomes. This study aimed to test the efficacy of the Azalee guided wave ultrasound in estimating the mechanical properties of the cortical bone in the neck of the femur. The aim was achieved by correlating peripheral ultrasound measurements from the tibia and radius with mechanical data from the clinically relevant femoral neck including the elastic modulus, ultimate strength, and density of the cortical bone. The results show weak to moderate correlations between the ultrasound measurements at the peripheral sites and the mechanical properties of cortical bone at the neck of the femur (Table 1). Further, correlation coefficients (R 2 ) and the level of statistical significance (p) between Azalee ultrasound measurements and the mechanical properties were not consistent across the four measurement sites in a single patient.

Mechanical properties at the hip

Although there were some weak to moderate correlations and weak to moderate predictions calculated with multiple regression analysis, we can conclude Azalee ultrasound scans at the peripheral sites cannot estimate the mechanical properties at the osteoarthritic neck of the femur, due to the strength and inconsistency of the correlations.

One explanation for weak and inconsistent correlations found across the neck of femur and peripheral sites could be regional variation in bone properties around the body. The regional variation could be attributed to several different causes. Firstly, localised bone remodelling as a consequence of osteoarthritis may have caused regional changes in bone properties.

Osteoarthritis is typically associated with wear and tear of articular cartilage, however, there J o u r n a l P r e -p r o o f is evidence to suggest that localised bone remodelling also plays a role [START_REF] Donell | Subchondral Bone Remodelling in Osteoarthritis[END_REF]. It is well documented that there are changes in subchondral bone, which affect both the subchondral cortical bone and the subchondral cancellous bone [START_REF] Zhu | Subchondral Bone Remodeling: A Therapeutic Target for Osteoarthritis[END_REF], and these changes can be found near to and even distal from the joint [START_REF] Auger | Spatial Assessment of Femoral Neck Bone Density and Microstructure in Hip Osteoarthritis[END_REF]. Therefore, the bone properties may exhibit distinct characteristics at the neck of femur compared to the peripheral sites, which may have given rise to weak correlations. However, the full spatial picture of remodelling is not well documented. Therefore, we cannot be sure of the extent to which this bone remodelling affected the cortical bone samples, and thus the resulting correlations.

Other factors may also contribute to regional variation in bone properties across the body, such as different loads experienced across a patient's lifetime. There is evidence to show that site-specific bone remodelling occurs as a response to strains experienced by the bone [START_REF] Carter | Mechanical Loading Histories and Cortical Bone Remodeling[END_REF]. Therefore, over time bone properties across different anatomical regions may vary with different loads experienced, potentially leading to the differences in the mechanical properties observed between the clinically relevant and the peripheral sites. One study has shown that the bone mineral density, as measured by DEXA, from the femur and distal radius showed significant differences between the two sites [START_REF] Wilson | Variation in Bone Mineral Density by Anatomical Site in Patients with Proximal Humeral Fractures[END_REF]). Additionally, handedness has been shown to cause asymmetries in the thickness and BMD of bone between the left and right sides of the body due to habitual loads [START_REF] Troy | Dominant and Nondominant Distal Radius Microstructure: Predictors of Asymmetry and Effects of a Unilateral Mechanical Loading Intervention[END_REF] This could explain the internal inconsistencies observed between correlations of the dominant and nondominant sides.

Verification against literature

J o u r n a l P r e -p r o o f

Experimental results in this study were verified by comparing against values found in the literature. The density measurements in this study varied from 1.4 kg/m 3 to 2.0 kg/m 3 (Table 1). This falls within the range of values reported for femoral cortical bone density in the literature [START_REF] Rho | An Ultrasonic Method for Measuring the Elastic Properties of Human Tibial Cortical and Cancellous Bone[END_REF]. The compressive elastic modulus of cortical bone at the neck of the femur calculated in this study was between 0.9 GPa and 6.3 GPa, which is in line with the reported values in the literature with a similar testing protocol [START_REF] Boughton | Measuring Bone Stiffness Using Spherical Indentation[END_REF].

The model that was utilised for estimating the thickness and porosity of the cortical bone was previously validated (J.-G. [START_REF] Minonzio | Bone Cortical Thickness and Porosity Assessment Using Ultrasound Guided Waves: An Ex Vivo Validation Study[END_REF]. In this study, the thickness of the cortical bone at a point on the distal one-third of the radius and at the midpoint of the tibia was estimated as 3.17±0.24 mm and 5.23±0.89 mm, respectively. The porosity of the cortical bone at the one-third distal radius and at half tibia was estimated as 9.3±2.3% and 11.7±2.3% respectively (Table 1). These values were all in good agreement with previously reported values (J.-G. [START_REF] Minonzio | Bone Cortical Thickness and Porosity Assessment Using Ultrasound Guided Waves: An Ex Vivo Validation Study[END_REF][START_REF] Schneider | In Vivo Measurements of Cortical Thickness and Porosity at the Proximal Third of the Tibia Using Guided Waves: Comparison with Site-Matched Peripheral Quantitative Computed Tomography and Distal High-Resolution Peripheral Quantitative Computed Tomography[END_REF][START_REF] J O U R N A L P R E -P R O O F | In Vivo Characterization of Cortical Bone Using Guided Waves Measured by Axial Transmission[END_REF].

We can conclude that all values found from the mechanical testing and ultrasound data were in valid ranges as prescribed in the literature.

Limitations

The displacement of bone in mechanical testing was measured using the displacement of the plates. However, utilising techniques that directly measures the displacement of bone such as laser would be more accurate. Moreover, the ultrasound signal processing algorithm relies heavily on several assumptions which could account for some of the variation reported in the ultrasound estimates. Firstly, the shape of the bone and waveguide physics were simplified into a 2D isotropic plate model (J.-G. [START_REF] Minonzio | Bone Cortical Thickness and Porosity Assessment Using Ultrasound Guided Waves: An Ex Vivo Validation Study[END_REF]. Secondly, the waveguide was simplified such that it was parameterised only by the thickness and porosity, by assuming the J o u r n a l P r e -p r o o f bone microstructure consisted of equally sized and equally distributed pores. Due to these assumptions, non-uniformity and irregularities within the bone may cause scattering or disruption in the guided wave propagation. Additionally, the device did not work as well on patients with a high BMI, due to a thicker layer of surrounding soft tissue. Soft tissue can attenuate the signal or introduce additional unwanted modes, disrupting the signal.

Future directions

The ultrasound measurements at peripheral sites cannot predict the mechanical properties at the neck of femur in osteoarthritic patients. Previous research has also found that peripheral measurements using QUS and QCT cannot predict fracture cases at clinically relevant sites [START_REF] Clowes | The Discriminative Ability of Peripheral and Axial Bone Measurements to Identify Proximal Femoral, Vertebral, Distal Forearm and Proximal Humeral Fractures: A Case Control Study[END_REF]. This shows the need for a device which can directly measure and characterise properties at the clinically relevant site. In-vivo characterisation at the neck of the femur would be more accurate as we scan the region of interest, and it would give surgeons the ability to characterise before the surgery; and therefore, plan the surgery in normal clinical consultations. This may reduce the time of the surgery as well.

There have been some attempts to characterise the bone at the femur, using ultrasound however, the method was focused on trabecular bone. Furthermore, the method required scans to be performed in a water bath which would be difficult to perform and establish in clinics [START_REF] Andre | Measurement of the Velocity of Ultrasound in the Human Femur in Vivo[END_REF][START_REF] Barkmann | A Device for in Vivo Measurements of Quantitative Ultrasound Variables at the Human Proximal Femur[END_REF]. There is a need to refer to physics principles to address the challenges of ultrasound scanning of the femoral neck, due to the complexity of its geometry and the thickness of the surrounding soft tissue.

To make more informed decisions in hip arthroplasty, assessment of the mechanical properties of bone at the neck of the femur can also be performed in-vitro or in-situ. However, J o u r n a l P r e -p r o o f the mechanical characterisation of bone must be intraoperative for in-vitro and in-situ methods. In-vitro ultrasound scans of cortical bone at the neck of the femur correlated strongly with the mechanical properties of cortical bone (Boughton et al. 2019;[START_REF] Fan | Cortical Bone Viscoelastic Damping Assessed with Resonant Ultrasound Spectroscopy Reflects Porosity and Mineral Content[END_REF].

Furthermore, the speed of the sound (SOS) of the first arrival signal in cadaveric bone showed a significant difference between fractured and non-fractured patients [START_REF] J O U R N A L P R E -P R O O F Fu | Fragility Fracture Discriminative Ability of Radius Quantitative Ultrasound: A Systematic Review and Meta-Analysis[END_REF].

However, the likelihood of surgical site infection (SSI), increase with increasing operative time [START_REF] Cheng | Prolonged Operative Duration Increases Risk of Surgical Site Infections: A Systematic Review[END_REF]. In-vitro ultrasound scanning of bone due to the preparation of the bone is lengthy and therefore the method is not practical intraoperatively.

In-situ testing of cortical bone at the neck of femur could be a potential way forward. In-situ testing of bone during surgery can enable us to measure scan bone directly without any surrounding soft tissue. In-situ testing also increases the time of surgery, though sample preparation is not needed in this process which makes it quicker compared to in-vitro testing.

While Azalee ultrasound cannot be used for estimating mechanical properties of bone at the neck of the femur in osteoarthritic patients indirectly, there is still potential for further studies to investigate the use of the device in monitoring systematic changes in bone properties at peripheral sites over time [START_REF] Horváth | FRI0547 Osteoporosis and Bone Metabolism in Systemic Sclerosis[END_REF]. Osteoporosis is a systemic disease, in which the bone composition changes over time (J. [START_REF] Minonzio | Bi-Directional Axial Transmission Measurements Is an Easy To Apply Methodology Allowing Risk Assessment of Fracture in Elderly[END_REF][START_REF] Minonzio | Ultrasound -Based Estimates of Cortical Bone Thickness and Porosity Are Associated With Nontraumatic Fractures in Postmenopausal Women : A J o u r n a l P r e -p r o o f Pilot Study[END_REF].

Therefore, Azalee may have the potential to be used in monitoring elderly populations who are vulnerable to osteoporosis or estimating changes in bone health under osteoporosis medications.

Conclusion

Azalee ultrasound scans of the peripheral sites were not estimated of the mechanical properties at the neck of the femur and were not internally consistent across various J o u r n a l P r e -p r o o f peripheral sites (e.g., upper versus lower limb or left versus right). Persisting with ultrasound scans at the peripheral sites, rather than directly at clinically relevant sites like the neck of femur, does not seem to be the way forward. There is, however, potential for a new device based on a better knowledge of ultrasound physics in the bone to cope with a thicker soft tissue to make direct ultrasound measurements.

Figures and Tables:

Figure 1 :

 1 Figure 1: Bone sample preparation and mechanical testing. a) sectioning b) compressive mechanical testing, c) plotting typical stress-strain graph, and d) a stress-strain graph of one of the bone samples.

Figure 2 :

 2 Figure 2: a) ultrasound scanning the radius and tibia of one of the participants and, b) the resulting ultrasound signal.

Figure 3 :

 3 Figure 3: Contour plots were used to depict the solution to the inverse problem, where the global maximum represents the value of thickness (Ct.Th) and porosity (Ct.Po) which best fits the measured with the theoretical dispersion curves. In this acquisition, the cortical thickness (Ct.Th) was estimated at 2.82±0.01mm and the cortical porosity (Ct.Po) was estimated at 9.1±0.2%.

Figure 4 :

 4 Figure 4: Bone properties at the femoral neck including (a) density and (b) elasticity are moderately negatively correlated with the porosity of the tibia and radius, respectively.

Figure 5 :

 5 Figure 5: Predictions from multiple linear regression models using (a) the match side porosity and match-side v fas at the tibia to predict the density of the femoral neck, and (b)the match side porosity of the radius the tibia to predict the elasticity.

Figure 5 :

 5 Figure 5: Predictions from multiple linear regression models using (a) the match side porosity and match side v fas at the tibia to predict the density of the femoral neck, and (b)the match side porosity of the radius the tibia to predict the elasticity.

  r n a l P r e -p r o o f

Figure 5 :

 5 Figure 5: Predictions from multiple linear regression models using (a) the match side porosity and match side vfas at the tibia to predict the density of the femoral neck, and (b) the match side porosity of the radius the tibia to predict the elasticity

Table 3 :

 3 Multiple linear regression models and predictive equations. Abbreviations of ultrasound variables: PT = porosity measured at tibia, PR = porosity measured at radius, and VT = v FAS measured at tibia

	Mechanical Testing Variable	Ultrasound Variables	R 2	p value	Equation	RMSE	NRMSE (%)
	Density	PT, VT	0.48	0.01	Density = 7526 -45 * PorosityTibia -1.3 * VfasTibia	137.0	8.0
	Elasticity	PR, PT	0.41	0.04	Elasticity = 9939 -319 * PorosityRadius -274 * PorosityTibia	1323.6	36.7

Acknowledgment:

This study was funded by Wellcome Trust and Sackler Charity fund. Jean-Gabriel Minonzio is supported by Grant ANID Fondecyt Regular 1201311 and ANID ECOS 200061. J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f

J o u r n a l P r e -p r o o f

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: