- Email: [email protected]
Shape-topology optimization of stokes flow via variational level set method
Available online at www.sciencedirect.com
Applied Mathematics and Computation 202 (2008) 200–209 www.elsevier.com/locate/amc
Shape-topology optimization of stokes flow via variational level set method q Xian-Bao Duan a,*, Yi-Chen Ma a, Rui Zhang b a
b
School of Science, Xi’an Jiaotong University, Xi’an 710049, PR China Institute of Applied Mathematics, Shandong University of Technology, Zibo 255049, PR China
Abstract In this paper, the conventional level set method is extended as an efficient approach for shape optimization. The motion of the dynamic interface is governed by Stokes equations. Although this method is not specifically designed for topology optimization, it can handle topology changes easily. Furthermore, a relatively smooth evolution can be maintained without re-initialization. The cost of this method is moderate. The accuracy and efficiency of the proposed algorithm are illustrated by several numerical examples. Ó 2008 Elsevier Inc. All rights reserved. Keywords: Shape optimization; Topology optimization; Stokes flow; Variational level set method; Inverse problem; Domain identification
1. Introduction Shape-topology optimization of the fluid flow is a very important and popular field and the applications of the optimal shape-topology design are uncountable [1]. The classical method of shape sensitivity, which has been much studied [2], is a very general method that can handle any type of objective functional and fluid models. But it has at least two main drawbacks: its computational cost and its prone to fail into local minima far away from global ones [3]. Firstly devised by Osher and Sethian [3–5], level set method is a collection of algorithms for solving a particular class of partial differential equations (PDE) for numerically tracking fronts and free boundaries [6]. The level set method is easy to implement and versatile for computing and analyzing the motion of an interface in two or three dimensions. Since it can easily develop sharp corners, break apart, and merge together in a robust and stable way, the level set method has wide range of successful applications, such as fields motion by mean curvature, graphics, image processing, computational fluid dynamics (CFD), materials science and so on.
q *
This work is supported by the National Natural Science Foundation of China, Grant Nos. 10671153 and 10371096. Corresponding author. E-mail address: [email protected] (X.-B. Duan).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.02.014
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
201
The level set method has been recently introduced in the field of shape optimization, enabling a smooth representation of the boundaries on a fixed mesh and therefore leading a fast numerical algorithms [3,7– 13]. The basic idea of level set method for optimization problems is to choose the direction (velocity) in such a way that a decrease of the objective functional is achieved, which assembles the classical speed method in shape optimization. Moreover, the weak formulation via level set method allows for more general evolutions and in particular for topological changes such as splitting or merging of domains. In conventional level set method, the moving fronts or interface are represented by the zero level set. But this may not maintain during the evolution. To solve this problem, re-initializing the level set function as a signed distance function during the evolution procedure has been extensively used as a remedial measure. But if the level set function is not smooth or is much steeper on one side of the interface than the other, the zero level set of the resulting function may be moved incorrectly from that of the origin function [4,14]. Furthermore, if the level set function is far away from a signed distance function, this remedy may not preserve to a signed distance function. In order to deal with this difficulty, Chan and Vese [15] and Li et al. [16] proposed a new variational formulation of level set method that needs no re-initialization at all. In the present paper, we study the shape-topology optimization of the Stokes flow based on the method proposed by Chan and Li. The rest of the paper is organized as follows: In Section 2, we give the fundamental mathematical formulations that we shall use in this paper. Section 3 is dedicated to the modal and formulas of our method, and in Section 4 several numerical examples are provided. The conclusions are given in the last section. 2. Mathematical formulations In this section, we shall give the mathematical formulations of the Stokes problem and level set method that we shall use in this paper, following mainly [4,14,17,18]. 2.1. The stokes problems Assume X is a bounded open set in R2 with a Lipschitz boundary C :¼ oX, and f 2 L2 ðXÞ is a given vector function. We denote by Q the cylinder X ð0; T Þ, where T is positive and fixed. We are looking for a vector function u ¼ ðu1 ; . . . ; un Þ and a scalar function p, representing, respectively, the velocity and the pressure of the fluid, which are defined in Q and satisfy the following equations and boundary/initial conditions: 8 ou mDu þ ðu rÞu þ rp ¼ f ; in Q; > ot > > < div u ¼ 0; in Q; ð1Þ > u ¼ g; on C ½0; T ; > > : uðx; 0Þ ¼ u0 ðxÞ; in X; where m is the constant kinematic viscosity and f is the body force. For slow and steady flow, the convection term ðu rÞu can be neglected and the time derivative ou=ot vanishes. Eq. (1) then reduce to the Stokes equations 8 > < mDu þ rp ¼ f ; in X; div u ¼ 0; in X; ð2Þ > : u ¼ g; on C: In general, the velocity u is prescribed to some known function g on C. By integrating the incompressibility condition div u ¼ 0 over X and using the divergence theorem it is easy to see that the function g must satisfy the following compatibility condition: Z g n dC ¼ 0; ð3Þ C
where n denotes the outward unit normal.
202
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
If f ; u and p are smooth functions and satisfy Eq. (2), then take the scalar product of the first equation of (2) with a test function v 2 V, where V ¼ fv 2 DðXÞ; div v ¼ 0g;
ð4Þ 1
and DðXÞ represents the space of C functions with compact support contained in X, we obtain the variational form of Eq. (2) Z Z Z m ru rv dX p div v dX ¼ f v dX; ð5Þ X
X
X
after having used the Green’s theorem. Using the bilinear form Z aðu; vÞ ¼ m ru rv dX;
ð6Þ
X
and the functionals Z bðv; pÞ ¼ p div v dX; ZX f v dX; hf ; viX ¼
ð7Þ ð8Þ
X
the equation of virtual power (5) can be compactly written as aðu; vÞ þ bðv; pÞ ¼ hf ; viX : The second equation of (2) multiplied by any pressure function q 2 bðu; qÞ ¼ 0:
ð9Þ L20 ðXÞ
and integrated over X is ð10Þ
The weak form of Stokes problem (2) can be formulated as follows: Find ðu; pÞ 2 H10 ðXÞ L20 ðXÞ such that ( aðu; vÞ þ bðv; pÞ ¼ hf ; viX 8v 2 H10 ðXÞ; ð11Þ bðu; qÞ ¼ 0 8q 2 L20 ðXÞ; R where L20 ðXÞ ¼ fq 2 L2 ðXÞj X q dX ¼ 0g. The Stokes problem (11) is uniquely solvable [19]. 2.2. The conventional level set method The main idea of the level set method is to represent curves or surfaces not via parameterizations, but in an implicit form as the zero level set or isosurface of a high dimension continuous function /, the so-called level set function, and then trace the deformation of the curves or surfaces by means of the deformation of this embedding function. So the implicit function can be used for both representing and evolving the interfaces or boundaries. For a given domain X with smooth boundary C :¼ oX, we assume there exists an implicit function / satisfies 8 > < /ðxÞ < 0 8x 2 X; /ðxÞ ¼ 0 8x 2 C; ð12Þ > : /ðxÞ > 0; others: The local unit normal n to the surface is given by r/ n¼ : ð13Þ jr/j In the level set method, it is convenient to use the Heaviside distribution H ð/Þ and the Dirac d-distribution dð/Þ defined as 1; / P 0; H ð/Þ ¼ ð14Þ 0; / < 0;
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
203
and dð/Þ ¼
DH ð/Þ ¼ dð/Þjr/j; D/
ð15Þ
respectively. Then, the interior and the boundary C of a shape can be described in terms of the level set function /, respectively, as X ¼ fx : H ð/Þ ¼ 1g
ð16Þ
C ¼ fx : dð/Þ > 0g:
ð17Þ
and For a function F ðxÞ, the volume integration is defined as Z F ðxÞH ð/ÞdX:
ð18Þ
X
If F ðxÞ ¼ 1, Eq. (18) yields the volume of the domain X as follows: Z V ðXÞ ¼ H ð/ÞdX:
ð19Þ
X
The surface integration of function F ðxÞ can also be defined as Z Z F ðxÞH ð/ÞdC ¼ F ðxÞdð/Þjr/jdX: C
ð20Þ
X
In order to let the level set function dynamically changes in time, a continuous velocity field V, which is a function of position x, is introduced and the evolution can be described as the following Cauchy problem: ( o/ þ V n r/ ¼ 0; ot ð21Þ /ðx; 0Þ ¼ /0 ðxÞ; where Vn ¼V
r/ jr/j
ð22Þ
is the normal velocity of the moving boundary. In the conventional level set method, it is numerically necessary to keep the evolving level set function / close to a signed distance function [4]. Re-initializing the level set function to a signed distance function during the evolution has been extensively used as a numerical remedy for maintaining stable curve evolution and ensuring desired results. Indeed, the re-initialization process is crucial and cannot be avoided in the conventional level set method [4,14]. But from the practical point of view, the re-initialization process is quite complicated, expensive, and have subtle side effects. Moreover, most of the level set methods are fraught with their own problems, such as when and how to re-initialize the level set function to a signed distance function [16,20]. In order to circumvent this difficulty, Chan and Vese [15] and Li et al. [16] proposed a new variational formulation of level set method without re-initialization. The algorithm we propose in this paper is mainly based on the ideas of Chan and Li. 3. Variational formulation of level set method In this section, we will develop the shape-topology optimization algorithm for a domain immersed in an incompressible fluid described by the Stokes equations. The shape optimization of a domain (obstacle) immersed in a stationary, viscous, incompressible fluid can be formulated as find a velocity field u such that the cost functional Z 1 JðX; uÞ :¼ ðu u Þ2 dX ð23Þ 2 X
204
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
is minimized, subject to the Stokes equations 8 > < mDu þ rp ¼ f ; in X; div u ¼ 0; in X; > : u ¼ u0 ; on C:
ð24Þ
The goal of the optimization problem is to identify the shape where the target velocity field u is reached. As mentioned above, it is crucial to keep the level set function as a distance function during the evolution. It is well known that a signed distance function must satisfy the property jr/ j¼ 1, and on the contrary, if any function / satisfies jr/j ¼ 1, then / is a signed distance function plus a constant [21]. So we will use the following functional as a metric to characterize how close a function / is to a signed distance function in X 2 R2 [16], Z 1 ðjr/j 1Þ2 dX: EI ¼ ð25Þ X 2 This metric plays an important role in our variational formulation. The energy functional can be defined as E ¼ k1 EI þ k0 JðX; uÞ:
ð26Þ
Throughout this paper, ki ði ¼ 0; 1; 2; 3; 4Þ will denote various positive constants. In our proposed algorithm, the energy functional (26) will be minimized, and some regularizing terms, such as the length of the curve (boundary) or the area of the region will be added. The regularizing terms are defined by Z Z ER :¼ k3 G dC þ k4 G dX; ð27Þ oX
X
where G is the edge indicator function G which is defined by [22] G¼
1 1 þ jrGr uj2
ð28Þ
;
and Gr u, a smoother version of u, is the convolution of the data u with the Gaussian kernel Gr ðx; yÞ ¼ r1=2 expfjx2 þ y 2 j=4rg. Therefore, the final total energy functional reads: ð29Þ E ¼ k1 EI þ k2 EF þ ER ; R where EF :¼ X Gðu u Þ2 dX and k2 ¼ k0 =2. Using Eqs. (18) and (20), the energy functional E can be rewritten as Z Z Z 1 ðjr/j 1Þ2 dX þ k2 Gðu u Þ2 H ð/ÞdX þ k3 Gdð/Þjr/j dX þ k4 Eð/Þ ¼ k1 X 2 X X Z GH ð/ÞdX: ð30Þ X
1.2
ε = 0.1 ε = 0.05
1
ε = 0.8 ε =1
1.4 1.2
0.8
1
0.6
0.8 0.6
0.4
0.4 0.2 0.2 0
0
—0.2 —1.5
—0.2 —1.5
—1
—0.5
0
0.5
1
1.5
—1
—0.5
0
0.5
Fig. 1. The regularization of the Heaviside function and d-function.
1
1.5
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
205
The variational formulation (30) will be applied in our shape optimization algorithm. In order to obtain the Euler–Lagrange equation of the above functional with respect to the unknown function /, the regularized version of the Heaviside distribution H ð/Þ and d-distribution dð/Þ, denoted, respectively, as H e ð/Þ and de ð/Þ ¼ H 0e ð/Þ are employed. Then we get the final energy functional Z Z 1 2 2 ðjr/j 1Þ dX þ k2 Gðu u Þ H e ð/ÞdX Ee ð/Þ ¼ k1 2 X Z ZX þ k3 Gde ð/Þjr/j dX þ k4 GH e ð/ÞdX: ð31Þ X
X
Fig. 2. The domain and the boundary conditions.
Fig. 3. The evolution of the level set functions in the first example.
206
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
Minimizing Ee ð/Þ with respect to /, we can deduce the corresponding Euler–Lagrange equation, which can be used to evolve the level set function / by an artificial time t P 0: 8 h i o/ r/ Gr/ 2 > > < ot ¼ k1 D/ div jr/j þ k2 Gde ð/Þðu u Þ þ k3 de ð/Þdiv jr/j þ k4 Gde ð/Þ; in X; ð32Þ /ð0; xÞ ¼ /0 ðxÞ; in X; > > : o/ ¼ 0; on oX; ox where /ð0; xÞ ¼ /0 ðxÞ is the initial contour and n is the exterior normal to the boundary oX. In our proposed algorithm, the level set function will be evolved by (32). Because the first term of functional (31) automatically maintained / as a signed distance function during the level set evolution, the re-initialization procedure can be completely eliminated. 4. Numerical examples In this section, we will illustrate the efficiency and applicability of the developed shape-topology optimization algorithm by presenting several numerical examples in two dimensions. From the computational point of view, the smeared-out Heaviside and Dirac function are preferred during optimization process. We will use the following regularized Heaviside and Dirac distribution (as shown in Fig. 1) in our numerical examples 8 if x > e; < 1; if x < e; H e ðxÞ ¼ 0; ð33Þ px :1 x 1 ; if jxj 6 e: 1 þ þ sin e p e 2 0; if jxj > e; px de ðxÞ ¼ 1 ð34Þ ; if jxj 6 e: 1 þ cos 2e e
Fig. 4. The evolution of the level set functions in the second example.
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
207
We demonstrate the application of shape optimization algorithm in a Stokes flow using the model as Fig. 2. A fluid flows from left to right in a rectangle region surrounding the domain (obstacle) X0 we want to reshape. No slip and impermeability boundary conditions have been used on all walls except at inlet and outlet. The domain and the boundary conditions are as below. At the top and bottom wall u ¼ 0; v ¼ 0; at the inlet u ¼ 1; v ¼ 0 and at the outlet p ¼ 0. Unless stated otherwise, in all the experimental results shown in this section the following parameters are assumed: k1 ¼ k2 ¼ k3 ¼ 1; k4 ¼ 0. The time step s can be chosen much larger than conventional level set method. But in order to ensure the stability of the evolution, we use s 6 10 in all our examples. We apply a Taylor–Hood element (P 2 P 1 ) discretization of the Stokes problem [23] and the update of the level set function by finite difference method [24]. The present algorithm will be terminated when the difference between two successive objective functional values is less than the prescribed error or when the given maximum number of iterations has been reached. Fig. 3 displays the contours of the level set functions and the velocity fields of the first numerical example by the proposed level set method. The zero-contour of the initial level set function /0 , which is assumed to be a rectangle that encloses the target domain, is shown in the top-left figure. The top-right and the bottom-left figures show the zero-contours after 50 iterations and 100 iterations, respectively, and the bottom-right is the final zero-contour. It can be seen that the final optimal shape is quite satisfiable. As for the second numerical example, the evolution of the level sets of / and the target velocity field which was generated for the obstacle that we want to reshape are shown in Fig. 4. The figure on the top-left shows the initial shape, a rectangle included in the target optimal shape. The figures on the top-right and bottom-left show the shape obtained after 20 iterations and 40 iterations, respectively. It can be seen that we are able to closely identify the shape that generated the target velocity field despite the initial shape and position.
Fig. 5. The evolution of the level set functions in the third example.
208
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
In the third numerical example, the initial shape is two circles. The typical evolution shapes are shown in Fig. 5. Due to the level set implementation, the algorithm allows automatical change of the topology. This example shows the proposed algorithm can handle topology change automatically during the optimization process. It can be concluded from the numerical examples, the presented algorithm is successful in accuracy, convergence speed and insensitivity to initial shape or topology. 5. Conclusions In this paper, we have proposed a variational level set method that completely eliminates the re-initialization procedure. This method has the following advantages over the conventional level set method: the re-initialization procedure can be completely eliminated and a significantly large artificial time can be used for numerically solving the evolutionary PDE, therefore speeding up the evolution procedure. The level set function can be initialized with any distance function with any shape or topology that are more efficient to construct and easier to use in practice. The level set method in our formulation can be easily implemented and is computationally more efficient. Acknowledgement The authors would like to thank the editor and referees for their valuable comments and suggestions that help us to improve this paper. References [1] B. Mohammadi, O. Pironneau, Shape optimization in fluid mechanics, Annual Review of Fluid Mechanics 36 (2004) 255–279. [2] J. Sokolowski, J.P. Zolesio, Introduction to Shape Optimization: Shape Sensitivity Analysis, Springer Series in Computational Mathematics, vol. 10, Springer, Berlin, 1992. [3] G. Allaire, F. Jouve, A.M. Toader, A level-set method for shape optimization, CR Academic Sciences Paris, Series I 334 (2002) 1125– 1130. [4] S. Osher, J.A. Sethian, Level Set Methods and Dynamic Implicit Surface, Springer-Verlag, 2002. [5] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, Journal of Computational Physics 79 (1988) 12–49. [6] M. Sussmann, P. Smereka, S. Osher, A level set approach for computing solution to incompressible two-phase flow, Journal of Computational Physics 114 (1994) 114–146. [7] M.Y. Wang, X.M. Wang, D.M. Guo, A level set method for structural topology optimization, Computer Methods in Applied Mechanics and Engineering 192 (2003) 227–246. [8] F. Santosa, A level-set approach for inverse problems involving obstacles, ESAIM: Control, Optimization, and Calculus of Variation 1 (1996) 17–33. [9] X.B. Duan, Y.C. Ma, R. Zhang, Optimal shape control of fluid flow using variational level set method, Physics Letters A 372 (9) (2008) 1374–1379. [10] X.B. Duan, Y.C. Ma, R. Zhang, Shape-topology optimization for Navier–Stokes problem using variational level set method, Journal of Computational and Applied Mathematics (2007), doi:10.1016/j.cam.2007.11.016. [11] A. Litman, D. Lesselier, F. Santosa, Reconstruction of 2-D binary obstacle by controlled evolution of a level-set, Inverse Problems 14 (1998) 685–706. [12] M. Burge, S. Osher, A survey on level set methods for inverse problems and optimal design, European Journal of Applied Mathematics 16 (2005) 263–301. [13] S. Amstutz, H. Andra¨, A new algorithm for topology optimization using a level-set method, Journal of Computational Physics 216 (2) (2006) 573–588. [14] J.A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science, Cambridge University Press, 1999. [15] T. Chan, L. Vese, Active contours without edges, IEEE Transactions on Image Processing 10 (2001) 266–277. [16] C.M. Li, C.Y. Xu, C.F. Gui, M.D. Fox, Level set evolution without re-initialization: a new variational formulation, in: CVPR’ 05. [17] R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, AMS Chelsea Publishing, 2000. [18] M.E. Gurtin, An Introduction to Continuum Mechanics, Academic Press, San Diego, 1981. [19] V. Girault, P. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, Heidelberg, 1986.
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
209
[20] J. Gomes, O. Faugeras, Reconciling distance functions and level sets, Journal of Visual Communication and Image Representation 114 (1994) 146–150. [21] V.I. Arnold, Geometrical Methods in the Theory of Ordinary Differential Equations, Springer-Verlag, New York, 1983. [22] M. Kass, A. Witkin, D. Terzopoulos, Snakes: active contour models, International Journal of Computer Vision 1 (1998) 321–331. [23] V. Girault, P. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, Heidelberg, 1986. [24] H.K. Zhao, T. Chan, B. Merriman, S. Osher, A variational level set approach to multiphase motion, Journal of Computational Physics 127 (1996) 179–195.
Applied Mathematics and Computation 202 (2008) 200–209 www.elsevier.com/locate/amc
Shape-topology optimization of stokes flow via variational level set method q Xian-Bao Duan a,*, Yi-Chen Ma a, Rui Zhang b a
b
School of Science, Xi’an Jiaotong University, Xi’an 710049, PR China Institute of Applied Mathematics, Shandong University of Technology, Zibo 255049, PR China
Abstract In this paper, the conventional level set method is extended as an efficient approach for shape optimization. The motion of the dynamic interface is governed by Stokes equations. Although this method is not specifically designed for topology optimization, it can handle topology changes easily. Furthermore, a relatively smooth evolution can be maintained without re-initialization. The cost of this method is moderate. The accuracy and efficiency of the proposed algorithm are illustrated by several numerical examples. Ó 2008 Elsevier Inc. All rights reserved. Keywords: Shape optimization; Topology optimization; Stokes flow; Variational level set method; Inverse problem; Domain identification
1. Introduction Shape-topology optimization of the fluid flow is a very important and popular field and the applications of the optimal shape-topology design are uncountable [1]. The classical method of shape sensitivity, which has been much studied [2], is a very general method that can handle any type of objective functional and fluid models. But it has at least two main drawbacks: its computational cost and its prone to fail into local minima far away from global ones [3]. Firstly devised by Osher and Sethian [3–5], level set method is a collection of algorithms for solving a particular class of partial differential equations (PDE) for numerically tracking fronts and free boundaries [6]. The level set method is easy to implement and versatile for computing and analyzing the motion of an interface in two or three dimensions. Since it can easily develop sharp corners, break apart, and merge together in a robust and stable way, the level set method has wide range of successful applications, such as fields motion by mean curvature, graphics, image processing, computational fluid dynamics (CFD), materials science and so on.
q *
This work is supported by the National Natural Science Foundation of China, Grant Nos. 10671153 and 10371096. Corresponding author. E-mail address: [email protected] (X.-B. Duan).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.02.014
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
201
The level set method has been recently introduced in the field of shape optimization, enabling a smooth representation of the boundaries on a fixed mesh and therefore leading a fast numerical algorithms [3,7– 13]. The basic idea of level set method for optimization problems is to choose the direction (velocity) in such a way that a decrease of the objective functional is achieved, which assembles the classical speed method in shape optimization. Moreover, the weak formulation via level set method allows for more general evolutions and in particular for topological changes such as splitting or merging of domains. In conventional level set method, the moving fronts or interface are represented by the zero level set. But this may not maintain during the evolution. To solve this problem, re-initializing the level set function as a signed distance function during the evolution procedure has been extensively used as a remedial measure. But if the level set function is not smooth or is much steeper on one side of the interface than the other, the zero level set of the resulting function may be moved incorrectly from that of the origin function [4,14]. Furthermore, if the level set function is far away from a signed distance function, this remedy may not preserve to a signed distance function. In order to deal with this difficulty, Chan and Vese [15] and Li et al. [16] proposed a new variational formulation of level set method that needs no re-initialization at all. In the present paper, we study the shape-topology optimization of the Stokes flow based on the method proposed by Chan and Li. The rest of the paper is organized as follows: In Section 2, we give the fundamental mathematical formulations that we shall use in this paper. Section 3 is dedicated to the modal and formulas of our method, and in Section 4 several numerical examples are provided. The conclusions are given in the last section. 2. Mathematical formulations In this section, we shall give the mathematical formulations of the Stokes problem and level set method that we shall use in this paper, following mainly [4,14,17,18]. 2.1. The stokes problems Assume X is a bounded open set in R2 with a Lipschitz boundary C :¼ oX, and f 2 L2 ðXÞ is a given vector function. We denote by Q the cylinder X ð0; T Þ, where T is positive and fixed. We are looking for a vector function u ¼ ðu1 ; . . . ; un Þ and a scalar function p, representing, respectively, the velocity and the pressure of the fluid, which are defined in Q and satisfy the following equations and boundary/initial conditions: 8 ou mDu þ ðu rÞu þ rp ¼ f ; in Q; > ot > > < div u ¼ 0; in Q; ð1Þ > u ¼ g; on C ½0; T ; > > : uðx; 0Þ ¼ u0 ðxÞ; in X; where m is the constant kinematic viscosity and f is the body force. For slow and steady flow, the convection term ðu rÞu can be neglected and the time derivative ou=ot vanishes. Eq. (1) then reduce to the Stokes equations 8 > < mDu þ rp ¼ f ; in X; div u ¼ 0; in X; ð2Þ > : u ¼ g; on C: In general, the velocity u is prescribed to some known function g on C. By integrating the incompressibility condition div u ¼ 0 over X and using the divergence theorem it is easy to see that the function g must satisfy the following compatibility condition: Z g n dC ¼ 0; ð3Þ C
where n denotes the outward unit normal.
202
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
If f ; u and p are smooth functions and satisfy Eq. (2), then take the scalar product of the first equation of (2) with a test function v 2 V, where V ¼ fv 2 DðXÞ; div v ¼ 0g;
ð4Þ 1
and DðXÞ represents the space of C functions with compact support contained in X, we obtain the variational form of Eq. (2) Z Z Z m ru rv dX p div v dX ¼ f v dX; ð5Þ X
X
X
after having used the Green’s theorem. Using the bilinear form Z aðu; vÞ ¼ m ru rv dX;
ð6Þ
X
and the functionals Z bðv; pÞ ¼ p div v dX; ZX f v dX; hf ; viX ¼
ð7Þ ð8Þ
X
the equation of virtual power (5) can be compactly written as aðu; vÞ þ bðv; pÞ ¼ hf ; viX : The second equation of (2) multiplied by any pressure function q 2 bðu; qÞ ¼ 0:
ð9Þ L20 ðXÞ
and integrated over X is ð10Þ
The weak form of Stokes problem (2) can be formulated as follows: Find ðu; pÞ 2 H10 ðXÞ L20 ðXÞ such that ( aðu; vÞ þ bðv; pÞ ¼ hf ; viX 8v 2 H10 ðXÞ; ð11Þ bðu; qÞ ¼ 0 8q 2 L20 ðXÞ; R where L20 ðXÞ ¼ fq 2 L2 ðXÞj X q dX ¼ 0g. The Stokes problem (11) is uniquely solvable [19]. 2.2. The conventional level set method The main idea of the level set method is to represent curves or surfaces not via parameterizations, but in an implicit form as the zero level set or isosurface of a high dimension continuous function /, the so-called level set function, and then trace the deformation of the curves or surfaces by means of the deformation of this embedding function. So the implicit function can be used for both representing and evolving the interfaces or boundaries. For a given domain X with smooth boundary C :¼ oX, we assume there exists an implicit function / satisfies 8 > < /ðxÞ < 0 8x 2 X; /ðxÞ ¼ 0 8x 2 C; ð12Þ > : /ðxÞ > 0; others: The local unit normal n to the surface is given by r/ n¼ : ð13Þ jr/j In the level set method, it is convenient to use the Heaviside distribution H ð/Þ and the Dirac d-distribution dð/Þ defined as 1; / P 0; H ð/Þ ¼ ð14Þ 0; / < 0;
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
203
and dð/Þ ¼
DH ð/Þ ¼ dð/Þjr/j; D/
ð15Þ
respectively. Then, the interior and the boundary C of a shape can be described in terms of the level set function /, respectively, as X ¼ fx : H ð/Þ ¼ 1g
ð16Þ
C ¼ fx : dð/Þ > 0g:
ð17Þ
and For a function F ðxÞ, the volume integration is defined as Z F ðxÞH ð/ÞdX:
ð18Þ
X
If F ðxÞ ¼ 1, Eq. (18) yields the volume of the domain X as follows: Z V ðXÞ ¼ H ð/ÞdX:
ð19Þ
X
The surface integration of function F ðxÞ can also be defined as Z Z F ðxÞH ð/ÞdC ¼ F ðxÞdð/Þjr/jdX: C
ð20Þ
X
In order to let the level set function dynamically changes in time, a continuous velocity field V, which is a function of position x, is introduced and the evolution can be described as the following Cauchy problem: ( o/ þ V n r/ ¼ 0; ot ð21Þ /ðx; 0Þ ¼ /0 ðxÞ; where Vn ¼V
r/ jr/j
ð22Þ
is the normal velocity of the moving boundary. In the conventional level set method, it is numerically necessary to keep the evolving level set function / close to a signed distance function [4]. Re-initializing the level set function to a signed distance function during the evolution has been extensively used as a numerical remedy for maintaining stable curve evolution and ensuring desired results. Indeed, the re-initialization process is crucial and cannot be avoided in the conventional level set method [4,14]. But from the practical point of view, the re-initialization process is quite complicated, expensive, and have subtle side effects. Moreover, most of the level set methods are fraught with their own problems, such as when and how to re-initialize the level set function to a signed distance function [16,20]. In order to circumvent this difficulty, Chan and Vese [15] and Li et al. [16] proposed a new variational formulation of level set method without re-initialization. The algorithm we propose in this paper is mainly based on the ideas of Chan and Li. 3. Variational formulation of level set method In this section, we will develop the shape-topology optimization algorithm for a domain immersed in an incompressible fluid described by the Stokes equations. The shape optimization of a domain (obstacle) immersed in a stationary, viscous, incompressible fluid can be formulated as find a velocity field u such that the cost functional Z 1 JðX; uÞ :¼ ðu u Þ2 dX ð23Þ 2 X
204
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
is minimized, subject to the Stokes equations 8 > < mDu þ rp ¼ f ; in X; div u ¼ 0; in X; > : u ¼ u0 ; on C:
ð24Þ
The goal of the optimization problem is to identify the shape where the target velocity field u is reached. As mentioned above, it is crucial to keep the level set function as a distance function during the evolution. It is well known that a signed distance function must satisfy the property jr/ j¼ 1, and on the contrary, if any function / satisfies jr/j ¼ 1, then / is a signed distance function plus a constant [21]. So we will use the following functional as a metric to characterize how close a function / is to a signed distance function in X 2 R2 [16], Z 1 ðjr/j 1Þ2 dX: EI ¼ ð25Þ X 2 This metric plays an important role in our variational formulation. The energy functional can be defined as E ¼ k1 EI þ k0 JðX; uÞ:
ð26Þ
Throughout this paper, ki ði ¼ 0; 1; 2; 3; 4Þ will denote various positive constants. In our proposed algorithm, the energy functional (26) will be minimized, and some regularizing terms, such as the length of the curve (boundary) or the area of the region will be added. The regularizing terms are defined by Z Z ER :¼ k3 G dC þ k4 G dX; ð27Þ oX
X
where G is the edge indicator function G which is defined by [22] G¼
1 1 þ jrGr uj2
ð28Þ
;
and Gr u, a smoother version of u, is the convolution of the data u with the Gaussian kernel Gr ðx; yÞ ¼ r1=2 expfjx2 þ y 2 j=4rg. Therefore, the final total energy functional reads: ð29Þ E ¼ k1 EI þ k2 EF þ ER ; R where EF :¼ X Gðu u Þ2 dX and k2 ¼ k0 =2. Using Eqs. (18) and (20), the energy functional E can be rewritten as Z Z Z 1 ðjr/j 1Þ2 dX þ k2 Gðu u Þ2 H ð/ÞdX þ k3 Gdð/Þjr/j dX þ k4 Eð/Þ ¼ k1 X 2 X X Z GH ð/ÞdX: ð30Þ X
1.2
ε = 0.1 ε = 0.05
1
ε = 0.8 ε =1
1.4 1.2
0.8
1
0.6
0.8 0.6
0.4
0.4 0.2 0.2 0
0
—0.2 —1.5
—0.2 —1.5
—1
—0.5
0
0.5
1
1.5
—1
—0.5
0
0.5
Fig. 1. The regularization of the Heaviside function and d-function.
1
1.5
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
205
The variational formulation (30) will be applied in our shape optimization algorithm. In order to obtain the Euler–Lagrange equation of the above functional with respect to the unknown function /, the regularized version of the Heaviside distribution H ð/Þ and d-distribution dð/Þ, denoted, respectively, as H e ð/Þ and de ð/Þ ¼ H 0e ð/Þ are employed. Then we get the final energy functional Z Z 1 2 2 ðjr/j 1Þ dX þ k2 Gðu u Þ H e ð/ÞdX Ee ð/Þ ¼ k1 2 X Z ZX þ k3 Gde ð/Þjr/j dX þ k4 GH e ð/ÞdX: ð31Þ X
X
Fig. 2. The domain and the boundary conditions.
Fig. 3. The evolution of the level set functions in the first example.
206
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
Minimizing Ee ð/Þ with respect to /, we can deduce the corresponding Euler–Lagrange equation, which can be used to evolve the level set function / by an artificial time t P 0: 8 h i o/ r/ Gr/ 2 > > < ot ¼ k1 D/ div jr/j þ k2 Gde ð/Þðu u Þ þ k3 de ð/Þdiv jr/j þ k4 Gde ð/Þ; in X; ð32Þ /ð0; xÞ ¼ /0 ðxÞ; in X; > > : o/ ¼ 0; on oX; ox where /ð0; xÞ ¼ /0 ðxÞ is the initial contour and n is the exterior normal to the boundary oX. In our proposed algorithm, the level set function will be evolved by (32). Because the first term of functional (31) automatically maintained / as a signed distance function during the level set evolution, the re-initialization procedure can be completely eliminated. 4. Numerical examples In this section, we will illustrate the efficiency and applicability of the developed shape-topology optimization algorithm by presenting several numerical examples in two dimensions. From the computational point of view, the smeared-out Heaviside and Dirac function are preferred during optimization process. We will use the following regularized Heaviside and Dirac distribution (as shown in Fig. 1) in our numerical examples 8 if x > e; < 1; if x < e; H e ðxÞ ¼ 0; ð33Þ px :1 x 1 ; if jxj 6 e: 1 þ þ sin e p e 2 0; if jxj > e; px de ðxÞ ¼ 1 ð34Þ ; if jxj 6 e: 1 þ cos 2e e
Fig. 4. The evolution of the level set functions in the second example.
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
207
We demonstrate the application of shape optimization algorithm in a Stokes flow using the model as Fig. 2. A fluid flows from left to right in a rectangle region surrounding the domain (obstacle) X0 we want to reshape. No slip and impermeability boundary conditions have been used on all walls except at inlet and outlet. The domain and the boundary conditions are as below. At the top and bottom wall u ¼ 0; v ¼ 0; at the inlet u ¼ 1; v ¼ 0 and at the outlet p ¼ 0. Unless stated otherwise, in all the experimental results shown in this section the following parameters are assumed: k1 ¼ k2 ¼ k3 ¼ 1; k4 ¼ 0. The time step s can be chosen much larger than conventional level set method. But in order to ensure the stability of the evolution, we use s 6 10 in all our examples. We apply a Taylor–Hood element (P 2 P 1 ) discretization of the Stokes problem [23] and the update of the level set function by finite difference method [24]. The present algorithm will be terminated when the difference between two successive objective functional values is less than the prescribed error or when the given maximum number of iterations has been reached. Fig. 3 displays the contours of the level set functions and the velocity fields of the first numerical example by the proposed level set method. The zero-contour of the initial level set function /0 , which is assumed to be a rectangle that encloses the target domain, is shown in the top-left figure. The top-right and the bottom-left figures show the zero-contours after 50 iterations and 100 iterations, respectively, and the bottom-right is the final zero-contour. It can be seen that the final optimal shape is quite satisfiable. As for the second numerical example, the evolution of the level sets of / and the target velocity field which was generated for the obstacle that we want to reshape are shown in Fig. 4. The figure on the top-left shows the initial shape, a rectangle included in the target optimal shape. The figures on the top-right and bottom-left show the shape obtained after 20 iterations and 40 iterations, respectively. It can be seen that we are able to closely identify the shape that generated the target velocity field despite the initial shape and position.
Fig. 5. The evolution of the level set functions in the third example.
208
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
In the third numerical example, the initial shape is two circles. The typical evolution shapes are shown in Fig. 5. Due to the level set implementation, the algorithm allows automatical change of the topology. This example shows the proposed algorithm can handle topology change automatically during the optimization process. It can be concluded from the numerical examples, the presented algorithm is successful in accuracy, convergence speed and insensitivity to initial shape or topology. 5. Conclusions In this paper, we have proposed a variational level set method that completely eliminates the re-initialization procedure. This method has the following advantages over the conventional level set method: the re-initialization procedure can be completely eliminated and a significantly large artificial time can be used for numerically solving the evolutionary PDE, therefore speeding up the evolution procedure. The level set function can be initialized with any distance function with any shape or topology that are more efficient to construct and easier to use in practice. The level set method in our formulation can be easily implemented and is computationally more efficient. Acknowledgement The authors would like to thank the editor and referees for their valuable comments and suggestions that help us to improve this paper. References [1] B. Mohammadi, O. Pironneau, Shape optimization in fluid mechanics, Annual Review of Fluid Mechanics 36 (2004) 255–279. [2] J. Sokolowski, J.P. Zolesio, Introduction to Shape Optimization: Shape Sensitivity Analysis, Springer Series in Computational Mathematics, vol. 10, Springer, Berlin, 1992. [3] G. Allaire, F. Jouve, A.M. Toader, A level-set method for shape optimization, CR Academic Sciences Paris, Series I 334 (2002) 1125– 1130. [4] S. Osher, J.A. Sethian, Level Set Methods and Dynamic Implicit Surface, Springer-Verlag, 2002. [5] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, Journal of Computational Physics 79 (1988) 12–49. [6] M. Sussmann, P. Smereka, S. Osher, A level set approach for computing solution to incompressible two-phase flow, Journal of Computational Physics 114 (1994) 114–146. [7] M.Y. Wang, X.M. Wang, D.M. Guo, A level set method for structural topology optimization, Computer Methods in Applied Mechanics and Engineering 192 (2003) 227–246. [8] F. Santosa, A level-set approach for inverse problems involving obstacles, ESAIM: Control, Optimization, and Calculus of Variation 1 (1996) 17–33. [9] X.B. Duan, Y.C. Ma, R. Zhang, Optimal shape control of fluid flow using variational level set method, Physics Letters A 372 (9) (2008) 1374–1379. [10] X.B. Duan, Y.C. Ma, R. Zhang, Shape-topology optimization for Navier–Stokes problem using variational level set method, Journal of Computational and Applied Mathematics (2007), doi:10.1016/j.cam.2007.11.016. [11] A. Litman, D. Lesselier, F. Santosa, Reconstruction of 2-D binary obstacle by controlled evolution of a level-set, Inverse Problems 14 (1998) 685–706. [12] M. Burge, S. Osher, A survey on level set methods for inverse problems and optimal design, European Journal of Applied Mathematics 16 (2005) 263–301. [13] S. Amstutz, H. Andra¨, A new algorithm for topology optimization using a level-set method, Journal of Computational Physics 216 (2) (2006) 573–588. [14] J.A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science, Cambridge University Press, 1999. [15] T. Chan, L. Vese, Active contours without edges, IEEE Transactions on Image Processing 10 (2001) 266–277. [16] C.M. Li, C.Y. Xu, C.F. Gui, M.D. Fox, Level set evolution without re-initialization: a new variational formulation, in: CVPR’ 05. [17] R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, AMS Chelsea Publishing, 2000. [18] M.E. Gurtin, An Introduction to Continuum Mechanics, Academic Press, San Diego, 1981. [19] V. Girault, P. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, Heidelberg, 1986.
X.-B. Duan et al. / Applied Mathematics and Computation 202 (2008) 200–209
209
[20] J. Gomes, O. Faugeras, Reconciling distance functions and level sets, Journal of Visual Communication and Image Representation 114 (1994) 146–150. [21] V.I. Arnold, Geometrical Methods in the Theory of Ordinary Differential Equations, Springer-Verlag, New York, 1983. [22] M. Kass, A. Witkin, D. Terzopoulos, Snakes: active contour models, International Journal of Computer Vision 1 (1998) 321–331. [23] V. Girault, P. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, Heidelberg, 1986. [24] H.K. Zhao, T. Chan, B. Merriman, S. Osher, A variational level set approach to multiphase motion, Journal of Computational Physics 127 (1996) 179–195.