Jérémie Labroquère
email: jeremie.labroquere@inria.fr

Régis Duvigneau
email: regis.duvigneau@inria.fr

OPTIMIZATION METHODS APPLIED TO AERODYNAMIC FLOW CONTROL

Keywords: flow control, optimization, aerodynamics

This study deals with the use of optimization algorithms to determine efficient parameters of flow control devices. To improve the performance of systems characterized by detached flows and vortex shedding, the use of flow control devices such as oscillatory jets, are intensively studied, using numerical as well as experimental methods. However, the determination of efficient control parameters is still a bottleneck for industrial problems. Therefore, we propose to couple a global optimization algorithm with an unsteady flow simulation to derive efficient flow control rules.

We consider as testcase the turbulent flow over a backward facing step, including a synthetic jet actuator. The aim is to reduce the time-averaged recirculation length behind the step by optimizing the jet blowing/suction amplitude and frequency. The Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations are solved within a Mixed finite-Element/finite-Volume (MEV) framework using the near-wall low-Reynolds number one-equation Spalart-Allmaras turbulence closure. The steady flow simulation without control is first validated by comparison with experimental and numerical data. Then, the optimization method EGO (Efficient Global Optimization), based on the construction of a Gaussian surrogate model, is coupled with the solver and applied to the unsteady flow with actuation. It is shown that the timeaveraged recirculation length can be shortened when suitable control parameters are used.

INTRODUCTION

Shape optimization methods have significantly matured for the last years and automated design loops based on CFD software are now commonly used in industry. However, shape optimization may not be efficient to improve the performance, when the flow is characterized by a strong unsteadiness related to a massive detachement. This is typically the case for the flow around a car or a wing in stall condition. To overcome this difficulty, flow control strategies have been recently developed, that aim at manipulating vortex dynamics by introducing some active actuators, such as periodic blowing/suction jets. In this context, the choice of the control parameters (location, amplitude, frequency) is critical and not straightforward. The design of flow control devices is usually carried out using systematic parametric studies, either on the basis of numerical simulations or experimental campaigns.

Therefore, the objective of this work is to extend some methodologies, developed successfully for shape optimization, to improve the efficiency of flow control devices. The main difficulty is related to the fact that the flows are strongly unsteady. As consequences, the computational cost of each simulation is huge, the objective functional can be highly multimodal and the computation of functional gradients is cumbersome using an adjoint approach.

To overcome these difficulties, a deterministic global optimization method based on surrogate models is used in this study: we extend the approach proposed in [START_REF] Duvigneau | Kriging-based optimization applied to flow control[END_REF] for laminar flows, to an unsteady turbulent flow control problem.

NUMERICAL FRAMEWORK

We describe in this section the numerical framework used to solve compressible turbulent flow equations. The multidisciplinary platform NUM3SIS [START_REF] Labroquère | Interactive computation and visualization towards a virtual wind tunnel[END_REF], developed at INRIA, is used for the implementation. We present the numerical methods employed, and more particularly the implementation of the synthetic jet boundary conditions.

Mean turbulent flow modeling

The equation solved for the turbulent compressible flow simulations are the Unsteady Favre-Averaged Navier-Stokes equations. For any variable x, with x the mean part and x the fluctuations of x, the Reynolds average or filter is denoted by x = x + x with x = 0. The density weighted average, called Favre averaging or filtering, is denoted by x = ρx/ρ with the related fluctuation part x = x -x.

In conservative form, the averaged turbulent compressible equations can be written as:

∂ t W c + ∇ • F (W c ) = ∇ • N (W c ) + S(W c ) (1)
with the averaged conservative variables W c = (ρ, ρ u, ρ v, ρ w, ρ E) T , the inviscid flux

F (W c ) = (F (W c ), G(W c ), H(W c )) T , the viscous flux N (W c ) = (R(W c ), S(W c ), T (W c ))
T and the source term S(W c ).

The components of the inviscid flux in the global frame R 0 ( x, y, z) are:

F (W c ) =       ρ u ρ u 2 + p ρ u v ρ u w u(ρ E + p)       , G(W c ) =       ρ v ρ v u ρ v 2 + p ρ v w v(ρ E + p)       H(W c ) =       ρ w ρ w u ρ w v ρ w 2 + p w(ρ E + p)      
The components of the viscous flux are:

R(W c ) =       0 τ xx τ xy τ xz uτ xx + vτ xy + wτ xz -q x       , S(W c ) =       0 τ yx τ yy τ yz uτ yx + vτ yy + wτ yz -q y       T (W c ) =       0 τ zx τ zy τ zz uτ zx + vτ zy + wτ zz -q z      
By averaging the Navier-Stokes equations, other quantities as the molecular diffusion

τ lam • V , the turbulent transport term ρ(V • V )V and the turbulent kinetic energy k = 1 2 ρV •V ρ
appears but are not shown here. Theses quantities are neglected in this study.

The averaged shear stress tensor τ is made of a laminar and turbulent tensor: τ = τ lam + τ turb . The laminar part is modeled considering the flow as beeing Newtonian:

τ lam = µ l ∇ V + ∇ V T - 2 3 (∇ • V )I n
The turbulent part uses the Boussinesq eddy-viscosity assumption :

τ turb = -ρV ⊗ V ≈ µ t ∇ V + ∇ V T - 2 3 (∇ • V )I n
As for the shear stress tensor, the heat flux q is split in two parts: q = q lam + q turb . Both the laminar and turbulent parts are modeled by using a Fick or Fourier law:

q lam ( , κ) = -γ µ l Pr ∇ q turb ( , κ) = ρV h ≈ -γ µ t Pr t ∇
with = p (γ-1)ρ and the specific enthalpy h = c p T where c p is the specific coefficient at constant pressure and T the temperature.

Finally, to close the equations, the total energy is modeled by the perfect gas state law E = + 1 2 V • V , the adiabatic index is set to γ = 1.4, the Prandtl number is set to Pr = 0.72, the turbulent Prandtl number is set to Pr t = 0.9 and the laminar viscosity µ l is modeled by the Sutherland law. The turbulent viscosity µ t is found by solving the Spalart-Allmaras equation.

For simplicity purpose, the average notations . and . are ommited in the following.

Spalart-Allmaras turbulent viscosity modeling

The turbulence model used in this study is the basic version of the one-equation Spalart-Allmaras model, where the turbulent viscosity is a function of the variable ν:

µ t = f v1 ρmax(ν, 0) with f v1 = χ 3 χ 3 +C 3 v1 , χ = ν µ l /ρ
and C v1 a model constant defined below. As for the mean equations, the Spalart-Allmaras equation in conservative form can be written as:

∂ t W c t + ∇ • F t (W c , W c t ) = ∇ • N t (W c , W c t ) + S t (W c , W c t ) (2)
with the Spalart-Allamars conservative variable W c t = (ρν), the inviscid turbulence flux F t (W c , W c t ) = (ρνu, ρνv, ρνw) T , the viscous flux N t (W c , W c t ) = µ t e ∇ν with µ t e = µ l +ρν σm , σ m a model constant and the source flux S t (W c , W c t ). The components of the turbulence source flux in the global frame R 0 ( x, y, z) are:

S t (W c , W c t ) = P -D + d
with P, D and d respectively the production, destruction and diffusion terms defined by:

