A continuous Lagrangian sensitivity equation method for incompressible flow

Journal of Computational Physics 231 (2012) 5989–6011

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

A continuous Lagrangian sensitivity equation method for incompressible flow L. Charlot ⇑, S. Etienne, D. Pelletier École Polytechnique de Montréal, Mechanical Engineering Department, C.P. 6079, succ. Centre-ville, Montéal, Québec, Canada H3C 3A7

a r t i c l e

i n f o

Article history: Received 25 November 2010 Received in revised form 15 December 2011 Accepted 20 February 2012 Available online 5 April 2012 Keywords: Sensitivity Lagrangian point of view Mesh deformation Finite element method Navier Stokes equations Incompressible flow simulations

a b s t r a c t A continuous Lagrangian sensitivity equation method (CLSEM) is presented as a cost effective alternative to the continuous (Eulerian) sensitivity equation method (CESEM) in the case of shape parameters. Boundary conditions for the CLSEM are simpler than those of the CESEM. However a mapping must be introduced to relate the undeformed and deformed configurations thus making the PDEs more complicated. We propose the use of pseudo-elasticity equations to provide a general framework to generate this mapping for unstructured meshes on complex geometries. The methodology is presented in details for the incompressible Navier–Stokes and sensitivity equations in variational form. The PDEs are solved with an adaptive FEM. Sensitivity data obtained with both approaches for a flow around a NACA 4512 are used to obtain estimates of flows around nearby geometries. Results indicate that the CLSEM produces significant improvements in terms of both accuracy and CPU time. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction In a design procedure, the CFD engineer has to deal with many parameters defining the flow problem and complex numerical simulations. An optimization algorithm is used to find an optimal configuration by minimizing a cost function. For the gradient type algorithms, the derivative of the cost function with respect to the parameter is needed. This is often achieved by computing an adjoint state. An alternative consists in computing the derivatives of the flow variables (velocity, pressure, etc.) referred as the flow sensitivities. For an airfoil, @u/@ a is the sensitivity of the velocity with respect to the angle of attack a around its nominal value. While adjoint method have proven more effective for many design problems, sensitivities find additional non optimization uses: they can be used to obtain quick and inexpensive estimates of the flow for nearby values of the design parameters without resorting to a full blown reanalysis. This is achieved using Taylor series in parameter space and is especially useful to answer what if questions. They can also serve to cascade input data uncertainty through CFD code to obtain uncertainty estimates of the flow response or to sort the numerous parameters from the most important to the less one. Hence they give to the engineer helpful informations to understand the physical phenomena involved in the problem. Generally speaking parameters can be classified in two categories. Value parameters (boundary condition, a physical property) do not affect the geometry of the domain and lead to a simple treatment. The case of a shape parameter is more challenging because a perturbation of its value changes the geometry of the computational domain so that differentiation of a flow variables necessitates implicit differentiation of the geometry of the domain, a non trivial issue.

⇑ Corresponding author. Tel.: +1 5143404711x5897; fax: +1 5143405917. E-mail address: [email protected] (L. Charlot). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.02.028

5990

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

In our previous work [1] we used the Eulerian sensitivity concept in developing a continuous (Eulerian) sensitivity equation method (CESEM). In this approach, computational challenges arise for shape parameters because boundary conditions along parameter dependent boundary segments necessitate the evaluation of first and/or second order partial derivatives of the flow solution at the boundary. While reconstruction techniques can yield improved accuracy of boundary derivatives [2], they prove costly in 2-D and impractical in 3-D. In this paper, we present an alternative: the continuous Lagrangian sensitivity method (CLSEM). Since this point of view follows the deformation of the boundary as the shape parameter is varied, the boundary conditions that were challenging with the CESEM become trivial and their treatment is straightforward. However, a different challenge arises: integrals over a parameter dependent domain must be differentiated. It is possible to map the deformed parameter dependent geometry to a referenced, parameter independent configuration. Differentiation of the equations can then be processed as usual in the reference plane and the result mapped back to the physical deformed configuration. Zolésio [3] and Céa [4] described this approach under the name ‘‘material derivative method’’ as early as in 1981. Delfour and Zolésio [5] set a general mathematical framework and developed general formulations. Existence and uniqueness of shape sensitivity of drag function are shown for compressible Navier–Stokes equations by Plotnikov et al. [6,7]. For the incompressible Navier–Stokes equations, the Fréchet differentiability of the flow solution with respect to a regular enough mapping has been proven by Bello et al. [8]. Similar Lagrangian point of views may be found in the literature. Tortorelli and Wang [9] developed a domain parameterization for shape sensitivity analysis. Many examples can be reported for heat conduction problems: Haftka [10] discusses it in the context of 1-dimensional problems using both sensitivity and adjoint techniques. Sluzalek and Kleiber [11] discuss the material derivative approach for non linear steady state heat conduction, Dems et al. [12] use it to obtain the first and second derivatives of an objective function by both sensitivity and adjoint method. In a set of two papers Dems and Rousselet [13,14] discuss sensitivity analysis for transient heat conduction and consider shape parameters for external boundaries as well as for internal interfaces between two materials. They use an analytical mapping. The Lagrangian point of view is discussed in structural mechanics with an adjoint method [15] or with the direct differentiation method to solve linear elasticity equations [16]. Bobaru and Mukherjee [17] use this approach with an element-free Galerkin method while Lee [18] uses a boundary integral equation formulation. Taroco [19] uses second order sensitivities in fracture mechanics. In fluid mechanics, Wang et al. [20] use the domain parameterization method for contraction design but only a few numerical results are presented. In most previous work, an explicit mapping is used in the differential equations or problem dependent or problem specific parameterizations are used. This means that a new transformation or parameterization must be developed, implemented and tested whenever a new geometry is encountered that cannot be accommodated with existing ones. This severely limits the general applicability of the overall methodology. Our Lagrangian sensitivity equation formulation is original in that the mapping is defined implicitly identifying it with the solution of linear elasticity equation. Thus we use the pseudo-solid approach of Sackinger et al. [21]. This technique offers the significant advantage that the sensitivity solver is as general as the finite element flow solver itself. Thus we propose in this paper a simple and robust implementation leading to accurate evaluation of shape sensitivities. The paper is organized as follows. Sections 2–4 present the flow equations, the continuous Eulerian sensitivity equation method (CESEM) and the continuous Lagrangian sensitivity equation method (CLSEM). Section 5 briefly discusses the Finite element method used to solve the flow and sensitivity equations. Section 6 presents detailed performance and accuracy comparisons of two variants of the CESEM and the proposed CLSEM using the method of manufactured solutions. The methodology is then applied to the evaluation of nearby flows around NACA airfoils of 4 digits series. 2. Flow equation We consider the laminar incompressible steady Navier–Stokes equations on the domain Xa. We use the superscript or subscript a to emphasize on the fact that geometry depends on the shape parameter a.

ru¼0 qðu  rÞu ¼ r  r þ f qcP u  rT ¼ r  ðkrTÞ þ qS

ð1Þ ð2Þ ð3Þ

where r is the stress tensor:

r ¼ pI þ s ¼ pI þ lðru þ rT uÞ with s the viscous stress tensor and I the identity matrix. Dirichlet boundary conditions are prescribed on Cau and CaT and Neumann boundary conditions on Cat and Caq for the velocity and the temperature. Obviously, Cau ; Cat ; CaT and Caq depend on a.

 on Cau u¼u tf ¼ r  n ¼ tf T¼T

ð4Þ a

on Ct

on CaT

 on Caq q ¼ krT  n ¼ q

ð5Þ ð6Þ ð7Þ

5991

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

 and t f are the prescribed values of the boundary velocity and tractions, T and q  are the prescribed values of the where u boundary temperature and heat flux and n is the outward normal of the domain. Let dp, du and dT be admissible test functions associated to the pressure, the flow velocity and the temperature. The weak form of the problem is:

Z

r  udp dX ¼ 0

ZXa

ð8Þ

½qðu  rÞu  du  pr  du þ s : rdu dX ¼

Xa

Z

ðqcP u  rTdT þ krT  rdT Þ dX ¼

