- Email: [email protected]
Shape derivative of discretized problems
Computer Methods in Applied Mechanics and Engineering North-Holland
108 (1993) 187-199
CMA 431
Shape derivative
of discretized
problems
M. Souli and J.P. Zolesio Universite de Nice, Institut Non LinPaire de Nice, UMR CNRS no. 129, Vniversid Part Valrose Nice 06034 Cedex. France
de Nice-Sophia
Antiplis,
Received 7 March 1989 Revised manuscript received 6 January 1993
A nodal derivative formulation is given for shape optimization problems. A review of the kinematics and field equations is presented. Since the major challenge in the free surface calculation lies in the mesh moving algorithm, a new technique is developed concerning the motion of the internal nodes, in order to preserve a regular mesh, and prevent a material distortion that leads to element entanglement. The procedures developed here are appropriate for modeling free surface problems.
1. Introduction Presently, one of the main areas for shape optimization in applied mathematics is the continuous gradient with respect to the domain, or to the boundary of the domain. This discipline has achieved a high degree of sophistication and success from a numerical and theoretical point of view. The continuous formulation of the gradient with respect to the domain introduces the curvature of the boundary. One generally uses cubic spline interpolation for boundary approximation; this method allows the curvature of the boundary to be approximated. Another aspect of shape optimization is the nodal derivative. This formulation is used due to its convenience to avoid the computation of the curvature of the boundary. For many free surface problems and problems of fluid-structure interaction, it is convenient to describe the problem as a shape optimization problem. Nevertheless, in some simple cases, the free surface may be computed using a fixed point method, which is less expensive and easier to implement. Shape optimization involves the minimization of a discrete function J(oh, uh), with respect to the discrete domain Q, described by the location of the boundary nodes; u,, is a finite element approximation on the discrete domain. The case where uh is approximated by a Pl Lagrange finite element is treated by Zolesio [l]. The aim of the paper is to extend the method to higher order triangular and quadrangular Lagrange finite elements. For these types of element, we show that the shape functions are convected, and thus their material derivatives vanish. This is not the case for Hermite type elements where the functions are no longer convected, and a contribution due to their. derivative has to be considered. It is not the goal of this paper to overcome this difficulty; a study of Hermite finite element analysis is in progress and will be discussed in a forthcoming paper. The aim of this paper is to introduce a sensitivity finite element analysis, and to give a description of the computation of the nodal derivatives, using higher order finite element approximation. In Section 1, we develop the general theory for the computation of the material derivative of a function that depends on the solution of a boundary value problem. A mathematical study for the material derivative and shape derivative of the boundary value problem is developed in (1) [2-41. In Section 2, a general theory is developed which serves as the basis of the material derivative of the variational formulation describing the flow velocity, as well as the material derivative of the constraint describing the free surface condition. In Section 3, the computer implementation of the nodal derivative is considered; it is shown that the Correspondence
to: M. Souli, Centric Engineering
System, 3801 East Bayshore Road, Palo Alto, CA 94303, USA. (Present
address.) 0045-7825/93/$06.00
0
1993 Elsevier Science Publishers B.V. All rights reserved
M. Souli, .I. P. Zolesio, Shape derivative of discretized problems
188
nodal derivative is a particular case of the general theory developed in shape optimization, using an appropriate choice of the nodal velocity. The domain 0 can be described by the position of the boundary nodes; in Section 4, we describe a way to move the internal nodes in order to maintain a regular mesh, and compute the contribution to the derivative of J, due to the motion of the internal nodes. Numerical computations for the second order problem using quadrilateral finite elements are illustrated for the numerical computation of the free surface of a steady inviscid flow. In order to focus on the free surface calculation, we consider a bidimensional potential and incompressible flow. The flow field is described by a stream function, and thus the free surface is a streamline; this gives a boundary condition on the free surface.
2. General theory For our analysis, we introduce a fictitious time t; 0 is to be thought of as the domain occupied at t = 0 by the material particles, 0, as the domain upon which the fluid takes place. Thus, the derivative with respect to the domain is analogous to a temporal derivative. The computation is realized by introducing the deformation T,, constructed from a velocity field V. The properties of T, are widely studied in [3]. If U, is the solution in the configuration 0, with boundary c, we denote by J(C?,, uI) the function to minimize with respect to the domain. The derivative of J at the reference domain, or at time t = 0, in the v-direction is denoted by dJ(0, u) = dJ(L&, IJ) /dt (at t = 0). We established the properties of its derivative in [3]. In particular, when r is a smooth boundary, dJ(Q, u) depends on the domain, only through its boundary. The computation of the gradient dJ(R, v) is done in general, via the calculation of the material derivative of the state equation. Material derivatives are indicated by a dot above the symbol, e.g., ti = d(u, 0 T,) /dt .
In practice, we are concerned with minimizing the function J(Q, u,), which in general is positive and represents the kinematical energy on the free surface, with respect to the domain a,, where U, is the solution of a boundary value problem defined in 4: a(u,, w,) = l(w,)
(7-l)
VW, E Iv, .
w, is the trial function defined in an appropriate Sobolev space on the domain 0,. For instance, let us minimize the function J with respect to the domain, defined by
I tu,- 4’ dx, R,
(2.2)
3
where U, is the solution of (2.1) and z is an H’(R2) function. The derivative t = 0, is given by
dJ in the V-direction,
dJ(R, v) = jn (U - z) - (u -Vz . u(x, 0)) dx + ; jr (u - z’) - v(x, 0). n ds . However,
introducing
the partial, with respect to t, derivative
U’ = ri -vu
* V(0) )
and combining
at
(2.3)
U’ of the state solution u, defined by (2.4)
(2.3) and (2.4), we obtain
dJ(fh
u>=IO (u -
z) * u’ dx + ; 1 (U - 2)‘. u(x, 0). n ds .
(2.5)
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
189
The computation of the gradient (2.5) requires the computation of the shape derivative u’, which is the solution of a boundary value problem deduced from deriving the variational formulation (2.1) with respect to t. For deriving the expression (2.9, we assume the regularity of the boundary r. This results from the structure theorems for the shape derivatives and gradients (see [l, 31) that for a smooth boundary, the value of the gradient df in (2.5) depends only on boundary integrals; thus
where g is a function of the solution U, the state adjoin p, and the curvature of the boundary J’; n is the outward normal to r. As long as we are concerned with the discretized approximations of U, we cannot use the shape derivative described by (2.6); this would make use of the regularity of the solution (which is true for the continuous soiution u but not necessary for its finite element approximation uh).
3. The deformation field We distinguish two types of deformation. The first is adapted to small displacements, and the second to large displacements, and to deformations governed by parameters. Let &?be an open domain in R” (we will take it = 2), and let V be the velocity field defined on the neighborhood R. The transformation T, is defined on a by T,:X+X+t*V(X),
T,(X) = X +
I
of
V(s, T,(X)) ds .
(3.1)
We are more concerned with the second deformations which, for small values of t, have the same properties as the first one, and preserve them for a large value of t. In addition, for shape optimization defined by parameters, it is easy to find the field u(x, t) which constructs the variation of the geometry as well as the variation of the parameters. So the derivative of the cost and states with respect to the parameters is simply a matter of a particular choice of the field V. For such situations, see {S]. The analogous situation is presented for nodal mesh displacements, however the coordinates of the nodes are the parameters of the problem. As for the continuous case [3,4], our purpose is to define a transformation T,, which maps the discrete domain fl into 0,. The transformation T, is solely defined through a velocity field V( *, t) and is the solution of the differential equation.
dT,(X) /dt = V(t, I;(x)) ,
TotX)= X,
(3.2)
where V(t, T,(X)) is the Eulerian velocity expressed in the configuration at time t and location x = T,(X). We denote the material points on the domain fl and a,, by X = (Xi, X,) and x = (x,, x,), respectively, so the solution of the differential equation (3.2) is given by (3.1). The transformation T is defined by its restriction to each element K. To avoid further notation, the algorithms are restricted to the higher order Lagrange quadrangular elements; the case of triangular elements is deduced in a straightforward way. From now on, we deal with discrete formulations, defined on finite dimensional spaces. We denote the discrete domain, the finite element solution and the discrete space to which the solution belongs by 0, u and W, respectively. The first aspect is to consider a regular mesh of the domain (see [6]) made up of quadrangles, established over the set 0, i.e. the set 0 is written as a finite union of finite elements K, each element is affine isoparametrically equivalent to a reference element g-, where the reference coordinates are denoted by (5, n), in the sense that there exists a mapping FE (Q,(i?))‘, polynomial functions of degree s, the degree of approximation, which maps the element R into the element K. But for the isoparametric Lagrange element, the mapping F degenerates and becomes of degree 1 on each variable, so F is written as
M. Souli, J. P. Zolesio, Shape derivative of discretired problems
190
(3.3) where fii is the shape function in the reference element, polynomial function of degree 1. The second aspect in developing the method is to construct a finite dimensional approximation of the space W (W = W, at t = 0) defined in (2.1). A function u in the discrete space W takes the form nNodes 2.4 =
C
ui ’ ei(xl,
$1
(3.4)
Y
i=l
where nNodes is the number of nodal points and e, is referred to a shape, basis or interpolation function associated with node i; its restriction to each quadrangular element K is also denoted by ei, and is defined by e, =
d,oF-’ ,
4E
c?,uG
(3.5)
7
where s is the degree of approximation. For our analysis, we use the polynomial function of degree 1 iVi for the coordinate interpolation, and the polynomial function of degree s ei, for the function interpolation. Let us notice that for s > 1, the polynomial function ei, is different from the polynomial function of degree 1 iVi. For the analysis, we use the function iVi to map the reference element R into the element K, and to define the nodal velocity. The transformation T, is completely determined on the discrete domain 0 from its restriction on each element K. For simplicity, let us describe the velocity field u(x, t), and hence the transformation T, from (3.1) on the element K, when only one vertex node M is moving in the x1 direction to the vertex node M,. Let K, denote the new element *defined from this motion. The element K, is isoparametrically equivalent to the reference element K through the mapping F,. At each time t, the transformation T,(X) and hence the Eulerian velocity u(x, t) must satisfy the following two properties: (a) The transformation T, maps the triangle or quadrangle K into the triangle or quadrangle K,. (b) The transformation T, transforms the finite element basis ei (i = 1, nNodes) defined by (3.5) into the function ef on the element K,. This property is very important from a numerical point of view in shape optimization problems. We show that the basis functions ei are convected; by this we mean that the functions e: are the basis function defined on the domain 0,. For the material derivatives of the discrete formulation of the variational problem (2.1), the nodal derivative of the basis functions ef vanish. This assumption considerably simplifies the computation. Through the transformation M to M,, the reference configuration 0 is transformed to the current configuration 4. The transformation T,, which maps 0 into a,, is invariant for all vertices Mj, j # i. A natural choice of the nodal velocity u(t, x), for the transformation T,, is the shape function N: , corresponding to the moving node M: defined on the current configuration, u(t, x) =
(Nf 0) .
(3-e)
)
The function N: is defined on each element K by N;=+F,‘, The transformation given by T,(X)
%&Q,(R). solution
T, of the differential
(3.7) equation
(3.2), defined on the element
= X + t . (N,(X), 0) .
Let us show that the transformation
T, defined by (3.8) satisfies the two previous properties
K, is then
(3.8) (a) and (b).
191
M. Souli, J.P. Zotesio, Shapederivative of discretizedp~o&Ie~
First 1; transforms Mi to M: and the element K to the element K,, second through the transformation the basis function ei on the element K is transformed to the basis function e: on K,.
T,
LEMMA 3.1. Let K be an element on the discrete domain 0, and K, be the element constructed on Q from K when a vertex Mi is moving to Mi, in the k-direction (k = 1,2); thus the transformation 1”, ~5 died by T,= FpF-‘.
(3.9)
PROOF. Without loss of generality, we consider the case k = 1. For t = 0, (3.9) is satisfied; let us show that dT,/dt
= d(F,o F-')fdt
(3.10)
.
Let us write the mapping F, in terms of the shape functions I’$ on the reference
element,
where nElemNodes is the number of vertex nodes per element, and thus nElemNodes = 4. Since K, is isoparamet~cally equivalent to the reference element 2. The transfo~ation T, satisfies Ts(Mi) = Mj for jZi; thus
x;= xi
i#j,
,
x; = xi + t )
y; = y, *
The mappings F, can be written as
@(f>q) + $
Ft( 5, n) = (xi + t, Yi>*
(xi, Yi). q( 6, n)
and F,(h?)=(t,O).~(~,~)~~~~(xj,Yi).~(i,~)=(t,O).~(~,~)+‘(~,~). Thus F, = F + t - (c,
0). Deriving the previous expression with respect to t, we obtain
dF,/‘dt = (g, 0). Let us derive the expression
(3.9) with respect to t:
d(F~~F-l)/dt=(dF*~dt)~F-l=(~,O)~F-l=(~i,O),
(3.11)
using (3.9), we finally obtain dT,/‘dt = (Ni, 0) .
The first property
(a) is then satisfied, and thus T, transforms
To satisfy the second property,
the element K to the element
K,.
q
we need the following lemma.
LEMMA 3.2. For each nodal point j, 1 -C - 1’ d nNodes, the basis function ei on the element K, satisfies eioT,kej.
(3.12)
192
PROOF.
M. Souli, J. P. Zolesio. Shape derivative of discretized problems
Using expression
(3.9)) we have (3.13)
e,foT,=e,foF,oF-‘.
Combining expressions (3.12). 0 COROLLARY
3.1.
(3.12)
using the relations
ej 0 F, = ii and Ei 0 F-
’=
ej, yields
the relation
For each node Mj, the basis function ef satisfies
[d(e:0 T,) /dtlc,,oj= 0 , which means the basis functions are convected.
The proof of the corollary independent.
is straightforward
from Lemma 3.2 since the function
ei 0 T, is time
4. Nodal derivatives In the general case, we are looking for partial derivatives with respect to the coordinates of the vertex nodes defining the discrete domain, for a functional J(f&, u,), where u, is a discretized solution of a variational problem defined on the discrete domain 4, which is called the state formulation. We show in this section how to derive bilinear forms defined in (2.1). To simplify the analysis, we only consider the two-dimensional case. 4.1. Nodal derivatives of bilinear forms Let u’ be the solution of the discrete variational formulation u, defined in 0 through the transformation T,; thus
in the current configuration
0,; let us set
Let us consider a bilinear form arising from a second order problem, on the domain a,, which is a finite union of elements K,; u’ is the solution in the discrete space W,;
(4.1) To compute the material or nodal derivatives of the bilinear form (4.1), let us define the integrals in which involves the Jacobian J, of the (4.1) on the element K, using variable transformation, transformation T,. Vu’ . Vei dx, =
(A(t)Vu,,
Ve,)) dx ,
where e, = e: 0 T, is the basis function on the element K for node i and A(t) is a 2 each element K by A(t) = J,(DT;‘)‘(DT,) .I, = det(DT,)
.
,
(4.2) x
2 matrix defined on
A(0) = Id, (4.3)
To derive the bilinear form (4.2), we first need to derive A(t) with respect to t, at t = 0, we need to derive the expression (4.3), which gives
193
M. Souli, J.P. Zolesio, Shape derivative of discretired problems
A’(0) = [dA(t) /dt] /t = 0 = div(V(0)) - (DV(0)
+ ‘DV(0))
.
Thus, using the fact that the material derivative of the basis function convected, the derivative of the bilinear form (4.2), at t = 0, is i
V(eio T,)) dx/(t = 0) = IK (A’(0)
jK (A(t)Vu,,
ei vanishes,
since they are
-Vu,, Ve,) + (Vzi, Ve,)) dx ,
$ JK F. (eio T,) dx/(t = 0) = IK div(F. u(x, O))e, dx .
(4.4)
For the two-dimensional problem, the matrix A’(0) is straightforwardly computed from (4.3). For x1 and xq coordinate motion of the nodal point Mi, respectively, the convective velocity u(x, 0) is represented by the basis functions (Ni, 0), and (0, Ni), and hence A’(0) is given by
(4.5) respectively. Combining expressions (4.4) and (4.5), the material derivative solution of the following variational problem:
+
u for u = (0, iVi) is the
d,(F- e,). Ni dx .
(4.6)
4.2. Nodal derivative of the functional J Let us consider the functional J to be minimized, I(&
defined by
u’) = n, (u’ - z)’ dx, ,
(4.7)
I
where U’ is a solution of the variational problem (4.2), and z is a real function defined on the current configuration 0,. Transforming the integral in (4.7) on the reference configuration 0, we obtain J(Q, u’) = ; I0 (u’ - z 0 T,)‘) . J, dx .
(4-g)
with respect to t at t = 0 in the u direction is given by
The derivative dJ(0,
u) = R (ti -Vz I
. u(x,
O))(u - z) + $((u - z)* div(u, 0)) dx .
A classical way to avoid the computation of the material derivative u in the solution of the variational problem (4.1) is to introduce an adjoin state,
I
R
Vp.Vwdx=
I
n(u-z).wdx
VWEW.
(4.9)
Using the expressions (4.6) and (4.9), where the state function p plays the role of the shape function w, the gradient dJ takes the form dJ(R,u)=;~KVpVtidx+/K Using expression
- Vz . u(x, O)(u - z) + 1 ((u - z)’ div(u(x, 0)) dx .
(4.6), we obtain the following form of the gradient:
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
194
dJ(f& u> =2
j- alei(-a,u
- a2Ni + alNi. a2u) + a2e,(dlNi.
a+ + a,N,u) dx
KK
+ c 1 KK
-Vz *v(x, O)(u - z) + $(u - z)’ div(u(x, 0)) dx .
(4.10)
The formulation (4.10) is valid for a general velocity field u(x, 0). To compute the partial derivative with respect to the x1 and x2 coordinates of the node Mj, we take u(x, 0) = (0, N,(x)) or u(x, 0) = (Ni , 0). This choice of velocity simplifies the expression (4.10). THEOREM. Let u and p be the solutions of the dkcrete variational problem (4.2) and (4.9), and let J be the functional defined in (4.8). The partial derivatives of J with respect to the x, and x2 coordinates of the node Mi, respectively, are a,u. a,p - a+ a2p) + a2Ni. (aiue alp + a, d,p) dx
(a#‘- Ni)* p - (Ni a,z(u - z) dx + 1 . (u - z)* a,Ni dx
(4.11)
and a,u. alp - a+ a2p) + a,Ni . (alue a2p + a2 alp) dx +c KK
I
(a,(F.N,)*p-Nia2z.(u-z)dx+
$.(u-z)*dZNidx.
(4.12)
These partial derivatives have already been computed by Pironneau [7]. However, the technique in this paper allows us to compute the derivatives (4.11) and (4.12) using the theory of the continuous gradient, and to select a velocity field for the motion of the nodes in such a way that the basis functions are convected, and thus the computation of (4.11) and (4.12) is simplified.
5. Motion of internal nodes The motion of the nodes on the free surface will also serve the purpose of defining the motion of the internal nodes. Each node on the free surface (for example node A), is considered to be the endpoint of a straight line segment, composed of element nodes, which begins at some base point B in space. To describe clearly our method of moving the internal nodes, let us consider a free surface flow for a two-dimensional problem. We assume henceforth that the reference domain is covered by a mesh of isoparametric finite elements, and we think of 0 as the location of the internal nodes. Let us assume that the line y = g(x) coincides with the free surface (Fig. 1). Furthermore, suppose that there is to be no horizontal motion of the nodes and that the vertical motion of the internal nodes is controlled by the
Fig. 1.
M. Souli, 3. P. Zolesio, Shape derivative of discretized problems
195
vertical motion of the nodes on the free surface. Then the following representation of 6 giving a linear distribution of the internal nodes on the segment AB is appropriate, let us define
where nl is the number of boundary vertex nodes, and n2 the number of internal vertex nodes. The application 8 defines the location of the internal nodes in terms of those on the boundary. To prevent the proliferation of notation, let us denote the coordinates of the node Mi by (xi, y,), @(y,, Y,, . * * 9 Y,*) = (Y,1+1,
Yn1+29 * . * 9 Ynl+n2).
(5.1)
Our purpose is to impose a motion on the internal nodes, and so define the application 8, in order to keep a regular mesh for the domain. Using (5. l), the functional J, to be minimized, can be expressed in terms of the coordinates of the vertex nodes n = nl + n2, Ju-4
u) = J(YI7
’ * . 7 YE17 Y,1+17
f * . ) Ynl+n2).
Using (5.1), where the coordinates of the internal nodes are expressed in terms of those on the boundary, the functional J is a function of the coordinates (y,, y,, . . . , y,,), J(y,,
. . . ) Y”l,
@(y,,
. . * 9 Y”1,
Ynl))
= J,(Yl,
Y27.
. *’
Yfzl) .
Using the new formulation Je, to minimize J, we only consider the partial derivative of Je with respect to the y-coordinates of the boundary vertex nodes of the mesh. These derivatives introduce new terms from the fact that the internal nodes depend on the motion of the boundary nodes through 0. Let us assume that 8 is a differentiable function; thus the Jacobian (De), gives the motion of the internal vertex nodes, 1 s i < nl and 1 s j s n2. A straightforward calculation of the partial derivative of J gives
The derivatives aJ/ay, and dJ,/ay, are computed from (4.11) and from the fact that 0 is known explicitly. Let us give an example for the calculation of 8 Je lay,,, . Consider the domain in Fig. 1, n = {(X, y) st a S x ab
and O
The domain is subdivided into quadrangular elements, each line y = i - h, h = (b - a) /n and i = 1 * * , n is divided into m - 1 equal segments. Let us denote by (xi, y,), 1 d i G N, the coordinates of the N vertex nodes, and number of boundary vertex nodes from 1 to nl. However, for an internal node Mj, j > n, (j can be written as j = p. n, + d where d < n,), we have yj = p*yd/(m - 1). This defines the transformation 0 by O( yj) = p*y,/(m - 1). From this expression, we can easily compute the partial derivatives of J using the chain rule.
6. Application
to problems
in hydrodynamics
6.1. The continuous case In this section, the nodal derivative formulation is applied to two different engineering problems. The first one is a numerical model of a wave problem generated by a potential flow around an obstacle resting on the bottom. A critical parameter for this problem is the Froude number F, which is a
196
M. Souli. J.P. Zolesio, Shape derivative of discretized problems
nondimensional parameter defined by F = (C/g * /z)~‘~), w h ere C is the fluid velocity upstream, g is the acceleration of gravity and h is the depth of the fluid. The amplitude wave analysis for the two-dimensional problem is studied in 181. For values of the Froude number less than one, a propagation of a solitary wave over a constant depth is generated by the obstacle. For Froude number greater than one, the wave amplitude decreases away from the obstacle, and the flow becomes uniform downstream. The second problem concerns a flow around an obstacle on the free surface. The difficulty with this problem comes from the singularities on the free surface at the two contact points on the obstacle. This problem is very important from an engineering point of view, but few analytical and numerical results are available. and deduce the equations and boundary conditions Let us explain the physical phenomena, governing the flow. In what follows, we consider the potential and incompressible flow. We define a stream function !P, such that the velocity u = (aP/ay
)
-N/ax).
61.1.
Problem with an obstacle at the bottom For simplicity, we consider a uniform flow upstream and downstream to infinity (see [8,9]) and write the equations concerning a flow around an obstacle resting at the bottom, in a nondimensional stream function formulation:
A?P=O,
0,
F=O,
B,
a?Pl&z=O,
S,,S,,
EW/an=a(y),
S,
(6.1)
where a(y) = (1 - (2 /F*F)*( y - 1))“’ . Let us explain the physical meaning of (6.1). (a) Equation in J2: since the flow is potential, the vorticity vanishes, so the stream function is harmonic. (b) Condition at the bottom B: the bottom B is a streamline; this implies that ly is constant on B and we can fix an arbitrary value of zero. (4 Conditions on lateral boundaries: the flow is uniform upstream and downstream to infinity. If the domain is large enough compared to the obstacle, then the homogeneous Neuman boundary conditions seem reasonable. This boundary condition is valid only for the Froude number greater than 1. In this case, the problem (6.1) is variational and symmetrical. For F < 1, (6.1) is no longer variational and non-symmetrical. Thus, we obtain more complicated boundary conditions on the lateral sides. This is described in detail in (12). (d) Conditions on the free surface: along the free surface S, we dispose of two boundary conditions. The first one is deduced from the Bernoulli equation, and is related to the continuity of the pressure through the free surface; this is given in (6.1), where F is the Froude number and y(x) the nondimensional equation of S. The second condition explains the fact that S is a streamline, so P is constant on S. The constant is fixed by the conservation of the mass flux across a lateral section. In nondimensional parameters, this gives W--l,
s.
(6.2)
The two boundary conditions (6.1), (6.2) on S express the fact that S is an unknown. We formulate the problems (6.1), (6.2) as a minimization problem. The functional to be minimized is J(S) = (l/2)*
k (W - l)* ds ,
where !P is the solution of (6.1).
(6.3)
M. Souli, J. P. Zolesio, Shape derivative
of discretized
197
problem
The shape of the domain fl is defined by the boundary nodes. We only consider the motion of the nodes on S. To simplify the calculations, only the Y-direction of these nodes is computed. Using the technique developed above, the partial derivative of J with respect to the boundary nodes in the Y-direction V = (0, e,) is
dJ(G vi) =
-I\ (A.V'P,
VP) dx + h a( y)*(?P - y) -t l)*B
J(df(du(y)/dy)*(~-
+ s
+ (l/2)*
.y
ds
y) - a(y)*P+ vl: ds
IS (!P - 1)2*B * v ds ,
The matrix A is defined in (3.5) and B = div V- (DV. n, n). The adjoin state P is solution of hP=O,
fi,
P=O,
B,
aP/an=(*-l),
S,
aPjan=o,
s,, s,.
(6.5)
6.1.2. Problem with an obstacle piercing the surface From a mathematical viewpoint, this problem has the same formulation, except that in this case, a part of the free surface C is known. C is the region where the obstacle pierces the free surface S. Thus, the problem is formulated as follows: Av’=O,
0,
?If=O,
B,
aPlafl=O,
S1,S2,
a~/an=~(y), SK, V=O,
C,
with the additional condition P=l,
s‘\c .
6.2. Numerical results The domain 0 in Fig. 1 is subdivided into quadrangles. A Q1 Lagrange finite element approximation is used to solve the variational problems deduced from (6.1) and (6.5) for the state u and the adjoint state I? To check the validity of the gradient, we compare it to the gradient obtained using a finite difference method. Thus, for each node i = 1, . . . , N on the surface S, the finite difference gradient in the Y-direction is
The perturbation of the initial free surface for the computation of the derivative dJi is generally used for an approximation of the Jacobian matrix in nonlinear Newton iterations. The values of dJi are stable with respect to E for E c 10w4. For the first numerical simulation, we consider the initial free surface y = 1. The flow is uniform upstream, and the obstacle is sinusoidal (see Fig. 2). The Froude number F = 2.55, the length of the
198
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
0.690 I 0.460
s
n
%
It
4.c Fmude I 2.55
homqenous
neumann
Fig. 2.
domain L = 1, and the mesh size h = 0.1. The values of the gradient dJ for the nodes on the free surface are compared in Table 1. We give numerical results for F < 1, for which the mathematical formulation is more complicated because of the boundary conditions on the lateral sides S, and S,; these are no longer the homogeneous Neuman boundary conditions (6.1) but satisfy the continuity of the mass flux through S, and S,. Outside the domain 0, the linearized problem for ly has an analytical solution. The gradient which obviously has a more complicated form is given in [lo]. To search for the free surface, or water wave, we use the minimization algorithm in (1)) where an approximation of the Hessian matrix of the gradient is used in order to obtain the optimal slope. To distinguish the water wave from the initial wave, y = 1, we take a long domain L = 4 (Fig. 3). The convergence is fast, because the problem (6.1) is variational; we reach the steady state after 4 iterations, and the cost decreases from J, = 1.23 to J4 = 0.65 x 10V4. This case is also treated in [19] where a fixed point method is used. Figure 4 corresponds to F = 0.5, as described in [9]. One observes a solitary wave propagating to the right of the obstacle; this wave is due to the special boundary condition on the lateral sides. Figure 5 corresponds to a flow around an obstacle piercing the surface. The results are analogous to those given in [9] with the fixed point method. The hxed point method in [9] is more efficient and faster than the gradient method, when it converges. But it fails for certain cases, whereas the gradient method converges. This case corresponds to Fig. 4 with F = 0.4 and the amplitude of the sinusoidal obstacle is 0.4. The gradient method converges after 40 iterations, and the cost decreases from J, = 0.270 to J40 = 0.12 X 10P4. For the search of the water wave in naval hydrodynamics problems, or in general free surface problems, we should first test the fixed point method, because it is faster and less expensive. The convergence is geometrical; for
Table 1 Gradient in (6.6)
Gradient in (6.4)
0.7453E-00 -O.l129E-00 -0.2367E-00 -0.2835E-00 -0.2577E-00 -O.l892E-00 -O.l105E-01 -0.4493E-01 -0.3621E-02 O.l433E-01 O.l116E-01
0.7455E-00 -O.l13OE-00 -0.2368E-00 -0.2839E-00 -0.2573E-00 -O.l892E-00 -O.l108E-01 -0.4498E-01 -0.3634E-01 O.l438E-01 O.l12OE-01
0
1.0
2.0
Fmtide
- 0.50
Fig. 3.
3.0
obstacle - 0.2
4.0
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
0
0.2
0.4
Frovde=OSO
Fig. 4.
0.6
0.8
obstacle.O.t,
1.0
0
0.5
199
1.0
1.5
2.0
Frwde -0.64 obstacle -0.05
Fig. 5
this method; at each iteration we only need to solve a linear system corresponding to the state problem. In case this method fails, the gradient method may converge, at least for some cases. The gradient method is more expensive; at each iteration, we solve two linear systems, the state and the adjoin state, and the convergence as it is known, is not geometrical.
Acknowledgment
This research was supported by DRET contract with the collaboration of ENSTA at Palaiseau. The authors would like to thank M. Lenoir and J. Cahouet at ENSTA for helpful discussions. The first author wishes to thank Charlotte Van Campen at Centric engineering for the help in making the figures.
References [l] J.P. ZolCsio, Les derives par rapport aux noeuds des triangularisations et leurs utilisation en identification de domaine, Ann. Sci. Math. Quebec 8 (1984) 97-120. [2] J. Cea, Problems of shape optimal design, in: J. Cea and Haug, eds, Optimization of Distributed Parameters Structures, Vol 12 (Sijthoff, Alphen aan den Rijn, 1981) 1005-1088. [3] J.P. Zolbsio, Identification de domaines par deformations, Thesis, Nice, 1979. [4] M. Souli and J.P. Zoltsio, Gradient with respect to a non smooth domain in shape optimization, Optim. Control Applic. Control, in press. [5] M. Delfour, M. Payre and J.P. Zoltsio, Optimal and suboptimal design of thermal diffusers for communications satellites, Comptes rendus du 3ieme Symposium on Control of Distributed Parameters Systems (1982). [6] P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, Amsterdam, 1978). [7] 0. Pironneau, Optimal Shape Design for Elliptic Systems (Springer, Berlin, 1986). [8] J. Cahouet and M. Lenoir, Resolution numerique du problbme non lineaire de la resistance de vague non lineaire, CRAS 297, 1985. [9] J. Cahouet, Etude numCrique et expi%mentale du problbme de la resistance de vague non lineaire, Thesis, Paris, 1984. [lo] M. Souli and J.P. Zolbsio, Le controle de domaine dans la resistance de vague non lintaire en hydrodynamique, Ann Sci. Quebec 1 (1990) 121-132. [ll] F. Neittamaki, Finite Element Approximation for Optimal Shape Design. Theory and Application (Wiley, New York, 1987). [12] A. Buckeley, An alternative implementation of Goldfard’s minimization algorithms, Math. Prog. 8 (1975) 207-231.
108 (1993) 187-199
CMA 431
Shape derivative
of discretized
problems
M. Souli and J.P. Zolesio Universite de Nice, Institut Non LinPaire de Nice, UMR CNRS no. 129, Vniversid Part Valrose Nice 06034 Cedex. France
de Nice-Sophia
Antiplis,
Received 7 March 1989 Revised manuscript received 6 January 1993
A nodal derivative formulation is given for shape optimization problems. A review of the kinematics and field equations is presented. Since the major challenge in the free surface calculation lies in the mesh moving algorithm, a new technique is developed concerning the motion of the internal nodes, in order to preserve a regular mesh, and prevent a material distortion that leads to element entanglement. The procedures developed here are appropriate for modeling free surface problems.
1. Introduction Presently, one of the main areas for shape optimization in applied mathematics is the continuous gradient with respect to the domain, or to the boundary of the domain. This discipline has achieved a high degree of sophistication and success from a numerical and theoretical point of view. The continuous formulation of the gradient with respect to the domain introduces the curvature of the boundary. One generally uses cubic spline interpolation for boundary approximation; this method allows the curvature of the boundary to be approximated. Another aspect of shape optimization is the nodal derivative. This formulation is used due to its convenience to avoid the computation of the curvature of the boundary. For many free surface problems and problems of fluid-structure interaction, it is convenient to describe the problem as a shape optimization problem. Nevertheless, in some simple cases, the free surface may be computed using a fixed point method, which is less expensive and easier to implement. Shape optimization involves the minimization of a discrete function J(oh, uh), with respect to the discrete domain Q, described by the location of the boundary nodes; u,, is a finite element approximation on the discrete domain. The case where uh is approximated by a Pl Lagrange finite element is treated by Zolesio [l]. The aim of the paper is to extend the method to higher order triangular and quadrangular Lagrange finite elements. For these types of element, we show that the shape functions are convected, and thus their material derivatives vanish. This is not the case for Hermite type elements where the functions are no longer convected, and a contribution due to their. derivative has to be considered. It is not the goal of this paper to overcome this difficulty; a study of Hermite finite element analysis is in progress and will be discussed in a forthcoming paper. The aim of this paper is to introduce a sensitivity finite element analysis, and to give a description of the computation of the nodal derivatives, using higher order finite element approximation. In Section 1, we develop the general theory for the computation of the material derivative of a function that depends on the solution of a boundary value problem. A mathematical study for the material derivative and shape derivative of the boundary value problem is developed in (1) [2-41. In Section 2, a general theory is developed which serves as the basis of the material derivative of the variational formulation describing the flow velocity, as well as the material derivative of the constraint describing the free surface condition. In Section 3, the computer implementation of the nodal derivative is considered; it is shown that the Correspondence
to: M. Souli, Centric Engineering
System, 3801 East Bayshore Road, Palo Alto, CA 94303, USA. (Present
address.) 0045-7825/93/$06.00
0
1993 Elsevier Science Publishers B.V. All rights reserved
M. Souli, .I. P. Zolesio, Shape derivative of discretized problems
188
nodal derivative is a particular case of the general theory developed in shape optimization, using an appropriate choice of the nodal velocity. The domain 0 can be described by the position of the boundary nodes; in Section 4, we describe a way to move the internal nodes in order to maintain a regular mesh, and compute the contribution to the derivative of J, due to the motion of the internal nodes. Numerical computations for the second order problem using quadrilateral finite elements are illustrated for the numerical computation of the free surface of a steady inviscid flow. In order to focus on the free surface calculation, we consider a bidimensional potential and incompressible flow. The flow field is described by a stream function, and thus the free surface is a streamline; this gives a boundary condition on the free surface.
2. General theory For our analysis, we introduce a fictitious time t; 0 is to be thought of as the domain occupied at t = 0 by the material particles, 0, as the domain upon which the fluid takes place. Thus, the derivative with respect to the domain is analogous to a temporal derivative. The computation is realized by introducing the deformation T,, constructed from a velocity field V. The properties of T, are widely studied in [3]. If U, is the solution in the configuration 0, with boundary c, we denote by J(C?,, uI) the function to minimize with respect to the domain. The derivative of J at the reference domain, or at time t = 0, in the v-direction is denoted by dJ(0, u) = dJ(L&, IJ) /dt (at t = 0). We established the properties of its derivative in [3]. In particular, when r is a smooth boundary, dJ(Q, u) depends on the domain, only through its boundary. The computation of the gradient dJ(R, v) is done in general, via the calculation of the material derivative of the state equation. Material derivatives are indicated by a dot above the symbol, e.g., ti = d(u, 0 T,) /dt .
In practice, we are concerned with minimizing the function J(Q, u,), which in general is positive and represents the kinematical energy on the free surface, with respect to the domain a,, where U, is the solution of a boundary value problem defined in 4: a(u,, w,) = l(w,)
(7-l)
VW, E Iv, .
w, is the trial function defined in an appropriate Sobolev space on the domain 0,. For instance, let us minimize the function J with respect to the domain, defined by
I tu,- 4’ dx, R,
(2.2)
3
where U, is the solution of (2.1) and z is an H’(R2) function. The derivative t = 0, is given by
dJ in the V-direction,
dJ(R, v) = jn (U - z) - (u -Vz . u(x, 0)) dx + ; jr (u - z’) - v(x, 0). n ds . However,
introducing
the partial, with respect to t, derivative
U’ = ri -vu
* V(0) )
and combining
at
(2.3)
U’ of the state solution u, defined by (2.4)
(2.3) and (2.4), we obtain
dJ(fh
u>=IO (u -
z) * u’ dx + ; 1 (U - 2)‘. u(x, 0). n ds .
(2.5)
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
189
The computation of the gradient (2.5) requires the computation of the shape derivative u’, which is the solution of a boundary value problem deduced from deriving the variational formulation (2.1) with respect to t. For deriving the expression (2.9, we assume the regularity of the boundary r. This results from the structure theorems for the shape derivatives and gradients (see [l, 31) that for a smooth boundary, the value of the gradient df in (2.5) depends only on boundary integrals; thus
where g is a function of the solution U, the state adjoin p, and the curvature of the boundary J’; n is the outward normal to r. As long as we are concerned with the discretized approximations of U, we cannot use the shape derivative described by (2.6); this would make use of the regularity of the solution (which is true for the continuous soiution u but not necessary for its finite element approximation uh).
3. The deformation field We distinguish two types of deformation. The first is adapted to small displacements, and the second to large displacements, and to deformations governed by parameters. Let &?be an open domain in R” (we will take it = 2), and let V be the velocity field defined on the neighborhood R. The transformation T, is defined on a by T,:X+X+t*V(X),
T,(X) = X +
I
of
V(s, T,(X)) ds .
(3.1)
We are more concerned with the second deformations which, for small values of t, have the same properties as the first one, and preserve them for a large value of t. In addition, for shape optimization defined by parameters, it is easy to find the field u(x, t) which constructs the variation of the geometry as well as the variation of the parameters. So the derivative of the cost and states with respect to the parameters is simply a matter of a particular choice of the field V. For such situations, see {S]. The analogous situation is presented for nodal mesh displacements, however the coordinates of the nodes are the parameters of the problem. As for the continuous case [3,4], our purpose is to define a transformation T,, which maps the discrete domain fl into 0,. The transformation T, is solely defined through a velocity field V( *, t) and is the solution of the differential equation.
dT,(X) /dt = V(t, I;(x)) ,
TotX)= X,
(3.2)
where V(t, T,(X)) is the Eulerian velocity expressed in the configuration at time t and location x = T,(X). We denote the material points on the domain fl and a,, by X = (Xi, X,) and x = (x,, x,), respectively, so the solution of the differential equation (3.2) is given by (3.1). The transformation T is defined by its restriction to each element K. To avoid further notation, the algorithms are restricted to the higher order Lagrange quadrangular elements; the case of triangular elements is deduced in a straightforward way. From now on, we deal with discrete formulations, defined on finite dimensional spaces. We denote the discrete domain, the finite element solution and the discrete space to which the solution belongs by 0, u and W, respectively. The first aspect is to consider a regular mesh of the domain (see [6]) made up of quadrangles, established over the set 0, i.e. the set 0 is written as a finite union of finite elements K, each element is affine isoparametrically equivalent to a reference element g-, where the reference coordinates are denoted by (5, n), in the sense that there exists a mapping FE (Q,(i?))‘, polynomial functions of degree s, the degree of approximation, which maps the element R into the element K. But for the isoparametric Lagrange element, the mapping F degenerates and becomes of degree 1 on each variable, so F is written as
M. Souli, J. P. Zolesio, Shape derivative of discretired problems
190
(3.3) where fii is the shape function in the reference element, polynomial function of degree 1. The second aspect in developing the method is to construct a finite dimensional approximation of the space W (W = W, at t = 0) defined in (2.1). A function u in the discrete space W takes the form nNodes 2.4 =
C
ui ’ ei(xl,
$1
(3.4)
Y
i=l
where nNodes is the number of nodal points and e, is referred to a shape, basis or interpolation function associated with node i; its restriction to each quadrangular element K is also denoted by ei, and is defined by e, =
d,oF-’ ,
4E
c?,uG
(3.5)
7
where s is the degree of approximation. For our analysis, we use the polynomial function of degree 1 iVi for the coordinate interpolation, and the polynomial function of degree s ei, for the function interpolation. Let us notice that for s > 1, the polynomial function ei, is different from the polynomial function of degree 1 iVi. For the analysis, we use the function iVi to map the reference element R into the element K, and to define the nodal velocity. The transformation T, is completely determined on the discrete domain 0 from its restriction on each element K. For simplicity, let us describe the velocity field u(x, t), and hence the transformation T, from (3.1) on the element K, when only one vertex node M is moving in the x1 direction to the vertex node M,. Let K, denote the new element *defined from this motion. The element K, is isoparametrically equivalent to the reference element K through the mapping F,. At each time t, the transformation T,(X) and hence the Eulerian velocity u(x, t) must satisfy the following two properties: (a) The transformation T, maps the triangle or quadrangle K into the triangle or quadrangle K,. (b) The transformation T, transforms the finite element basis ei (i = 1, nNodes) defined by (3.5) into the function ef on the element K,. This property is very important from a numerical point of view in shape optimization problems. We show that the basis functions ei are convected; by this we mean that the functions e: are the basis function defined on the domain 0,. For the material derivatives of the discrete formulation of the variational problem (2.1), the nodal derivative of the basis functions ef vanish. This assumption considerably simplifies the computation. Through the transformation M to M,, the reference configuration 0 is transformed to the current configuration 4. The transformation T,, which maps 0 into a,, is invariant for all vertices Mj, j # i. A natural choice of the nodal velocity u(t, x), for the transformation T,, is the shape function N: , corresponding to the moving node M: defined on the current configuration, u(t, x) =
(Nf 0) .
(3-e)
)
The function N: is defined on each element K by N;=+F,‘, The transformation given by T,(X)
%&Q,(R). solution
T, of the differential
(3.7) equation
(3.2), defined on the element
= X + t . (N,(X), 0) .
Let us show that the transformation
T, defined by (3.8) satisfies the two previous properties
K, is then
(3.8) (a) and (b).
191
M. Souli, J.P. Zotesio, Shapederivative of discretizedp~o&Ie~
First 1; transforms Mi to M: and the element K to the element K,, second through the transformation the basis function ei on the element K is transformed to the basis function e: on K,.
T,
LEMMA 3.1. Let K be an element on the discrete domain 0, and K, be the element constructed on Q from K when a vertex Mi is moving to Mi, in the k-direction (k = 1,2); thus the transformation 1”, ~5 died by T,= FpF-‘.
(3.9)
PROOF. Without loss of generality, we consider the case k = 1. For t = 0, (3.9) is satisfied; let us show that dT,/dt
= d(F,o F-')fdt
(3.10)
.
Let us write the mapping F, in terms of the shape functions I’$ on the reference
element,
where nElemNodes is the number of vertex nodes per element, and thus nElemNodes = 4. Since K, is isoparamet~cally equivalent to the reference element 2. The transfo~ation T, satisfies Ts(Mi) = Mj for jZi; thus
x;= xi
i#j,
,
x; = xi + t )
y; = y, *
The mappings F, can be written as
@(f>q) + $
Ft( 5, n) = (xi + t, Yi>*
(xi, Yi). q( 6, n)
and F,(h?)=(t,O).~(~,~)~~~~(xj,Yi).~(i,~)=(t,O).~(~,~)+‘(~,~). Thus F, = F + t - (c,
0). Deriving the previous expression with respect to t, we obtain
dF,/‘dt = (g, 0). Let us derive the expression
(3.9) with respect to t:
d(F~~F-l)/dt=(dF*~dt)~F-l=(~,O)~F-l=(~i,O),
(3.11)
using (3.9), we finally obtain dT,/‘dt = (Ni, 0) .
The first property
(a) is then satisfied, and thus T, transforms
To satisfy the second property,
the element K to the element
K,.
q
we need the following lemma.
LEMMA 3.2. For each nodal point j, 1 -C - 1’ d nNodes, the basis function ei on the element K, satisfies eioT,kej.
(3.12)
192
PROOF.
M. Souli, J. P. Zolesio. Shape derivative of discretized problems
Using expression
(3.9)) we have (3.13)
e,foT,=e,foF,oF-‘.
Combining expressions (3.12). 0 COROLLARY
3.1.
(3.12)
using the relations
ej 0 F, = ii and Ei 0 F-
’=
ej, yields
the relation
For each node Mj, the basis function ef satisfies
[d(e:0 T,) /dtlc,,oj= 0 , which means the basis functions are convected.
The proof of the corollary independent.
is straightforward
from Lemma 3.2 since the function
ei 0 T, is time
4. Nodal derivatives In the general case, we are looking for partial derivatives with respect to the coordinates of the vertex nodes defining the discrete domain, for a functional J(f&, u,), where u, is a discretized solution of a variational problem defined on the discrete domain 4, which is called the state formulation. We show in this section how to derive bilinear forms defined in (2.1). To simplify the analysis, we only consider the two-dimensional case. 4.1. Nodal derivatives of bilinear forms Let u’ be the solution of the discrete variational formulation u, defined in 0 through the transformation T,; thus
in the current configuration
0,; let us set
Let us consider a bilinear form arising from a second order problem, on the domain a,, which is a finite union of elements K,; u’ is the solution in the discrete space W,;
(4.1) To compute the material or nodal derivatives of the bilinear form (4.1), let us define the integrals in which involves the Jacobian J, of the (4.1) on the element K, using variable transformation, transformation T,. Vu’ . Vei dx, =
(A(t)Vu,,
Ve,)) dx ,
where e, = e: 0 T, is the basis function on the element K for node i and A(t) is a 2 each element K by A(t) = J,(DT;‘)‘(DT,) .I, = det(DT,)
.
,
(4.2) x
2 matrix defined on
A(0) = Id, (4.3)
To derive the bilinear form (4.2), we first need to derive A(t) with respect to t, at t = 0, we need to derive the expression (4.3), which gives
193
M. Souli, J.P. Zolesio, Shape derivative of discretired problems
A’(0) = [dA(t) /dt] /t = 0 = div(V(0)) - (DV(0)
+ ‘DV(0))
.
Thus, using the fact that the material derivative of the basis function convected, the derivative of the bilinear form (4.2), at t = 0, is i
V(eio T,)) dx/(t = 0) = IK (A’(0)
jK (A(t)Vu,,
ei vanishes,
since they are
-Vu,, Ve,) + (Vzi, Ve,)) dx ,
$ JK F. (eio T,) dx/(t = 0) = IK div(F. u(x, O))e, dx .
(4.4)
For the two-dimensional problem, the matrix A’(0) is straightforwardly computed from (4.3). For x1 and xq coordinate motion of the nodal point Mi, respectively, the convective velocity u(x, 0) is represented by the basis functions (Ni, 0), and (0, Ni), and hence A’(0) is given by
(4.5) respectively. Combining expressions (4.4) and (4.5), the material derivative solution of the following variational problem:
+
u for u = (0, iVi) is the
d,(F- e,). Ni dx .
(4.6)
4.2. Nodal derivative of the functional J Let us consider the functional J to be minimized, I(&
defined by
u’) = n, (u’ - z)’ dx, ,
(4.7)
I
where U’ is a solution of the variational problem (4.2), and z is a real function defined on the current configuration 0,. Transforming the integral in (4.7) on the reference configuration 0, we obtain J(Q, u’) = ; I0 (u’ - z 0 T,)‘) . J, dx .
(4-g)
with respect to t at t = 0 in the u direction is given by
The derivative dJ(0,
u) = R (ti -Vz I
. u(x,
O))(u - z) + $((u - z)* div(u, 0)) dx .
A classical way to avoid the computation of the material derivative u in the solution of the variational problem (4.1) is to introduce an adjoin state,
I
R
Vp.Vwdx=
I
n(u-z).wdx
VWEW.
(4.9)
Using the expressions (4.6) and (4.9), where the state function p plays the role of the shape function w, the gradient dJ takes the form dJ(R,u)=;~KVpVtidx+/K Using expression
- Vz . u(x, O)(u - z) + 1 ((u - z)’ div(u(x, 0)) dx .
(4.6), we obtain the following form of the gradient:
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
194
dJ(f& u> =2
j- alei(-a,u
- a2Ni + alNi. a2u) + a2e,(dlNi.
a+ + a,N,u) dx
KK
+ c 1 KK
-Vz *v(x, O)(u - z) + $(u - z)’ div(u(x, 0)) dx .
(4.10)
The formulation (4.10) is valid for a general velocity field u(x, 0). To compute the partial derivative with respect to the x1 and x2 coordinates of the node Mj, we take u(x, 0) = (0, N,(x)) or u(x, 0) = (Ni , 0). This choice of velocity simplifies the expression (4.10). THEOREM. Let u and p be the solutions of the dkcrete variational problem (4.2) and (4.9), and let J be the functional defined in (4.8). The partial derivatives of J with respect to the x, and x2 coordinates of the node Mi, respectively, are a,u. a,p - a+ a2p) + a2Ni. (aiue alp + a, d,p) dx
(a#‘- Ni)* p - (Ni a,z(u - z) dx + 1 . (u - z)* a,Ni dx
(4.11)
and a,u. alp - a+ a2p) + a,Ni . (alue a2p + a2 alp) dx +c KK
I
(a,(F.N,)*p-Nia2z.(u-z)dx+
$.(u-z)*dZNidx.
(4.12)
These partial derivatives have already been computed by Pironneau [7]. However, the technique in this paper allows us to compute the derivatives (4.11) and (4.12) using the theory of the continuous gradient, and to select a velocity field for the motion of the nodes in such a way that the basis functions are convected, and thus the computation of (4.11) and (4.12) is simplified.
5. Motion of internal nodes The motion of the nodes on the free surface will also serve the purpose of defining the motion of the internal nodes. Each node on the free surface (for example node A), is considered to be the endpoint of a straight line segment, composed of element nodes, which begins at some base point B in space. To describe clearly our method of moving the internal nodes, let us consider a free surface flow for a two-dimensional problem. We assume henceforth that the reference domain is covered by a mesh of isoparametric finite elements, and we think of 0 as the location of the internal nodes. Let us assume that the line y = g(x) coincides with the free surface (Fig. 1). Furthermore, suppose that there is to be no horizontal motion of the nodes and that the vertical motion of the internal nodes is controlled by the
Fig. 1.
M. Souli, 3. P. Zolesio, Shape derivative of discretized problems
195
vertical motion of the nodes on the free surface. Then the following representation of 6 giving a linear distribution of the internal nodes on the segment AB is appropriate, let us define
where nl is the number of boundary vertex nodes, and n2 the number of internal vertex nodes. The application 8 defines the location of the internal nodes in terms of those on the boundary. To prevent the proliferation of notation, let us denote the coordinates of the node Mi by (xi, y,), @(y,, Y,, . * * 9 Y,*) = (Y,1+1,
Yn1+29 * . * 9 Ynl+n2).
(5.1)
Our purpose is to impose a motion on the internal nodes, and so define the application 8, in order to keep a regular mesh for the domain. Using (5. l), the functional J, to be minimized, can be expressed in terms of the coordinates of the vertex nodes n = nl + n2, Ju-4
u) = J(YI7
’ * . 7 YE17 Y,1+17
f * . ) Ynl+n2).
Using (5.1), where the coordinates of the internal nodes are expressed in terms of those on the boundary, the functional J is a function of the coordinates (y,, y,, . . . , y,,), J(y,,
. . . ) Y”l,
@(y,,
. . * 9 Y”1,
Ynl))
= J,(Yl,
Y27.
. *’
Yfzl) .
Using the new formulation Je, to minimize J, we only consider the partial derivative of Je with respect to the y-coordinates of the boundary vertex nodes of the mesh. These derivatives introduce new terms from the fact that the internal nodes depend on the motion of the boundary nodes through 0. Let us assume that 8 is a differentiable function; thus the Jacobian (De), gives the motion of the internal vertex nodes, 1 s i < nl and 1 s j s n2. A straightforward calculation of the partial derivative of J gives
The derivatives aJ/ay, and dJ,/ay, are computed from (4.11) and from the fact that 0 is known explicitly. Let us give an example for the calculation of 8 Je lay,,, . Consider the domain in Fig. 1, n = {(X, y) st a S x ab
and O
The domain is subdivided into quadrangular elements, each line y = i - h, h = (b - a) /n and i = 1 * * , n is divided into m - 1 equal segments. Let us denote by (xi, y,), 1 d i G N, the coordinates of the N vertex nodes, and number of boundary vertex nodes from 1 to nl. However, for an internal node Mj, j > n, (j can be written as j = p. n, + d where d < n,), we have yj = p*yd/(m - 1). This defines the transformation 0 by O( yj) = p*y,/(m - 1). From this expression, we can easily compute the partial derivatives of J using the chain rule.
6. Application
to problems
in hydrodynamics
6.1. The continuous case In this section, the nodal derivative formulation is applied to two different engineering problems. The first one is a numerical model of a wave problem generated by a potential flow around an obstacle resting on the bottom. A critical parameter for this problem is the Froude number F, which is a
196
M. Souli. J.P. Zolesio, Shape derivative of discretized problems
nondimensional parameter defined by F = (C/g * /z)~‘~), w h ere C is the fluid velocity upstream, g is the acceleration of gravity and h is the depth of the fluid. The amplitude wave analysis for the two-dimensional problem is studied in 181. For values of the Froude number less than one, a propagation of a solitary wave over a constant depth is generated by the obstacle. For Froude number greater than one, the wave amplitude decreases away from the obstacle, and the flow becomes uniform downstream. The second problem concerns a flow around an obstacle on the free surface. The difficulty with this problem comes from the singularities on the free surface at the two contact points on the obstacle. This problem is very important from an engineering point of view, but few analytical and numerical results are available. and deduce the equations and boundary conditions Let us explain the physical phenomena, governing the flow. In what follows, we consider the potential and incompressible flow. We define a stream function !P, such that the velocity u = (aP/ay
)
-N/ax).
61.1.
Problem with an obstacle at the bottom For simplicity, we consider a uniform flow upstream and downstream to infinity (see [8,9]) and write the equations concerning a flow around an obstacle resting at the bottom, in a nondimensional stream function formulation:
A?P=O,
0,
F=O,
B,
a?Pl&z=O,
S,,S,,
EW/an=a(y),
S,
(6.1)
where a(y) = (1 - (2 /F*F)*( y - 1))“’ . Let us explain the physical meaning of (6.1). (a) Equation in J2: since the flow is potential, the vorticity vanishes, so the stream function is harmonic. (b) Condition at the bottom B: the bottom B is a streamline; this implies that ly is constant on B and we can fix an arbitrary value of zero. (4 Conditions on lateral boundaries: the flow is uniform upstream and downstream to infinity. If the domain is large enough compared to the obstacle, then the homogeneous Neuman boundary conditions seem reasonable. This boundary condition is valid only for the Froude number greater than 1. In this case, the problem (6.1) is variational and symmetrical. For F < 1, (6.1) is no longer variational and non-symmetrical. Thus, we obtain more complicated boundary conditions on the lateral sides. This is described in detail in (12). (d) Conditions on the free surface: along the free surface S, we dispose of two boundary conditions. The first one is deduced from the Bernoulli equation, and is related to the continuity of the pressure through the free surface; this is given in (6.1), where F is the Froude number and y(x) the nondimensional equation of S. The second condition explains the fact that S is a streamline, so P is constant on S. The constant is fixed by the conservation of the mass flux across a lateral section. In nondimensional parameters, this gives W--l,
s.
(6.2)
The two boundary conditions (6.1), (6.2) on S express the fact that S is an unknown. We formulate the problems (6.1), (6.2) as a minimization problem. The functional to be minimized is J(S) = (l/2)*
k (W - l)* ds ,
where !P is the solution of (6.1).
(6.3)
M. Souli, J. P. Zolesio, Shape derivative
of discretized
197
problem
The shape of the domain fl is defined by the boundary nodes. We only consider the motion of the nodes on S. To simplify the calculations, only the Y-direction of these nodes is computed. Using the technique developed above, the partial derivative of J with respect to the boundary nodes in the Y-direction V = (0, e,) is
dJ(G vi) =
-I\ (A.V'P,
VP) dx + h a( y)*(?P - y) -t l)*B
J(df(du(y)/dy)*(~-
+ s
+ (l/2)*
.y
ds
y) - a(y)*P+ vl: ds
IS (!P - 1)2*B * v ds ,
The matrix A is defined in (3.5) and B = div V- (DV. n, n). The adjoin state P is solution of hP=O,
fi,
P=O,
B,
aP/an=(*-l),
S,
aPjan=o,
s,, s,.
(6.5)
6.1.2. Problem with an obstacle piercing the surface From a mathematical viewpoint, this problem has the same formulation, except that in this case, a part of the free surface C is known. C is the region where the obstacle pierces the free surface S. Thus, the problem is formulated as follows: Av’=O,
0,
?If=O,
B,
aPlafl=O,
S1,S2,
a~/an=~(y), SK, V=O,
C,
with the additional condition P=l,
s‘\c .
6.2. Numerical results The domain 0 in Fig. 1 is subdivided into quadrangles. A Q1 Lagrange finite element approximation is used to solve the variational problems deduced from (6.1) and (6.5) for the state u and the adjoint state I? To check the validity of the gradient, we compare it to the gradient obtained using a finite difference method. Thus, for each node i = 1, . . . , N on the surface S, the finite difference gradient in the Y-direction is
The perturbation of the initial free surface for the computation of the derivative dJi is generally used for an approximation of the Jacobian matrix in nonlinear Newton iterations. The values of dJi are stable with respect to E for E c 10w4. For the first numerical simulation, we consider the initial free surface y = 1. The flow is uniform upstream, and the obstacle is sinusoidal (see Fig. 2). The Froude number F = 2.55, the length of the
198
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
0.690 I 0.460
s
n
%
It
4.c Fmude I 2.55
homqenous
neumann
Fig. 2.
domain L = 1, and the mesh size h = 0.1. The values of the gradient dJ for the nodes on the free surface are compared in Table 1. We give numerical results for F < 1, for which the mathematical formulation is more complicated because of the boundary conditions on the lateral sides S, and S,; these are no longer the homogeneous Neuman boundary conditions (6.1) but satisfy the continuity of the mass flux through S, and S,. Outside the domain 0, the linearized problem for ly has an analytical solution. The gradient which obviously has a more complicated form is given in [lo]. To search for the free surface, or water wave, we use the minimization algorithm in (1)) where an approximation of the Hessian matrix of the gradient is used in order to obtain the optimal slope. To distinguish the water wave from the initial wave, y = 1, we take a long domain L = 4 (Fig. 3). The convergence is fast, because the problem (6.1) is variational; we reach the steady state after 4 iterations, and the cost decreases from J, = 1.23 to J4 = 0.65 x 10V4. This case is also treated in [19] where a fixed point method is used. Figure 4 corresponds to F = 0.5, as described in [9]. One observes a solitary wave propagating to the right of the obstacle; this wave is due to the special boundary condition on the lateral sides. Figure 5 corresponds to a flow around an obstacle piercing the surface. The results are analogous to those given in [9] with the fixed point method. The hxed point method in [9] is more efficient and faster than the gradient method, when it converges. But it fails for certain cases, whereas the gradient method converges. This case corresponds to Fig. 4 with F = 0.4 and the amplitude of the sinusoidal obstacle is 0.4. The gradient method converges after 40 iterations, and the cost decreases from J, = 0.270 to J40 = 0.12 X 10P4. For the search of the water wave in naval hydrodynamics problems, or in general free surface problems, we should first test the fixed point method, because it is faster and less expensive. The convergence is geometrical; for
Table 1 Gradient in (6.6)
Gradient in (6.4)
0.7453E-00 -O.l129E-00 -0.2367E-00 -0.2835E-00 -0.2577E-00 -O.l892E-00 -O.l105E-01 -0.4493E-01 -0.3621E-02 O.l433E-01 O.l116E-01
0.7455E-00 -O.l13OE-00 -0.2368E-00 -0.2839E-00 -0.2573E-00 -O.l892E-00 -O.l108E-01 -0.4498E-01 -0.3634E-01 O.l438E-01 O.l12OE-01
0
1.0
2.0
Fmtide
- 0.50
Fig. 3.
3.0
obstacle - 0.2
4.0
M. Souli, J. P. Zolesio, Shape derivative of discretized problems
0
0.2
0.4
Frovde=OSO
Fig. 4.
0.6
0.8
obstacle.O.t,
1.0
0
0.5
199
1.0
1.5
2.0
Frwde -0.64 obstacle -0.05
Fig. 5
this method; at each iteration we only need to solve a linear system corresponding to the state problem. In case this method fails, the gradient method may converge, at least for some cases. The gradient method is more expensive; at each iteration, we solve two linear systems, the state and the adjoin state, and the convergence as it is known, is not geometrical.
Acknowledgment
This research was supported by DRET contract with the collaboration of ENSTA at Palaiseau. The authors would like to thank M. Lenoir and J. Cahouet at ENSTA for helpful discussions. The first author wishes to thank Charlotte Van Campen at Centric engineering for the help in making the figures.
References [l] J.P. ZolCsio, Les derives par rapport aux noeuds des triangularisations et leurs utilisation en identification de domaine, Ann. Sci. Math. Quebec 8 (1984) 97-120. [2] J. Cea, Problems of shape optimal design, in: J. Cea and Haug, eds, Optimization of Distributed Parameters Structures, Vol 12 (Sijthoff, Alphen aan den Rijn, 1981) 1005-1088. [3] J.P. Zolbsio, Identification de domaines par deformations, Thesis, Nice, 1979. [4] M. Souli and J.P. Zoltsio, Gradient with respect to a non smooth domain in shape optimization, Optim. Control Applic. Control, in press. [5] M. Delfour, M. Payre and J.P. Zoltsio, Optimal and suboptimal design of thermal diffusers for communications satellites, Comptes rendus du 3ieme Symposium on Control of Distributed Parameters Systems (1982). [6] P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, Amsterdam, 1978). [7] 0. Pironneau, Optimal Shape Design for Elliptic Systems (Springer, Berlin, 1986). [8] J. Cahouet and M. Lenoir, Resolution numerique du problbme non lineaire de la resistance de vague non lineaire, CRAS 297, 1985. [9] J. Cahouet, Etude numCrique et expi%mentale du problbme de la resistance de vague non lineaire, Thesis, Paris, 1984. [lo] M. Souli and J.P. Zolbsio, Le controle de domaine dans la resistance de vague non lintaire en hydrodynamique, Ann Sci. Quebec 1 (1990) 121-132. [ll] F. Neittamaki, Finite Element Approximation for Optimal Shape Design. Theory and Application (Wiley, New York, 1987). [12] A. Buckeley, An alternative implementation of Goldfard’s minimization algorithms, Math. Prog. 8 (1975) 207-231.