P = c b1 ρν S, D = c w1 f w ρ ν d wall 2 , d = c b2 σ m ρ ∇ νt 2 2
with S = w 2 + ν κ 2 d 2 wall f v2 limited to be positive, d wall the distance to the wall and the vorticity w defined as w = ∇ × V .

The missing model variables are set to:

f v1 = χ 3 χ 3 + C 3 v1 , f v2 = 1 - χ 1 + χf v1 , f w = g 1 + c 6 w3 g 6 + c 6 w3 1/6 g = r + c w2 (r 6 -r), r = min ν Sκ 2 d 2 wall , 10 , χ = ν µ l /ρ κ = 0.41, c b1 = 0.1355, σ m = 2/3, c b2 = 0.622, κ = 0.41, c w2 = 0.3, c w3 = 2, c v1 = 7.1, c w1 = c b1 κ 2 + 1+c b2 σm .

Spatial discretization

A Mixed finite-Element/finite-Volume (MEV) discretization is used with a vertex centered approach. The inviscid fluxes and source fluxes are discretized with a finite-volume approach while the viscous fluxes are discretized with a finite-element approach [START_REF] Fezoui | Resolution numerique des equations de Navier-Stokes pour un fluide compressible en maillage triangulaire[END_REF][START_REF] Francescatto | Méthodes multigrilles par agglomération dirrectionnelle pour le calcul d'écoulements turbulents[END_REF].

A polygonal bounded domain Ω ⊂ R n is considered with a boundary Γ, sub-divided into a tetrahedrization or triangulation T h with elements T i . Around each vertex s i a finite-volume control cell C i of a measure m(C i ) is constructed. The set of vertices which are joined to the vertex s i is denoted by N (s i ). The subset of all the highest topological dimension polygons sharing the vertex s i is denoted by T (s i ).

The inviscid fluxes are computed on the dual control cells C i while the viscous fluxes are computed on the elements T i . A weak formulation of the Navier-Stokes equations can be expressed using a Galerkin approach.

By integrating the equation (1) over a control cell S i (dual control cell or element) against a test function ϕ i , the weak formulation can be written as:

S i (∂ t W c + ∇ • F (W c ))ϕ i dΩ = S i (∇ • N (W c ))ϕ i dΩ + S i (S(W c ))ϕ i dΩ (3)

Convective fluxes discretization

The finite-volume method can be interpreted as a Galerking method using the control cell S i = C i and a test function defined as:

ϕ C i (x) = 1 if x is in C i 0 else
The variables W c are supposed to be constant on each control cell C i and denoted W c i (see Fig. 1(a By using the Green-Ostrogradski theorem, the left-hand side of equation ( 3) yields:

C i (∂ t W c + ∇ • F (W c ))ϕ C i dΩ = m(C i )∂ t W c i + ∂C i (F (W c ) • η))dσ (4)
with ∂C i the boundary of the cell C i and η the outward unit normal vector. Furthermore, as the domain is discretized,

∂C i (F (W c ) • η))dσ = s j ∈N (s i ) ∂C i ∩∂C j (F (W c ) • η))dσ ≈ s j ∈N (s i ) F (W c ) |ij • ∂C i ∩∂C j ηdσ = s j ∈N (s i ) F (W c ) |ij • η ij
The F (W c ) |ij term is estimated by using an approximate Riemann solver for which the associated numerical flux is

Φ ij = Φ(W c i , W c j , η ij ).
In the current work, for the mean flow equations, the HLLC [START_REF] Ef Toro | Restoration of the contact surface in the hll-riemann solver[END_REF] numerical flux is used for the internal domain and a modified version of the Steger-Warming flux is used for boundaries [START_REF] Fezoui | Resolution numerique des equations de Navier-Stokes pour un fluide compressible en maillage triangulaire[END_REF]. For the estimation of the inviscid flux of the turbulence equation, the density flux computed for the mean flow equations is used [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF].

To obtain a high-order approximation in space, a MUSCL reconstruction technique is used. The reconstructed primitive state at the common interface of the cells C i and C j is denoted by

W p ij = W p i + 1 2 α ij (∇W p i ) • IJ . The slope (∇W p i ) • IJ is approximated by ∆W p i = 2/3∆ |N W p i + 1/3∆ |C W p ij .
The nodal slope ∆ |N W p i is the weighted average of the P 1 -Galerkin gradients computed on each element T ∈ T (s i ). The slope ∆ |C W p ij corresponds to the centered slope W p j -W p i . We denote by α ij (∆ |N W p i , ∆ |C W p ij ) the limiting coefficient. Finally, the high-order numerical flux is computed using extrapolated states

Φ ij High order = Φ(W c ij , W c ji , η ij ).

Viscous terms discretization

The finite-element method is retrieved by using the control cell S i = T i . The polygon T i is defined as a Lagrange P 1 finite-element with the canonical basis functions denoted ϕ T i for each vertex s i . Thus, ϕ T i is equal to unity at the vertex s i and zero at the other vertices of the element. The P 1 approximation of any function f on a polygon T i can be obtained by :

f (X) ≈ f T i h (X) = Ns∈T i i=1 f i ϕ T i (X)
where f i is the value of f at the vertex s i .

By integrating by part and using the Green-Ostrogradski theorem, the diffusive term of equation (3) becomes:

T i (∇ • N (W c ))ϕ T i dΩ = - T i N (W c ) • ∇ϕ T i dΩ + ∂T i (N (W c ) • n)ϕ T i dσ (5)
The velocity gradient N and ∇ϕ T i being assumed to be constant by element :

T i N (W c ) • ∇ϕ T i dΩ = N (T i ) • ∇ϕ T i T i dΩ = N (T i ) • ∇ϕ T i m(T i )
By denoting η T i = -∇ϕ T i m(T i ), the viscous numerical flux is finally defined as:

Υ T i = N (T i ) • η T i
Thus, the equation ( 5) reads:

T i (∇ • N (W c ))ϕ T i dΩ = Υ T i + ∂T i (N (W c ) • n)ϕ T i dσ (6)
The term ∂T i (N (W c ) • n)ϕ T i dσ concerns the boundary conditions.

Source terms discretization

The source terms are discretized in a finite volume way, by considering the variables are constant on the control cell C i . The source term of 3 becomes:

C i (S(W c ))ϕ C i dΩ = m(C i )S(W c i )
If the source terms require gradient values, a weighted average of the P 1 -Galerkin gradients is used, as for the inviscid fluxes MUSCL reconstruction.

Time integration

An implicit second-order time discretization is obtained by using a dual time step approach [START_REF] Jameson | Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings[END_REF] and a backward time integration. We introduce the computed residual R i for the control cell i. By denoting δ 2 1 Λ = Λ 2 -Λ 1 the variation of the variable Λ from the state 1 to the state 2, we obtain as second-order implicit time integration scheme:

m(C i ) 3W c i n+1 -4W c i n + W c i n-1 2δ n+1 n t + R i (W cn+1 ) = 0 (7)
To solve this problem, we introduce a sub-iteration state of index k, such as:

W c i n+1 = lim k→∞ W c i (k+1)
The linearization of equation ( 7) around the state of index k, with a local time step δ k+1 k t i yields:

m(C i ) δ k+1 k t i + 3m(C i ) 2δ n+1 n t I n + J * (W c(k) ) δ k+1 k W c i = -R i (W c(k) ) + m(C i ) δ n n-1 W c i -3δ k n W c i 2δ n+1 n t with J * (W c(k)
) an approximate Jacobian of the numerical residuals. The approximate Jacobian is composed of three parts arising from the inviscid fluxes, the viscous terms and the source terms. The Jacobian of the inviscid fluxes is based on the firstorder Rusanov flux. The use of the Rusanov flux, or spectral radius Jacobian approximation, is usually used in matrix-free approaches [START_REF] Kloczko | Développement d'une méthode implicite sans matrice pour la simulation 2D-3D des écoulements compressibles et faiblement compressibles en maillages nonstructurés[END_REF]. The Jacobian of the viscous terms is based on an exact linearization, as in [START_REF] Fezoui | Resolution numerique des equations de Navier-Stokes pour un fluide compressible en maillage triangulaire[END_REF]. Concerning the source terms, for the turbulence variables, the full source Jacobian is computed.

The resulting Jacobian matrix is inversed with the Lower-Upper Symmetric Gauss-Seidel (LU-SGS) iterative algorithm. The use of a delta form for the implicit formulation allows to employ an approximated Jacobian without loosing the second order accuracy in space reached with the MUSCL extrapolation.

The time integration of the mean flow equations and the Spalart-Allmaras equation is carried out using a weak coupling. Thus, the turbulence variable is frozen during the resolution of the system for mean flow variables, whereas the system for the turbulence variable is solved with frozen mean flow variables.

W c b = W c b (W c i , W c o )
. In this study, the computation of the fictive state is carried out using imposed density and velocity with extrapolation of the pressure for subsonic inlets, imposed pressure with extrapolation of characteristics for subsonic outlets (characteristic-based methods), and symmetry conditions for slipping walls. Furthermore, for all the inlet/outlet boundaries, a non-reflexive version of the Steger-Warming flux is used [START_REF] Fezoui | Resolution numerique des equations de Navier-Stokes pour un fluide compressible en maillage triangulaire[END_REF], allowing to correctly filter the entering or exiting characteristics.

Strong approach for non-slip condition

The strong boundary condition approach consists in imposing the condition at the boundary and modifying the Jacobian and the residuals in such a way that the boundary condition is verified implicitly when solving the linear system. This boundary condition is used when the value of one particular conservative variable is completely defined at the boundary. The interior state W c i is set to be equal to the fictive state W c b . The resulting flux computed at the boundary is a physical flux as shown on the Fig. 3. Thus, the non-slip boundary condition is imposed by setting at the boundary a zero momentum for the mean flow equations and a zero turbulent variable for the Spalart-Allmaras equation.

Iterative approach for jet condition

The iterative boundary conditions are implemented in the cases where it is cumbersome to use strong boundary conditions. This can occur for instance when the user would like to impose a primitive variable while the Jacobian and residuals are expressed using conservative variables. We use the dual time step approach to implement iteratively the boundary condition, in order to reach a target W c i target at the step n + 1, at the boundary :

m(C i ) δ n+1 n t i δ n+1 n W c i + R i n+1 (W c i n ) = 0 lim k→∞ δ n,k+1 n W c i = δ n+1 n W c i = δ target n W c i
In other words, this boundary condition consists in imposing at each dual time step subiterations the targeted value of the variable. As described in [START_REF] Ashcroft | Numerical modelling of wake-jet interaction with application to active noise control in turbomachinery[END_REF], the internal fields and the boundary conditions converge towards the solution with the targeted boundary condition at the next physical step.

The jet boundary conditions are implemented using a mixed iterative/weak approach. In particular, the velocity profile is set with the iterative boundary condition while the other variables are imposed using a weak form.

The jet velocity profile is modeled by:

V j = A(x) sin(ωt + φ) d j
with A(x) a profile function, ω the angular frequency of the jet, φ the phase and d j the direction of the jet. In this study, A(x) is a sine squared function. The other variables are extrapolated from the interior of the domain.

OPTIMISATION METHODS

Kriging model of the flow response

We use a kriging model to describe the variation of the objective function value, for instance the drag coefficient or the recalculation length, with respect to control parameters. Kriging models (also called Gaussian process models) belong to response surface methods, that allow to predict a function value f at a given point x, on the basis of a set of known function values

F N = {f 1 , f 2 , . . . , f N } at some points X N = {x 1 , x 2 , . . . , x N } ⊂ R d , that are stored in a database.
Kriging models [START_REF] Mackay | Introduction to gaussian processes[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF] treat the response of some experiments as if it were a realization of a stochastic process, because the model relies on a finite number of observations and is therefore subject to uncertainties. In the following, we adopt the Bayesian viewpoint of Gaussian processes [START_REF] Williams | Prediction with gaussian processes: From linear regression to linear prediction and beyond[END_REF]. The vector of known function values F N is assumed to be one realization of a multivariate Gaussian process with joint probability density:

p(F N ) = exp -1 2 F N C -1 N F N (2π) N det(C N ) , (8)
where C N is the N × N covariance matrix. The element C mn of the covariance matrix gives the correlation between the function values f m and f n obtained respectively at points x m and

x n . We assume that these values are correlated, since they correspond to underlying physical phenomena. This is expressed in terms of a correlation function c, i.e., C mn = c(x m , x n ). Now, we suppose that we would like to evaluate the function value at a new point x N +1 . When adding a new point x N +1 , the resulting vector of function values F N +1 is assumed to be a realization of the (N + 1)-variable Gaussian process with joint probability density:

p(F N +1 ) = exp -1 2 F N +1 C -1 N +1 F N +1 (2π) N +1 det(C N +1 ) . ( 9
)
Using the rule of conditional probabilities p(A|B) = p(A, B)/p(B), we can write the probability density for the unknown function value f N +1 , given the data (F N ) as:

p(f N +1 |F N ) = p(F N +1 ) p(F N ) . ( 10
)
In order to simplify this expression we notice that the (N + 1)-variable covariance matrix C N +1 can be written as

C N +1 = C N k k κ,
where

k = [c(x 1 , x N +1 ), c(x 2 , x N +1 ), . . . , c(x N , x N +1 )] and κ = c(x N +1 , x N +1 ).
This allows us to write the inverse of correlation matrix as

C -1 N +1 = M m m µ, where M = C -1 N + 1 µ mm , m = -µC -1 N k, µ = (κ -k C -1 N k) -1 .
Then, we obtain that the probability density for the function value at the new point is

p(f N +1 |F N ) ∝ exp - (f N +1 -fN+1 ) 2 2σ 2 f N +1
,

where fN+1 = k C -1 N F N σ 2 f N +1 = κ -k C -1 N k. ( 11
)
Thus, the probability density for the function value at the new point x N +1 is also Gaussian with mean fN+1 and standard deviation σ f N +1 . Therefore, the most likely value at the new point x N +1 is fN+1 . This value will be considered as the prediction of the kriging model. The variance σ f N +1 can be interpreted as a measure of uncertainty in the value prediction. The function value can be expected to vary in some range like [ f -3σ, f + 3σ].

The covariance function must reflect the characteristics of the output of the computer code since, it represents the way the different point values are correlated. In the absence of any knowledge regarding the unknown function, the most commonly used correlation function is an exponential; in this work we take the correlation function of the form:

c(x, y) = θ 1 exp - 1 2 d i=1 (x i -y i ) 2 r 2 i + θ 2 ,
where Θ = (θ 1 , θ 2 , r 1 , r 2 , . . . , r d ) are some parameters to be determined. The first term is a distance-dependent correlation between the function values at two data points; if their distance is small compared to the length scales r i , the exponential term is close to one while it decays exponentially as their distance increases. Practical studies have shown that the use of a nonisotropic scaling is beneficial for the problems considered here. The parameter θ 1 scales this correlation. In the second term, θ 2 gives an offset of the function values from zero. It now remains to find the parameters Θ in the correlation function. These parameters are determined by maximizing the joint probability density p(F N ). Indeed, we want the constructed model to be most consistent with the observed data. This is equivalent to minimizing the loglikelihood function given by:

L = F N C -1 N F N + log det(C N )
This function is known to be multi-modal; hence an evolution strategy is employed for this work, that has the capability to avoid local minima. There are many practical issues that must be taken care of so that all the computations are stable; we follow [START_REF] Buche | Accelerating evolutionary algorithms with gaussian process fitness function models[END_REF] in this respect. The hyperparameters must be non-negative; hence it is better to work with the logarithm of the parameters so that they always remain positive. The correlation matrix can be ill-conditioned in which case the computation of its inverse will not be accurate. If the condition number is above a specified tolerance, then the log-likelihood is taken to be a large positive value. If the condition number is within the specified tolerance, an LU decomposition of C N is first performed; then

C -1 N F N and C -1
N k are computed using the LU decomposition. The logarithm of the determinant of C N is also computed using the LU decomposition

log det(C N ) = N n=1 log L nn + N n=1 log U nn
This avoids the under-flow problem associated with multiplying a large number of small numbers which would be the case if the matrix C N is badly scaled. An upper limit of 1/ is imposed on the condition number of C N , where is the machine precision.

Again following [START_REF] Buche | Accelerating evolutionary algorithms with gaussian process fitness function models[END_REF] we scale the coordinates and function values. Moreover, the parameters in the covariance function are constrained by upper and lower bounds:

θ 1 ∈ [10 -3 , 1] θ 2 ∈ [10 -3 , 1] r i ∈ [10 -2 , 10], i = 1, . . . , d

Optimization strategy based on kriging models

The optimization strategy used in this study is based on the iterative construction of a kriging model (see [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF] for a review paper on these methods). The use of such a model for optimization must be an iterative process since it is not possible to build a model accurate enough to find the optimal parameters in only one step. The model should be updated with the results from flow simulations until some convergence criterion is fulfilled. Therefore, the algorithm is organized in two phases. During the first one, an initial a priori database is constructed, that gathers the flow responses (objective function values) corresponding to different control parameter values. The control parameters are chosen in order to explore uniformly the search space, according to a Design Of Experiments (DOE) method. During the second phase, a kriging model is constructed on the basis of available data and is used to determine which flow simulations should be carried out and added into the database. This second phase is then repeated until convergence of the algorithm.

Using such an algorithm, the database and the model are iteratively enriched and reconstructed. The choice of the optimization algorithm used to solve the minimization problem in step 3 is not critical, since this is an inexpensive task. An evolution strategy is used in practice since it is robust and not sensitive to local optima.

The robustness and efficiency of this algorithm depends critically on the choice of the new simulations to be carried out during the step 2. To obtain satisfactory results, the capability of a kriging model to predict the uncertainty related to value predictions is used. Indeed, for a given design vector x, a kriging model provides not only a prediction of the function value f (x), but also an estimation of the uncertainty of this prediction σ(x). Therefore, one can use this knowledge to determine not only the point that minimizes the model, but also areas for which the model is uncertain, i.e. of poor quality. Evaluating such points promotes the minimization of the cost function as well as the improvement of the kriging model.

Thus, we define a merit function that can be interpreted as a statistical lower bound of the model, by aggregating the function value and the standard deviation value [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF]:

f M (x) = f (x) -ρσ(x), ( 12
)
where ρ is a positive parameter that allows to balance the two terms. The algorithm can finally be summarized as:

1. Build an a priori database 2. Construct a kriging model 3. For i = 1 to p : Find the point x i that minimizes the merit function f (x) -ρ i σ(x) (with ρ i = i -1 in general)

4. Evaluate the p points (x i ) i=1,...,p and add them in the database

Return to step 2 until convergence

The use of ρ 1 = 0 yields the minimization of the kriging model. Moreover, the use of other values, e.g. ρ 2 = 1, ρ 3 = 2, promotes the exploration of areas where the standard deviation is high and where possible interesting points can be found. These extra points will be useful since they allow to improve the model quality. This approach has been found robust and applied to various problems [START_REF] Buche | Accelerating evolutionary algorithms with gaussian process fitness function models[END_REF][START_REF] Duchaine | Computational-fluid-dynamics-based kriging optimization tool for aeronautical combustion chambers[END_REF][START_REF] Praveen | Study of some strategies for global optimization using Gaussian process models with application to aerodynamic design[END_REF].

The optimization process is stopped according to three possible criteria : the user can set a maximum number of simulations to perform or a target for the cost function value. Moreover, there is a third criterion, related to the iterative construction of the kriging model: the process is stopped when the minimization of the merit functions yields to points already in the database (for a given accuracy).

Nevertheless, this algorithm exhibits some limitations: first, the kriging model can be sensitive to the error originating from the evaluation process of the cost function. Typically, if the cost function is evaluated through a numerical simulation, a low convergence of the solver can lead the search to spurious local optima. Secondly, this approach is limited to problems of rather low dimension. The number of parameters should be typically lower than twenty in order to have an efficient search, according to the literature results. However, this is usually the case for flow control applications. In a previous work [START_REF] Praveen | Study of some strategies for global optimization using Gaussian process models with application to aerodynamic design[END_REF], this algorithm was used to solve aerodynamic shape optimization problems on the basis of inviscid flow models, with an increasing number of design parameters. In this context, it has been found that the proposed algorithm provides satisfactory results up to 32 variables.

APPLICATION TO FLOW CONTROL ON A BACKWARD FACING STEP

The backward facing step has been selected to test the optimization method on a flow control problem. This test-case has been widely used for code and turbulence models validation [START_REF] Rumsey | Application of fun3d and cfl3d to the third workshop on cfd uncertainty analysis[END_REF]. Despite the fact that the turbulence model and the jet model used in this study are quite simple, this problem will demonstrate the capability to optimize control parameters for detached flows.

Computational conditions

The backward facing step configuration and flow conditions correspond to the case presented by Driver and Seegmiller in [START_REF] Driver | Features of a reattaching turbulent shear layer in divergent channel flow[END_REF]. The boundary conditions used for the computation are shown on the Fig. 4. The inlet flow Mach number is set to M = 0.128 and the velocity magnitude is set to U ref = 44.2 m/s. The outlet pressure is set to P = 1.0015 P ∞ . The Reynolds number based on the boundary layer thickness is 5000. In our computations, the top wall is replaced by a symmetry boundary condition. The first point of the synthetic jet is located at -h/50 from the step corner. Its diameter is set to h/10 and its direction is perpendicular to the wall. The amplitude and frequency ranges are defined in section 4.1.2. The unsteady computation is initialized using a steady solution, corresponding to the same case with an inactive jet. Calculation is performed with a time step of t * = tV ref /h = 0.004. The dual time step sub-iterations are stopped when they have reached 15 iterations or when the residual has been reduced of 5 orders.