Xa

Z

Z

f  du dX þ

Xa

qS dT dX þ

Z

Xa

Caq

Z Cat

t f  du dC

ð9Þ

krT  ndT dC

ð10Þ

3. Eulerian sensitivity equation method We consider now that the flow variables are functions of the space coordinates and of the shape parameter. Thus for the flow velocity we write u = u(x, a). The Eulerian sensitivity equations method solves differential equations for the partial derivatives of the flow variables with respect to the shape parameter a:

su ¼ ðsu ; sv Þ ¼ lim

da!0

uðx; a þ daÞ  uðx; aÞ @u ¼ da @a

3.1. General formulation of the CESEM The CESEM equations are obtained by formal implicit differentiation of the flow Eqs. (1) and (2) with respect to the design parameter a. For the continuity equation, this yields

@ r  u ¼ 0 ) r  su ¼ 0 @a

ð11Þ

Similarly for the momentum and energy equations we have

q0 u  ru þ qsu  ru þ qu  rsu ¼ rsp þ r  ½l0 ðru þ rT uÞ þ lðrsu þ rT su Þ þ f

0

ð12Þ



 q0 cp þ qc0p u  rT þ qcp ðsu  rT þ u  rsT Þ ¼ r  ðk0 rT þ krsT Þ þ q0S

ð13Þ

The derivatives with respect to a of fluid properties and other flow coefficients are denoted with a (0 ). This leads to the following weak form of the Eulerian sensitivity equations:

Z ZX X

r  su dp dX ¼ 0 ½q0 u  ru þ qsu  ru þ qu  rsu   du dX  Z Z ¼ f 0  du dX þ t 0f  du dC

Z h X

ð14Þ

X

Z

sp ra  du dX þ X

Z

½l0 ðru þ rT uÞ þ lðrsu þ rT su Þ : rdu dX

X

ð15Þ

Ct

Z Z Z i  q0 cp þ qc0p u  rT þ qcp ðsu  rT þ u  rsT Þ dT dX þ ðk0 rT þ krsT Þ  rdT dX ¼ q0S dX þ X

X

q0 dT dC

ð16Þ

Cq

with

  t 0f ¼ ½l0 ru þ rT u þ lðrsu þ rT su Þ  n  sP n and q0 ¼ ðk0 rT þ krsT Þ  n

ð17Þ

Contrary to the flow equations, the sensitivity equations are linear. 3.2. Boundary conditions Boundary conditions for su for the Dirichlet boundary conditions and for the traction of the sensitivity t 0f for the Neumann boundary conditions are obtained with a similar approach. However, since the boundary shape depends on the parameters of interest, see Fig. 1, one must rely on the material derivative of the flow because the partial derivatives of the flow cannot be expressed directly on this boundary segment. Hence taking material derivative of Eqs. (4) and (7), we have

5992

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

Γ(α+δα)

Xf(ξ,α+δα)

Γ(α)

Xf(ξ,α)

Fig. 1. Deformation of the boundary with a variation of the shape parameter.

 Du Du ¼ on Cau Da Da Dt f Dt f ¼ on Cat Da Da DT DT ¼ on CaT Da Da  Dq Dq ¼ on Caq Da Da

ð18Þ ð19Þ ð20Þ ð21Þ

where the overbar denotes an imposed (specified) boundary quantity. This material derivatives are expanded to make terms in su ; sT ; t 0f and q0 appear and to obtain their boundary conditions. For example, for the Dirichlet boundary conditions on the velocity, we have:

 @u @u Dxf @u Dyf Du ¼ þ þ Da @ a @x Da @ a Da  DX f Du ) su ¼  ru  on Cau Da Da

ð22Þ ð23Þ

and for the temperature:

sT ¼

DX f DT  rT  Da Da

on CaT

ð24Þ

The flow boundary conditions (5) and (7) are differentiated in a similar manner to provide Neumann boundary conditions for the sensitivity. Thus, we have:

    Dt f @ s Dxf @ s Dyf @p Dxf @p Dyf Dn Dn n ns  þ þ þp on Cat @x Da @y Da Da Da Da @x Da @y Da " ! ! #    Dq @k Dxf @k Dyf @ 2 T Dxf @ 2 T Dyf @ 2 T Dxf @ 2 T Dyf 0 nx þ ny q ¼  þ rT  n  k þ þ Da @x Da @y Da @x2 Da @x@y Da @x@y Da @y2 Da

t 0f ¼

 krT 

Dn Da

on Caq

ð25Þ

ð26Þ

Here Xf = t[xf(a), yf(a)] are the coordinates of the points on the boundary. For a value parameter, the shape sensitivity of the boundary DXf/Da vanishes. However for a shape parameter, Eqs. (23), (25), (24) and (26) indicate that the boundary conditions require the accurate evaluation of the first derivatives of the flow at the boundary for the Dirichlet conditions and the second derivatives for the Neumann conditions through the terms @ s/@x and @ s/@y in (24). This is a challenging task because the accuracy of the FEM derivatives decreases near the boundary exactly where Eqs. (23) and (25) require highly accurate values. This introduces inaccuracies that degrade grid convergence of the sensitivity solution [22]. Techniques for extracting accurate flow gradients at the boundary have been developed [2] but they are complicated and costly in CPU time for 2D and the situation is exacerbated in 3D. 4. Lagrangian sensitivity equation method The Lagrangian point of view leads to simpler boundary conditions because this method uses the material derivative of the flow as an unknown. Hence expressions (18)–(21) directly provide the appropriate boundary conditions for the material/ Lagrangian sensitivities. Because the material derivatives follow the geometric changes of the domain, the solution u exhibits both explicit and implicit dependences on the shape parameter a. We introduce a mapping between a reference undeformed domain and the deformed configuration, so as to make clear the dependence of the flow solution on a. Thus we write

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

5993

u = u(x(a), a), where x(a) represents the implicit dependence on a through the mapping. We denote the Lagrangian sensitivity by:

S u ¼ lim

da!0

uðxða þ daÞ; a þ daÞ  uðxðaÞ; aÞ Du ¼ da Da

ð27Þ

Note that we use lowercase symbols (su, sp, sT) to denote Eulerian sensitivities and uppercase letters (Su, Sp, ST) for the Lagrangian sensitivities. We will suppose for the Lagrangian approach that the mapping and the domain are regular enough so that we can use the existence results of Bello et al. [8]. 4.1. General formulation of CLSEM The CLSEM equations are derived formally by implicit differentiation of the weak form of the Navier–Stokes Eqs. (8) and (9). We choose this approach rather than differentiate the PDEs directly. Differentiation of the PDEs would require the introduction of an explicit expression of the mapping from the physical parameter dependent domain to a reference fixed domain which would not depend on the parameter of interest. Thus the parameter dependence would now be born by the mapping itself. Applying the implicit function theorem on the PDEs requires an explicit algebraic expression of the mapping which would limit the scope of the resulting CLSEM. Differentiation of the weak form of the Navier–Stokes while not trivial provides room for using a mapping defined by the solution of linear elasticity equations. Such an implicit approach is very general and imposes no restriction to the development of a CLSEM. In fact such CLSEM formulation has the same scop and generality as that of the flow solver. 4.1.1. Mapping b between the reference, fixed and undeformed configuration and the actual parameter dependent geomThe mapping U etry is illustrated on Fig. 2. The Lagrangian sensitivity equations are obtained by mapping the weak form of the flow equations onto the reference fixed domain X0 which does not depend on the shape parameter a. The operation of differentiation with respect to a can be brought inside the integral on X0, so that we can proceed with the differentiation of the integrand with respect to a. Finally, the resulting weak form of the Lagrangian sensitivities is mapped back to the deformed parameter dependent domain Xa. We use superscript a to label the coordinates, domain and operators expressed on the deformed configuration. The mathematical proofs and development are standard and may be found elsewhere [5,9]. We simply summarize some of the notations associated with the mapping as follows:  The deformation tensor:

" @xa ^ aÞ ¼ Fðx; aÞ ¼ r/ðx;

@x

@xa @y

@ya @x

@ya @y

# ð28Þ

