Shape optimization for the Stokes system with threshold leak boundary conditions

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Shape optimization for the Stokes system with threshold leak boundary conditions © 2024 The Author(s). Published by Elsevier B.V. on behalf of International Association for Mathematics and Computers in Simulation (IMACS). Published version Haslinger, Jaroslav; Mäkinen, Raino A. E. Haslinger, J., & Mäkinen, R. A. E. (2024). Shape optimization for the Stokes system with threshold leak boundary conditions. Mathematics and Computers in Simulation, 221, 180-196. https://doi.org/10.1016/j.matcom.2024.03.002 2024 Mathematics and Computers in Simulation 221 (2024) 180–196 Available online 4 March 2024 0378-4754/© 2024 The Author(s). Published by Elsevier B.V. on behalf of International Association for Mathematics and Computers in Simulation (IMACS). This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Mathematics and Computers in Simulation journal homepage: www.elsevier.com/locate/matcom Original articles Shape optimization for the Stokes system with threshold leak boundary conditions Jaroslav Haslingera, Raino A.E. Mäkinenb,∗ aFaculty of Mechanical Engineering, VŠB-Technical University of Ostrava, 17 listopadu 2172/15, 708 00 Ostrava-Poruba, Czech Republic bFaculty of Information Technology, University of Jyväskylä, P.O. Box 35, FIN-40014, Jyväskylä, Finland ARTICLE INFO Keywords: Shape optimization Stokes problem Threshold leak boundary condition Variational inequality Finite element method ABSTRACT This paper discusses the process of optimizing the shape of systems that are controlled by the Stokes flow with threshold leak boundary conditions. In the theoretical part it focuses on studying the stability of solutions to the state problem in relation to a specific set of domains. In order to facilitate computation, the slip term and impermeability condition are regulated. In the computational part, the optimized portion of the boundary is defined using Bézier polynomials, in order to create a finite dimensional optimization problem. The paper also includes numerical examples to demonstrate the computational efficiency of this approach. 1. Introduction Control and optimization of fluid mechanics models including shape optimization is nowadays well established discipline with many practical applications, see [12,19] and references therein. Typically the behavior of the controlled system is governed by generally nonlinear partial differential equations comprising appropriate boundary conditions. Their solutions are usually smooth functions of control parameters. Some thirty years ago, mathematicians introduced into fluid models the so-called threshold boundary conditions which are well-known in contact mechanics of solids as unilateral and friction conditions. Fujita in his pioneering paper [8] studied two types of such conditions in the Stokes and Navier–Stokes model, namely slip and leak boundary conditions of Tresca type, when slip, leak on the boundary may occur only if the shear, and normal stress, respectively, attains a threshold bound given a-priori. A possible way how to express these conditions is to write them in the form of inclusions involving multivalued mappings which represent the subdifferential of appropriate nonsmooth convex functions. The whole mathematical model then leads to an inequality type problem whose complexity depends partly on the flow model and partly on the choice of the slip/leak law see [2,3,18], e.g. Optimization of systems governed by nonsmooth state relations gives rise to possible nonsmoothness of the whole optimization problem. This fact creates some difficulties from the computational point of view. If we use the original nonsmooth formulation then (to be correct) discretized models should be solved by methods which are tailored just for this type of problems [21]. But their successful application needs some elementary knowledge of tools of nonsmooth analysis. On the other hand, classical gradient type methods when used for solving nonsmooth problems usually fail or give unsatisfactory results. One of ways how to overcome these difficulties is to replace the original state problem by a sequence of smooth ones and to use them as the new state relation in optimization. The resulting problem becomes smooth (provided that the cost function is smooth, too) and so it can be solved by standard methods. Just this way is used in this paper. The present paper deals with a class of 2D shape optimization problems governed by the Stokes equations with threshold leak boundary conditions of Tresca–Navier type which are prescribed on an optimized part of the boundary. It extends the previous ∗Corresponding author. E-mail addresses: [email protected] (J. Haslinger), [email protected] (R.A.E. Mäkinen). https://doi.org/10.1016/j.matcom.2024.03.002 Received 31 May 2023; Received in revised form 20 December 2023; Accepted 2 March 2024 Mathematics and Computers in Simulation 221 (2024) 180–196 181 J. Haslinger and R.A.E. Mäkinen papers [15,16] which are devoted to the Stokes system but with the threshold slip conditions and, in addition it improves several results obtained there. To simplify our presentation we shall consider a very simple geometry of admissible domains. Moreover, the optimized part of the boundary will be represented by the graph of 𝐶1,1functions which will play the role of the design variables. The velocity formulation of the state relation leads to a variational inequality of the 2nd kind using the terminology from [10] due to the presence of the nonsmooth leak term 𝑗. To regularize the problem, 𝑗is replaced by an appropriate sequence of smooth functionals 𝑗𝜀,𝜀→0+. There is yet another troublesome thing from the computational point of view: namely the zero tangential velocity condition 𝑢𝜏= 0 prescribed on the optimized part of the boundary. This condition is realized in computations by a smooth penalty technique. Thus we use simultaneously a penalty and regularization approach for solving the state problem. The paper is organized as follows: in Section 2, the state and shape optimization problems in their original, i.e. nonsmooth form, are defined together with the assumptions guaranteeing the existence of a solution. The most important result needed in the existence analysis is the proof of a stability of solutions with respect to domains. i.e. to show that the solutions to the state problem considered as a function of domains depend continuously (in an appropriate sense) on domain variations. To this end one needs another very important property: to show that any test function used in the weak formulation to the Stokes equations on any admissible domain can be approximated by functions which can be used as test functions on close domains. In [16] this property has been proven for functions satisfying the impermeability condition 𝑣𝜈= 0 on the slip part of the boundary. In the present paper this result is extended to the more general boundary condition of the form 𝒗⋅𝒔= 0 prescribed on the optimized part of the boundary, where 𝒔is a sufficiently smooth vector field depending continuously on the boundary variations. Shape optimization problems with the penalized/regularized state equations are introduced in Section 3. Their solutions now depend on the regularization/penalization parameter 𝜀. It is shown that if 𝜀→0+, they tend on subsequences to a solution of the original nonsmooth optimization problem. Also this convergence result is stronger than these ones in [15,16]. Section 5deals with computational aspects. Optimized part of the boundary with the leak conditions is parametrized by Bézier polynomials, while the regularized-penalized state problem is discretized by stable P1-bubble/P1 elements. The gradient of the cost function is evaluated using the algebraic adjoint state approach. Finally, Section 6presents computational results for two model problems. The paper uses the following notation. If 𝑄is a bounded domain in R𝑛,𝑛= 1,2then 𝐻𝑘(𝑄), 𝑘 ≥0integer, denotes the standard Sobolev space of functions defined in 𝑄which are together with their derivatives up to order 𝑘square integrable in 𝑄. We set 𝐻0(𝑄) = 𝐿2(𝑄). The norm in 𝐻𝑘(𝑄)will be denoted by ‖⋅‖𝑘,𝑄 and the scalar product by (,)𝑘,𝑄. If 𝑋is an ordered vector space then 𝑋+stands for the cone of its non-negative elements. Algebraic vectors and vector functions will be denoted by bold characters. If 𝐚,𝐛are two vectors from 𝑅𝑑, 𝑑 = 1,2,…their scalar product is denoted by 𝐚⋅𝐛. If 𝐀=(𝑎𝑖𝑗 ),𝐁=(𝑏𝑖𝑗)are two 𝑛×𝑛matrices then 𝐀∶𝐁∶= 𝑎𝑖𝑗𝑏𝑖𝑗 (the summation convention is used). The symbol 𝑐stands for a generic positive constant, which may take different values at different places of its occurrence. 2. State problem Let 𝛺 ⊂ R2be a bounded domain with the Lipschitz boundary 𝜕𝛺 =𝛤∪𝛤N∪𝑆, where 𝛤,𝛤N, and 𝑆are non-empty, disjoint parts open in 𝜕𝛺. The classical formulation of the state problem reads as follows: find the velocity vector 𝒖∶𝛺→R2and the pressure 𝑝∶𝛺→Rsuch that ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ −2𝜇div(D𝒖)+∇𝑝=𝒇in 𝛺, div𝒖= 0 in 𝛺, 𝒖=𝟎on 𝛤, 𝝈𝝂 =𝝈𝑁on 𝛤N, 𝑢𝜏= 0 on 𝑆, |𝜎𝜈+𝜅𝑢𝜈|≤𝑔, (𝜎𝜈+𝜅𝑢𝜈)𝑢𝜈+𝑔|𝑢𝜈|= 0 on 𝑆. (2.1) Here 𝜇 > 0is the dynamic viscosity of the fluid, 𝒇∈ (𝐿2(𝛺))2,𝝈𝑁∈ (𝐿2(𝛤N))2,𝑔, 𝜅 ∶𝑆→R+denote an external force, a given value of the stress vector, a non-negative leak threshold, and leak coefficient, respectively. Further D𝒖=1 2(∇𝒖+ (∇𝒖)T)is the symmetric part of the gradient of 𝒖,𝝂,𝝉are the unit normal, and tangential vector, respectively, to 𝜕𝛺. Finally, 𝑣𝜈=𝒗⋅𝝂,𝑣𝜏=𝒗⋅𝝉are the normal, and tangential components of a vector 𝒗∈R2on 𝜕𝛺, respectively, 𝝈= 2𝜇(D𝒖) − 𝑝𝑰is the stress tensor, and 𝜎𝜈=𝝈𝝂⋅𝝂is the normal component of the stress vector 𝝈𝝂 on 𝜕𝛺. From (2.1)6it follows: ∙if 𝑢𝜈(𝑥)=0then |𝜎𝜈(𝑥)|≤𝑔(𝑥), 𝑥 ∈𝑆, ∙if 𝑢𝜈(𝑥)≠0then 𝜎𝜈(𝑥)=−𝜅(𝑥)𝑢𝜈(𝑥) − 𝑔(𝑥)sign𝑢𝜈(𝑥), 𝑥 ∈𝑆.}(2.2) The relation between 𝜎𝜈and −𝑢𝜈is depicted in Fig. 1. Thus a leak at 𝑥∈𝑆occurs only if |𝜎𝜈(𝑥) + 𝜅(𝑥)𝑢𝜈(𝑥)|=𝑔(𝑥). The weak formulation of (2.1) reads: ⎧ ⎪ ⎨ ⎪ ⎩ Find 𝒖∈V(𝛺), 𝑝 ∈𝐿2(𝛺)such that 𝑎(𝒖,𝒗−𝒖) + 𝑏(𝒗−𝒖, 𝑝) + 𝑗(𝑣𝜈, 𝑢𝜈) − 𝑗(𝑢𝜈, 𝑢𝜈)≥𝐿(𝒗−𝒖) ∀𝒗∈V(𝛺) 𝑏(𝒖, 𝑞) = 0 ∀𝑞∈𝐿2(𝛺), () Mathematics and Computers in Simulation 221 (2024) 180–196 182 J. Haslinger and R.A.E. Mäkinen Fig. 1. Relation between normal velocity and normal stress on 𝑆. where V(𝛺)={𝒗∈ (𝐻1(𝛺))2∣𝒗=𝟎on 𝛤, 𝑣𝜏= 0 on 𝑆}, 𝑎(𝒖,𝒗)=2𝜇∫𝛺 D𝒖∶D𝒗𝑑𝑥, 𝒖,𝒗∈ (𝐻1(𝛺))2, 𝑏(𝒗, 𝑞)=−∫𝛺 𝑞div 𝒗𝑑𝑥, 𝒗∈ (𝐻1(𝛺))2, 𝑞 ∈𝐿2(𝛺), 𝐿(𝒗) = ∫𝛺 𝒇⋅𝒗𝑑𝑥 +∫𝛤N 𝝈𝑁⋅𝒗𝑑𝑠, 𝒇∈ (𝐿2(𝛺))2,𝝈𝑁∈ (𝐿2(𝛤N))2,𝒗∈ (𝐻1(𝛺))2, 𝑗(𝑣𝜈, 𝑢𝜈) = ∫𝑆 (𝑔|𝑣𝜈|+𝜅 𝑣𝜈𝑢𝜈)𝑑𝑠, 𝑔, 𝜅 ∈𝐿∞ +(𝑆),𝒖,𝒗∈ (𝐻1(𝛺))2. Problem ()has been studied in [8,9] provided that 𝛤N= ∅ and 𝜅≡0on 𝑆. Fujita proved that 𝒖is unique, whereas 𝑝is determined up to an additive constant which is subject to appropriate constraints arising from the leak conditions (2.1)6. In our case the pressure is unique since the boundary condition of 𝛤Nfixes the value of 𝑝on 𝑆. In what follows we shall suppose that 𝜇=1 2. Concerning the existence and uniqueness of the solution to ()we have the following Theorem 2.1. Problem ()has a unique solution (𝒖, 𝑝)for any 𝒇∈ (𝐿2(𝛺))2,𝝈𝑁∈ (𝐿2(𝛤N))2, and 𝑔, 𝜅 ∈𝐿∞ +(𝑆). In addition, ‖𝒖‖1,𝛺 ≤1 𝑐𝐾‖𝐿‖∗≤1 𝑐𝐾(‖𝒇‖0,𝛺 +𝑐𝑡𝑟‖𝝈𝑁‖0,𝛤N)(2.3) and 𝛽‖𝑝‖0,𝛺 ≤‖𝑎‖‖𝒖‖1,𝛺 +𝑐𝑡𝑟‖𝑔‖∞,𝑆 |length 𝑆|1∕2 +𝑐2 𝑡𝑟‖𝜅‖∞,𝑆 ‖𝒖‖1,𝛺,(2.4) where 𝑐𝐾>0is the constant in Korn’s inequality, 𝑐𝑡𝑟 >0is the norm of the trace mapping 𝑡𝑟 ∶V(𝛺)→(𝐿2(𝜕𝛺))2,‖𝑎‖,‖𝐿‖∗is the norm of 𝑎, and 𝐿, respectively, and 𝛽 > 0is the constant in the inf-sup condition for 𝑏. Proof. The existence and uniqueness of the solution to ()follows from V(𝛺)-ellipticity of the bilinear form 𝑎which is a consequence of Korn’s inequality ∃𝑐𝐾=const. >0 ∶ ∫𝛺 D𝒗∶D𝒗𝑑𝑥 ≥𝑐𝐾‖𝒗‖2 1,𝛺 ∀𝒗∈V(𝛺)(2.5) and the inf-sup condition satisfied by the form 𝑏on V(𝛺) × 𝐿2(𝛺)[17]: ∃𝛽=const. >0 ∶ sup 𝒗∈V(𝛺)⧵{𝟎} 𝑏(𝒗, 𝑞) ‖𝒗‖1,𝛺 ≥𝛽‖𝑞‖0,𝛺 ∀𝑞∈𝐿2(𝛺).(2.6) Inserting 𝒗=𝟎,2𝒖into ()1, we obtain: 𝑎(𝒖,𝒖) + 𝑗(𝑢𝜈, 𝑢𝜈) = 𝐿(𝒖),(2.7) 𝑎(𝒖,𝒗) + 𝑏(𝒗, 𝑝) + 𝑗(𝑢𝜈, 𝑣𝜈)≥𝐿(𝒗) ∀𝒗∈V(𝛺).(2.8) From (2.7) and (2.5) we have: 𝑐𝐾‖𝒖‖2 1,𝛺 ≤𝑎(𝒖,𝒖) + 𝑗(𝑢𝜈, 𝑢𝜈)≤‖𝐿‖∗‖𝒖‖1,𝛺.(2.9) Mathematics and Computers in Simulation 221 (2024) 180–196 183 J. Haslinger and R.A.E. Mäkinen Fig. 2. Decomposition of the boundary of 𝛺(𝛼). It is readily seen that ‖𝐿‖∗≤‖𝒇‖0,𝛺 +𝑐𝑡𝑟‖𝝈𝑁‖0,𝛤N. From this and (2.9) the estimate (2.3) follows. To prove (2.4) we use the inf-sup condition (2.6) and (2.8). It is easy to see that (div 𝒗, 𝑞)0,𝛺 ‖𝒗‖1,𝛺 ≤‖𝑎‖‖𝒖‖1,𝛺 +𝑐𝑡𝑟‖𝑔‖∞,𝑆 |length 𝑆|1∕2 +𝑐2 𝑡𝑟‖𝜅‖∞,𝑆 ‖𝒖‖1,𝛺 holds for any 𝒗∈V(𝛺),𝒗≠𝟎. From this and (2.6) the estimate (2.4) follows. □ 3. Optimal shape design problem: definition and existence analysis The aim of this section is to present and analyze a class of shape optimization problems with the state problem introduced in Section 2. To this end we use the following system of admissible domains: = {𝛺(𝛼) ∣ 𝛼∈𝑎𝑑} where 𝛺(𝛼) = {(𝑥1, 𝑥2) ∣ 𝑥2∈ (0,1), 𝛼(𝑥1)< 𝑥2< 𝛾} and 𝑎𝑑 = {𝛼∈𝐶1,1([0,1]) ∣ −𝛾+𝛥<𝛼min ≤𝛼≤𝛼max < 𝛾 in [0,1], |𝛼(𝑗)|≤𝐶𝑗a.e. in [0,1], 𝑗 = 1,2}.(3.1) Here 𝛼min is a real constant and 𝛼max, 𝛾, 𝐶1, 𝐶2, 𝛥 are positive constants such that 𝑎𝑑 ≠∅. By  𝛺= (0,1) × (−𝛾, 𝛾)we denote the hold-all domain, i.e. 𝛺(𝛼)⊆ 𝛺∀𝛼∈𝑎𝑑. The boundary of any 𝛺(𝛼)will be decomposed as follows: 𝜕𝛺(𝛼) = 𝛤∪𝛤N(𝛼) ∪ 𝑆(𝛼), where 𝛤= (0,1) × {1}, 𝑆(𝛼) = graph of 𝛼, 𝛤N(𝛼) = 𝜕𝛺(𝛼)⧵(𝛤∪𝑆(𝛼)) (see Fig. 2). Since the state problem on 𝛺(𝛼)will be defined for variable 𝛼∈𝑎𝑑 we shall suppose that 𝒇∈ (𝐿2( 𝛺))2and 𝝈𝑁∈ (𝐿2( 𝛤N))2, where  𝛤Nis the union of the vertical sides of  𝛺. To define the functions 𝑔, 𝜅 appearing in the leak term 𝑗on any 𝑆(𝛼), 𝛼 ∈𝑎𝑑 we use functions 𝑔, 𝜅 ∈𝐿∞ +((0,1)) and set 𝑔(𝑥1, 𝑥2) = 𝑔(𝑥1), 𝜅(𝑥1, 𝑥2) = 𝜅(𝑥1) ∀(𝑥1, 𝑥2) ∈  𝛺. Hence ‖𝑔‖∞,𝑆(𝛼)=‖𝑔‖∞,(0,1),‖𝜅‖∞,𝑆(𝛼)=‖𝜅‖∞,(0,1) ∀𝛼∈𝑎𝑑 .(3.2) On any 𝛺(𝛼), 𝛼 ∈𝑎𝑑 we consider the following state problem: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Find (𝒖(𝛼), 𝑝(𝛼)) ∈ V(𝛺(𝛼)) × 𝐿2(𝛺(𝛼)) such that 𝑎𝛼(𝒖(𝛼),𝒗−𝒖(𝛼)) + 𝑏𝛼(𝒗−𝒖(𝛼), 𝑝(𝛼)) + 𝑗𝛼(𝒗⋅𝝂𝛼,𝒖(𝛼)⋅𝝂𝛼) −𝑗𝛼(𝒖(𝛼)⋅𝝂𝛼,𝒖(𝛼)⋅𝝂𝛼)≥𝐿𝛼(𝒗−𝒖(𝛼)) ∀𝒗∈V(𝛺(𝛼)) 𝑏𝛼(𝒖(𝛼), 𝑞) = 0 ∀𝑞∈𝐿2(𝛺(𝛼)), ((𝛼)) Mathematics and Computers in Simulation 221 (2024) 180–196 184 J. Haslinger and R.A.E. Mäkinen where V(𝛺(𝛼)) is the space V(𝛺)defined in Section 2with 𝛺∶= 𝛺(𝛼),𝑆∶= 𝑆(𝛼), and 𝛤N∶= 𝛤N(𝛼). To point out that the forms 𝑎, 𝑏, 𝐿 and the leak term 𝑗depend on 𝛼∈𝑎𝑑, we use notation 𝑎𝛼, 𝑏𝛼, 𝐿𝛼, and 𝑗𝛼, respectively in what follows. The same convention holds for the vectors 𝝂𝛼,𝝉𝛼. Finally, let 𝐽∶𝑎𝑑 × (𝐻1( 𝛺))2×𝐿2( 𝛺)→Rbe a cost functional. The optimal shape design problem we shall study reads as follows: {Find 𝛼∗∈𝑎𝑑 such that 𝐽(𝛼∗,𝒖(𝛼∗), 𝑝(𝛼∗)) ≤𝐽(𝛼, 𝒖(𝛼), 𝑝(𝛼)) ∀𝛼∈𝑎𝑑 ,(P) where (𝒖(𝛼), 𝑝(𝛼)) is the solution to ((𝛼)). Our aim is to show that under appropriate assumptions on 𝐽, problem (P)has at least one solution. We start with Lemma 3.1. Solutions to ((𝛼)) are uniformly bounded with respect to 𝛼∈𝑎𝑑: there exists a positive constant 𝑐, which does not depend on 𝛼∈𝑎𝑑 such that ‖𝒖(𝛼)‖1,𝛺(𝛼)+‖𝑝(𝛼)‖0,𝛺(𝛼)≤𝑐∀𝛼∈𝑎𝑑.(3.3) Proof. From (2.3) and the assumptions on 𝒇and 𝝈𝑁it follows: ‖𝒖(𝛼)‖1,𝛺(𝛼)≤1 𝑐𝐾(‖𝒇‖0, 𝛺+𝑐𝑡𝑟𝛼‖𝝈𝑁‖0, 𝛤N).(3.4) The constant 𝑐𝐾of Korn’s inequality can be chosen to be independent of 𝛼∈𝑎𝑑 (see [20]). Further 𝑡𝑟𝛼stands for the norm of the trace mapping 𝑡𝑟𝛼∶V(𝛺(𝛼)) →(𝐿2(𝜕𝛺(𝛼)))2. It is readily seen that 𝑡𝑟𝛼can be chosen to be independent of 𝛼∈𝑎𝑑 [14]. From this and (3.4) uniform boundedness of ‖𝒖(𝛼)‖1,𝛺(𝛼)follows. The prove the same for ‖𝑝(𝛼)‖0,𝛺(𝛼)we use (2.4): ‖𝑝(𝛼)‖0,𝛺(𝛼)≤1 𝛽(‖𝑎𝛼‖‖𝒖(𝛼)‖1,𝛺(𝛼)+𝑐𝑡𝑟𝛼‖𝑔‖∞,𝑆(𝛼)|lenght 𝑆(𝛼)|1∕2 +𝑐2 𝑡𝑟𝛼‖𝜅‖∞,𝑆(𝛼)‖𝒖(𝛼)‖1,𝛺(𝛼)),(3.5) The constant 𝛽 > 0of the inf-sup condition can be chosen again to be independent of 𝛼∈𝑎𝑑 [4] and the same holds for ‖𝑎𝛼‖. Finally |length𝑆(𝛼)|1∕2 ≤√1 + 𝐶2 1∀𝛼∈𝑎𝑑 as follows from the definition of 𝑎𝑑. Taking into account all these facts together with (3.2),(3.4),(3.5), we obtain uniform boundedness of ‖𝑝(𝛼)‖0,𝛺(𝛼)with respect to 𝛼∈𝑎𝑑 .□ The solution (𝒖(𝛼), 𝑝(𝛼)) to ((𝛼)),𝛼∈𝑎𝑑 will be extended from 𝛺(𝛼)on the hold-all domain  𝛺and denoted as ( 𝒖(𝛼), 𝑝(𝛼)) ∈ (𝐻1( 𝛺))2×𝐿2( 𝛺)in what follows. We can use any extension mapping which preserves the uniform boundedness property of ( 𝒖(𝛼), 𝑝(𝛼)) with respect to 𝛼∈𝑎𝑑 : ∃𝑐 > 0 ∶ ‖ 𝒖(𝛼)‖1, 𝛺+‖𝑝(𝛼)‖0, 𝛺≤𝑐(‖𝒖(𝛼)‖1,𝛺(𝛼)+‖𝑝(𝛼)‖0,𝛺(𝛼))(3.3) ≤𝑐∀𝛼∈𝑎𝑑,(3.6) where 𝑐is a positive constant which does not depend on 𝛼∈𝑎𝑑 . For the pressure 𝑝(𝛼) ∈ 𝐿2(𝛺(𝛼)) we simply use the extension by zero on  𝛺⧵𝛺(𝛼). The extension of 𝒖(𝛼) ∈ (𝐻1(𝛺(𝛼)))2is more involved. One can use either a general result from [5] on the uniform extension property of domains satisfying the uniform cone property or to construct himself an extension mapping taking advantage of a simple shape of 𝛺(𝛼) ∈ and the condition |𝛼′|≤𝐶1 in [0,1]. The key role in the existence analysis plays Theorem 3.1. For any sequence {(𝛼𝑛,𝒖𝑛, 𝑝𝑛)}, where 𝛼𝑛∈𝑎𝑑 and (𝒖𝑛, 𝑝𝑛) ∶= (𝒖(𝛼𝑛), 𝑝(𝛼𝑛)) ∈ V(𝛺(𝛼𝑛)) × 𝐿2(𝛺(𝛼𝑛)) solves ((𝛼𝑛)), 𝑛→∞there exist: its subsequence (denoted by the same symbol) and functions 𝛼∈𝑎𝑑 ,(𝒖, 𝑝)∈(𝐻1( 𝛺))2×𝐿2( 𝛺)such that ⎧ ⎪ ⎨ ⎪ ⎩ 𝛼𝑛→𝛼in 𝐶1([0,1]),  𝒖𝑛⇀𝒖(weakly) in (𝐻1( 𝛺))2, 𝑝𝑛⇀𝑝in 𝐿2( 𝛺), 𝑛 →∞. (3.7) In addition, (𝒖, 𝑝)|𝛺(𝛼)= (𝒖(𝛼), 𝑝(𝛼)) solves ((𝛼)). Proof. The existence of a subsequence satisfying (3.7) results from compactness of 𝑎𝑑 in 𝐶1([0,1]) and (3.6). To prove that (𝒖, 𝑝)|𝛺(𝛼) solves ((𝛼)) we first verify that 𝒖|𝛺(𝛼)∈V(𝛺(𝛼)). To this end it is sufficient to show that 𝒖|𝑆(𝛼)⋅𝝉𝛼= 0 on 𝑆(𝛼).(3.8) From Lemma 2.21 in [14] we know that  𝒖𝑛◦𝛼𝑛∶=  𝒖𝑛(𝑥1, 𝛼𝑛(𝑥1)) →𝒖◦𝛼in (𝐿2((0,1)))2, 𝑛 →∞ and also 0 =  𝒖𝑛◦𝛼𝑛⋅𝝉𝛼𝑛◦𝛼𝑛→𝒖◦𝛼⋅𝝉𝛼◦𝛼in 𝐿2((0,1)) making use of (3.7)1. Thus (3.8) holds and so 𝒖|𝛺(𝛼)∈V(𝛺(𝛼)). Mathematics and Computers in Simulation 221 (2024) 180–196 185 J. Haslinger and R.A.E. Mäkinen Let 𝜒𝑛, 𝜒𝛼be the characteristic functions of 𝛺(𝛼𝑛), and 𝛺(𝛼), respectively and 𝜒𝑛, 𝜒𝛼∈𝐿2( 𝛺)their extensions by zero on  𝛺. From (3.7)1it easily follows that 𝜒𝑛→𝜒𝛼in 𝐿2( 𝛺), 𝑛 →∞.(3.9) The definition of ((𝛼𝑛)),(3.7)2and (3.9) yield: 0 = 𝑏𝑛(𝒖𝑛, 𝑞) = 𝑏 𝛺( 𝒖𝑛, 𝜒𝑛𝑞)→𝑏 𝛺(𝒖, 𝜒𝛼𝑞) = 𝑏𝛼(𝒖, 𝑞) ∀𝑞∈𝐿2( 𝛺),(3.10) where for brevity of notation 𝑏𝑛∶= 𝑏𝛼𝑛and similarly for other forms in the sequel. Hence div 𝒖|𝛺(𝛼)= 0. To accomplish the proof it remains to show that the couple (𝒖, 𝑝)|𝛺(𝛼)satisfies the inequality in ((𝛼)). Let 𝒗∈V(𝛺(𝛼)) be given. Then accordingly to Theorem A.1 and Remark A.2 from Appendix there exist: a sequence {𝒗𝑘}, 𝒗𝑘∈ (𝐻1( 𝛺))2and a function 𝒗∈ (𝐻1( 𝛺))2such that 𝒗|𝛺(𝛼)=𝒗and 𝒗𝑘→𝒗in (𝐻1( 𝛺))2, 𝑘 →∞.(3.11) Moreover, for any 𝑘∈Nthere exists 𝑛𝑘∈Nsuch that 𝒗𝑘|𝛺(𝛼𝑛𝑘)∈V(𝛺(𝛼𝑛𝑘)) (3.12) and consequently 𝒗𝑘|𝛺(𝛼𝑛𝑘)can be used as a test function in ((𝛼𝑛𝑘)): 𝑎𝑛𝑘(𝒖𝑛𝑘,𝒗𝑘−𝒖𝑛𝑘) + 𝑏𝑛𝑘(𝒗𝑘−𝒖𝑛𝑘, 𝑝𝑛𝑘) + 𝑗𝑛𝑘(𝒗𝑘⋅𝝂𝑛𝑘,𝒖𝑛𝑘⋅𝝂𝑛𝑘) − 𝑗𝑛𝑘(𝒖𝑛𝑘⋅𝝂𝑛𝑘,𝒖𝑛𝑘⋅𝝂𝑛𝑘)≥𝐿𝑛𝑘(𝒗𝑘−𝒖𝑛𝑘)(3.13) holds for any 𝑘∈N, where 𝑎𝑛𝑘∶= 𝑎𝛼𝑛𝑘,𝝂𝑛𝑘∶= 𝝂𝛼𝑛𝑘, etc. Next, we pass to the limit with 𝑘→∞in (3.13). From (3.7)2,(3.9) and (3.11) we obtain: lim sup 𝑘→∞ 𝑎𝑛𝑘(𝒖𝑛𝑘,𝒗𝑘−𝒖𝑛𝑘) = lim sup 𝑘→∞∫ 𝛺 𝜒𝑛𝑘D 𝒖𝑛𝑘∶D(𝒗𝑘− 𝒖𝑛𝑘)𝑑𝑥 ≤∫ 𝛺 𝜒𝛼D𝒖∶D(𝒗−𝒖)𝑑𝑥 =𝑎𝛼(𝒖,𝒗−𝒖)(3.14) using weak lower semicontinuity of 𝑎 𝛺and the fact that 𝒗|𝛺(𝛼)=𝒗. Similarly lim 𝑘→∞𝑏𝑛𝑘(𝒗𝑘−𝒖𝑛𝑘, 𝑝𝑛𝑘) = lim 𝑘→∞𝑏𝑛𝑘(𝒗𝑘, 𝑝𝑛𝑘) = 𝑏𝛼(𝒗−𝒖, 𝑝)(3.15) as follows from (3.10) and lim 𝑘→∞𝐿𝑛𝑘(𝒗𝑘−𝒖𝑛𝑘) = 𝐿𝛼(𝒗−𝒖).(3.16) Finally, the leak term: 𝑗𝑛𝑘(𝒗𝑘⋅𝝂𝑛𝑘,𝒖𝑛𝑘⋅𝝂𝑛𝑘) = ∫𝑆(𝛼𝑛𝑘) 𝜅(𝒗𝑘⋅𝝂𝑛𝑘)(𝒖𝑛𝑘⋅𝝂𝑛𝑘)𝑑𝑠 +∫𝑆(𝛼𝑛𝑘) 𝑔|𝒗𝑘⋅𝝂𝑛𝑘|𝑑𝑠 =∫1 0 𝜅(𝒗𝑘⋅𝝂𝑛𝑘)◦𝛼𝑛𝑘(𝒖𝑛𝑘⋅𝝂𝑛𝑘)◦𝛼𝑛𝑘√1+(𝛼′ 𝑛𝑘)2𝑑𝑥1 +∫1 0 𝑔|𝒗𝑘⋅𝝂𝑛𝑘|◦𝛼𝑛𝑘√1+(𝛼′ 𝑛𝑘)2𝑑𝑥1→𝑗𝛼(𝒗⋅𝝂𝛼,𝒖⋅𝝂𝛼)(3.17) using (3.7)1,2,(3.11) and convergence of {𝝂𝑛𝑘◦𝛼𝑛𝑘}to 𝝂𝛼◦𝛼in (𝐿2([0,1]))2. Similarly for the second leak term. From (3.13)–(3.17) we arrive at the assertion of the theorem. □ Remark 3.1. Besides (3.7)2one can prove strong convergence of 𝒖𝑛to 𝒖(𝛼)in the 𝐻1 𝑙𝑜𝑐 (𝛺(𝛼))-norm. Indeed, from (2.7) it follows that ‖𝜒𝑛D 𝒖𝑛∶D 𝒖𝑛‖0, 𝛺→‖𝜒𝛼D𝒖∶D𝒖‖0, 𝛺. From this, (3.7)2, and the fact that we already know that 𝒖|𝛺(𝛼)=𝒖(𝛼), we have ‖𝒖(𝛼) − 𝒖𝑛‖1,𝐷 →0, 𝑛 →∞,(3.18) that holds for any subdomain 𝐷 ⊂ 𝛺(𝛼)such that dist(𝐷, 𝑆(𝛼)) >0. The same result has been proven in [16, Remark 3] for the Stokes system with the threshold slip boundary condition of Tresca type. Next we show that the sequence {𝑝𝑛}tends strongly to 𝑝(𝛼)in the 𝐿2 𝑙𝑜𝑐 (𝛺(𝛼))-norm. To this end we introduce the set 𝐺𝛿(𝛼) = 𝛺(𝛼)⧵BL𝛿(𝛼), where BL𝛿(𝛼) = {(𝑥1, 𝑥2) ∈ 𝛺(𝛼) ∣ 𝑥1∈ (0,1), 𝛼(𝑥1)< 𝑥2< 𝛼(𝑥2) + 𝛿} Mathematics and Computers in Simulation 221 (2024) 180–196 186 J. Haslinger and R.A.E. Mäkinen is the boundary layer along 𝑆(𝛼)and 𝛿 > 0is an arbitrary but sufficiently small. Let such 𝛿 > 0be fixed. On any 𝐺𝛿(𝛼)we consider the spaces V0(𝐺𝛿(𝛼)) = {𝒗∈ (𝐻1(𝐺𝛿(𝛼)))2∣𝒗=𝟎on 𝛤∪𝑆𝛿(𝛼), 𝑆𝛿(𝛼) = 𝑆(𝛼) + 𝛿} and  V0(𝐺𝛿(𝛼)) = {  𝒗∈ (𝐻1( 𝛺))2∣ 𝒗|𝐺𝛿(𝛼)=𝒗∈V0(𝐺𝛿(𝛼)), 𝒗=𝟎in  𝛺⧵𝐺𝛿(𝛼)}. Owing to (3.7)1there exists 𝑛1∶= 𝑛1(𝛿)such that ‖𝛼𝑛−𝛼‖𝐶([0,1]) < 𝛿∕2 ∀𝑛≥𝑛1. Hence any function from  V0(𝐺𝛿(𝛼)) can be used as a test function in ((𝛼)), and ((𝛼𝑛)),𝑛≥𝑛1. The inequality (2.8) corresponding to ((𝛼𝑛)) with test functions  𝒗∈ V0(𝐺𝛿(𝛼)) changes into the equation 𝑎 𝛺( 𝒖𝑛, 𝒗) + 𝑏 𝛺( 𝒗, 𝑝𝑛) = 𝐿 𝛺( 𝒗) 𝒗∈ V0(𝐺𝛿(𝛼)).(3.19) Since 𝒗=𝟎on 𝑆(𝛼𝑛),𝑛≥𝑛1, the leak term 𝑗𝛼𝑛disappears. From the definition of  V0(𝐺𝛿(𝛼)) we see that (3.19) is equivalent to 𝑎𝐺𝛿(𝛼)(𝒖𝑛,𝒗) + 𝑏𝐺𝛿(𝛼)(𝒗, 𝑝𝑛) = 𝐿𝐺𝛿(𝛼)(𝒗) ∀𝒗∈V0(𝐺𝛿(𝛼)), 𝑛 ≥𝑛1. The same holds for the solution (𝒖(𝛼), 𝑝(𝛼)) to ((𝛼)): 𝑎𝐺𝛿(𝛼)(𝒖(𝛼),𝒗) + 𝑏𝐺𝛿(𝛼)(𝒗, 𝑝(𝛼)) = 𝐿𝐺𝛿(𝛼)(𝒗) ∀𝒗∈V0(𝐺𝛿(𝛼)). Subtracting the second equation from the first one we obtain ∫𝐺𝛿(𝛼) div 𝒗(𝑝𝑛−𝑝(𝛼)) 𝑑𝑥 =𝑎𝐺𝛿(𝛼)(𝒖𝑛−𝒖(𝛼),𝒗) ∀𝒗∈V0(𝐺𝛿(𝛼)). Finally from this and the inf-sup condition we obtain 𝛽‖𝑝𝑛−𝑝(𝛼)‖0,𝐺𝛿(𝛼)≤𝑐‖𝒖𝑛−𝒖(𝛼)‖1,𝐺𝛿(𝛼) (3.18) ⟶0, where 𝑐=const. >0which does not depend on 𝑛.□ To guarantee the existence of a minimizer of 𝐽in the optimal shape design problem (P), we shall suppose that 𝐽is lower semicontinuous in the following sense: for any sequence {(𝛼𝑛,𝒚𝑛, 𝑧𝑛)},𝛼𝑛∈𝑎𝑑 ,𝒚𝑛∈ (𝐻1( 𝛺))2and 𝑧𝑛∈𝐿2( 𝛺)such that 𝛼𝑛→𝛼in 𝐶1([0,1]), 𝒚𝑛⇀𝒚in (𝐻1( 𝛺))2, 𝑧𝑛⇀𝑧in 𝐿2( 𝛺), 𝑛 →∞, it holds that lim inf 𝑛→∞𝐽(𝛼𝑛,𝒚𝑛|𝛺(𝛼𝑛), 𝑧𝑛|𝛺(𝛼𝑛))≥𝐽(𝛼, 𝒚|𝛺(𝛼), 𝑧|𝛺(𝛼)).(3.20) Theorem 3.2. Problem (P)has a solution. Proof. The result follows from (3.20) using compactness arguments stated in Theorem 3.1.□ Remark 3.2. In the next computational section we shall use two cost functionals: 𝐽1(𝛼) = 1 2∫1 0 (𝑢𝜈(𝛼)◦𝛼−𝑢)2𝑑𝑥1, 𝑢 ∈𝐿∞((0,1)) given, and 𝐽2(𝛼) = 1 2𝑎𝛼(𝒖(𝛼),𝒖(𝛼)), where 𝒖(𝛼)is the velocity component of the solution to ((𝛼)). It is easy to see that on the basis of Theorem 3.1 both cost functionals satisfy (3.20) (𝐽1is in fact even continuous). 4. Shape optimization with penalized/regularized state problem Problem (P)studied in the previous section posseses two inconveniences from the computational point of view. First of all, the problem is nonsmooth since the state relation is represented by the variational inequality ((𝛼)). This fact restricts the use of numerical minimization methods. Secondly, the tangential no-slip condition 𝑢𝜏= 0 is prescribed on the designed part 𝑆(𝛼). To avoid these drawbacks we use a penalization to release this condition on 𝑆(𝛼)and a regularization of the nonsmooth leak term. To simplify the presentation, the smooth part of the leak functional 𝑗∶ (𝐻1(𝛺(𝛼)))2× (𝐻1(𝛺(𝛼)))2→R+defined in Section 2will be added to the bilinear form 𝑎𝛼, 𝛼 ∈𝑎𝑑 and the new form will be denoted as 𝑎𝜅,𝛼(𝒖,𝒗) ∶= 𝑎𝛼(𝒖,𝒗)+(𝜅𝒖⋅𝝂𝛼,𝒗⋅𝝂𝛼)0,𝑆(𝛼),𝒖,𝒗∈ (𝐻1(𝛺(𝛼)))2(4.1) Mathematics and Computers in Simulation 221 (2024) 180–196 187 J. Haslinger and R.A.E. Mäkinen and set 𝑗𝛼(𝒗⋅𝝂𝛼) ∶= ∫𝑆(𝛼) 𝑔|𝒗⋅𝝂𝛼|𝑑𝑠, 𝒗∈ (𝐻1(𝛺(𝛼)))2.(4.2) We use the simplest penalty functional 𝑡𝛼 𝜀(𝒗⋅𝝉𝛼) = 1 2𝜀‖𝒗⋅𝝉𝛼‖2 0,𝑆(𝛼),𝒗∈ (𝐻1(𝛺(𝛼)))2.(4.3) Regularization of 𝑗𝛼defined by (4.2) consists in its approximation by an appropriate sequence {𝑗𝜀 𝛼}, 𝜀 →0+ of smooth functionals 𝑗𝜀 𝛼. We do not specify their particular choice at the moment, only summarize their properties which will be needed in what follows: ∙𝑗𝜀 𝛼∶𝐿2(𝑆(𝛼)) →R+are convex, 𝐶2-functionals ∀𝜀 > 0, 𝛼 ∈𝑎𝑑,(4.4) ∙𝛼𝑛→𝛼in 𝐶1([0,1]), 𝛼𝑛, 𝛼 ∈𝑎𝑑 𝒗𝑛⇀𝒗in (𝐻1( 𝛺))2, 𝑛 →∞}⟹𝑗𝜀𝑛 𝛼𝑛(𝒗𝑛⋅𝝂𝛼𝑛)⟶ 𝜀𝑛→0+ 𝑗𝛼(𝒗⋅𝝂𝛼),(4.5) ∙ ∃𝑐0>0 ∃𝜀0>0 ∶ 𝑗𝜀 𝛼(0) ≤𝑐0∀𝜀∈ [0, 𝜀0], 𝛼 ∈𝑎𝑑.(4.6) On any 𝛺(𝛼), 𝛼 ∈𝑎𝑑 and 𝜀 > 0we define the following penalized/regularized state problem: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Find (𝒖𝜀(𝛼), 𝑝𝜀(𝛼)) ∈  V(𝛺(𝛼)) × 𝐿2(𝛺(𝛼)) such that 𝑎𝜅,𝛼(𝒖𝜀(𝛼),𝒗) + 𝑏𝛼(𝒗, 𝑝𝜀(𝛼)) + (∇𝑗𝜀 𝛼(𝒖𝜀(𝛼)⋅𝝂𝛼),𝒗⋅𝝂𝛼)0,𝑆(𝛼) +1 𝜀(𝒖𝜀(𝛼)⋅𝝉𝛼,𝒗⋅𝝉𝛼)0,𝑆(𝛼)=𝐿𝛼(𝒗) ∀𝒗∈ V(𝛺(𝛼)) 𝑏(𝒖𝜀(𝛼), 𝑞) = 0 ∀𝑞∈𝐿2(𝛺(𝛼)), (𝜀(𝛼)) where  V(𝛺(𝛼)) = {𝒗∈ (𝐻1(𝛺(𝛼)))2∣𝒗=𝟎on 𝛤}. From (4.4) and (4.5) it follows that 𝜀(𝛼)has a unique solution for any 𝜀 > 0and 𝛼∈𝑎𝑑 .1 We define the new shape optimization problem in which the state equation (𝜀(𝛼)) instead of ((𝛼)) is used: given 𝜀 > 0, {Find 𝛼∗ 𝜀∈𝑎𝑑 such that 𝐽(𝛼∗ 𝜀,𝒖𝜀(𝛼∗ 𝜀), 𝑝𝜀(𝛼∗ 𝜀)) ≤𝐽(𝛼, 𝒖𝜀(𝛼), 𝑝𝜀(𝛼)) ∀𝛼∈𝑎𝑑 ,(P𝜀) where 𝐽is the same cost functional as in (P)and (𝒖𝜀(𝛼), 𝑝𝜀(𝛼)) solves (𝜀(𝛼)). In what follows we shall study if there is a relation between (P)and (P𝜀)as 𝜀→0+. We start with Lemma 4.1. Solutions to (𝜀(𝛼)) are uniformly bounded with respect to 𝜀 > 0and 𝛼∈𝑎𝑑 : ∃𝑐= const. > 0 ∶ ‖𝒖𝜀(𝛼)‖1,𝛺(𝛼)+‖𝑝𝜀(𝛼)‖0,𝛺(𝛼)+1 2𝜀‖𝒖𝜀(𝛼)⋅𝝉𝛼‖2 0,𝑆(𝛼)≤𝑐, (4.7) where 𝑐does not depend on 𝜀 > 0and 𝛼∈𝑎𝑑. Proof. It is well-known that the velocity component 𝒖𝜀(𝛼)solves the following minimization problem: ⎧ ⎪ ⎨ ⎪ ⎩ Find 𝒖𝜀(𝛼) ∈  Vdiv(𝛺(𝛼)) such that 𝜀 𝛼(𝒖𝜀(𝛼)) = min 𝒗∈ Vdiv(𝛺(𝛼)) 𝜀 𝛼(𝒗)(4.8) where 𝜀 𝛼(𝒗) = 1 2𝑎𝜅,𝛼(𝒗,𝒗) + 𝑗𝜀 𝛼(𝒗⋅𝝂𝛼) + 1 2𝜀‖𝒗⋅𝝉𝛼‖2 0,𝑆(𝛼)−𝐿𝛼(𝒗)(4.9) and  Vdiv(𝛺(𝛼)) = {𝒗∈ V(𝛺(𝛼)) ∣ div 𝒗= 0 in 𝛺(𝛼)}.(4.10) From (4.8),(4.9), nonnegativeness of 𝑗𝜀 𝛼,𝑡𝛼, Korn’s inequality and (4.6) we have: 1 2𝑎𝜅,𝛼(𝒖𝜀(𝛼),𝒖𝜀(𝛼)) + 1 2𝜀‖𝒖𝜀(𝛼)⋅𝝉𝛼‖2 0,𝑆(𝛼)≤𝜀 𝛼(𝒖𝜀(𝛼)) + 𝐿𝛼(𝒖𝜀(𝛼)) ≤𝜀 𝛼(𝟎) + (‖𝒇‖0, 𝛺+𝑐𝑡𝑟𝛼‖𝝈𝑁‖0, 𝛤N)‖𝒖𝜀(𝛼)‖1,𝛺(𝛼)≤𝑐, (4.11) where the meaning of  𝛺,𝑐𝑡𝑟𝛼, and  𝛤Nis the same as in Section 3. 1In fact, the existence and uniqueness of the solution can be established under weaker assumptions (see [10]). The stronger assumptions (4.4)–(4.6) are needed because 𝑗𝜀 𝛼depends also on 𝛼∈𝑎𝑑 . Mathematics and Computers in Simulation 221 (2024) 180–196 194 J. Haslinger and R.A.E. Mäkinen Fig. 8. Normal velocity 𝑢𝜈and normal stress 𝜎𝜈. Appendix The aim of this part is to justify (3.11) and (3.12) in Theorem 3.1 which plays the key role in the existence analysis. Roughly speaking, we want to prove that any function from the function space on the limit domain which is used in the weak formulation can be approximated by functions from the same type of spaces on close domains. Before we start, let us recall some notation which will be used in the sequel. If 𝝋,𝝃∶𝑄↦R2,𝑄 ⊆ R𝑑,𝑑= 1,2, are two vector functions, then ‖𝝋‖∶= ‖𝝋‖∞,𝑄 =ess sup 𝑥∈𝑄‖𝝋(𝑥)‖, where ‖𝝋(𝑥)‖denotes the Euclidian norm of 𝝋(𝑥) ∈ R2and 𝝋⋅𝝃∶𝑄→R1,(𝝋⋅𝝃)(𝑥) = 𝝋(𝑥)⋅𝝃(𝑥) ∀𝑥∈𝑄. The system of admissible domains is exactly the same as in Section 3. Unlike domains 𝛺(𝛼) ∈ considered in Section 3the boundaries of which are decomposed into 𝛤, 𝛤𝑁(𝛼), 𝑆(𝛼), boundaries of domains considered in this appendix are split into two parts: 𝜕𝛺(𝛼) = 𝛤(𝛼) ∪ 𝑆(𝛼),𝑆(𝛼) = graph of 𝛼. On any 𝛺(𝛼), 𝛼 ∈𝑎𝑑 we define the space V(𝛼)={𝒗∈ (𝐻1(𝛺(𝛼)))2∣𝒗=𝟎on 𝛤(𝛼),𝒗⋅𝒔𝛼= 0 on 𝑆(𝛼)}, where 𝒔𝛼∶ [0,1] →R2is a given unit vector field defined on 𝑆(𝛼): 𝒔𝛼(𝑥1) ∶= 𝒔𝛼(𝑥1, 𝛼(𝑥1)) ‖𝒔𝛼(𝑥1)‖= 1 }∀𝑥1∈ (0,1) ∀𝛼∈𝑎𝑑. Let 𝒓𝛼∶ [0,1] →R2be another unit vector field on 𝑆(𝛼)which is perpendicular to 𝒔𝛼at any point of 𝑆(𝛼): 𝒔𝛼⋅𝒓𝛼= 0 ‖𝒓𝛼(𝑥1)‖= 1}∀𝑥1∈ (0,1) ∀𝛼∈𝑎𝑑.(A.1) Both these vector fields will be extended from 𝑆(𝛼)on the hold-all domain  𝛺= (0,1) × (−𝛾, 𝛾)as follows: 𝒔𝛼(𝑥1, 𝑥2) = 𝒔𝛼(𝑥1),𝒓𝛼(𝑥1, 𝑥2) = 𝒓𝛼(𝑥1) ∀𝑥= (𝑥1, 𝑥2) ∈  𝛺. (A.2) Convention: from now on the symbols 𝒔𝛼and 𝒓𝛼will denote the vector fields defined by (A.2) in the whole  𝛺. Since the pair (𝒔𝛼,𝒓𝛼)is the orthonormal basis at any 𝑥∈ 𝛺as follows from (A.1), any function 𝒗∶ 𝛺→R2can be written in the form 𝒗=𝒗⋅𝒔𝛼𝒔𝛼+𝒗𝒓𝛼=𝒗⋅𝒔𝛼𝒔𝛼+𝒗⋅𝒓𝛼𝒓𝛼in  𝛺. (A.3) In the next theorem we shall need the following assumptions imposed on 𝒔𝛼and 𝒓𝛼: ∙𝒔𝛼,𝒓𝛼∈ (𝐶0,1( 𝛺))2∀𝛼∈𝑎𝑑,(A.4) ∙ ∃𝐶3=const. >0 ∶ ‖∇𝒔𝛼‖∞, 𝛺≤𝐶3∀𝛼∈𝑎𝑑,(A.5) ∙𝛼𝑛→𝛼in 𝐶1([0,1]), 𝛼𝑛, 𝛼 ∈𝑎𝑑 ⟹𝒔𝛼𝑛→𝒔𝛼in (𝐶( 𝛺))2.(A.6) Mathematics and Computers in Simulation 221 (2024) 180–196 195 J. Haslinger and R.A.E. Mäkinen Theorem A.1. Let (A.4)–(A.6) be satisfied and 𝛼𝑛, 𝛼 ∈𝑎𝑑 be such that 𝛼𝑛→𝛼in 𝐶1([0,1]). Then for any 𝒗∈V(𝛼)there exists a sequence {𝒗𝑘},𝒗𝑘∈ (𝐻1( 𝛺))2and a function 𝒗∈ (𝐻1( 𝛺))2such that 𝒗|𝛺(𝛼)=𝒗and 𝒗𝑘→𝒗in (𝐻1( 𝛺))2, 𝑘 →∞.(A.7) In addition, for any 𝑘∈Nthere exists 𝑛𝑘∈Nsuch that 𝒗𝑘|𝛺(𝛼𝑛𝑘)∈V(𝛼𝑛𝑘).(A.8) Proof. Let 𝒗∈V(𝛼), 𝛼 ∈𝑎𝑑 be fixed and 𝛼𝑛→𝛼in 𝐶1([0,1]), as 𝑛→∞. We denote 𝜑∶= 𝒗⋅𝒔𝛼,𝝍∶= 𝒗𝒓𝛼=𝒗⋅𝒓𝛼𝒓𝛼in 𝛺(𝛼).(A.9) From the definition of V(𝛼)it follows that 𝜑∈𝐻1 0(𝛺(𝛼)) and 𝝍∈ (𝐻1(𝛺(𝛼)))2,𝝍=𝟎on 𝛤(𝛼). Using the density arguments we know that there exist sequences {𝜑𝑘},𝜑𝑘∈𝐶∞ 0(𝛺(𝛼)) and {𝝍𝑘},𝝍𝑘∈ (𝐶∞(𝛺(𝛼)))2such that dist(supp𝝍𝑘, 𝛤(𝛼)) >0for all 𝑘∈Nand {𝜑𝑘→𝜑in 𝐻1 0(𝛺(𝛼)), 𝝍𝑘→𝝍in (𝐻1(𝛺(𝛼)))2.(A.10) Therefore3 {𝜑𝑘→𝜑 in 𝐻1 0( 𝛺),  𝝍𝑘→ 𝝍in (𝐻1( 𝛺))2.(A.11) Moreover we may suppose that dist(supp  𝝍𝑘, 𝛤)>0 ∀𝑘∈N, where  𝛤= {0} × (−𝛾, 𝛾) ∪ {1} × (−𝛾, 𝛾). To construct the sequence {𝒗𝑘}which satisfies (A.7) and (A.8), let us suppose for the moment that for any 𝑘∈Nthere exist: 𝑛𝑘∈N such that ∙𝑆(𝛼𝑛) ∩ supp 𝜑𝑘= ∅ ∀𝑛≥𝑛𝑘(A.12) and a function 𝐍𝑛𝑘∈ (𝐶0,1( 𝛺))2satisfying: ∙𝐍𝑛𝑘|𝑆(𝛼𝑛𝑘)=𝒔𝛼𝑛𝑘|𝑆(𝛼𝑛𝑘),(A.13) ∙𝐍𝑛𝑘 →𝒔𝛼in (𝐻1( 𝛺))2as 𝑘→∞,(A.14) ∙ ∃𝐶4=const. >0 ∶ ‖𝐍𝑛𝑘‖∞, 𝛺+‖∇𝐍𝑛𝑘‖∞, 𝛺≤𝐶4∀𝑘∈N.(A.15) Then the sequence {𝒗𝑘}is defined as follows: 𝒗𝑘=𝜑𝑘𝐍𝑛𝑘+ 𝝍𝑘− 𝝍𝑘⋅𝐍𝑛𝑘𝐍𝑛𝑘.(A.16) Clearly 𝒗𝑘∈ (𝐻1( 𝛺))2,𝒗𝑘=𝟎on 𝛤(𝛼𝑛𝑘)and (𝒗𝑘⋅𝒔𝛼𝑛𝑘)|𝑆(𝛼𝑛𝑘)= ( 𝜑𝑘𝐍𝑛𝑘⋅𝒔𝛼𝑛𝑘)|𝑆(𝛼𝑛𝑘)+ (  𝝍𝑘⋅𝒔𝛼𝑛𝑘)|𝑆(𝛼𝑛𝑘) −(  𝝍𝑘⋅𝐍𝑛𝑘)|𝑆(𝛼𝑛𝑘)(𝐍𝑛𝑘⋅𝒔𝛼𝑛𝑘)|𝑆(𝛼𝑛𝑘) = (  𝝍𝑘⋅𝒔𝛼𝑛𝑘)|𝑆(𝛼𝑛𝑘)− (  𝝍𝑘⋅𝐍𝑛𝑘)|𝑆(𝛼𝑛𝑘)= 0 making use of (A.12) and (A.13). From (A.16),(A.11),(A.14), and (A.15) we see that 𝒗𝑘⟶ 𝑘→∞𝜑𝒔𝛼+ 𝝍− 𝝍⋅𝒔𝛼𝒔𝛼=∶ 𝒗in (𝐻1( 𝛺))2. Finally from (A.3) and (A.9) it follows that 𝒗|𝛺(𝛼)=𝒗. It remains to construct the functions 𝐍𝑛𝑘and the sequence {𝑛𝑘},𝑘→∞satisfying (A.13)–(A.15). Let 𝜉𝑘∈𝐶∞([0,∞)),𝑘→∞be functions such that 0≤𝜉𝑘≤1in [0,∞),𝜉𝑘|[0,1∕(2𝑘)] = 1,𝜉𝑘|[1∕𝑘,∞) = 0 ∀𝑘∈N. For any 𝑘, 𝑛 ∈Nwe define 𝐍𝑛,𝑘(𝑥) = 𝜉𝑘(|𝑥2−𝛼(𝑥1)|)(𝒔𝛼𝑛−𝒔𝛼) + 𝒔𝛼in  𝛺. (A.17) It is readily seen that 𝐍𝑛𝑘∈ (𝐶0,1( 𝛺))2. Further ‖𝐍𝑛,𝑘‖∞, 𝛺≤3 ∀𝑘, 𝑛 ∈N(A.18) and from (A.17) and (A.6) ‖𝐍𝑛,𝑘 −𝒔𝛼‖0, 𝛺≤‖𝒔𝛼𝑛−𝒔𝛼‖0, 𝛺→0, 𝑛 →∞(A.19) 3Let us observe that for functions from 𝐻1 0(𝛺(𝛼)) the symbol ‘‘ ’’ above the function means its extension by zero from 𝛺(𝛼)on  𝛺. Mathematics and Computers in Simulation 221 (2024) 180–196 196 J. Haslinger and R.A.E. Mäkinen uniformly with respect to 𝑘. Let 𝑘∈Nbe fixed. Since 𝛼𝑛→𝛼in 𝐶1([0,1]), there exists an index 𝑛0∶= 𝑛0(𝑘)such that (A.12) holds for any 𝑛≥𝑛0and so 𝐍𝑛,𝑘|𝑆(𝛼𝑛)=𝒔𝛼𝑛|𝑆(𝛼𝑛)∀𝑛≥𝑛0, i.e. (A.13) is satisfied. To estimate ‖∇𝐍𝑛,𝑘‖∞, 𝛺we only need to estimate the term max 𝑥∈ 𝛺(‖‖‖∇𝜉𝑘(|𝑥2−𝛼(𝑥1)|)‖‖‖)‖𝒔𝛼𝑛−𝒔𝛼‖∞, 𝛺‖𝒔𝛼‖∞, 𝛺.(A.20) Since 𝜉′ 𝑘is unbounded on [1∕(2𝑘),1∕𝑘]as 𝑘→∞, one has to compensate this fact by (A.6). Thus for 𝑘fixed, there exists 𝑛1∶= 𝑛1(𝑘) such that the expression (A.20) is bounded by (say) 2 for any 𝑛≥𝑛1. The remaining terms appearing in ∇𝐍𝑛,𝑘 are uniformly bounded due to (A.5). This, together with (A.18) proves (A.15). To verify (A.14) it remains to estimate ‖∇(𝐍𝑛,𝑘 −𝒔𝛼)‖0, 𝛺. From (A.17) and the definition of 𝜉𝑘is follows: ‖∇(𝐍𝑛,𝑘 −𝒔𝛼)‖0, 𝛺≤ max 𝑥∈ 𝛺‖∇𝜉𝑘(|𝑥2−𝛼(𝑥1)|)‖‖𝒔𝛼𝑛−𝒔𝛼‖0, 𝛺+‖∇(𝒔𝛼𝑛−𝒔𝛼)‖0,{|𝑥2−𝛼(𝑥1)|<1∕𝑘} ≤√1 + 𝐶2 1‖𝜉′ 𝑘‖∞,[0,∞)‖𝒔𝛼𝑛−𝒔𝛼‖0, 𝛺+(1∕𝑘), where 𝐶1, 𝐶3are the constants from (3.1) and (A.5). Then we proceed in the same way as in the estimation of (A.20). One can find 𝑛2∶= 𝑛2(𝑘) ∈ Nsuch that ‖∇(𝐍𝑛,𝑘 −𝒔𝛼)‖0, 𝛺=(1∕𝑘)for any 𝑛≥𝑛2. This, together with (A.19) proves (A.14). The function 𝐍𝑛𝑘∶= 𝐍𝑛𝑘,𝑘 having the required properties is defined by (A.17) with 𝑛𝑘= max{𝑛0, 𝑛1, 𝑛2}.□ Remark A.1. It is easy to show that the assertion of Theorem A.1 remains valid also for the space V(𝛺(𝛼)),𝛺(𝛼) ∈ introduced in Section 3when 𝛤𝑁≠∅. Remark A.2. In the previous part of the paper we use Theorem A.1 with 𝒔𝛼∶= 𝝉𝛼, and 𝒓𝛼∶= 𝝂𝛼, where 𝝉𝛼,𝝂𝛼are the unit tangential, and outward normal vectors at points of 𝑆(𝛼),𝛼∈𝑎𝑑 , respectively. Since 𝝉𝛼=(1 √1+(𝛼′)2 ,𝛼′ √1+(𝛼′)2),𝝂𝛼=(𝛼′ √1+(𝛼′)2 ,−1 √1+(𝛼′)2), it is effortless to show that (A.4)–(A.6) are satisfied. This justifies the use of Theorem A.1 to this particular choice of 𝒔𝛼,𝒓𝛼. References [1] D. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337–344. [2] L. Balilescu, J.S. Martín, T. Takahashi, On the Navier-Stokes system with the Coulomb friction law boundary condition, Math. Phys. 68 (3) (2017) 1–25. [3] M. Boukrouche, L. Paoli, Global existence for a 3D non-stationary Stokes flow with Coulomb’s type friction boundary conditions, Appl. Anal. (2017) 1–31. [4] D. Bucur, E. Feireisl, Š. Nečasová, Influence of wall roughness on the slip behavior of viscous fluids, Proc. R. Soc. Edinb., Sect. A, Math. 138 (2008) 957–973. [5] D. Chenais, On the existence of a solution in a domain identification problem, J. Math. Anal. Appl. 52 (1975) 189–219. [6] G. Farin, Curves and Surfaces for CAGD, Fifth Ed., Morgan Kaufmann, 2002. [7] F.N. Fritsch, J. Butland, A method for constructing local monotone piecewise cubic interpolants, SIAM J. Sci. Stat. Comput. 5 (2) (1984) 300–304. [8] H. Fujita, A mathematical analysis of motions of viscous incompressible fluid under leak and or slip boundary conditions, Res. Inst. Math. Sci. Kokyuroku 888 (1994) 199–216. [9] H. Fujita, A coherent analysis of Stokes flows under boundary conditions of friction type, J. Comput. Appl. Math. 149 (2002) 57–69. [10] R. Glowinski, Numerical methods for nonlinear variational problems, Springer Series in Computational Physics, Springer-Verlag, New York, Berlin, Heidelberg, Tokyo, 1984. [11] A. Griewank, A. Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Second ed., Society for Industrial and Applied Mathematics, Philadelphia, PA, 2008. [12] M. Gunzburger, Perspectives in flow control and optimization, in: Advances in design and control, 2002. [13] J. Haslinger, R. Kučera, K. Motyčková, V. Šátek, Numerical modeling of the leak through semipermeable walls for 2D/3D Stokes flow: Experimental scalability of dual algorithms, Mathematics 9 (22) (2021) DOI: 10.3390/math9222906. [14] J. Haslinger, R.A.E. Mäkinen, Introduction to Shape Optimization: Theory, Approximation, and Computation, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2003. [15] J. Haslinger, R.A.E. Mäkinen, J. Stebel, Shape optimization for Stokes problem with threshold slip boundary conditions, Discrete Continuous Dyn. Syst.: Ser. S 10 (6) (2017) 1281–1301. [16] J. Haslinger, J. Stebel, T. Sassi, Shape optimization for Stokes problem with threshold slip, Appl. Math. 59 (6) (2014) 631–652. [17] J.S. Howel, N.J. Walkington, Inf-sup conditions for two-fold Saddle point problems, Numer. Math. 118 (2011) 663–693. [18] C. Le Roux, A. Tani, Steady solutions of the Navier-Stokes equations with threshold slip boundary conditions, Math. Methods Appl. Sci. 30 (2007) 595–624. [19] B. Mohammadi, O. Pironneau, Applied shape optimization for fluids, Second ed., 2009. [20] J.A. Nitsche, On Korn’s second inequality, RAIRO Anal. Numer. 15 (1981) 237–248. [21] J. Outrata, M. Kočvara, J. Zowe, Nonsmooth approach to optimization problems with equilibrium constraints, Nonconvex Optimization and its Applications, vol. 28, Kluwer Academic Publishers, Dordrecht, Boston, London, 1998. [22] The MathWorks Inc, MATLAB version: 9.13.0 (R2022b), The MathWorks Inc, Natick, Massachusetts, USA, 2022, URL https://www.mathworks.com.