Computational grid

All the grid generated are unstructured grids on which control cells are constructed following Barth method [START_REF] Barth | Aspects of unstructured grids and finite-volume solvers for the euler and navierstokes equations[END_REF][START_REF] Viozat | On vertex-centered unstructured finitevolume methods for stretched anisotropic triangulations[END_REF].

The steady computation is done on three different grids. From coarse to fine, the grids contains respectively 15946, 29606 and 61344 vertices. For each grid, the distance between the first interior node and the wall is 5.10 -6 . For all the cases, this restriction on the grid to be simulated and included into the database. Unfortunately, these computations are still in progress.

Controlled case

Although the optimization process is not finished, the kriging model has already determined a parameter region for which the actuation allows to reduce the time-averaged length of the recirculation bubble. A suitable set of parameters to achieve this result is typically :

• Large blowing amplitude: U jet /U ref ≈ 0.8

• Low frequency actuation: F jet h/U ref ≈ 0.05 With this configuration, the figure 8 provides some snapshots of the flow along an actuation period. As can be seen, a large recirculation bubble is generated by the actuation, which makes the reattachment point moving. Figure 9 shows a comparison between the stationary flow without actuation and the time-averaged flow with actuation. In particular, one can observe a reduction of the time-averaged recirculation length of about 5%. One can also notice the change of the shape of the reticulation bubble, whose thickness is higher in the case of actuation. Better results might be obtained by moving the jet or changing the jet direction.

Conclusion and prospects

The kriging-based optimization algorithm has been coupled with an unsteady RANS solver. This approach has been applied to optimize the amplitude and frequency of a synthetic jet located just before a backward facing step, in order to shorten the time-averaged recirculation length. This study has shown the capability of the kriging model to determine efficient control parameters, although the complete optimization process has not been achieved. These results should be then carefully validated by quantifying the influence of numerical parameters, such as grid size and time step, as well as numerical models and in particular by assessing the role of the turbulence closure.

  ) and Fig.1(b)).

  (a) A control cell in 2D. (b) A control cell in 3D.

Figure 1 :

 1 Figure 1: Illustration of control cells.

2. 5

 5 Boundary conditions 2.5.1 Weak approach for inlet and outlet conditions Some boundary conditions are implemented in a weak form through numerical fluxes, as illustrated in the Fig.2.

Figure 2 :

 2 Figure 2: Illustration of the implementation of weak boundary conditions.

Figure 3 :

 3 Figure 3: Illustration of the implementation of strong boundary conditions.

Figure 4 :

 4 Figure 4: Illustration of the backward facing step configuration.

  (a) Surrogate model of the function of interest. (b) Database and new points to be evaluated.

Figure 7 :

 7 Figure 7: Surrogate model of recirculation length.

  (a) t = t 0 (b) t ≈ t 0 + T /8

Figure 8 :Figure 9 :

 89 Figure 8: Streamlines over a cycle for the large blowing amplitude and low frequency actuation.

size ensures to fulfill the condition y+ < 1. Presently, the unsteady computations and the optimization exercise are only done on the coarse grid.

Design of experiment

The optimization parameters considered in this study are restricted to the amplitude and frequency of the synthetic jet. To construct the a priori database used to explore the design space, these jet parameters are supposed to vary in the following ranges:

Then, an initial sampling of 16 flows is generated using an optimized Latin Hypercube Sampling (LHS) , as illustrated in figure 5.

Results and discussion

Validation without actuation

We first validate our computations by considering the case without actuation. The case is the same as the backward facing step configuration presented in 4.1, but with a wall instead of the jet boundary condition. The friction coefficient, the pressure and velocity distibution profiles obtained with the NUM3SIS platform are compared to experimental data and numerical results from the CFL3D code developed at NASA [START_REF] Driver | Features of a reattaching turbulent shear layer in divergent channel flow[END_REF][START_REF] Rumsey | Application of fun3d and cfl3d to the third workshop on cfd uncertainty analysis[END_REF]. As seen in figure 6, the results converge to those obtained by the CLF3D code with a very fine grid (about 300 000 nodes). However, one can notice that the Spalart-Allmaras model is not able to reproduce correctly the flow in the recirculation area. Nevertheless, the length of the recirculation bubble and the reattachment point are quite well predicted, even on the coarsest grid. Therefore, the coarsest grid is used in the following for the unsteady computations.

Surrogate model construction

Based on the LHS, 16 different actuations are simulated using NUM3SIS. The time-averaged recirculation lenght l/h computed at the wall is selected to be the function of interest. A kriging model is constructed and can be seen in figure 7(a). Using such a model, one can easily identify some parameter regions for which the length of the recirculation bubble is reduced. According to the algorithm described above, three points are then selected by minimization of the three merit functions. These new candidates can be seen in figure 7(b). Then, these new points have


First A Zerbinati
email: adrien.zerbinati@inria.fr

Second J A Désidéri

Third R Duvigneau
email: regis.duvigneau@inria.fr

APPLICATION OF METAMODEL-ASSISTED MULTIPLE-GRADIENT DESCENT ALGORITHM (MGDA) TO AIR-COOLING DUCT SHAPE OPTIMIZATION

Keywords: optimization, multi-objective, gradient descent direction, Pareto front, metamodel

. Using MGDA in a multi-objective optimization problem requires the evaluation of a substantial number of points with regard to criteria, and their gradients. In industrial test cases, in which computing the objective functions is CPU demanding, a variant of the method was to be found. Here, a metamodel-assisted MGDA is proposed and tested. The MGDA is assisted by a Kriging surrogate model construction. A first database is computed as an Latin Hypercube Sampling (LHS) distribution in the admissible design space, which is problem-dependent. Then, MGDA leads each database point to a non dominated set of the surrogate model. In this way, each function computation is made on the surrogate model at a negligible computational cost.

INTRODUCTION

MGDA stands for Multiple-Gradient Descent Algorithm [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF]. It is a generalization of the classical steepest-descent method [START_REF] Gill | Practical optimization[END_REF] that applies to cases in which an arbitrary number of criteria, of known gradients, are to minimized, or simply reduced. It is based on the observation that if ω is the minimal norm element in the convex hull of these gradients, then -ω is a descent direction for all criteria simultaneously. If ω = 0, the considered point belongs to the Pareto set 1 . MGDA was tested successfully in [START_REF] Zerbinati | Comparison between MGDA and PAES for multi objective optimization[END_REF] on several analytical test cases proposed in [START_REF] Deb | Scalable Test Problems for Evolutionary Multi-Objective Optimization[END_REF] to asses evolution strategies. One of then is reported here. Using MGDA in a multi-objective optimization problem requires the evaluation of a substantial number of points with regard to criteria, and their gradients. In the particular case of Computational Fluid Dynamics (CFD) problems, each point evaluation is very costly, in term of CPU. Thus, here, we propose to alleviate this difficulty by constructing metamodels and calculating approximate gradients. As a test case for numerical experiments, we consider the problem of shape optimization of an automobile air-cooling duct. It consists of the simultaneous minimization of two criteria associated with a compressible Navier-Stokes flow inside the duct.

MGDA

MGDA principles