The determinant of F is the Jacobian of the transformation and is defined as follows: a a a a ^ aÞ ¼ @x @y  @y @x Jðx; aÞ ¼ det½Fðx; aÞ ¼ det½r/ðx; @x @y @x @y

ð29Þ

 The deformation velocity or material velocity is the derivative of the mapping with respect to the shape parameter a.

V M ðx; aÞ ¼

^ aÞ @ / ^ a ½/ðx; ^ aÞ @ /ðx; ¼ ¼ V aM ðxa ; aÞ @a @a

It leads to the following tensor equality:

^

Φ

(x α, y α)

^−1 Φ

(x 0, y 0)

Ωα

Ω0 Fig. 2. Domain deformation.

5994

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

2 @V Mx DF @F @x ¼ ¼ rV M ¼ 4 @V My Da @ a

@V Mx @y @V My @y

@x

3 5

ð30Þ

We now present a few useful formula for obtaining the equations of the Lagrangian sensitivities (more details are given by Delfour and Zolesio [5] or by Tortorelli and Wang [9]).  Derivative of the Jacobian of the mapping: it is written on the reference domain as follows:

@J ¼ rV M : JF 1 @a

ð31Þ

It takes the following form on the deformed domain:

DJa ¼ J a ra  V aM Da

ð32Þ

 Gradient and divergence operators: These operators involve space derivatives, therefore, their expressions depend on which domain they are expressed. Consider a given vector field defined on the deformed Xa by the representation fa(xa, a) and its corresponding representation f(x, a) on the reference domain X0. The divergence of fa and f are related by

ra  f a ¼ F 1 : rf

ð33Þ

ra f a ¼ rf  F 1

ð34Þ

 Nanson formula [5]: It is used to handle integrals along a boundary

Z

a

f ðxa ; ya Þ na dCa ¼

Ca

Z

fðx; yÞ JF T  n0 dC0

ð35Þ

C0

where f is a vector field or a scalar function.  Other useful formulas:

DF 1 ¼ ra V aM Da DJ ¼ ra  V aM J 1 Da r  ðJF 1 Þ ¼ 0 F

ð36Þ ð37Þ ð38Þ

4.1.2. Example of the continuity equation We illustrate the differentiation and mapping processes on the weak form of the continuity equation. Let dp be an admissible test function associated to the pressure. The weak form of the continuity equation is

Z

ra  udp dXa ¼ 0

Xa

1. We map this expression onto the reference domain X0 which does not depend on a.

Z

F 1 : rudp JdX0 ¼ 0

X0

2. We evaluate the material derivative with respect to a.

D Da

Z

F 1 : rudp JdX0 ¼ X0

Z X0

D 1 ðF : rudp JÞdX0 ¼ 0 Da

ð39Þ

which yields

Z

" F

X0

1

# ! Du DF 1 DJ 1 dX0 ¼ 0 :r þ : ru dpJ þ F : rudp Da Da Da

3. We map back to the deformed domain and obtain:

Z Xa





ra  S u  ra V aM : ra u þ ra  ura  V aM dp dXa ¼ 0

ð40Þ

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

5995

b a which is the conThis operation introduces the so-called domain deformation velocity: V M ¼ ð@xa =@ a; @ya =@ aÞ ¼ @ U=@ tinuous analogue of mesh sensitivities encountered in discrete approaches. 4.1.3. Lagrangian sensitivity equations The method is applied in the same manner to the momentum and energy equations. The detailed expressions are given in Appendix A. The deformation velocity also appears in the momentum and energy sensitivity equations. Hence we need to know VM everywhere in the domain to solve the sensitivity equations. However it is known only on the parameter dependent part of the boundary. Thus the boundary deformation velocity must be extended to the whole domain to provide the required geometric sensitivities. There are many ways to achieve this extension. As mentioned by Plotnikov et al. [7], the shape sensitivity solution depends of the definition of the mapping only at the boundary. Hence the definition of VM over Xa is not unique. To ensure convergence and regularity of the solution, the derivatives of this velocity must be continuous and its evaluation must be compatible with our mesh adaptation. We use the pseudo-solid approach which generates the domain deformation by solving linear elasticity equations [21]:

r

  1 kps trðrV M þ rT V M Þ þ lps ðrV M þ rT V M Þ ¼ 0 in X 2

subject to the boundary conditions

VM ¼ VM

on @ X

where kps and lps are the Lamé coefficient of the pseudo-solid and V M is a given deformation velocity. The pseudo-solid approach is very general and allows to efficiently treat complex geometries (in fact as complex and complicated as can be treated by the most advanced FEM). 4.2. Boundary conditions The CLSEM leads to simple boundary conditions by direct differentiation of the flow boundary conditions.

 Du Du ¼ Da Da DT DT ¼ ST ¼ Da Da

Su ¼

on Cau

ð41Þ

on CaT

ð42Þ

Similarly for Neumann boundary conditions we have:

 Dt f þ r  V M t f  t f ðrV M  nÞ  n  du dCa Da Ca Cat  Z t Z   D Dq a q ðrV M  nÞ  n dT dCa krT  n dST dCa ¼ þ r  VMq Da Caq Caq Da D Da

Z

ðpna þ s  na Þ  du dCa ¼

Z



ð43Þ ð44Þ

These boundary conditions are simpler than those of the Eulerian formulation because the material derivatives follow the deformation of the boundary segments and they do not require evaluation of flow gradient at the boundary. From this perspective, we expect the Lagrangian formulation to be free of the convergence and accuracy problems encountered in the CESEM in the case of shape parameters. Recall that the CESEM requires very fine meshes near parameter dependent boundary segment to ensure accuracy of the flow gradient at the boundary. We will show that the CLSEM performs better and more accurately on coarser meshes than the CESEM on finer grids. 4.3. CESEM vs CLSEM Now we compare the CLSEM with the CESEM. Both methods provide enough sensitivity information to evaluate the total derivatives of an output functional with respect to a shape parameter. For the CESEM, the weak form is easily obtained while more tedious mathematical developments are required for the CLSEM. Furthermore the use of elasticity equations to define the deformation velocity field increases the size of the discrete system. In the Eulerian approach, one must solve for the Eulerian sensitivities of the velocity (su, sv) = (@u/@ a, @v/ @ a) and pressure sp in the domain X, whereas in the Lagrangian method, the unknowns are the Lagrangian sensitivities of the velocity (Su, Sv) = (Du/Da, Dv/da), of the pressure Sp and of the deformation velocity field (VMx, VMy). While the boundary conditions required to solve the Lagrangian sensitivity problem are exact and trivial, those for the Eulerian sensitivities are approximate in the sense that they necessitate the evaluation of the flow derivatives at the boundary precisely where the FEM derivatives are least accurate. In other words, the big issue for the Lagrangian formulation is the evaluation of the deformation velocity while accurate evaluation of the flow gradients at the boundary is the main challenge for the Eulerian formulation. We will now determine if the additional numerical cost of the CLSEM is counterparted by improvement of the sensitivity results.

5996

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

5. Finite element formulation and adaptive remeshing The flow and sensitivity equations are solved by a finite element method using the 6 nodes Taylor–Hood element. The discretization of the velocity, temperature and deformation fields is quadratic and that of the pressure is linear. The sensitivity fields are discretized in the same manner as the flow field. The resulting system of non-linear algebraic equations for the flow is linearized by Newton’s method and solved using a sparse direct solver. The flow data are then used to build the sensitivity system. Note that this system is linear with respect to the sensitivity unknowns. Mesh independent solutions of the equations are obtained with an adaptive finite element algorithm [23]. The adaptive methodology is set to reduce the error on each variable by a factor of two at each cycle of adaptation. On any given element of the current mesh, error estimates are obtained by a local least-squares reconstruction [24] of the derivatives of the flow and their sensitivities. Ideal element size is determined for each variable and the smallest value is assigned to this element. This process is repeated for all elements of the grid to define the grid density function for the next adaptive mesh. Using this information a new improved mesh is produced by an advancing front mesh generator. 6. Numerical results : verification The code is first verified with the method of manufactured solutions [25]. The verification aims at checking that the implemented equations are correctly solved. This is achieved by comparing an analytical solution to the finite element one. The error can be evaluated and the convergence rate determined through a grid refinement study. The analytical solution is first choosen. It must be divergence free. Then an appropriate source term is defined in the momentum equations to ensure equilibrium and satisfaction the Navier–Stokes equations augmented with the source term. Such a manufactured solution serves as a closed form solution. A similar approach is used for the shape sensitivities We study and compare three variants of the sensitivity equation method:  The classical Eulerian formulation that we term CESEM-1.  The Eulerian formulation using the technique proposed by Duvigneau [2] to reconstruct more accurate values of the boundary conditions for shape parameters. We refer to this technique as CESEM-2.  The present Lagrangian sensitivity equation method (CLSEM). 6.1. Analytical solution We choose a flow solution complex enough to ensure that all terms in the Navier–Stokes and sensitivity equations are exercized. Our solution is:

uðx; y; aÞ ¼ 4x2 expðax2 yÞð1  2ax2 yÞ

v ðx; y; aÞ ¼ 8xy expðax2 yÞð1  2ax2 yÞ Dð1  yÞ @v

ð45Þ

þ 2l 1 þ Cðx2 þ y2 Þ @y   x Tðx; y; aÞ ¼ exp c dð1 þ byÞ

pðx; y; aÞ ¼

where a, C, c, d and b are arbitrary coefficients. The domain and its deformation are presented in Fig. 3. Points on the curved boundary satisfy the implicit equation 2ax2y = 1 in which a is our shape parameter set to a baseline value of 50. Neumann conditions are prescribed on the curved boundary, Dirichlet conditions elsewhere. Our manufactured solution is substituted into the Navier–Stokes Eqs. (1)–(3) to determine the source terms f and qs that ensure equilibrium. Manufactured solution fields for the Eulerian sensitivities are easily obtained by direct partial differentiation of the flow solution (45) with respect to a. The partial derivative of f in Eq. (2) yields the required source term f0 in Eq. (12). The manufactured Lagrangian sensitivity, that is the analytical solution for the Lagrangian sensitivity problem, fields are evaluated with the chain rule:

Su ¼

@u þ ru  V M ; @a

Sp ¼

@p þ rp  V M @a

and ST ¼

@T þ rT  V M @a

ð46Þ

with VM sufficiently smooth. Since the Lagrangian sensitivities depend on the deformation velocity which itself depends on the pseudo-solid solution, we must either find an appropriate manufactured solution for the sensitivities of pseudo-solid displacement to get the required deformation velocity or hard code a known admissible deformation velocity VM into the solver. For simplicity we opt to determine VM analytically and to impose it to ensure that the computed Lagrangian flow sensitivity converge to the manufactured solution with mesh refinement. The required source term in the sensitivity equations is computed in the code using f0 , @f/@x and @f/@y, the computed VM and the chain rule relation

Df @f @f 0 ¼ f þ V Mx þ V My Da @x @y

ð47Þ

5997

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

u(x,y) Su(x,y)

(−0.1 ; 1)

Outlet

(0 ; 1)

Curve 2ax²y=1

u(x,y) Su(x,y)

Wall t(x,y) St(x,y)

Ω

Ω

(−0.316 ; 0.1) u(x,y) Su(x,y)

Inlet (−0.316 ; 0)

Symmetry

O

u(x,y) Su(x,y)

Fig. 3. Domain and deformation.

6.2. Grid refinement study A grid refinement study is carried out to evaluate the performance of each formulation in terms of convergence rate. At each cycle, the errors are estimated and the mesh is refined where needed. We will compute and compare two types of error for each field.  The true error: the finite element solution, denoted with the superscript h is compared to the exact solution field denoted h exa with the superscriptexa. For example for the Lagrangian sensitivity of the velocity field eexa S u ¼ jS u  S u j  The error estimate: in general, the exact solution is not available. It is estimated with a projection reconstruction [26] h  denoted with the superscript ⁄. For the Lagrangian sensitivity of the velocity field eest S u ¼ jS u  S u j Checking and verifying the accuracy and reliability of the error estimator against the true error is important in practice since for real cases, the true error is not available and its estimator is the sole tool that we have to drive the adaptive process and simultaneously assess the mesh convergence of the solution. Error estimates are evaluated in energy norm for both the velocity and its sensitivity fields and in H1 semi norm for the pressure and temperature and their sensitivities:

Z    keu k2energie ¼ reu þ rT eu : reu þ reTu dX X Z keP k2H1 ¼ reP  reP dX ZX reT  reT dX keT k2H1 ¼

ð48Þ ð49Þ ð50Þ

X

where eu, eP and eT are the local error or its estimation for the velocity, pressure and temperature. Error norms for the sensitivity solution are computed in the same manner. The expected rates of convergence are given by:

      1 1 1 jesu jenergie ¼ O jeSu jenergie ¼ O jeu jenergie ¼ O Nnd N nd N nd

5998

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

1=2

 jeT jH1 ¼ O

!

1

jep jH1 ¼ O

Nnd

1 Nnd

jesp jH1 ¼ O

 jesT jH1

1

!

1=2

jeSp jH1 ¼ O

1

!

1=2

N Nnd  nd   1 1 jeST jH1 ¼ O ¼O Nnd Nnd

where Nnd is the number of nodes in the mesh. The convergence rate should be the same with the estimate or true errors. 6.3. Results The flow and the sensitivity equations are solved for the three variants of the sensitivity equation method (CESEM-1, CESEM-2, CLSEM). In each case, the adaptive mesh refinement is driven by the velocity and pressure error estimates in energy norm and H1P norm for the flow and the sensitivities. The efficiency index (ratio of estimate to true error) is also computed. Figs. 4 and 5 present a complete portrait of the convergence behavior of the three approaches to compute flow sensitivities for shape parameters. In addition to the true errors for the pressure, velocity and temperature sensitivity, Figs. 4 and 5 include the convergence histories of the error estimators for all variables and the efficiency indexes. We can make the following observations:  For the CESEM-1 (Figs. 4(a) and 5(a)), the pressure sensitivity error does not converge with the mesh refinement and the velocity error converges at a suboptimal rate of 1.14 instead of 2. The estimators never asymptotically approach the true error. This illustrates how approximate boundary conditions can cause inaccuracies of the sensitivity solution in the whole domain. The flow gradients evaluated at the boundary are not accurate enough to ensure the convergence of the solution.  For the CESEM-2 (Figs. 4(b) and 5(b)) very fine meshes are needed to achieve the required accuracy. This is due to the fact that the reconstruction of the gradient requires very fine mesh at the boundary to compute sufficiently accurate boundary flow gradients. Therefore, for further applications, only fine meshes should be used.  The CLSEM (Figs. 4(c) and 5(c)) behaves according to theory: the convergence rate is 2 and the efficiency index tends toward 1, the error estimator is accurate. Note that the Eulerian and Lagrangian sensitivities are two different fields that cannot be compared directly: we cannot determine which is more accurate by simply looking at the absolute error. Fig. 6 represents the norm of the exact relative h exa error defined for a field u as erel j=juexa j, for the H1P and energy norm. Note that we do not present the relative u ¼ ju  u error for the temperature sensitivity field because the analytical Eulerian sensitivity of the temperature is zero over the domain while the Lagrangian one is not, the comparison would not be relevant. For the velocity and pressure sensitivity fields, the relative error is smaller with the Lagrangian formulation: the error for CESEM-1 on a mesh with 100,000 nodes is the same as the error for the CLSEM at 1000 nodes. CESEM-2 yields at the same error levels than the CLSEM but only for fine meshes. The CLSEM is sufficiently accurate that sensitivity data obtained on meshes that would be considered much too coarse for the CESEM can be used safely and reliably. The implementation of the Lagrangian sensitivity equation is verified and the comparison between the two methods indicates that the Lagrangian approach yields better accuracy for the computed sensitivity. As previously mentioned, the Lagrangian and Eulerian sensitivity fields are different. The interpretation will also be different. We will now validate the Lagrangian approach. 7. Application: flows on nearby geometries We consider the flow around airfoils at a Reynolds number of 1000 at an incidence of 5°. We use the sensitivity data to evaluate flows on nearby geometries as a test case to compare the accuracy of the CESEM and CLSEM. We consider NACA 4digit airfoils parametrized with the three standard parameters: maximum thickness (e), vertical displacement at the point of maximum camber (m) and maximum camber location along the chord the vertical displacement at the point of maximum camber (p). The shape of the mean line is expressed by two tangent parabolic arcs:

m ð2px  x2 Þ 0 6 x 6 xp p2 m ¼ ð1  2p þ 2px  x2 Þ xp 6 x 6 c ð1  pÞ2

yml ¼ yml

where c is the chord , m  c the maximum mean line ordinate and xp = p  c the position of the maximum mean line ordinate The thickness is given by:

pffiffiffi yth ¼ 5eð0:2969 x  0:1260x  0:3537x2 þ 0:2843x3  0:1015x4 Þ Our baseline airfoil is the NACA 4512. Sensitivities are computed with respect to the airfoils 3 parameters, e, m and p. The nearby geometries are shown Fig. 7 and the shape parameters are listed in Table 1.

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

1

5999

True H1P Estimate H1P True energy Estimate energy True H1T Estimate H1T

0.1 0.01 0.001 0.0001 1e-05 1e-06 1e-07 10

100

1000

10000

100000

1e+06

number of nodes

100 10 1 0.1 0.01 0.001 0.0001 1e-05 1e-06 1e-07 1e-08

True H1P Estimate H1P True energy Estimate energy True H1T Estimate H1T

10

100

1000

10000

100000

number of nodes

1

True H1P Estimate H1P True energy Estimate energy True H1T Estimate H1T

0.1 0.01 0.001 0.0001 1e-05 1e-06 10

100

1000

10000

100000

1e+06

number of nodes

Fig. 4. Estimate and exact energy norm error and H1P norm error for Eulerian and Lagrangian sensitivity.

The NACA 4515, 5512 and 4412 may be viewed as neighboring geometries of the NACA 4512 because only one parameter changes. However, the NACA 9714 and 6314 are not. A first order Taylor series in parameter space yields estimates of the flow around nearby geometries. Let a0 denote the baseline value of a = (e, m or p), Da its perturbation and w will be either a local surface quantity (pressure or skin friction distribution), an integral quantity (Drag, lift and moment coefficient) or a dependent variable. The extrapolated field is given by:

wða0 þ DaÞ  wða0 Þ þ

X Dw Dai ða Þ D ai 0i i

ð51Þ

Consistent Taylor series are constructed using sensitivity information from CESEM-1, CESEM-2 and the proposed CLSEM. The total derivatives Dw/Dai arise naturally in the Lagrangian approach because they are the Lagrangian sensitivities, whereas if

6000

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

0.6

Efficiency-H1P Efficiency-energy Efficiency-H1T

0.5 0.4 0.3 0.2 0.1 0 10

100

1000

10000

100000

1e+06

number of nodes

1

Efficiency-H1P Efficiency-energy Efficiency-H1T

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 10

100

1000 10000 number of nodes

100000

1.2

Efficiency-H1P Efficiency-energy Efficiency-H1T

1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 10

100

1000

10000

100000

1e+06

number of nodes

Fig. 5. Efficiency index for Eulerian and Lagrangian sensitivity.

one is using the CESEM, one must reconstruct the deformation velocity on the parameter dependent boundary @Xf/@ ai to obtain the required material derivatives.

Dw @w @w @xf @w @yf ¼ þ þ Dai @ a @x @ a @y @ a

7.1. Sensitivities fields We compute the sensitivities with respect to each shape parameter e, m and p on the baseline configuration NACA 4512. Fig. 8 shows the contours of the horizontal velocity and of its Eulerian and Lagrangian sensitivity with respect to e. As can be seen, the Eulerian and Lagrangian fields are very different. The contour of the Eulerian sensitivity are clustered highly close to the airfoil. The Eulerian sensitivity of the velocity is the partial derivative of the velocity with respect to the shape parameter.

6001

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

10

1

CLSEM CESEM CESEM-improved BC

CLSEM CESEM CESEM-improved BC

0.1

1 0.01 0.1 0.001 0.01

0.0001 10

100

1000 10000 100000 1e+06 number of nodes

10

100

1000 10000 100000 1e+06 number of nodes

Fig. 6. True relative error for Eulerian and Lagrangian sensitivity.

Fig. 7. Geometry of NACA airfoils.

Table 1 Parameters of NACA airfoils. NACA

4512

4515

5512

4412

9714

6314

e m p

0.12 0.04 0.50

0.15 0.04 0.50

0.12 0.05 0.50

0.12 0.04 0.40

0.14 0.09 0.70

0.14 0.06 0.30

Its value is negative close to the airfoil. If the thickness of the airfoil increases, points close to the profile with non zero velocity will be closer, so their velocity will significantly decrease. The approach with the Lagrangian sensitivity must be different: the deformation must be taken into account. The Lagrangian sensitivity of the velocity is zero on the airfoil because the material point will stay on the profile whatever the deformation is and its velocity will still be null. The Lagrangian sensitivity field is more uniform and provide information all over the domain. Thus, the negative sensitivity near the leading edge indicates that if the thickness of the profile increase, the blocking effect in this area will be higher. The negative sensitivity in the wake of the profile indicates that if the thickness increases the recirculation zone will be larger. Such information can be extracted from the Eulerian sensitivity but not so clearly.

7.2. Grid refinement study A grid refinement study is carried out in order to check the convergence and to compare the performance of the three formulations. We will focus on the error norm of the sensitivity fields and then on the aerodynamic coefficient.

7.2.1. Sensitivity field Mesh adaptation is driven by the energy norm for the flow and sensitivity fields. Fig. 9 shows the evolution of the error in energy norm with mesh refinement for the flow and sensitivities. For the CESEM, see Fig. 9(a), the convergence of the error estimate appears good. However we have seen in the previous section that for the CESEM, the error estimator was not accurate so the true error might be higher. The behavior of the CESEM-2 Fig. 9(b) is not smooth and presents irregularities, so that it is difficult to determine whether the solutions are in the asymptotic range of the discretization. Moreover, the convergence rate degrades significantly for grids with more than 100,000 nodes. For the CLSEM, convergence is as good as observed with

6002

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

SU_E -21.66

-16.991

-12.322

-7.6534

-2.9844

1.6845

SU_E -6.0049

-4.9734

-3.942

-2.9106

-1.8791

-0.84766

0.18378

1.2152

Fig. 8. Contours of horizontal velocity and sensitivity with respect to e – NACA4512.

the manufactured solution. However serious degradation occurs for the convergence of sensitivities with respect to m and p in the last few cycles of mesh adaptation. An analysis of the estimated elemental error norms reveal that the maximum error occurs at the trailing edge. Furthermore the deformation velocity exhibits a singular point. Fig. 10(a) presents the projection on the vertical axis of the deformation velocity for the parameter e, which is smooth and for the parameter m which is singular. The spatial derivatives of the sensitivities are not continuous at the trailing edge which causes the convergence degradation observed in the Lagrangian sensitivities. This phenomenon is intrinsic to the sharp trailing edge of the airfoil. The fluid domain and the mapping do not meet the regularity criteria to ensure existence of the sensitivity solution [8,7]. The effect of the singularity can be mitigated somewhat by stiffening the pseudo-solid (increasing the Young modulus E) near the airfoil, so that the pseudo-solid deforms less. Convergence results with this approach are shown on Fig. 9(d). The results with a stiffened pseudo-solid show improvements compared to the approach with a constant Young modulus. It is quite probable that we have only postponed the onset of the degradation to finer grids on which we would be faced with the same problem. Fig. 10(b) shows that the deformation velocity for m is still singular. An alternative could be to find another smoother deformation velocity so as not to affect grid convergence of the sensitivity fields. 7.2.2. Sensitivities of the lift and drag coefficient In many instances one is interested in output functionals such as drag and lift coefficients CD and CL. 7.2.2.1. Influence of the mesh adaptation. We consider the following approaches:     

