Optimal internal stabilization of the linear system of elasticity

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Optimal internal stabilization of the linear system of elasticity A. M¨ unch1, P. Pedregal2, F. Periago3 1Laboratoire de Math´ematiques de Besan¸con, Universit´e de Franche-Comt´e, UMR CNRS 6623, 16, route de Gray 25030 Besan¸con, France. E-mail: [email protected]. 2Dpto. de Matem´aticas, ETSI Industriales, Universidad de Castilla-La Mancha, 13071 Ciudad Real. E-mail: [email protected]. 3Dpto. de Matem´atica Aplicada y Estad´ıstica, ETSI Industriales, Universidad Polit´ecnica de Cartagena, 30203 Cartagena. E-mail: [email protected]. Palabras clave: Relaxation, stabilization, system of elasticity, nonlinear optimal design, gradient descent method. Resumen We consider the nonlinear optimal design problem which consists in finding the best position and shape of the internal viscous damping set for the stabilization of the linear system of elasticity. Since non-existence of classical designs is usual in this context, a relaxation of the original problem is proposed. Then, the relaxed problem is solved numerically. Finally, a penalization technique to recover quasi-optimal classical designs from the relaxed ones and the over-damping phenomenon are analyzed in several numerical experiments. 1. Problem formulation Consider the following nonlinear optimal design problem: (P) ´ınf ω∈ΩL J(Xω) = 1 2ZT 0ZΩ³¯¯u0¯¯ 2+σ(u) : ε(u)´dxdt (1) where for a fixed 0 < L < 1, ΩL={ω⊂Ω : |ω|=L|Ω|} , 1 A. M¨unch, P. Pedregal, F. Periago |ω|and |Ω|being the Lebesgue measure of ωand Ω,respectively, and uis the solution of the elasticity system        u00 − ∇x·σ+a(x)Xω(x)u0= 0 in (0, T)×Ω, u= 0 on (0, T)×Γ0, σ·n= 0 on (0, T)×Γ1, u(0,·) = u0,u0(0,·) = u1in Ω. (2) As usual, ε=ε(u) = 1 2³∇xu+ (∇xu)T´(3) is the linearized strain tensor and σ=σ(u) = (σij =aijklεkl) (4) the stress tensor. The coefficients aijkl ∈W1,∞(Ω), i, j, k, l = 1,· · · , N, are such that aijkl =aklij =ajikl and aijklεijεkl ≥αεijεij in Ω (5) for some fixed α > 0. Moreover, Xωis the characteristic function of ω, ∇x·is the divergence operator considered with respect to the spatial variable x, n= (n1,· · · , nN) is the outward unit normal vector to Γ1,0< T ≤ ∞,and a=a(x)∈L∞(Ω) is a damping potential satisfying a(x)≥a0>0 a. e. x∈ω. As for the physical meaning of system (2), the dissipative term a(x)Xω(x)u0is usually referred in the literature as a viscous damping because it is caused for the viscosity of the medium in which the vibrations of the system take place. From an engineering viewpoint, this term may also be seen as a feedback control mechanism which measures the velocity of vibrations, through use of sensors, and acts on the system accordingly to these measures by means of actuators. In this sense, Xωindicates the place and shape of sensors and actuators. It is then natural and very important in practise to analyze the question of determining the best position and shape of sensors and actuators that minimize the energy of the system over a time interval. This is our problem (P). 2. Relaxation As shown recently by the authors for the case of the wave equation [4, 5] (see also [2, 3] for some related works), (P) is usually ill-posed in the sense that there is no minimizer in the class of characteristic functions. Then, a full relaxation of the original problem was carried out by using a suitable representation of divergence-free vector fields which enables to transform the original problem into a non-convex, vector variational one. The aim of this work is to extend the results in [4, 5] to the case of the system of linear elasticity. However, our approach here does not require the introduction of auxiliary potentials associated with divergence-free vector fields. Instead of those, we use div-curl Young measures as given by [7]. This makes the treatment much more direct and dimensionindependent. Precisely, consider the relaxed problem (RP) ´ınf s∈L∞(Ω) J(s) = ZT 0ZΩ³¯¯u0¯¯ 2+σ(u) : ε(u)´dxdt (6) 2 Optimal internal stabilization of the elasticity system where usolves the new system        u00 − ∇x·σ+a(x)s(x)u0= 0 in (0, T)×Ω, u= 0 on (0, T)×Γ0, σ·n= 0 on (0, T)×Γ1, u(0,·) = u0,u0(0,·) = u1in Ω, (7) and now the competing functions ssatisfy the point-wise and volume constraints 0≤s(x)≤1 and ZΩ s(x)dx =L|Ω|.(8) Our main theoretical result follows. Theorem 1 Assume that the initial data of system (2) have the regularity (u0,u1)∈³¡H2(Ω)¢N∩V0´×V0.(9) Then (RP)is a full relaxation of (P)in the sense: (i) There are optimal solutions for (RP ). (ii) The minimum of (RP )equals the infimum of (P). (iii) Minimizing sequences for (P)are recovered by first-order laminates with any normal in space and independent of time. For a full proof of this result we refer to [6]. Nevertheless, a few comments on it are now in order. First, the regularity condition (9) on the initial data is a sufficient condition in order to avoid concentration of energy phenomena and therefore this enables us to use the Young measure theory to compute the cost limit of a minimizing sequence for problem (P). Second, we will show later on that for some values of the damping potential athere is a numerical evidence that problem (P) is ill-posed. This justifies the relaxation stated in points (i) and (ii) above. Finally, in what concerns point (iii), it tells us how the microstructure of the optimal damping designs looks like. This information is codified by the optimal Young measure associated with the relaxed problem (RP). In fact, from this it can be proved that if Xωjis a minimizing sequence for (P), then the associated displacements ujconverge to the optimal displacement ufor the relaxed problem (RP) in a strong sense. 3. Numerical resolution of problems (P)and (RP ) 3.1. Algorithm of minimization We propose a first-order gradient descent method to solve (RP). Precisely, we define the descent direction s1(x) = −µa(x)ZT 0 u0(t, x)·p(t, x)dt +γ¶,∀x∈Ω,(10) 3 A. M¨unch, P. Pedregal, F. Periago where the multiplier γis determined so that for any function η∈L∞(Ω,R+), with η6= 0, and ||s+ηs1||L1(Ω) =L|Ω|we have Jγ(s+ηs1)≤Jγ(s) for Jγ(s) = J(s) + γ||s||L1(Ω).(11) This leads to take γ=(RΩs(x)dx −L|Ω|)−RΩη(x)a(x)RT 0u0(t, x)·p(t, x)dtdx RΩη(x)dx .(12) As for p, this is the solution of the adjoint problem        p00 − ∇x·σ(p)−a(x)s(x)p0=u00 +∇x·σ(u),in (0, T)×Ω, p= 0,on (0, T)×Γ0, p·n= 0,on (0, T)×Γ1, p(T, ·) = 0,p0(T, ·) = u0(T, ·) in Ω. (13) At last, the function ηis chosen so that s(x) + η(x)s1(x)∈[0,1], for all x∈Ω. A simple and efficient choice consists in taking η(x) = ²s(x)(1 −s(x)) for all x∈Ω with ²a small real positive. Consequently, the descent algorithm to solve numerically the relaxed problem (RP) may be structured as follows : let Ω ⊂RN, (u0,u1)∈((H2(Ω))N∩V0)×V0,L∈(0,1), T > 0, and ² < 1, ²1<< 1 be given : Initialization of the density function s0∈L∞(Ω; ]0,1[); For k≥0, iteration until convergence (i.e. |J(sk+1)−J(sk)| ≤ ²1|J(s0)|) as follows : •Computation of the solution uskof (7) and then the solution pskof (13), both corresponding to s=sk. •Computation of the descent direction sk 1defined by (10) where the multiplier γkis defined by (12). •Update the density function in Ω: sk+1 =sk+²sk(1 −sk)sk 1(14) with ε∈R+small enough in order to ensure the decrease of the cost function and sk+1 ∈L∞(Ω,[0,1]). 3.2. Numerical experiments Next, we present some numerical simulations for N= 2 and the unit square Ω = (0,1)2. For simplicity, we consider the case Γ0=∂Ω and assume that Ω is composed of an isotropic homogeneous material for which aijkl =λδijδkl +µ(δikδjl +δilδjk). λ > 0 and µ > 0 are the Lam´e coefficients and δdesignates the Kronecker symbol. The stress tensor becomes simply: σ(u) = λtr(∇x·u)IN×N+ 2µε(u). 4 Optimal internal stabilization of the elasticity system Systems (7) and (13) are solved in space using a C0-finite element method and the time discretization is performed in a standard way using centered finite differences of order two. Without loss of generality, we consider a constant damping function a(x) = aXΩ(x) in Ω since the dependence in xis contained in the density s. We treat the simple conditions : u0= (sin(πx1) sin(πx2),sin(πx1) sin(πx2)),u1= (0,0).(15) Results are obtained with h= 10−2,²1= 10−5,L= 10−1,T= 1, s0(x) = Lon Ω and ²= 10−2(see the algorithm). 3.2.1. Influence of the damping constant value a Numerical simulations exhibit a bifurcation phenomenon with respect to the value of the damping constant a. When this value is small enough, the optimal density is always a characteristic function which suggests that the original problem (P) is well-posed. On the other hand, for alarge enough it appears that the optimal density takes values strictly in (0,1). In this case, the ill-posedness is related to the over-damping phenomenon : when m´ınωa(x)goes to infinity, the damping term a(x)Xωacts as penalization term and enforces the solution uto be constant in time in ω: at the limit, there is no more dissipation in ω (and so in Ω) and the energy is constant (see also [1]). In order to avoid this phenomenon (which thus appears if as(x) is too large), the density smust take (at least locally) values lower than 1 in order to compensate a. As a result, (P) is no more well-posed. This justifies the introduction of the relaxed problem (RP ). For (λ, µ) = (1/2,1), Figure 1 depicts the iso-values of the optimal density sopt - obtained at the convergence of the algorithm - for several values of a. 3.2.2. From optimal relaxed designs to quasi-optimal classical designs In the case where the optimal density sopt is not in L∞((0, T)×Ω; {0,1}), one may associate to sopt a characteristic function spen ∈L∞((0, T )×Ω; {0,1}) whose cost J(spen) is arbitrarily near to J(sopt). Following [5], one may proceed as follows: we first decompose the domain (0,1)×(0,1) into M×Ncells such that Ω = ∪i=1,M [xi, xi+1]×∪j=1,N [yj, yj+1] where {xi}(i=1,M+1) and {yj}(j=1,N+1) designate two uniform subdivisions of the interval (0,1). Then, we associate to each cell the mean value mi,j ∈[0,1] defined by mi,j =1 (xi+1 −xi)(yj+1 −yj)Zxi+1 xiZyj+1 yj sopt(x, y)dx dy (16) At last, we define the function spen M,N in L∞(Ω,{0,1}) by spen M,N (x, y) = M X i=1 N X j=1 X[xi,(1−√mi,j )xi+√mi,j xi+1]×[yj,(1−√mi,j )yj+√mi,j yj+1](x, y).(17) We easily check that ||spen M,N ||L1(Ω) =||sopt||L1(Ω), for all M, N > 0. Thus, the characteristic function spen M,N takes advantage of the information codified in the density sopt. Let us illustrate this point with the optimal density obtained for (λ, µ) = (1/2,1) and a= 50 (see Figure 1 bottom right). The corresponding value of the cost function is J(sopt)≈2,0883. Table 1 collects the value of J(spen M,N ) for several values of M=Nand suggests the convergence of J(spen M,N ) toward J(sopt)≈2,0883 as Nincreases. 5 A. M¨unch, P. Pedregal, F. Periago 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x1 x2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x1 x2 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x1 x2 0.25 0.2 0.15 0.1 0.05 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x1 x2 Figura 1: T= 1, (λ, µ) = (1/2,1), a(x) = aXΩ(x) - Iso-value of the optimal density sopt on Ω for a= 5 (top left), a= 10 (top right), a= 25 (bottom left) and a= 50 (bottom right). 6 Optimal internal stabilization of the elasticity system 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x1 x2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x1 x2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x1 x2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 8 9 t E(t) N=1 N=2 N=5 N=16 Figura 2: T= 1, (λ, µ) = (1/2,1), a(x) = 50XΩ(x) - Iso-value of the penalized density for N= 2, N= 5 and N= 16 - Bottom right:E(t) vs. tassociated to spen N,N . 7 A. M¨unch, P. Pedregal, F. Periago N1 2 3 4 5 6 7 16 N→ ∞ J(spen N,N ) 3.824 3.261 2.446 2.238 2.161 2.137 2.109 2.096 2.0883 Tabla 1: (λ, µ) = (1/2,1), T= 1 - a= 50- Value of the cost function for the penalized characteristic density. Agradecimientos The second author is supported by projects MTM2004-07114 from Ministerio de Educaci´on y Ciencia (Spain), and PAI05-029 from JCCM (Castilla-La Mancha). The third author is supported by projects MTM2004-07114 from Ministerio de Educaci´on y Ciencia (Spain) and 00675/PI/04 from Fundaci´on S´eneca (Gobierno Regional de Murcia). Referencias [1] C. Castro and S.J. Cox, Achieving arbitrarily large decay in the damped wave equation, SIAM J. Control Optim. 39(6), (2001) 1748-1755 [2] S.J. Cox, Designing for optimal energy absorption II, The damped wave equation, International series of numerical mathematics 126, (1998) 103–109. [3] P. Hebrard and A. Henrot, Optimal shape and position of the actuators for the stabilization of a string, Systems and control letters 48, (2003) 199–209. [4] A. M¨unch, P. Pedregal and F. Periago, A variational approach to a shape design problem for the wave equation, C. R. Acad. Sci. Paris, Ser. I 343 (2006) 371–376. [5] A. M¨unch, P. Pedregal and F. Periago, Optimal design of the damping set for the stabilization of the wave equation, J. Differential Equations 231 (2006) 331–358. [6] A. M¨unch, P. Pedregal and F. Periago, Optimal internal stabilization of the linear system of elasticity. Preprint available at http://matematicas.uclm.es/omeva . [7] P. Pedregal, Div-Curl Young measures and optimal design in any dimension, Revista Mat. Complutense 20 (2007), no1, 239–255. 8