MGDA is a gradient based method of multi-objective optimization. This method could be applied to cases in which an arbitrary number of criteria, of known gradients, are to be minimized. A particular point in the design space is called Pareto-optimal when it is non dominated every other acceptable point. Thus consider n smooth criteria J i (Y ) (Y : design vector ; here in R N ). For particular reasons explained in [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF], n ≤ N and cost functions are assumed to be C 2 in some working ball in the design space R N . Let Y 0 be a Pareto-optimal point of the smooth criteria (J i (Y 0 )) and define the gradient vectors u 0 i = ∇J i (Y 0 ) in which ∇ denotes the gradient operator. There exists a convex combination of the gradient-vectors that is equal to zero. Inversely, if the smooth criteria J i (Y ) are not Pareto-optimal at a given point Y , descent directions common to all criteria exist. Let U (Equation 1) be the convex hull of the gradient vectors :

U = w ∈ R N , w = n i=1 α i u i , , α i ≥ 0, n i=1 α i = 1 . (1)
U is a closed and convex set. This implies existence and uniqueness of the element ω of minimum norm in U (Equation 2) :

ω = min w∈U w . (2)
Then, -ω is a descent direction for all criteria simultaneously. If ω = 0, the point is on the Pareto set. The particular case of two criteria is illustrated in Figure 1. Assume Y 0 ∈ R N , a particular design vector and J i (Y 0 ) the i th criteria value computed at this particular point. For each criterion, at Y 0 , the gradient value been computed. Then, the minimal norm vector in the convex hull of the gradients is denoted ω. If ω > ε, (ε is a user fixed parameter), a new point in -ω direction which dominates (in Pareto sense) Y 0 exists. In this particular direction, an optimal step size ρ is computed as into by fitting a local problem explained in [START_REF] Zerbinati | Comparison between MGDA and PAES for multi objective optimization[END_REF]. Then, the new point Y 0 -ρω is used as starting point by the MGDA. This method is summarized in Figure 2

Analytical validation

In this section, MGDA is illustrated on a problem proposed by Fonseca [START_REF] Deb | Scalable Test Problems for Evolutionary Multi-Objective Optimization[END_REF]. This test case corresponds to the two-objective unconstrained minimization of the following functions :

               f 1 (Y ) = 1 -exp - 3 i=1 y i - 1 √ 3 2 f 2 (Y ) = 1 -exp - 3 i=1 y i + 1 √ 3 2 , Y = (y 1 , y 2 , y 3 ). (3)
The design variable is Y = (y 1 , y 2 , y 3 ) ∈ R 3 . This test case is known to yield a continuous but non convex Pareto set in function space. The Pareto front was identified by Deb using the well-known genetic algorithm NSGA-II [START_REF] Deb | Scalable Test Problems for Evolutionary Multi-Objective Optimization[END_REF]. Furthermore, the Pareto front can be calculated analytically. From a starting point, MGDA converges rapidly (6 steps in the example chosen) and provides an accurately defined point on the Pareto set Figure 3. A set of initial design points produces an accurate set of non dominated points. To obtain an accurate representation of the Pareto set by MGDA, we have applied the method starting from a set of 60 initial design points distributed on a sphere in the design space, around the known non dominated set on Figure 4. The results are similar to those obtained by Deb using NSGA2.

In summary, MGDA converges rapidly from an initial design point to a non dominated point. At each iteration, the cost functions are evaluated twice at least in order to determine the step size. This remains possible for fast function evaluations, in term of computational time. In case of functions with an expensive computational cost, the strategy must be adapted.

3 Metamodel assisted MGDA

Metamodel assisted MGDA principles

Using MGDA in a multi-objective optimization problem requires the evaluation of a substantial number of points with regard to the criteria of interest along with their gradients. In the particular case of non-linear problems, such as CFD or non-linear mechanics problems, each point evaluation is very costly. Thus, we propose to overcome this difficulty by estimating function values on a surrogate model corresponding to the response of the functions within the design search space.

An initial set of design points is generated by using of a LHS method in R N . As illustrated in Figure 5, the databases are also used to provide starting points to initiate the MGDA iteration, which is conducted until convergence using gradients that are calculated on the basis of the metamodels. Each converged point belongs to the Pareto set associated with the two-criterion Figure 5: MGDA with surrogate model scheme. A surrogate model based on an initial database is trained. Then, MGDA [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF] is applied from each database point. Thus a non dominated set on the surrogate model is obtained. Unless a non dominated point obtained is too close to an initial database element, its exact performance is computed and used to increase the database.

problem related to the metamodels. This point is then reevaluated by a flow computation, and added to the database. A filtering method is used to remove this point if it is too close to an already existing point. At completion of this database enrichment process, the metamodel is updated, which describes the cycle. For more robustness and to reduce the antagonism the conflict between parameters of different nature (angle, length, ...) and also different scales, during the optimization, all parameter intervals are scaled to be in [-1, 1].

Numerical experimentation

As a test case for numerical experiments, we consider the problem of shape optimization of an air-cooling duct in which an compressible Navier-Stokes flow is considered. The inlet and outlet sections are fixed. Thus, only the shape of the elbow is optimized. The elbow between these two sections is defined by 8 parameters, 3 angles (X 2 , X 5 and X 6 ), and 5 lengths (X 1 , X 3 , X 4 , X 7 and X 8 ). See Figure 6. The shape is represented by 8 B-spline functions joining the inlet duct and the middle duct of the elbow. Eight more B-splines maintain a smooth junction between the middle part and the outlet. Each parameter can vary in the interval [-1, 1]. The proposed approach is used here to solve a two-criterion optimization problem in subsonic conditions M ∞ = 0.1 and transition flow (Re ≈ 2800). For each evaluation, the compressible Navier-Stokes equations are solved with our in-house CFD code Num3sis. The CAD and the mesh are made with a script style software GMSH [START_REF] Geuzaine | a three-dimensional finite element mesh generator with builtin pre-and post-processing facilities[END_REF]. It consists in minimizing the velocity variance as well as the pressure loss computed at a particular duct section, simultaneously, with bound constraints. The two objectives are computed as depiced in Figure 7 In practice, an initial database from R N is considered. Then, a Kriging metamodel is generated using this database, for the pressure loss and the velocity variance values. MGDA is applied using each database element as a starting point. Both pressure loss and velocity variance computations are evaluated on the surrogate model. In the same way, gradients are evaluated on the surrogate model. Thus, a non dominated set (with regard to the surrogate model) is obtained. Until a non dominated point is too close to an element of the database (in the search space), its

Velocity variance

Pressure loss

• Um = 1 Vtot d(cel,Po)≤ε u(cel) * V ol(cel) • σ 2 velx = 1 Vtot d(cel,Po)≤ε (Umx -ux(cel)) 2 * V ol(cel) • σ 2 = σ 2 velx + σ 2 vely + σ 2 velz • p k = P k m • u k = U k m • ∆p = p i -po + ρ i u 2 i 2 - ρou 2 o 2
Figure 7: Objective cost functions computation. The velocity variance is computed in a particular section of the outlet duct, designed in the last line of these tabular. U m (resp. P m ) is the weighted mean velocity (resp. pressure) in relation to cell volume V ol(cel) and the considered cells total volume V tot . The pressure loss is computed between two particular sections, on of the inlet duct and the second of the outlet duct.