CESEM-1 with mesh adaptation on the energy norm of the flow field only. CESEM-1 with mesh adaptation on the energy norm of the flow and sensitivities fields. CESEM-2 with mesh adaptation on the energy norm of the flow and sensitivities fields. CLSEM with mesh adaptation on the energy norm of the flow field only. CLSEM with a variable Young modulus with mesh adaptation on the energy norm of the flow and sensitivities fields.

Fig. 11 indicates that, for the 3 parameters, grid convergence of the sensitivities of drag and lift coefficient is better for the CLSEM than for the CESEMs. In fact, for the Lagrangian formulation, DCD/Dai and DCL/Dai on the grids containing 10,000 nodes differ by less than 1% from the values obtained on grids of 100,000. The difference are approximately 5 times larger

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

0.1

Flow e m p

0.01

1

6003

Flow e m p

0.1 0.01

0.001 0.001 0.0001 1e-05 100

0.0001 1000 10000 100000 1e+06 number of nodes

0.1

Flow e m p

0.01

1e-05 100

0.1

0.001

0.0001

0.0001

1000

10000 100000 1e+06

number of nodes

Flow e m p

0.01

0.001

1e-05 100

1000 10000 100000 1e+06 number of nodes

1e-05 100

1000

10000 100000 1e+06 number of nodes

Fig. 9. Convergence of the flow and sensitivities for the NACA airfoil in energy norm (absolute error).

Fig. 10. Vertical deformation velocity at x = 1 for the parameter e and m.

for the CESEM-1. Note also that the CESEM-2 yields inaccurate results on meshes of less than 60,000 nodes. In general, for the same level of accuracy, the CLSEM will require coarser meshes than CESEMs. For the CESEM, grid adaptation on the flow and its sensitivities leads in general to faster grid convergence, see for example Fig. 11(a). This is to be expected: the boundary conditions for the sensitivity depend on the boundary gradient of the velocity. Thus finer grids at the boundary improve

6004

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

CESEM (flow) CESEM (flow+sens) CESEM with improved BC (flow+sens) CLSEM (flow) CLSEM−variable E (flow+sens)

0.215 0.21 0.205

−2.4

−2.5 −2.55

0.2

−2.6

0.195

−2.65

0.19

−2.7

0.185

−2.75

0.18

−2.8

0.175 0.17

CESEM (flow) CESEM (flow+sens) CESEM with improved BC (flow+sens) CLSEM (flow) CLSEM−variable E (flow+sens)

−2.45

−2.85 4

10

5

number of nodes

CESEM (flow) CESEM (flow+sens) CESEM with improved BC (flow+sens) CLSEM (flow) CLSEM−variable E (flow+sens)

0.28

4

10

0.27

10

1.5

5

number of nodes

10

CESEM (flow) CESEM (flow+sens) CESEM with improved BC (flow+sens) CLSEM (flow) CLSEM−variable E (flow+sens)

1.45 1.4 1.35 1.3

0.26

1.25 1.2

0.25

1.15 1.1

0.24

1.05 0.23

4

10

−0.01

5

number of nodes

10

CESEM (flow) CESEM (flow+sens) CESEM with improved BC (flow+sens) CLSEM (flow) CLSEM−variable E (flow+sens)

−0.0102 −0.0104 −0.0106

1

4

10

0.06

5

number of nodes

10

CESEM (flow) CESEM (flow+sens) CESEM with improved BC (flow+sens) CLSEM (flow) CLSEM−variable E (flow+sens)

0.055 0.05 0.045

−0.0108 −0.011

0.04

−0.0112

0.035

−0.0114 0.03 −0.0116 0.025

−0.0118 −0.012

4

10

5

number of nodes

10

0.02

4

10

5

number of nodes

10

Fig. 11. Convergence of the sensitivities of the drag and lift coefficients.

accuracy but increase cost. The CLSEM does not exhibit this behavior. Indeed it appears to produce grid converged sensitivities for the drag and lift coefficient with approximately 5000–10,000 nodes while the CESEM require 100,000 nodes.

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

6005

We now look at the effects of adaptive strategy on the convergence of each formulation. We can now analyze the behavior of each formulation with the mesh adaptation. For the CESEM, the mesh adaptation has consequences on the convergence. The sensitivities of the drag and lift coefficients seem to convergence faster when mesh adaptation is driven by both sensitivity and flow error estimates. This is obvious on Fig. 11(a). The Eulerian sensitivities require careful treatment of the mesh. However the CLSEM appears to yield nearly grid converged sensitivities on meshes with as few as 5000 nodes whatever the mesh adaptation strategy. To illustrate this, Fig. 12 shows details of the adaptive grids obtained at the last cycles for the computation with mesh adaptation on flow only, on flow and Lagrangian sensitivities and finally on the flow and the Eulerian sensitivities. Fig. 12(a) and (b) are very similar, indicating that the grid requirements of the CLSEM are close to those of the flow alone whereas Fig. 12(a) and (c) show that the CESEM required finer mesh close to the boundary. This is due to the fact that the sensitivity gradients are largest near the airfoil because Eulerian sensitivities are always large near the parameter dependent boundary. On the one hand, the rapid variations of the Eulerian sensitivities near the airfoil force the adaptive procedure to cluster the grid near the airfoil surface. On the other hand, the Lagrangian formulation requires fewer grid points in the area so that the mesh can be refined elsewhere in the domain. It also explains why the Lagrangian sensitivities are more accurate. 7.2.2.2. Effect of the deformation velocity for the CLSEM. Fig. 13 compares the effects of a variable Young modulus on the singular deformation velocity and on the convergence of sensitivities of the drag and lift coefficients. Grid convergence is assessed with the Least-Squares Richardson extrapolation to zero mesh size due to Eça and Hoekstra [27]. Results on Fig. 13 are obtained by least squares on the five finest grids (grids with 5000 to 130,000 nodes). Note that changing the distribution of the Young modulus (constant or variable values) will lead to different meshes. Here the constant Young modulus approach yields oscillatory grid convergence in some cases so that we could not compute the Richardson extrapolation. The Table 2 summarizes the results and includes the numerical uncertainty estimates for the finer mesh, the sensitivities at a zero size mesh and the convergence rate obtained by the Least-Squares procedure. It indicates that the variable Young modulus technique is an effective measure to control the detrimental effects arising from the sharp trailing edge which acts as a rentrant corner in the flow domain. For the parameter e, the values of the sensitivity of the drag and lift coefficients differ on

Fig. 12. Grids.

6006

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011 0.21

variable E (extrapolation) variable E fixed E

0.208

−2.62

0.206

−2.66

0.204

−2.68

0.202

−2.7

0.2

−2.72

0.198

−2.74

0.196

−2.76

0.194

4

10

0.2635

5

number of nodes

10

variable E (extrapolation) variable E fixed E

0.263

variable E (extrapolation) variable E fixed E (extrapolation) fixed E

−2.64

−2.78

4

10

1.31

5

number of nodes

10

variable E (extrapolation) variable E fixed E (extrapolation) fixed E

1.3

0.2625 1.29

0.262 0.2615

1.28

0.261

1.27

0.2605

1.26

0.26 1.25

0.2595 0.259

4

10

5

number of nodes

10

−0.0107

1.24

4

10

5

number of nodes

10

0.05

−0.0107 0.045

−0.0108 −0.0108

0.04

−0.0109 −0.0109

0.035

−0.011 0.03

−0.011 variable E (extrapolation) variable E fixed E (extrapolation) fixed E

−0.0111 −0.0111 4

10

5

number of nodes

10

variable E (extrapolation) variable E fixed E (extrapolation) fixed E

0.025

0.02

4

10

5

number of nodes

10

Fig. 13. Effect of the Young modulus on the convergence of the sensitivities of the drag ans lift coefficients.

the coarsest meshes depending on whether a constant or variable Young modulus is used. However as the grid is refined, this two distributions of the Young modulus lead to identic values of the sensitivities. We see from Table 2 that the difference between the 2 extrapolated values is smaller than the estimated uncertainties for the lift coefficient. This shows that the final value of the sensitivity of the drag and lift coefficients do not depend on the deformation velocity. But for the parameters m and p, the results are different. For the fixed Young modulus the singular point in the deformation velocity degrades the

