Un método de elementos finitos mixtos para un problema de interacción sólido–fluido

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Un m´etodo de elementos finitos mixtos para un problema de interacci´on s´olido–fluido S. Meddahi1, G.N. Gatica2, A. M´ arquez3 1Dpto. de Matem´aticas, Universidad de Oviedo, Calvo Sotelo s/n, Oviedo. E-mail: [email protected]. 2Dpto. de Ingenier´ıa Matem´atica, Universidad de Concepci´on, Casilla 160-C, Concepci´on, Chile. E-mail: [email protected]. 3Dpto. de Construcci´on e Ingenier´ıa de Fabricaci´on, Universidad Oviedo, Campus de Viesques, Gij´on. E-mail: [email protected]. Palabras clave: elementos finitos mixtos, ecuaci´on de Helmholtz, ecuaci´on elastodin´amica Resumen En este trabajo consideramos un s´olido el´astico lineal e is´otropo, rodeado de un fluido perfecto compresible, sobre el que incide una onda ac´ustica arm´onica. Nuestro prop´osito es presentar un esquema num´erico para determinar tanto la respuesta en el s´olido como la distribuci´on de ondas ac´usticas en el fluido linealizado. En el s´olido utilizamos una formulaci´on variacional mixta de la que, posteriormente, eliminamos el campo de desplazamientos. As´ı, las ´unicas inc´ognitas en el s´olido son los campos de tensiones y de rotaciones. Esta formulaci´on mixta se acopla, mediante dos condiciones de transmisi´on (una de equilibrio y otra de continuidad en desplazamientos) sobre la frontera h´umeda, con la ecuaci´on de Helmholtz que satisface la presi´on sobre el medio ac´ustico. Para definir el correspondiente esquema discreto utilizamos elementos PEERS en el s´olido y elementos finitos de Lagrange de primer orden en el dominio ac´ustico. Finalmente, ilustramos las propiedades de convergencia del esquema propuesto con algunos experimentos num´ericos. 1. Introducci´on En este trabajo presentamos un m´etodo de elementos finitos mixto–primal para resolver un problema plano de interacci´on s´olido–fluido arm´onico en tiempo. Consideramos un s´olido el´astico Ωssobre el que incide una onda ac´ustica. Suponemos que el medio ac´ustico ocupa una regi´on anular Ωfcuya frontera exterior Γ est´a situada lejos del obst´aculo (el s´olido) e imponemos sobre esta curva artificial cerrada una condici´on de contorno que reproduce el comportamiento del campo ac´ustico reflejado en el infinito. De esta manera, 1 S. Meddahi, G.N. Gatica, A. M´arquez nuestro problema modelo queda planteado sobre una regi´on acotada. La mayor parte de los m´etodos num´ericos propuestos para resolver esta clase de problema de interacci´on s´olido– fluido utilizan una formulaci´on en desplazamientos para la elasticidad lineal (ver, por ejemplo, [2, 5, 6, 7] y las referencias citadas all´ı). En este trabajo, en cambio, empleamos una formulaci´on variacional mixta para la elasticidad sobre el s´olido y conservamos la formulaci´on primal habitual en el fluido linealizado. As´ı, el tensor de tensiones en Ωsy la presi´on en el fluido ac´ustico en Ωfser´an nuestras principales inc´ognitas. 2. Planteamiento del problema de interacci´on s´olido–fluido Consideramos un s´olido Ωsinmerso en un fluido ac´ustico sobre el que incide una onda. El contorno de Ωsse denomina frontera h´umeda y se denota mediante Σ. Para definir el dominio ac´ustico introducimos una circunferencia Γ de radio suficientemente grande, centrada en el origen, y denotamos mediante Ωfa la regi´on anular acotada entre Σ y Γ. Suponemos que la onda incidente y las fuerzas de volumen exhiben un comportamiento arm´onico en tiempo con frecuencia ωy amplitudes piyf, respectivamente, de manera que pisatisfaga la ecuaci´on de Helmholtz en Ωf. Suponemos, tambi´en, que el fluido es perfecto, compresible y homog´eneo, con densidad ρfy n´umero de onda κf:= ω v0, siendo v0la velocidad de sonido en el fluido linealizado. Adem´as, se asume un comportamiento el´astico lineal e is´otropo para el s´olido, con densidad ρsy constantes de Lam´e µyλ. Es decir, la ecuaci´on constitutiva del material s´olido, σ=Cε(u) en Ωs,(1) donde ε(u) := 1 2(∇u+ (∇u)t) es el tensor de peque˜nas deformaciones y Ces el tensor de constantes el´asticas, viene dada por la ley de Hooke: Cζ:= λtr(ζ)I+ 2 µζ,(2) siendo Ila matriz identidad y tr(ζ) := P2 i=1 ζii. Las inc´ognitas del problema son la amplitud σ: Ωs→C2×2del tensor de tensiones, la amplitud u: Ωs→C2del campo de desplazamientos y la amplitud de la presi´on global (incidente + reflejada) p: Ωf→C. Bajo la hip´otesis de peque˜nas oscilaciones en el s´olido y en el fluido, y dados f∈[L2(Ωs)]2y g∈H−1/2(Γ), consideramos el siguiente problema de interacci´on s´olido–fluido: Encontrar σ∈H(div; Ωs), u∈[L2(Ωs)]2yp∈H1(Ωf), tales que: σ=Cε(u) en Ωs, div(σ) + κ2 su=−fen Ωs, ∆p+κ2 fp= 0 en Ωf, σν =−pνsobre Σ , ρfω2u·ν=∂p ∂νsobre Σ , ∂p ∂ν−ι κfp=g:= ∂pi ∂ν−ι κfpisobre Γ , (3) donde κs:= √ρsωes el n´umero de onda en el s´olido, νes la normal exterior unitaria sobre Γ y utilizamos el s´ımbolo ıpara √−1. La segunda ecuaci´on de (3) es la ecuaci´on de 2 Un m´etodo de elementos finitos mixtos para un problema de interacci´on s´olido–fluido equilibrio de la elastodin´amica en r´egimen arm´onico, mientras que la tercera es la ecuaci´on de Helmholtz. La primera condici´on de transmisi´on en (3) representa el equilibrio de fuerzas sobre Σ y la segunda la continuidad en desplazamientos normales. 3. Formulaci´on variacional del problema Para tratar el sistema de ecuaciones (3) empleamos una formulaci´on primal en el fluido Ωfy otra mixta en el s´olido Ωs. As´ı, si multiplicamos la ecuaci´on ac´ustica por q∈H1(Ωf), integramos por partes y utilizamos la condici´on de contorno Robin obtenemos ZΩf∇p·∇q−κ2 fZΩf p q +h∂p ∂ν, qiΣ−ı κfZΓ p q =hg, qiΓ,(4) donde, dado S ∈ {Σ,Γ}, representamos mediante h·,·iSel producto de dualidad entre H−1/2(S) y H1/2(S) con respecto al L2(S)-producto escalar. A continuaci´on, utilizamos la condici´on de transmisi´on en desplazamientos y reemplazamos ∂p ∂νpor ρfω2u·νsobre Σ, introducimos la inc´ognita auxiliar ϕ:= u|Σ∈[H1/2(Σ)]2, y dividimos por ρfω2, para reescribir (4) como 1 ρfω2ZΩf∇p·∇q−κ2 f ρfω2ZΩf p q +hqν,ϕiΣ−ıκf ρfω2ZΓ p q =1 ρfω2hg, qiΓ,(5) donde h·,·iΣdenota, en lo que sigue, el producto de dualidad entre [H−1/2(Σ)]2y [H1/2(Σ)]2 respecto al [L2(Σ)]2-producto escalar. Por otra parte, para obtener una formulaci´on variacional mixta en el s´olido Ωs, seguimos la metodolog´ıa habitual (ver[1] y [9]) e introducimos la rotaci´on γ:= 1 2(∇u−(∇u)t)∈[L2(Ωs)]2×2 asym como inc´ognita adicional, donde [L2(Ωs)]2×2 asym denota el espacio de los tensores antisim´etricos con componenetes en L2(Ωs). De esta manera, la ley de comportamiento se puede reescribir como C−1σ=ε(u) = ∇u−γ.(6) A continuaci´on, multiplicamos (6) por τ∈H(div; Ωs) e integramos por partes para obtener ZΩsC−1σ:τ+ZΩs u·div(τ)− hτν,ϕiΣ+ZΩs τ:γ= 0 ,(7) donde σ:τ:= P2 i,j=1 σij τij. La ecuaci´on elastodin´amica nos proporciona la siguiente expresi´on para el campo de desplazamientos: u=−1 κ2 s¡f+div(σ)¢,(8) 3 S. Meddahi, G.N. Gatica, A. M´arquez que introducida en la ecuaci´on constitutiva nos permite escribir (7) como ZΩsC−1σ:τ−1 κ2 sZΩs div(σ)·div(τ)−hτν,ϕiΣ+ZΩs τ:γ=1 κ2 sZΩs f·div(τ).(9) Finalmente, la simetr´ıa de σy la condici´on de transmisi´on en fuerzas sobre Σ se imponen d´ebilmente mediante ZΩs σ:η= 0 ∀η∈[L2(Ωs)]2×2 asym (10) y hpν+σν,ψiΣ= 0 ∀ψ∈[H1/2(Σ)]2.(11) Consecuentemente, si sumamos (5) y (9), y sustraemos (11) de (10), obtenemos la siguiente formulaci´on variacional del problema (3): Encontrar ((σ, p),(ϕ,γ)) ∈H×Qtales que A((σ, p),(τ, q)) + B1((τ, q),(ϕ,γ)) = F(τ, q)∀(τ, q)∈H, B2((σ, p),(ψ,η)) = 0 ∀(ψ,η)∈Q, (12) donde HyQson los espacios producto H:= H(div; Ωs)×H1(Ωf),Q:= [H1/2(Σ)]2×[L2(Ωs)]2×2 asym ,(13) F:H→Ces el funcional lineal F(τ, q) := 1 κ2 sZΩs f·div(τ) + 1 ρfω2hg, qiΓ∀(τ, q)∈H,(14) yA:H×H→C,B1:H×Q→C, y B2:H×Q→Cson las formas bilineales definidas como A((ζ, r),(τ, q)) := ZΩsC−1ζ:τ−1 κ2 sZΩs div(ζ)·div(τ) + 1 ρfω2ZΩf∇r·∇q −κ2 f ρfω2ZΩf rq −ıκf ρfω2ZΓ rq ∀(ζ, r),(τ, q)∈H, (15) B1((τ, q),(ψ,η)) := hqν−τ ν,ψiΣ+ZΩs τ:η∀(τ, q)∈H,∀(ψ,η)∈Q,(16) y B2((τ, q),(ψ,η)) := −hqν+τ ν,ψiΣ+ZΩs τ:η∀(τ, q)∈H,∀(ψ,η)∈Q.(17) La demostraci´on del siguiente teorema se puede encontrar en [4]. Teorema 3.1 El problema (12) posee una ´unica soluci´on. 4 Un m´etodo de elementos finitos mixtos para un problema de interacci´on s´olido–fluido 4. Un m´etodo de elementos finitos mixto–primal En esta secci´on presentamos una aproximaci´on de Galerkin del problema (12). Sean {Th}h>0:= {Ths}hs>0∪{Thf}hf>0, donde {Ths}hs>0y{Thf}hf>0son familias regulares de triangulaciones de las regiones poligonales ¯ Ωsy¯ Ωf, respectivamente, mediante tri´angulos T de di´ametro hTcon tama˜nos de malla hs:= max{hT:T∈ Ths},hf:= max{hT:T∈ Thf}yh:= max{hs, hf}, y tal que los v´ertices de {Ths}hs>0y{Thf}hf>0coincidan sobre Σ. Tambi´en introducimos una partici´on independiente nˆ Σ1,ˆ Σ2,···,ˆ Σmode la frontera h´umeda Σ y denotamos ˆ h:= max{|ˆ Σj|:j∈ {1,··· , m}}. En estas condiciones definimos los subespacios de elementos finitos Hσ h,Hp h,Qϕ ˆ hyQγ hpara las inc´ognitas σ, p,ϕyγde (12), respectivamente, como: Hσ h:= ©τh∈H(div; Ωs) : τh,i|T∈RT0(T)t⊕P0(T)curltbT,∀T∈ Thsª,(18) Hp h:= ©qh∈C(¯ Ωf) : qh|T∈P1(T)∀T∈ Thfª,(19) Qϕ ˆ h:= nψˆ h∈[C(Σ)]2:ψˆ h|ˆ Σj∈[P1(ˆ Σj)]2∀j∈ {1,··· , m}o,(20) Qγ h:= ½µ 0ηh −ηh0¶:ηh∈C(¯ Ωs), ηh|T∈P1(T)∀T∈ Ths¾,(21) donde τh,i es la fila i-´esima de τh(i∈ {1,2}), RT0(T) es el espacio local de RaviartThomas de orden 0 (cf. [3], [8]), bTes la funci´on burbuja c´ubica habitual sobre T∈ Ths, curltbT:= (∂bT ∂x2,−∂bT ∂x1), C(·) es el espacio de las funciones continuas sobre el correspondiente dominio, y, dado un entero `≥0 y un subconjunto Kde R2,P`(K) denota el espacio de los polinomios definidos en Kde grado ≤`. Destacamos que si definimos Qu h:= ©vh∈[L2(Ωs)]2:vh|T∈[P0(T)]2∀T∈ Thsª,(22) el espacio Hσ h×Qu h×Qγ hconstituye el elemento PEERS introducido en [1] para una aproximaci´on con elementos finitos mixtos de la elasticidad lineal plana. Sean Hh:= Hσ h×Hp h,Qˆ h,h := Qϕ ˆ h×Qγ h.(23) Entonces, el esquema de elementos finitos mixto–primal asociado al problema de transmisi´on (12) es: Encontrar ((σh, ph),(ϕˆ h,γh)) ∈Hh×Qˆ h,h tales que A((σh, ph),(τh, qh)) + B1((τh, qh),(ϕˆ h,γh)) = F(τh, qh)∀(τh, qh)∈Hh, B2((σh, ph),(ψˆ h,ηh)) = 0 ∀(ψˆ h,ηh)∈Qˆ h,h . (24) La demostraci´on del siguiente teorema se puede encontrar en [4]. Teorema 4.1 El esquema de elementos finitos mixto–primal (24) poseee una ´unica soluci´on ((σh, ph),(ϕˆ h,γh)) ∈Hh×Qˆ h,h. Adem´as, si existe un δ∈(0,1] tal que σ∈ 5 S. Meddahi, G.N. Gatica, A. M´arquez [Hδ(Ωs)]2×2,div(σ)∈[Hδ(Ωs)]2,p∈H1+δ(Ωf),ϕ∈[H1/2+δ(Σ)]2yγ∈[Hδ(Ωs)]2×2, entonces se verifica k((σ, p),(ϕ,γ)) −((σh, ph),(ϕˆ h,γh))kH×Q≤Cˆ hδkϕk[H1/2+δ(Σ)]2 +C hδnkσk[Hδ(Ωs)]2×2+kdiv(σ)k[Hδ(Ωs)]2+kpkH1+δ(Ωf)+kγk[Hδ(Ωs)]2×2o, (25) con una constante C > 0independiente de hyˆ h. 5. Resultados num´ericos En esta secci´on presentamos un ejemplo que ilustra la convergencia del esquema de elementos finitos mixto–primal (24) sobre una sucesi´on de triangulaciones cuasi-uniformes del dominio. Denotamos mediante Nel n´umero de grados de libertad que define los subespacios de elementos finitos HhyQˆ h,h. Los errores individuales se definen como: e(σ) := kσ−σhkH(div; Ωs),e(p) := kp−phkH1(Ωf), e(ϕ) := kϕ−ϕˆ hk[H1/2(Σ)]2ye(γ) := kγ−γhk[L2(Ωs)]2×2. Tambi´en, sean r(σ), r(p), r(ϕ) y r(γ) los ratios de convergencia experimental: r(σ) := log(e(σ)/e0(σ)) log(h/h0), r(p) := log(e(p)/e0(p)) log(h/h0), r(ϕ) := log(e(ϕ)/e0(ϕ)) log(h/h0)yr(γ) := log(e(γ)/e0(γ)) log(h/h0), donde hyh0denotan dos tama˜nos de malla consecutivos con errores eye0. h N e(σ)r(σ)e(p)r(p)e(ϕ)r(ϕ)e(γ)r(γ) 0.0982 1297 2.20E-01 — 1.03E-01 — 2.53E-02 — 9.42E-03 — 0.0654 2763 1.45E-01 1.02 6.54E-02 1.12 1.21E-02 1.81 5.56E-03 1.30 0.0490 4885 1.09E-01 0.99 4.64E-02 1.18 7.02E-03 1.90 3.20E-03 1.91 0.0321 11506 6.93E-02 1.06 2.88E-02 1.11 2.86E-03 2.10 1.59E-03 1.63 0.0245 19616 5.30E-02 1.00 2.23E-02 0.95 1.97E-03 1.39 1.19E-03 1.09 0.0164 43140 3.57E-02 0.97 1.43E-02 1.10 9.81E-04 1.72 6.64E-04 1.44 0.0123 78271 2.59E-02 1.11 1.08E-02 0.97 5.92E-04 1.75 4.40E-04 1.43 0.0082 175630 1.73E-02 0.99 7.10E-03 1.03 2.96E-04 1.70 2.55E-04 1.34 0.0061 310084 1.31E-02 0.96 5.30E-03 1.01 1.87E-04 1.59 1.81E-04 1.18 Tabla 1: Tama˜nos de malla h, grados de libertad N, errores y convergencia. Consideramos los dominios Ωs:= ] −0,3,0,3[2y Ωf:= B(0,1) \Ωs, donde B(0,1) es el c´ırculo unidad, y elegimos los par´ametros ω= 10, ρs=ρf=λ=µ= 1, de donde κf= 1 y κs= 10. Por otra parte, sean K0,K1yK2las funciones de Bessel modificadas de segunda clase y ´ordenes 0, 1 y 2, respectivamente, y sea H(1) 0la funci´on de Hankel de 6 Un m´etodo de elementos finitos mixtos para un problema de interacci´on s´olido–fluido 0 0.5 1 1.5 2 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0 0.5 1 1.5 2 −0.02 −0.01 0 0.01 0.02 0.03 Figura 1: Partes real e imaginaria, exacta (azul) y aproximada (puntos en rojo) de ϕ1. primera clase y orden 0. Entonces, elegimos los datos fygde manera que la soluci´on exacta de (3) sea u(x) =      1 2πψ(x)−(x1−1)2 r2 1 χ(x) −(x1−1) x2 r2 1 χ(x)     ∀x∈Ωs,yp(x) = H(1) 0(ω|x|)∀x∈Ωf, donde r1:= p(x1−1)2+x2 2,ψ(x) := K0(ı ω r1) + 1 ı ω r1³K1(ı ω r1)−1 √3K1(ı ω r1 √3)´, yχ(x) := K2(ı ω r1)−1 3K2(ı ω r1 √3). La funci´on ues la soluci´on fundamental, centrada en (1,0), de la ecuaci´on elastodin´amica, lo cual supone f=0en Ωs, y pes la soluci´on fundamental, centrada en el origen, de la ecuaci´on de Helmholtz en Ωf. De esta forma, (u, p) es soluci´on de (3) con condiciones de transmisi´on no homog´eneas sobre Σ y condiciones de contorno adecuadas sobre Γ. En la Tabla 1 presentamos la convergencia del ejemplo elegido para una sucesi´on de triangulaciones cuasi-uniformes del dominio ¯ Ωs∪¯ Ωf. Podemos observar que el error dominante es e(σ), algo bastante frecuente en los esquemas de elementos finitos mixtos. Destacamos, tambi´en, que la convergencia O(h) que predice el Teorema 4.1 (cuando δ= 1) se alcanza para todas las inc´ognitas. Adem´as, en algunos casos la convergencia de e(ϕ) y e(γ) es aun m´as r´apida que O(h), lo cual puede indicar un fen´omeno de superconvergencia en estas inc´ognitas o simplemente una caracter´ıstica especial del ejemplo elegido. Finalmente, en la Figura 1 representamos las partes real e imaginaria de la componente ϕ1del multiplicador desplazamiento (para N= 43140), donde la l´ınea azul de trazo continuo es la soluci´on exacta y la soluci´on aproximada aparece con puntos rojos. 7 S. Meddahi, G.N. Gatica, A. M´arquez Referencias [1] D.N. Arnold, F. Brezzi y J. Douglas. PEERS: A new mixed finite element method for plane elasticity. Japan Journal of Applied Mathematics, vol. 1 (1984), 347–367. [2] J. Bielak y R.C. MacCamy. Symmetric finite element and boundary integral coupling methods for fluid-solid interaction. Quarterly of Applied Mathematics, vol. 49 (1991), 107–119. [3] F. Brezzi y M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Verlag, 1991. [4] G.N. Gatica, A. M´arquez y S. Meddahi. Analysis of the coupling of primal and dual-mixed finite element methods for a two-dimensional fluid-solid interaction problem. Aceptado en SIAM Journal on Numerical Analysis, (2007). [5] F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, Springer-Verlag, New York, 1998. [6] A. M´arquez, S. Meddahi y V. Selgas. A new BEM-FEM coupling strategy for two-dimensional fluidsolid interaction problems. Journal of Computational Physics, vol. 199 (2004), 205–220. [7] S. Meddahi y F.-J. Sayas. Analysis of a new BEM-FEM coupling for two dimensional fluid-solid interaction. Numerical Methods for Partial Differential Equations, vol. 21 (2005), 1017–1042. [8] J.E. Roberts y J.M. Thomas, Mixed and Hybrid Methods. En: Handbook of Numerical Analysis, editado por P.G. Ciarlet y J.L. Lions, vol. II, Finite Element Methods (Part 1), North-Holland, Amsterdam, 1991. [9] R. Stenberg. A family of mixed finite elements for the elasticity problem. Numerische Mathematik, vol. 53, 5 (1988), 513–538. 8