performance is computed using a solver simulation and the initial database is enhanced.

Numerical results

In the following example, an initial database of 30 points in R 8 is considered. The metamodel assisted MGDA cycle was made 9 times. After each step, the previous database is enhanced by the resulting points of the MGDA applied on the surrogate model, filtered and evaluated with the solver. Finally, a set of 223 points is considered with a non dominated set of 4 points. The total number of evaluations and the best values computed is resumed in Figure 8. magnitude on different section of the air-cooling duct. One section is located just after the elbow and the second close to the outlet. The second section indicates the velocity variance is computed. The nominal shape, the lowest velocity variance shape and the lowest pressure loss shape velocity magnitude are represented for comparison. The comparison between the best points computed flows and the nominal flow shows the increase of the longitudinal vortex during the optimization. This raise explains the reduction of velocity variance shown in the outlet of the duct for the two optimal duct shapes. At the same time, the cross-wise vortex is reduced which explain the pressure loss reduction. This is evident on Figure 10 The reduction of the objectives (minimization of variance velocity and pressure loss simultaneously) requires in both cases an narrower air-cooling duct shapes. This explains the weak difference between the two optimal shapes obtained by the optimization shown on Figure 12 with the nominal shape for comparison.

CONCLUSIONS

In this article, we have tested by numerical experiment a recently proposed gradient-based algorithm MGDA [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF] for multi-objective optimization. The convergence to Pareto-optimal solutions has been demonstrated in an analytical test case proposed by Fonseca [START_REF] Deb | Scalable Test Problems for Evolutionary Multi-Objective Optimization[END_REF] corresponding to a concave Pareto front.

A metamodel assisted MGDA has been developed and successfully tested on a two-objective industrial shape optimization test case. It consists in minimizing variance velocity and pressure loss simultaneously for an air-cooling duct for car using. The result obtained is good with regards to the number of solver evaluations.

In the future, the metamodel-assisted MGDA will be tested on optimization cases with more than two objectives. An ambitious target will be to solve multi-disciplinary optimization problem of concurrent engineering.

Figure 1 :

 1 Figure 1: Assume u and v two vectors. In each case, ω represent the minimal norm vector of the {u, v} convex hull.

Figure 2 :

 2 Figure 2: Modular scheme of the Multiple Gradient Descent Algorithm (MGDA)

Figure 3 :

 3 Figure 3: Convergence of MGDA from an initial design point to the non dominated set.

Figure 4 :

 4 Figure 4: Convergence of MGDA from initial design points around Pareto set, for a classical test case proposed by Fonseca, in design space (left), in function space (right). In the functional space, the analytical Pareto front is plotted to compare.

Figure 6 :

 6 Figure 6: Air-cooling duct for car use. Inlet and outlet ducts are fixed. The elbow shape is defined by 8 parameters called X i , i = 1, . . . , 8 and B-splines.

Figure 9 Figure 8 :Figure 9 :

 989 Figure 8: Initial and final statistics (after 9 metamodel-assisted MGDA cycles).

Figure 10 :

 10 Figure 10: Sections of the air-cooling duct colored according to velocity magnitude. Just after the elbow on the left side and close to the outlet on the right side. The nominal, lowest velocity variance and lowest pressure loss are represented for comparison.

Figure 12 :

 12 Figure 12: Shapes corresponding to cases 1, 2 and 3.

Actually a slightly lees stringent condition holds : Pareto-stationarity[START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF]


Gang Xu

Bernard Mourrain
email: bernard.mourrain@inria.fr

Régis Duvigneau
email: regis.duvigneau@inria.fr

André Galligo
email: galligo@unice.fr

A New Error Assessment Method in Isogeometric Analysis of 2D Heat Conduction Problems

Keywords: Isogeometric analysis, posteriori error estimation, k -refinement, heat conduction

In this paper, we propose a new error assessment method for isogeometric analysis of 2D heat conduction problems. The posteriori error estimation is obtained by resolving the isogeometric analysis problem with several k -refinement steps.

The main feature of the proposed method is that the resulted error estimation surface has a B-spline form, according to the main idea of isogeometric analysis. Though the error estimation method is expensive, it can be used as an error assessment method for isogeometric analysis. Two comparison examples are presented to show the efficiency of the proposed method.

I. INTRODUCTION

Usually, CAD modeling software relies on splines or NURBS representations, while the CAE software for CAD object uses mesh-based geometric descriptions (structured or unstructured). Therefore, in conventional approaches, several information transfers occur during the design phase, yielding approximations and non-linear transformations that can significantly deteriorate the overall efficiency of the design optimization procedure.

The isogeometric analysis (IGA for short) method proposed by Hughes et al. in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF] can be employed to overcome the gap between CAD and CAE. The approach uses the same type of mathematical representation (spline representation), both for the geometry and for the physical solutions, and thus avoids data transfers between the design and analysis phases. Moreover it reduces the number of parameters needed to describe the geometry, which is of particular interest for shape optimization. This framework allows to compute the analysis solution on the exact geometry (not a discrete geometry), obtain a more accurate solution (high-order approximation), reduce spurious numerical sources of noise that deteriorate convergence. Moreover, NURBS representation is naturally hierarchical and allows to perform refinement operations to improve the analysis result.

Since the concept of isogeometric analysis was proposed, many researchers in the fields of computational mechanical and geometric modeling were involved in this topic. The current work on isogeometric analysis can be classified into three categories: (1) application of IGA to various simulation and analysis problems [START_REF] Auricchio | A fully locking-free isogeometric approach for plane linear elasticity problems: A stream function formulation[END_REF] The topic of this paper belongs to the third field. As far as we know, there are few works on the error estimation method in isogeometric analysis. Bazilevs et al. studied the error estimation for h -refined mesh in isogeometric analysis [START_REF] Bazilevs | Isogeometric analysis: approximation, stability and error estimates for h -refined meshes[END_REF]; some estimates for h p k -refinement in isogeometric analysis are investigated in [START_REF] Beirao Da Veiga | Some estimates for h p k  -refinement in isogeometric analysis[END_REF]; Dörfel et al. proposed a posteriori error estimation for local h -refinement with T- splines [START_REF] Dörfel | Adaptive isogeometric analysis by local h-refinement with T-splines[END_REF]. In [START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF], an error assessment method is proposed based on h -refinement operation. In this paper, we propose a new error assessment method for isogeometric analysis of two-dimensional heat conduction problems, which is obtained by resolving the isogeometric analysis problem with several k -refinement steps.

The remainder of the paper is organized as follows.

Section II introduces the isogeometric analysis and k - refinement for two dimensional heat conduction problems.

Section III presents the error assessment method based on k - refinement for isogeometric analysis of two dimensional heat conduction problems. Some examples and comparisons with h -refinement method are presented in Section IV. Finally, we conclude this paper and outline future works in Section V.

II. ISOGEOMETRIC ANALYSIS OF HEAT CONDUCTION PROBLEM

Given a domain  with DN      , for ease of presentation, we consider the two dimensional second order elliptic PDE with homogeneous Dirichlet boundary condition as an illustrative model problem :

( ( )) ( ) ( ) ( ) ( ), D U d f d H                 x x x x
where ()  x are test functions. After integrating by parts and using boundary conditions, we obtain:

( ) ( ) ( ) ( ) . U d f d          x x x x (2)
According to the IGA paradigm, the temperature field is represented using B-spline basis functions. For a 2D problem, we have:

11 ˆ( ) ( , ) ( ) ( ) , j i j i n n p p i j ij ij U N N T        x 
where ˆi N functions are B-Spline basis functions and  u ( , )

  are domain parameters. Then, we define the test functions ()

 x in the physical domain such as: ˆˆ( ) ( , ) ( ( , )) ( , ) ( ) ( ). j i p p ij ij ij ij i j N N x y N N N N           x  The weak formulation Eq. (2) reads: 11 ( ) ( ) ( ) ( ) . kl nn kl kl ij ij kl T N N d f N d         x x x x
Finally, we obtain a linear system similar to that resulting from the classical finite-element methods, with a matrix and a right-hand side defined as: , ( ) ( )

( ) ( ) ( ) ( ) ( ) ij ij kl kl T kl kl M N N d N B N J d           u u x x u u u u u    
where J is the Jacobian of the transformation, T B is the transposed of the inverse of the Jacobian matrix. The above integrations are performed in the parameter space using classical Gauss quadrature rules.

Starting from a planar B-spline surface as computational domain, an isogeometric solver for thermal conduction problem (1) has been implemented in the AXEL 1 platform, yielding a B-spline surface as solution field. Gauss-Seidel algorithm is employed to solve the linear system. Fig. 1 shows an example of planar B-spline surface as computational domain and the corresponding isogeometric analysis results for two-dimensional heat conduction problems.

In order to improve the simulation results, refinement operation can be performed for two parametric directions. There are three kinds of refinement operations in isogeometric analysis: h -refinement by knot insertion, p -refinement by degree elevation operation, and k -refinement combining by knot insertion and degree elevation. The k -refinement operation is performed by elevating the degree of basis function to a desired order firstly followed by knot insertion thus obtaining the maximum available continuity. Compared with finite element analysis, the main advantage of refinement operations in isogeometric analysis is that the geometry of the computational domain can be kept while the degree of freedom increases. Fig. 2 presents an example to compare h - refinement, p -refinement and k -refinement in isogeometric analysis. Note that the number of control points increases during the refinement operation, and more degree of freedom can be achieved by k -refinement.

III. ERROR ASSESSMENT METHOD BASED ON

k -REFINEMENT Suppose that () U x is the exact solution, and () h U x is the approximation solution obtained by isogeometric method in Section II, then the discrete error e can be written as

.

h e U U  (3)
After performing Laplacian operator  on each side of (3), a posteriori error assessment can be obtained by resolving the following problem, in 0 on

h e f U e      Δ Δ (4)
From ( 4), the crucial point of a posteriori error estimation is the computation of () ,

hh h UU U xy JL JL L L JK                      where 22 2 2 2 2 , ( ) ( ) ,
,,

J x y x y K x y x y L x y L x y                       ( ) ( ) , L L y L y x L x L x y                ( ) ( ) . L L y L y x L x L x y               
Proof. The idea of isogeometric analysis is to use the same mathematical representation for the computational domain and solution field. Suppose that the computational domain  is parameterized by the following planar B-spline surface: 11 ( , ) ( ( , ), ( , )) ( ) ( ) ,

j i j i n n d d i j ij ij x y N N             c
In isogeometric analysis, the solution field of heat conduction problem (1) over the computational domain  has the following form, 11 ( ( , ), ( , )) ( , ) ( ) ( ) ,

j i j i n n d d h i j ij ij U x y N N T             Here ij
T are the unknown variables in isogeometric analysis to be solved from the boundary condition and Eq. ( 1). From ( ( , ), ( , )) ( , )

h U x y         ,

( ) ( ) ,

h h h h U U U U x x y y x x y y                             22 2 2 2 22 2 2 2 2 2 ( ) ( ) , h h h h U U U U x x y y x x y y                            
From above two equations, we have

( )] / , h h h hh U U U y x y x x y UU y x y K xy                        2 2 2 2 [ ( ) ( )] / , h h h hh U U U x x y y x y UU x x y K xy                        where 22 ( ) ( ) K x y x y     2 2 2 2 [ ( )
.

Hence, we can obtain

JL JL L L JK                                              where 2 2 2 2 ,, L x y L x y           ( ) ( ) , L L y L y x L x L x y                ( ) ( ) . L L y L y x L x L x y               
This completes the proof. 

The approximation error surface e from (4) also has a B- spline form. Different from the method in [START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF], we perform several k -refinement operations to achieve more accurate results. Though it is much more expensive, we can use it as an error assessment method for isogeometric simulation solutions. As a summary, the procedure of error assessment method for model problem (1) can be described as follows, From this example, we can find that k -refinement method can achieve a better approximation of the exact error surface than h -refinement method.

As an example with unknown exact solution, the error color map of simulation result in Fig. 1 obtained by h -refinement method and k -refinement method are shown in Fig. 4. Note that the original design of computational domain has significant impact on the simulation results. Hence, for different parameterization of computational domain, the resulted error surface obtained by the proposed method is also different.

V. CONCLUSIONS A new error assessment method for isogeometric analysis of two-dimensional heat conduction problems is proposed in this paper. The basic idea is to resolve the isogeometric analysis problem with several k -refinement steps. The main feature of the proposed method is that the resulted error estimation surface also has a B-spline representation. It can be used as an error assessment method for isogeometric analysis results. The efficiency of the proposed method is proved by several comparison examples.

In the future, we will generalize the proposed method to 3D case and employ it to validate the r -refinement method in 3D isogeometric analysis as in [START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF].

  [START_REF] Bazilevs | Isogeometric fluid structure interaction: Theory, algorithms, and computations[END_REF][9][11]; (2) application of various modeling tools in geometric computation to IGA [6][10] [7] [14];(3)error estimation, accuracy and efficiency improvement of IGA framework by reparameterization and refinement operations [1][4][3][8][13][15][16].

  the Cartesian coordinates,  is a Lipschitz domain with boundary  , on  and:

1

 1

Fig. 2

 2 Fig.1 An example of isogeometric analysis for twodimensional heat conduction problem: (a) computational domain with control points and iso-parametric curves; (b) isogeometric simulation result.

Fig. 3

 3 Fig.3 Comparison of error assessment method based on h -refinement and k -refinement: (a) isogoemetric solution surface with control points; (b) error surface obtained by three h -refinement steps; (c) error surface obtained by three k -refinement steps; (d) exact error color map; (e) color map of error surface in Fig. 3 (b); (f) colormap of error surface in Fig. 3 (c).

2 .Fig. 3 Fig. 4

 234 Fig. 3 illustrates an example over the computational domain [0,3] [0,3]  , which has exact solution for problem (1). The isogeometric solution surface with control points are shown in Fig. 3(a), the error surface obtained by three h -refinement steps is illustrated in Fig. 3 (b), Fig. 3 (c) shows the error

ACKNOWLEDGMENT The authors are partially supported by the 7th Framework Program of the European Union, project SCP8-218536 ``EXCITING''. The first author is partially supported by the National Nature Science Foundation of China (No.61004117), the Zhejiang Provincial Natural Science Foundation of China under Grant No.Y1090718, Defence Industrial Technology Development Program (No.A3920110002), and the Open Project Program of the State Key Lab of CAD &CG (A1105), Zhejiang University.