6007

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011 Table 2 Evaluation of the drag and lift coefficient. Variable Young modulus

Fixed Young modulus

DCd/De

Estimate Uncertainty Rate

0.195719724 6.6E05 3.7

DCl/De

Estimate Uncertainty Rate

2.77528769 4.8E04 4.1

DCd/Dm

Estimate Uncertainty Rate

0.261186139 9.9E07 1.4

DCl/Dm

Estimate Uncertainty Rate

1.26892727 2.8E04 2.9

1.27096733 1.8E03 1.3

DCd/Dp

Estimate Uncertainty Rate

0.0107968593 1.6E06 2.7

0.0107981736 2.0E06 2.4

DCl/Dp

Estimate Uncertainty rate

0.0469876808 4.6E06 4.4

0.0467386288 2.8E05 3.0

2.77530323 1.6E04 3.9

convergence of the sensitivity of CD and CL. Hence the uncertainties are larger and the convergence rate lower. Thus by controlling the Young modulus, we can find a deformation velocity that ensures convergence of the sensitivities of the output functional. 7.3. Fast evaluation of nearby flows We now apply the CESEM and CLSEM to fast evaluation of flows on nearby geometries. Here again, we will show other advantages of the Lagrangian formulation. 7.3.1. Extrapolation of the flow Because the deformation velocity is part of the solution and is computed consistently, the CLSEM provides Taylor series approximations of the flow around nearby airfoil geometries, the domain, the mesh and the flow field. In the following example, we use the flow around a NACA 4512 (Re = 1000, a = 5°) as a baseline flow to extrapolate this flow to obtain an approximation of the flow around a NACA 9714. Calculations were performed on an adaptive grid of 10,000 nodes as it was deemed to yield sufficient accuracy in the previous section. Fig. 14(a) and (b) show contours of the pressure and the u-velocity on the baseline airfoil. Fig. 14(c) and (d) present extrapolation of the baseline flow to the target airfoil and Fig. 14(e) and (f) are the solution obtained by reanalysis of the flow on the target airfoil. Contours have been set to the same values for the 3 flows to ease comparison. We note that the variations of the design parameter are significant (16% for e, 125% for m and 40% for p). As can be seen, the extrapolation from the NACA 4512 captures the flow features of the NACA 9714 well. However, comparisons of Fig. 14(c) to (e) and (d) to (f) show that there is a more pronounced difference between the extrapolated and true geometry than between the extrapolated and true contours of u on the NACA 9714. This probably indicates that while the flow field exhibit an overall fairly linear dependence on (e, m, p), the geometry has a stronger non linear dependence on p. As can be seen, the extrapolation of the geometry is not exact because it is not linear with respect to the parameter p, this can explain the difference for the contours of the solution. The perturbation of the design parameter may be too high to stay in a linear range. 7.3.2. Aerodynamics coefficients We now turn our attention to the evaluation of aerodynamic coefficients for a range of nearby airfoils. Six cycles of grid adaptation (on flow only) were carried out and yield a mesh with approximated 40,000 nodes. Drag and lift coefficients and their Lagrangian sensitivities are evaluated on the baseline airfoil. We then use formula (51) to compute the extrapolation. Flow reanalysis on the real geometry provides the true or exact values of CL and CD. Table 3 reports computed and extrapolated coefficients and the relative error. The agreement between extrapolation and reanalysis is very good for the neighboring airfoils. When the 3 parameters change, the extrapolated coefficient are less accurate. The relative error are very close to those obtained by Pelletier et al. [1] using CESEM-2 which required grids with more than 80,000 nodes. This confirms the cost-effectiveness of the CLSEM. An adaptive grid refinement study is carried out to compare the performances of the 3 formulations and assess both the grid convergence and accuracy of the extrapolation of drag and lift coefficients to nearby geometries. The mesh is adapted on the flow error only to ensure that all formulations are evaluated on the same grids. Fig. 15 presents the results from our

6008

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

Fig. 14. Pressure and velocity contours.

Table 3 Evaluation of the drag and lift coefficient. NACA

Computed coefficient

Extrapolated coefficient

Relative error (%)

0.14110 0.19569

0.14067 0.19694

0.30 0.64

CD CL

0.13766 0.29250

0.13740 0.29287

0.18 0.13

4412

CD CL

0.13595 0.27611

0.13588 0.27550

0.05 0.22

9714

CD CL

0.15197 0.30895

0.14960 0.29751

1.56 3.70

6314

CD CL

0.14856 0.22339

0.14609 0.24071

1.66 7.74

4512

CD CL

0.13479 0.28018

4515

CD CL

5512

study. The relative errors on CL and CD are plotted against the number of nodes in the adaptive grids. As can be seen the proposed CLSEM yields the most accurate results and does so on coarser meshes than any of the 2 CESEMs. We note that the relative errors of the extrapolations reach a non zero plateau as the mesh is refined. This is due to the fact that the error of the extrapolation has two contributions arising from the spatial discretization of the sensitivity PDEs

6009

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011 2.5

2

CLSEM CESEM CESEM with improved BC

4 3.5 Relative error (%)

Relative error (%)

4.5

CLSEM CESEM CESEM with improved BC

1.5

1

3 2.5 2 1.5 1

0.5

0.5 0 0

1

2

3

12

7

8

0

1

2

3

14

Relative error (%)

6

4

2

4 5 6 Number of nodes

7

8

9 4

x 10

CLSEM CESEM CESEM with improved BC

12

8

0

0

9 4

x 10

CLSEM CESEM CESEM with improved BC

10

Relative error (%)

4 5 6 Number of nodes

10

8

6

4

0

1

2

3

4 5 6 Number of nodes

7

8

9 4

x 10

2

0

1

2

3

4 5 6 Number of nodes

7

8

9 4

x 10

Fig. 15. Grid convergence: relative error between exact and extrapolated aerodynamic coefficient (%).

and arising also from the use of a Taylor series truncated to the first order. Thus the level of the plateau observed with the mesh refinement is a measure of accuracy of the Taylor series. For instance if CL depends linearly on a parameter a or if itnac exhibits a weak non linear dependence on a, the plateau should have a very small value whereas if the dependence were highly nonlinear, the value of the plateau would be higher. Generally speaking, the CLSEM delivers more accuracy at lower cost than the CESEMs.

8. Conclusion A continuous Lagrangian sensitivity equation method has been presented for incompressible viscous flows. The Lagrangian approach offers advantages over the Eulerian sensitivity equation method. First, the boundary conditions are simpler in the Lagrangian case than in the Eulerian: flow gradients, first and second derivatives are not required on the parameter depend boundary. The Lagrangian sensitivity equation method delivers more accurate shape sensitivities at a lower cost than the Eulerian technique. The parameterization of the geometry of the domain is defined by elasticity equations to manage the geometry and to provide a consistent deformation velocity over the fluid domain. This pseudo-solid approach is generic, general and applicable to very complex geometries and arbitrary parameterization. However, the CLSEM requires the solution of additional partial differential equations to obtain the deformation velocity, that is two additional unknowns in the finite element system. The CLSEM was described in details and compared to the CESEM. Code verification by the method of manufactured solutions showed that the finite element implementation of the approach is correct. It also highlights difference in behavior compared to the two CESEMs. It reveals that CLSEM yields better grid convergence than those of the CESEMs. It also indicates that the CLSEM yields more accuracy at a lower cost than CESEMs, even if more PDEs must be solved.

6010

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

The CLSEM was applied to fast evaluation of flows around airfoil geometries nearby that of the baseline flow. Results on the airfoils problem confirm the observation made on the manufactured solution. The advantages of the CLSEM apply equally well to extrapolation of the flow as well as the lift and drag coefficients (output functionals). On coarse meshes, the Lagrangian sensitivity data are accurate enough to be used in further applications such as fast evaluation of nearby flows, uncertainty analysis or optimization. Appendix A. Detailed Lagrangian sensitivity equations A.1. Momentum equations The momentum Eq. (2) has the following expression:

qðu  rÞu ¼ rp þ r  s þ f Let dSu be an admissible function. The weak form of this Eq. (9) is

Z



u  ra u  dS u dXa 

Z

Xa

pra  dS u dXa þ

Xa

Z

s : ra dS u dXa ¼

Xa

Z

f  dS u dXa dCa

Xa

We assume of simplicity that the border term vanishes. We can identity each term:

Convection þ Pressure þ Diffusion ¼ Force Thus, the weak form of the Lagrangian sensitivity equation is the following

D Convection D Pressure D Diffusion D Force þ þ ¼ Da Da Da Da Each term is derived using the method explained in Section 4. The resulting terms are A.1.1. Derivation of the convection term

D Da

Z





q u  ra u  dS u dXa ¼ Xa

  @q þ qra  V aM ðu  ra Þu  dS u dXa @a Xa Z  þ q S u  ra u  ra u  ra V aM  u þ ðu  ra ÞS u  dS u dXa Z

ðA:1Þ

Xa

A.1.2. Derivation of the pressure term

D Da

Z

pra  dS u dXa ¼  Xa

Z Xa



Sp þ pra  V aM ra  dS u  pra V aM : ra dS u dXa

ðA:2Þ

A.1.3. Derivation of the diffusion term

D Da

Z

  Z   Dl þ lra  V aM ðra u þ raT uÞ : ra dS u dXa þ l ra S u þ raT S u Da X X Z a  Z a     a a a a a a a T a l r ur V M þ r ur V M : r dS u dXa  l ra V aM  ra u þ raT u : r dS u dXa 

s : ra dS u dXa ¼ Xa

Z

Xa

Xa

: ra dS u dXa

ðA:3Þ

A.1.4. Derivation of the force term

D Da

Z

f  dS u dXa ¼ Xa

Z Xa



Df þ ra  V aM f Da



 dS u dXa

Putting back together these terms leads to the complete Lagrangian sensitivity equations of the momentum A.2. Energy equation The energy Eq. (3) is

qcP u  rT ¼ r  ðkrT Þ þ qS

ðA:4Þ

L. Charlot et al. / Journal of Computational Physics 231 (2012) 5989–6011

6011

Let dST be an admissible test function, the weak form of this equation is:

Z

ðqcP u  ra T dST þ kra T  ra dST Þ dXa ¼

Xa

Z

qS dST dXa

Xa

We differentiate the weak form with respect to the shape parameter a with the method presented in Section 4 and we obtain the Lagrangian sensitivity equation of the energy:

   DqcP þ qcP ra  V aM u  ra T þ qcP S u  ra T þ qcP u  ra ST  qcP u  ra T  ra V aM dST dXa Da Xa   Z    Dk þ þ kra  V aM ra T  ra dST þ kra ST  ra dST  kra T  ra V aM þ raT V aM  ra dST dXa Da Xa  Z  DqS ¼ þ qS ra  V aM dST dXa Da Xa

Z



ðA:5Þ

References [1] D. Pelletier, A. Hay, S. Etienne, J. Borggaard, The sensitivity equation method in fluid mechanics, REMN-EJCM 17 (2008) 1–32. [2] R. Duvigneau, D. Pelletier, On accurate boundary conditions for a shape sensitivity equation method, International Journal for Numerical Methods in Fluids 50 (2006) 147–164. [3] J.-P. Zolésio, The material derivative (or speed) method for shape optimization, Optimization of Distributed Parameter Structures, vol. 2, Sijhoff and Noordhoff, The Netherlands, 1981, pp. 1089–1151. [4] J. Céa, Problems of shape optimal design, Optimization of Distribiuted Parameter Structures, vol. 2, Sijhoff and Noordhoff, The Netherlands, 1981, pp. 1005–1048. [5] M.C. Delfour, J.-P. Zolesio, Shapes and Geometries: Analysis, Differential Calculus, and Optimization, Advances in Design and Control, SIAM, Philadelphia, 2001. [6] P.I. Plotnikov, E. Ruban, J. Sokolowski, Inhomogeneous boundary value problems for compressible Navier–Stokes equations, well-posedness and sensitivity analysis, SIAM Journal of Mathematics Analysis 50 (3) (2008) 1152–1200. [7] P.I. Plotnikov, J. Sokolowski, Shape derivative of drag functional, SIAM Journal of Control Optimizations 48 (7) (2010) 4680–4706. [8] J. Bello, E. Fernandez-Cara, J. Simon, Optimal shape design for Navier–Strokes flow, in: L. Davisson, A. MacFarlane, H. Kwakernaak, J. Massey, Y. Tsypkin, A. Viterbi, P. Kall (Eds.), System Modelling and Optimization, Lecture Notes in Control and Information Sciences, vol. 18, Springer, Berlin/ Heidelberg, 1992, pp. 481–489. [9] D.A. Tortorelli, Z.X. Wang, A systematic approach to shape sensitivity analysis, International Journal of Solids Structures 30 (1993) 1181–1212. [10] R.T. Haftka, Techniques for thermal sensitivity analysis, International Journal for Numerical Methods in Engineering 17 (1) (1981) 71–80. [11] A. Sluzalec, M. Kleiber, Shape sensitivity analysis for nonlinear steady-state heat conduction problems, International Journal of Heat and Mass Transfer 39 (12) (1996/08/) 2609–2613. [12] K. Dems, R. Korycki, B. Rousselet, Application of first- and second-order sensitivities in domain optimization for steady conduction problem, Journal of Thermal Stresses 20 (7) (1997/10) 697–728. [13] K. Dems, B. Rousselet, Sensitivity analysis for transient heat conduction in a solid body – Part I: External boundary modification, Structural Optimization 17 (1) (1999) 36–45. [14] K. Dems, B. Rousselet, Sensitivity analysis for transient heat conduction in a solid body – Part II: Interface modification, Structural Optimization 17 (1) (1999) 46–54. [15] E.J. Haug, K.K. Choi, V. Komkov, Design Sensitivity Analysis of Structural Systems, Academic Press, New York, 1986. [16] J.S. Arora, An exposition of the material derivative approach for structural shape sensitivity analysis, Computer Methods in Applied Mechanics and Engineering 105 (1993) 41–62. [17] F. Bobaru, S. Mukherjee, Shape sensitivity analysis and shape optimization in planar elasticity using element-free Galerkin method, Computer Methods in Applied Mechanics and Engineering 190 (2001) 4319–4337. [18] B.Y. Lee, Direct differentiation formulation for boundary element shape sensitivity analysis of axisymmetric elastic solid, International Journal of Solids Structures 34 (1997) 99–112. [19] E. Taroco, Shape sensitivity analysis in linear elastic fracture mechanics, Computer Methods in Applied Mechanics and Engineering 188 (2000) 697– 712. [20] Z.X. Wang, D.A. Tortorelli, J.A. Dantzig, Sensitivity analysis and optimization of coupled thermal and flow problems with applications to contraction design, International Journal for Numerical Methods in Fluids 23 (1996) 991–1020. [21] P.A. Sackinger, P. Schunk, R.R. Rao, A Newton–Raphson pseudo-solid domain mapping technique for free and moving boundary problems: a finite element implementation, Journal of Computational Physics 125 (1996) 83–103. [22] E. Turgeon, Méthode d’éléments finis adaptative pour la résolution des sensibilités d’écoulements, Ph.D. Thesis, École Polytechnique de Montréal, 2001. [23] D. Pelletier, Adaptive finite element computations of complex flows, International Journal for Numerical Methods in Fluids 31 (1999) 189–202. [24] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique, International Journal for Numerical Methods in Engineering 33 (1992) 1331–1364. [25] P.J. Roache, Verification and Validation in Computational Science and Engineering, Hermosa publishers, Albuquerque, New Mexico, 1998. [26] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity, International Journal for Numerical Methods in Engineering 33 (1992) 1365–1382. [27] L. Eça, M. Hoekstra, Testing uncertainty estimation and validation procedures in the flow around a backward facing step, in: Eça L., Hoekstra M. (Eds.) Proceedings of the 3rd Workshop on CFD Uncertainty Analysis, Lisbon, October 2008.