Estudo do problema de deformación dunha placa elástica rectangular usando o operador bilaplaciano

Traballo Fin de Grao Estudo do problema de deformación dunha placa elástica rectangular usando o operador bilaplaciano María del Carmen Mosquera Vázquez 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Estudo do problema de deformación dunha placa elástica rectangular usando o operador bilaplaciano María del Carmen Mosquera Vázquez 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Dedicado a: Meus pais e avós por apoiarme durante estes anos. José Antonio Álvarez Dios e María del Carmen Muñiz Castiñeira polo esforzo e axuda para realizar este traballo. Traballo proposto Área de Coñecemento: Matemática Aplicada Título: Estudo do problema de deformación dunha placa elástica rectangular usando o operador bilaplaciano Breve descrición do contido Estudarase o operador bilaplaciano para describir un modelo de deformación dunha placa cunhas determinadas condicións de borde. Dende un punto de vista práctico, abordarase a resolución de dito problema empregando o método de diferenzas nitas. Recomendacións Dominio da teoría clásica de derivación numérica, o método de diferenzas nitas para resolver ecuacións diferenciais en derivadas parciais, e unha linguaxe de programación (MatLab). Outras observacións III Índice xeral Resumo VII 1. Ecuacións 3D da elasticidade lineal 5 1.1. Modeladodoproblema.............................. 5 1.2. Ecuacións de equilibrio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3. Ecuacións constitutivas dos materiais elásticos . . . . . . . . . . . . . . . . . 7 1.4. Ecuacións da elasticidade non lineal . . . . . . . . . . . . . . . . . . . . . . 9 1.5. Ecuacións da elasticidade lineal . . . . . . . . . . . . . . . . . . . . . . . . . 9 2. Modelo de placa en exión 11 2.1. Principio dos traballos virtuais . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2. Relacións constitutivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3. Hipóteses de Kirchho-Love . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3. Estudo do problema discreto 21 3.1. Construción dunha malla . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2. Discretización do operador bilaplaciano . . . . . . . . . . . . . . . . . . . . . 23 3.3. Construción da matriz Lh ............................ 29 4. Exemplos numéricos 37 4.1. Exemplo analítico de comprobación . . . . . . . . . . . . . . . . . . . . . . . 37 4.2. Exemplosdeplacas................................ 41 A. Programas empregados 45 Bibliografía 53 V Capítulo 1 Obtención das ecuacións tridimensionais da elasticidade lineal 1.1. Modelado do problema O desprazamento e a tensión dun corpo elástico en resposta a unha determinada carga son predecibles como a solución dun sistema de ecuacións en derivadas parciais. O antedito sistema deriva de leis físicas en forma de dous conxuntos de ecuacións: as ecuacións de equilibrio e as ecuacións constitutivas. As ecuacións da elasticidade lineal obtéñense a partir das ecuacións da elasticidade non lineal, prescindindo dos termos non lineais no tensor de deformacións, como vemos a continuación. A obtención das ecuacións de equilibrio e as ecuacións constitutivas dadas neste capítulo pretende ser unha breve introducción ás mesmas, para unha mellor comprensión da memoria, non sendo a súa formulación un obxectivo da mesma. 1.2. Ecuacións de equilibrio Nesta sección, comezamos estudando o desprazamento e a tensión do corpo elástico en resposta ás forzas aplicadas. Consideramos un corpo que ocupa a clausura dun dominio Ω⊂R3 , sobre o que non actúan forzas e ao que chamamos conguración de referencia do corpo . Calquera outra conguración do corpo que poida tomar ao aplicarlle determinadas forzas queda establecida mediante unha deformación, denida por unha aplicación φ : Ω→R3 que conserva a orientación (é dicir, det ∇φ(x)>0 , ∀x∈Ω ) e inxectiva no aberto Ω . Á imaxe φ(Ω) chamarémoslle conguración deformada do corpo denida pola deformación φ . 5 6 CAPÍTULO 1. ECUACIÓNS 3D DA ELASTICIDADE LINEAL Denotamos a imaxe por φ de calquera subconxunto A⊂Ω mediante φ(A) = ˜ A . A diferenza entre a conguración deformada e a conguración de referencia do corpo denomínase desprazamento , que vén dado polo campo vectorial u:= φ−I. (1.1) O obxectivo desta sección é determinar a maior cantidade de información posible sobre a deformación do corpo, o cal se atopa en equilibrio estático baixo a acción das forzas aplicadas. As forzas aplicadas nunha conguración deformada represéntanse polas densidades ˜ f : ˜ Ω→R3 (forza volúmica), ˜g : ˜ Γ1→R3 (forza supercial), onde Γ1 é un subconxunto da fronteira Γ . Aplicando o axioma fundamental debido a Euler e Cauchy (ver Mardare [8]), postúlase que o equilibrio do corpo coas forzas aplicadas se reicte na existencia dunha forza que actúa na fronteira de calquera dominio ˜ A⊂˜ Ω , que depende só do punto ˜x e do vector normal ao corpo no dito punto, que denotamos por n(˜x) . Como consecuencia, existe un campo tensorial ˜ T:˜ Ω→M3×3 , tal que − div ˜ T(˜x) = ˜ f(˜x),∀˜x∈˜ Ω, (1.2) ˜ T(˜x)˜n(˜x) = ˜g(˜x),∀˜x∈˜ Γ1, (1.3) ˜ T(˜x)∈S3×3,∀˜x∈˜ Ω. (1.4) A estas ecuacións denomínaselles ecuacións de equilibrio da conguración deformada. As ditas ecuacións poden reformularse na conguración de referencia (coñecida) mediante un cambio de variable. Para este n denimos o primeiro tensor de tensións de Piola-Kirchho dado por T(x) := ˜ T(φ(x)) cof ∇φ(x),∀x∈Ω. (1.5) Tendo en conta a identidade (∇φ(x))−1T(x)=(∇φ(x))−1( det (∇φ(x)) ˜ T(φ(x)))(∇φ(x))−t, é inmediato vericar tendo en conta (1.4) que o segundo tensor de tensións de PiolaKirchho denido por Σ(x) := (∇φ(x))−1T(x) , ∀x∈Ω (1.6) é simétrico. Ademais, verifícanse as seguintes ecuacións de equilibrio na conguración de referencia: − div T(x) = f(x),∀x∈Ω, (1.7) 1.3. ECUACIÓNS CONSTITUTIVAS DOS MATERIAIS ELÁSTICOS 7 T(x)n(x) = g(x),∀x∈Γ1, (1.8) Σ(x)∈S3×3,∀x∈Ω, (1.9) onde f(x) = ˜ f(φ(x)) , g(x) = ˜g(φ(x)) . 1.3. Ecuacións constitutivas dos materiais elásticos O tensor de tensións depende da deformación inducida polas forzas aplicadas ao corpo. Esta dependencia reíctese nas ecuacións constitutivas do material. De (1.1) dedúcese que φ(x) = x+u(x) , ∀x= (x1, x2, x3)∈Ω, (1.10) polo tanto, (∂jφi) = (δij)+(∂jui), i, j = 1,2,3, ou, escrito matricialmente, ∇φ=I+∇u. (1.11) Calculemos agora canto se deforma o segmento [x, x +dx] mediante a deformación φ , supoñéndoa sucientemente regular. Denimos unha parametrización do segmento deformado como: [0,1] →[x, x +dx]→R3 t;s(t) = x+tdx ;φ(s(t)) := ˜s(t) de modo que a lonxitude do segmento deformado é: lφ=Z1 0 ||˜s0(t)||dt. Por outra parte, ˜s(t)=(φ1(s(t)), φ2(s(t)), φ3(s(t))) , entón ˜s0(t) = ∇φ(s(t))s0(t) = ∇φ    dx1 dx2 dx3    , 8 CAPÍTULO 1. ECUACIÓNS 3D DA ELASTICIDADE LINEAL de onde se deduce ||˜s(t)|| = (˜s0(t)) ·˜s0(t))1 2= ((˜s0(t))t˜s0(t))1 2= ((dx1dx2dx3)(∇φ(s(t)))t∇φ(s(t))     dx1 dx2 dx3    )1 2. Chamando C= (∇φ)t∇φ ao denominado tensor métrico tense lφ=Z1 0 ((dx1dx2dx3)C(s(t))     dx1 dx2 dx3    )1 2dt, polo que, a lonxitude da bra dependerá do tensor métrico C . Ademais, tendo en conta (1.11), obtemos C= (I+∇u)t(I+∇u)=(I+ (∇u)t)(I+∇u) = I+∇u+ (∇u)t+ (∇u)t∇u. Denindo o tensor de deformacións de Green-St Venant E =1 2(∇u+ (∇u)t+ (∇u)t∇u), (1.12) verifícase que C=I+ 2 E . (1.13) Por outra parte, defínese o tensor de deformacións de Green-St Venant linealizado como e(u) = 1 2(∇u(x)+(∇u(x))t). (1.14) Finalmente, dedúcese a partir do teorema de Rivlin-Ericksen que a ecuación constitutiva para un material elástico de Saint Venant-Kirchho é (ver Destuynder [3]): Σ(x) = λ( tr ( E (x))I+ 2µ E (x) , ∀x∈Ω, (1.15) sendo λ e µ as constantes de Lamé do material. A ecuación constitutiva (1.15) é invertible no sentido de que E pódese expresar como función de Σ da seguinte forma E (x) = 1 + ν EΣ(x)−ν E( tr (Σ(x))I , ∀x∈Ω, sendo ν o coeciente de Poisson e E o módulo de Young , os cales verican as desigualdades E > 0 e 0< ν < 1 2 . Podemos expresalos en termos das constantes de Lamé da seguinte forma ν=λ 2(λ+µ) , E=µ(3λ+ 2µ) λ+µ. (1.16) 1.4. ECUACIÓNS DA ELASTICIDADE NON LINEAL 9 Tamén podemos expresar os coecientes λ e µ en función do módulo de Young e o coeciente de Poisson. A partir da expresión de E , E=µ3λ+ 2µ λ+µ=µλ λ+µ+ 2. Agora empregamos a expresión (1.16) de ν e obtemos, E=µ(2ν+ 2) = 2µ(ν+ 1). Despexando µ , chegamos a: µ=E 2(ν+ 1). (1.17) Despexando λ na expresión (1.16) de ν e empregando (1.17) obtemos a expresión λ=Eν (1 + ν)(1 −2ν). 1.4. Ecuacións da elasticidade non lineal Combinando as ecuacións de equilibrio coas ecuacións constitutivas do material de Saint Venant-Kirchho, conclúese que a deformación orixinada no corpo como resposta á presenza das forzas aplicadas satisfai o problema de fronteira non lineal: − div T(x) = f(x),∀x∈Ω, (1.18) φ(x) = x, x ∈Γ0, (1.19) T(x)n(x) = g(x),∀x∈Γ1, (1.20) Σ(x) = λ( tr (E(x))I+ 2µE(x) , ∀x∈Ω, (1.21) sendo Γ0:= ∂Ω\Γ1. 1.5. Ecuacións da elasticidade lineal Estas ecuacións obtéñense a partir das ecuacións non lineais da elasticidade baixo a hipótese de pequenas deformacións, é dicir, φ(x) = x+u(x) con |∇u(x)|  1 , ∀x∈Ω. 10 CAPÍTULO 1. ECUACIÓNS 3D DA ELASTICIDADE LINEAL Un caso particular importante destas ecuacións formúlase cando o corpo está feito dun material elástico isótropo 1 e homoxéneo 2 . Tendo en conta o tensor de deformacións linealizado dado en (1.14) e a ecuación constitutiva dun material de St Venant-Kirchho (1.15), a parte lineal do tensor T(x) queda T(x) = ∇φ(x)Σ(x) = ∇φ(x)(λ( tr (e(x))I+ 2µe(x)) = = (I+∇u(x))(λ( tr (e(x))I+ 2µe(x)) = =λ( tr (e(x))I+ 2µe(x) + O(|∇u(x)|) =: σ(x) + O(|∇u(x)|). (1.22) Finalmente, as ecuacións da elasticidade lineal para un material homoxéneo e isótropo veñen dadas por − div σ(x) = f(x),∀x∈Ω, u(x)=0,∀x∈Γ0, σ(x)n(x) = g(x),∀x∈Γ1, (1.23) onde σ(x) = λ( tr (e(x))I+ 2µe(x) sendo e(x) = 1 2((∇u(x))t+∇u(x)). 1 O comportamento elástico do corpo non depende da dirección considerada. 2 Ten as mesmas propiedades en todos os puntos. Capítulo 2 Modelo de placa en exión Neste capítulo consideramos como referencia un sistema ortonormal. Supoñemos que ω é un aberto de R2 que se sitúa no plano x3= 0 , e que corresponde á sección media da placa. Sexa Ωε=ω×(−ε, ε) un sólido que tomamos como conguración de referencia, sendo a súa fronteira lateral Γε 0=∂ω ×(−ε, ε) que consideramos empotrada. Supoñemos que non actúan forzas superciais e as forzas volúmicas veñen dadas por f= (fi) con i= 1,2,3 . Recordamos a denición do tensor de deformacións de Green-St Venant linealizado (ver 1.14) eij(u) = 1 2(∂iuj+∂jui), i, j = 1,2,3. Denimos o espazo das deformacións admisibles: Vε=v= (vi)∈[H1(Ωε)]3;vi= 0 sobre Γε 0, i = 1,2,3. 11 12 CAPÍTULO 2. MODELO DE PLACA EN FLEXIÓN 2.1. Principio dos traballos virtuais Denotando por σ= (σij) o tensor de tensións (σij =σji) , entón dito tensor satisfai o principio dos traballos virtuais (ver Destuynder [3]): ZΩε σijeij(v) = ZΩε fivi , ∀v∈Vε, (2.1) onde se adopta o convenio de que os índices latinos se toman de 1 a 3 e os gregos de 1 a 2, ademais de que os índices repetidos súmanse, así σijeij(v) = σ11e11(v) + σ22e22(v) + σ33e33(v)+2σ12e12(v)+2σ13e13(v)+2σ23e23(v) . Nótese que empregando a fórmula de Green (ver Brezis [1]) e tomando v∈[C∞ c(Ω)ε]3 en (2.1), chégase ás seguintes ecuacións de equilibrio local: ∂iσij +fj= 0 en Ω, j = 1,2,3. (2.2) 2.2. Relacións constitutivas Cando o material do corpo é homoxéneo e isótropo, o modelo de elasticidade lineal proporciona as seguintes relacións constitutivas entre tensións e deformacións: σij =E 1 + νeij(u) + ν 1−2νepp(u)δij , (2.3) onde E é o módulo de Young e ν o coeciente de Poisson. A relación (2.3) pódese invertir para obter: eij(u) = 1 + ν Eσij −ν Eσppδij . (2.4) Tendo en conta (2.1), denimos o problema de elasticidade tridimensional da seguinte forma: ( Atopar uε∈Vε tal que: a(uε, v) = l(v) , ∀v∈Vε, (2.5) onde a(·,·) é a forma bilineal denida por a(u, v) = ZΩε E 1 + ν{eij(u)eij(v) + ν 1−2νepp(u)eqq(v)}, (2.6) e l(·) é a forma lineal denida por l(v) = ZΩε fivi. (2.7) 2.2. RELACIÓNS CONSTITUTIVAS 13 Empregando o Teorema de Lax-Milgram 1 (ver Brezis [1]), pódese probar que (2.5) ten solución única. Consideramos o espacio de tensións: Σε={τ= (τij)∈L2(Ωε)9 , τij =τji}, dotado da norma: kτkΣε=  3 X i,j=1 kτijk2 L2  1 2 . Se uε é a solución de (2.5), asociámoslle o tensor de tensións: σε ij =E 1 + ν{eij(uε) + ν 1−2νepp(uε)δij }. Entón (2.5) é equivalente á formulación mixta de Hellinger-Reissner:                      Atopar (σε, uε)∈Σε×Vε tal que ZΩε (1 + ν Eσε ijτij −ν Eσε ppτqq)−ZΩε τijeij(uε)=0 , ∀τ∈Σε, ZΩε σε ijeij(v) = ZΩε fivi , ∀v∈Vε. (2.8) 1 Sexa Vε un espazo de Hilbert e a:Vε×Vε→R unha forma bilineal continua e coercitiva, entón para toda l:Vε→R lineal e continua, existe un único v∈Vε tal que a(u, v) = l(v) , para todo v∈Vε . Capítulo 3 Estudo do problema discreto Para simplicar a notación poñemos neste capítulo, x1=x e x2=y ; así denotarase por ux e uy as derivadas parciais respecto de x e y , respectivamente; analogamente, as de segunda orde uxx , uxy , uyx e uyy , e así sucesivamente. Por outra parte recordamos a denición do operador bilaplaciano da función regular u no punto (a, b) , ∆2u(a, b) = uxxxx(a, b)+2uxxyy(a, b) + uyyyy(a, b). (3.1) 3.1. Construción dunha malla A partir de agora consideramos soamente o caso dunha placa rectangular. Sexa ω= (0, Lx)×(0, Ly) o dominio correspondente á sección media da placa. Consideramos unha malla uniforme con N+ 2 puntos na dirección x e M+ 2 puntos na dirección y de modo que as lonxitudes internodais son hx=Lx N+1 e hy=Ly M+1 respectivamente. Supoñemos que se verica que hx=hy=h , así podemos denir a malla como ωh={(hi, hj),1⩽i⩽N, 1⩽j⩽M}. 21 22 CAPÍTULO 3. ESTUDO DO PROBLEMA DISCRETO Figura 3.1: Malla uniforme ωh e conxunto fronteira ∂ωh . Empregando esta notación, un punto do interior do conxunto (i, j) ten por coordenadas (hi, hj) para i= 1...N, j = 1...M e un punto da fronteira pertence ao conxunto ∂ωh={(0, hj),(Lx, hj),0⩽j⩽M+ 1}[{(hi, 0),(hi, Ly),0⩽i⩽N+ 1}. En ∂ωh a solución é coñecida; entón para ter unha solución do problema é suciente con calcular a solución nos puntos de ωh . 3.2. DISCRETIZACIÓN DO OPERADOR BILAPLACIANO 23 3.2. Discretización do operador bilaplaciano Para a discretización do operador bilaplaciano empregaremos un desenvolvemento de Taylor de orde 6 para unha función no punto (i, j) , coa malla que xa construímos. Figura 3.2: Esquema de trece puntos. O desenvolvemento de Taylor de orde 6 proporciónanos un esquema de trece puntos. Para a discretización empregamos o seguinte resultado. Teorema 3.1. Sexa S un rectángulo aberto que conteña a P1, ... , P13 puntos do esquema representado na Figura 3.2 con P7= (a, b) , e sexa v:S→R unha función de clase 6 en S . Entón existen c1, ... , c13 ∈R tales que: 13 X i=1 civ(Pi) = h4(vxxxx(a, b)+2vxxyy(a, b) + vyyyy(a, b)) + O(h6). (3.2) Demostración. Escribimos primeiro o desenvolvemento de Taylor de orde 6 para unha 24 CAPÍTULO 3. ESTUDO DO PROBLEMA DISCRETO función v de dúas variables arredor de (a, b) avaliado no punto (a+h, b +k) , v(a+h, b +k) = v(a, b) + h vx(a, b) + k vy(a, b)+ +1 2!(h2vxx(a, b)+2hk vxy(a, b) + k2vyy(a, b)) +1 3!(h3vxxx(a, b)+3h2k vxxy(a, b)+3hk2vxyy(a, b) + k3vyyy(a, b)) +1 4!(h4vxxxx(a, b)+4h3k vxxxy(a, b)+6h2k2vxxyy(a, b)+4hk3vxyyy(a, b)+ +k4vyyyy(a, b)) + 1 5!(h5vxxxxx(a, b)+5h4k vxxxxy(a, b)+ + 10h3k2vxxxyy(a, b) + 10h2k3vxxyyy(a, b)+5hk4vxyyyy(a, b)+ +k5vyyyyy(a, b)) + O(hαk6−α), sendo α índice de suma, α= 0,1, ... , 6 no termo do resto. Agora aplicamos a expresión anterior a cada un dos puntos do esquema anteriormente denido de 13 puntos: v(P1) = v(a, b)−2h vy(a, b) + (−2h)2 2! vyy(a, b) + (−2h)3 3! vyyy(a, b) + (−2h)4 4! vyyyy(a, b)+ +(−2h)5 5! vyyyyy(a, b) + O(h6), v(P2) = v(a, b)−h(vx(a, b) + vy(a, b))+ +(−h)2 2! (vxx(a, b)+2vxy(a, b) + vyy(a, b))+ +(−h)3 3! (vxxx(a, b)+3vxxy(a, b)+3vxyy(a, b) + vyyy(a, b))+ +(−h)4 4! (vxxxx(a, b)+4vxxxy(a, b)+6vxxyy(a, b)+4vxyyy(a, b) + vyyyy(a, b))+ +(−h)5 5! (vxxxxx(a, b)+5vxxxxy(a, b) + 10 vxxxyy(a, b) + 10 vxxyyy(a, b)+ + 5 vxyyyy(a, b) + vyyyyy(a, b)) + O(h6), v(P3) = v(a, b)−h vy(a, b) + (−h)2 2! vyy(a, b) + (−h)3 3! vyyy(a, b) + (−h)4 4! vyyyy(a, b)+ +(−h)5 5! vyyyyy(a, b) + O(h6), 3.2. DISCRETIZACIÓN DO OPERADOR BILAPLACIANO 25 v(P4) = v(a, b) + h(vx(a, b)−vy(a, b))+ +h2 2! (vxx(a, b)−2vxy(a, b) + vyy(a, b))+ +h3 3! (vxxx(a, b)−3vxxy(a, b)+3vxyy(a, b)−vyyy(a, b))+ +h4 4! (vxxxx(a, b)−4vxxxy(a, b)+6vxxyy(a, b)−4vxyyy(a, b) + vyyyy(a, b))+ +h5 5! (vxxxxx(a, b)−5vxxxxy(a, b) + 10 vxxxyy(a, b)−10 vxxyyy(a, b)+ + 5 vxyyyy(a, b)−vyyyyy(a, b)) + O(h6), v(P5) = v(a, b)−2h vx(a, b) + (−2h)2 2! vxx(a, b) + (−2h)3 3! vxxx(a, b) + (−2h)4 4! vxxxx(a, b)+ +(−2h)5 5! vxxxxx(a, b) + O(h6), v(P6) = v(a, b)−h vx(a, b) + (−h)2 2! vxx(a, b) + (−h)3 3! vxxx(a, b) + (−h)4 4! vxxxx(a, b)+ +(−h)5 5! vxxxxx(a, b) + O(h6), v(P7) = v(a, b) , v(P8) = v(a, b) + h vx(a, b) + h2 2! vxx(a, b) + h3 3! vxxx(a, b) + h4 4! vxxxx(a, b) + h5 5! vxxxxx(a, b)+ +O(h6), v(P9) = v(a, b)+2h vx(a, b) + (2h)2 2! vxx(a, b) + (2h)3 3! vxxx(a, b) + (2h)4 4! vxxxx(a, b)+ +(2h)5 5! vxxxxx(a, b) + O(h6), v(P10) = v(a, b)−h(vx(a, b)−vy(a, b))+ +h2 2! (vxx(a, b)−2vxy(a, b) + vyy(a, b))+ +h3 3! (−vxxx(a, b)+3vxxy(a, b)−3vxyy(a, b) + vyyy(a, b))+ +h4 4! (vxxxx(a, b)−4vxxxy(a, b)+6vxxyy(a, b)−4vxyyy(a, b) + vyyyy(a, b))+ +h5 5! (−vxxxxx(a, b)+5vxxxxy(a, b)−10 vxxxyy(a, b) + 10 vxxyyy(a, b)−5vxyyyy(a, b)+ +vyyyyy(a, b)) + O(h6), v(P11) = v(a, b) + h vy(a, b) + h2 2! vyy(a, b) + h3 3! vyyy(a, b) + h4 4! vyyyy(a, b) + h5 5! vyyyyy(a, b)+ +O(h6), 26 CAPÍTULO 3. ESTUDO DO PROBLEMA DISCRETO v(P12) = v(a, b) + h(vx(a, b) + vy(a, b))+ +h2 2! (vxx(a, b)+2vxy(a, b) + vyy(a, b))+ +h3 3! (vxxx(a, b)+3vxxy(a, b)+3vxyy(a, b) + vyyy(a, b))+ +h4 4! (vxxxx(a, b)+4vxxxy(a, b)+6vxxyy(a, b)+4vxyyy(a, b) + vyyyy(a, b))+ +h5 5! (vxxxxx(a, b)+5vxxxxy(a, b) + 10 vxxxyy(a, b) + 10 vxxyyy(a, b)+5vxyyyy(a, b)+ +vyyyyy(a, b)) + O(h6), v(P13) = v(a, b)+2h vy(a, b) + (2h)2 2! vyy(a, b) + (2h)3 3! vyyy(a, b) + (2h)4 4! vyyyy(a, b)+ +(2h)5 5! vyyyyy(a, b) + O(h6). Tendo en conta as expresións para v(Pi) , i= 1 ... 13 anteriores, obtemos un sistema de ecuacións lineais sendo ci , i= 1 ... 13 as incógnitas. Identicamos en ambos membros da expresión (3.2) os coecientes das derivadas que aparecen ata orde 4:                                  1 ↓2 ↓3 ↓4 ↓5 ↓6 ↓7 ↓8 ↓9 ↓1 ↓0 1 ↓1 1 ↓2 1 ↓3 v→1 1 1 1 1 1 1 1 1 1 1 1 1 vx→0−h0h−2h−h0h2h−h0h0 vy→ −2h−h−h−h0 0 0 0 0 h h h 2h vxx →0h2 20h2 22h2h2 20h2 22h2h2 20h2 20 vxy →0h20−h20 0 0 0 0 −h20h20 vyy →2h2h2 2 h2 2 h2 20 0 0 0 0 h2 2 h2 2 h2 22h2 vxxx →0−h3 60h3 6 −4h3 3 −h3 60h3 6 4h3 3 −h3 60h3 60 vxxy →0−h3 20−h3 20 0 0 0 0 h3 20h3 20 vxyy →0−h3 20h3 20 0 0 0 0 −h3 20h3 20 vyyy →−4h3 3 −h3 6 −h3 6 −h3 60 0 0 0 0 h3 6 h3 6 h3 6 4h3 3 vxxxx →0h4 24 0h4 24 2h4 3 h4 24 0h4 24 2h4 3 h4 24 0h4 24 0 vxxxy →0h4 60−h4 60 0 0 0 0 −h4 60h4 60 vxxyy →0h4 40h4 40 0 0 0 0 h4 40h4 40 vxyyy →0h4 60−h4 60 0 0 0 0 −h4 60h4 60 vyyyy →2h4 3 h4 24 h4 24 h4 24 0 0 0 0 0 h4 24 h4 24 h4 24 2h4 3                                                              c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13                             =                                  0 0 0 0 0 0 0 0 0 0 h4 0 2h4 0 h4                                  Xa que as las 5, 12 e 14 da matriz ampliada do sistema son linealmente dependentes, podemos prescindir delas quedando un sistema de ecuacións lineal con 13 ecuacións e 13 3.2. DISCRETIZACIÓN DO OPERADOR BILAPLACIANO 27 incógnitas. Simplicando h coas súas correspondentes potencias, podemos escribir:                               1 ↓2 ↓3 ↓4 ↓5 ↓6 ↓7 ↓8 ↓9 ↓1 ↓0 1 ↓1 1 ↓2 1 ↓3 v→1 1 1 1 1 1 1 1 1 1 1 1 1 vx→0−1 0 1 −2−1 0 1 2 −1 0 1 0 vy→ −2−1−1−1 0 0 0 0 0 1 1 1 2 vxx →01 201 221 201 221 201 20 vxy →010−1 0 0 0 0 0 −1 0 1 0 vyy →21 2 1 2 1 20 0 0 0 0 1 2 1 2 1 22 vxxx →0−1 601 6 −4 3 −1 601 6 4 3 −1 601 60 vxxy →0−1 20−1 20 0 0 0 0 1 201 20 vxyy →0−1 201 20 0 0 0 0 −1 201 20 vyyy →−4 3 −1 6 −1 6 −1 60 0 0 0 0 1 6 1 6 1 6 4 3 vxxxx →01 24 01 24 2 3 1 24 01 24 2 3 1 24 01 24 0 vxxyy →01 401 40 0 0 0 0 1 401 40 vyyyy →2 3 1 24 1 24 1 24 0 0 0 0 0 1 24 1 24 1 24 2 3                                                            c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13                              =                              0 0 0 0 0 0 0 0 0 0 1 2 1                              O sistema anterior é compatible determinado xa que o determinante da matriz asociada é 1 3 e polo tanto distinto de 0 . Así obtemos os valores seguintes das constantes: c1= 1 c2= 2 c3=−8 c4= 2 c5= 1 c6=−8 c7= 20 c8=−8 c9= 1 c10 = 2 c11 =−8 c12 = 2 c13 = 1 28 CAPÍTULO 3. ESTUDO DO PROBLEMA DISCRETO Comprobamos que se anulan os térmos correspondentes ás derivadas de orde 5 tendo en conta as constantes que calculamos anteriormente. Coecientes de vxxxxx : c5 (−5h)5 5! +c6 (−h)5 5! +c8 h5 5! +c9 (2h)5 5! +c4 h5 5! +c10 (−h)5 5! +c12 h5 5! =h5 5! (−25c5−c6+ c8+ 25c9−c2+c4−c10 +c12) = 0 Coecientes de vxxxxy : c2 5h5 5! −c4 5h5 5! +c10 5h5 5! +c12 5h5 5! =5h5 5! (−c2−c4+c10 +c12)=0 Coecientes de vxxxyy : c2 −10h5 5! +c4 10h5 5! −c10 −10h5 5! +c12 10h5 5! =10h5 5! (−c2+c4−c10 +c12)=0 Coecientes de vxxyyy : c2 −10h5 5! +c4 10h5 5! −c10 −10h5 5! +c12 10h5 5! =10h5 5! (−c2+c4−c10 +c12)=0 Coecientes de vxyyyy : c2 −5h5 5! +c4 5h5 5! +c10 −5h5 5! +c12 5h5 5! =5h5 5! (−c2+c4−c10 +c12)=0 Coecientes de vyyyyy : c5 (−2h)5 5! +c6 (−h)5 5! +c8 h5 5! +c9 (2h)5 5! +c4 (−h)5 5! +c10 (−h)5 5! +c12 h5 5! =h5 5! (−25c5− c6+c8+ 25c9−c2+c4−c10 +c12)=0 Para toda función u denimos Uij =u(hi, hj), i = 1, ... N, j = 1, ... M, sendo (hi, hj), i = 1, ...N, j = 1, ...M puntos da malla ωh . Ademais dotámolos da orde seguinte: ˜u= [U11, U21, ... UN1, U12, ... , UN2, ... ... , UNM ]t. Coa devandita notación Uij é a k -ésima compoñente de ˜u , con k=i+ (j−1)N . Analogamente referímonos a calquera compoñente do vector de RN M con doble índice (i, j) ou con 3.3. CONSTRUCIÓN DA MATRIZ LH 29 índice k=i+ (j−1)N indistintamente. Denimios o operador lineal Lh:RN M →RN M como: (Lh˜u)k=1 h4(20Ui,j −8(Ui−1,j +Ui+1,j +Ui,j−1+Ui,j+1)+ + 2(Ui−1,j−1+Ui−1,j+1 +Ui+1,j−1+Ui+1,j+1)+ + (Ui−2,j +Ui+2,j +Ui,j−2+Ui,j+2), (3.3) cuxa matriz asociada denotamos tamén por Lh se non houbese risco de confusión. Construímos a dita matriz na seguinte sección e veremos que a matriz resultante é pentadiagonal por bloques grazas á orden dos nodos elixida anteriormente. Empregando (3.3) e (3.2) obtemos que (Lh˜u)k= ∆2u(hi, hj) + O(h2), para toda función u∈C6(ω) . Entón a discretización é consistente de orde dúas, é dicir, que o operador discreto aproxima ao continuo nos puntos da malla coa orde anterior (ver Lui [7]). 3.3. Construción da matriz Lh Figura 3.3: Malla ωh xunto cos nodos cticios. No noso problema xa coñecemos a solución nos puntos da fronteira así como a súa derivada normal nos mesmos, entón abonda con resolver o problema para os puntos do interior do Capítulo 4 Exemplos numéricos 4.1. Exemplo analítico de comprobación Neste capítulo empezamos comprobando, mediante un exemplo cuxa solución coñecemos de antemán, a bondade dos cálculos realizados nos capítulos anteriores. Esto permitiranos vericar se está ben aproximado o problema (2.23) e tamén se a orde de converxencia é a esperada. Consideramos ω= (0,1) ×(0,1) e a función u(x, y) = sen 2(πx) sen 2(πy), (4.1) sendo o seu bilaplaciano, ∆2u(x, y) = 8 π4(3 sen 2(πx) sen 2(πy) + cos 2(πx) cos 2(πy) −2 cos 2(πx) sen 2(πy)−2 cos 2(πy) sen 2(πx)). Para que a función (4.1) sexa solución do problema (2.23) debemos tomar, F3=2 3 ε3E 1−ν2∆2u(x, y). Por simplicidade, neste primeiro exemplo supoñemos que a constante que multiplica ao bilaplaciano é a unidade, é dicir, 2 3 ε3E 1−ν2= 1. Recordemos que para resolver o problema teórico impuxemos as condicións de que tanto a función como a súa derivada se anulan sobre ∂ω . Para resolver de forma numérica o problema consideramos a malla que construímos no capítulo anterior, coa notación (fh)k=F3(ih, jh) con k=i+ (j−1)N, i = 1, ... , N, j = 1, ... , M , e resolvemos o sistema para obter unha solución aproximada uh . 37 38 CAPÍTULO 4. EXEMPLOS NUMÉRICOS Para comprobar a converxencia do método, vemos a relación entre algúns valores de N e M e o erro cometido na aproximación. Este valor calculámolo mediante a norma euclídea e a discreta de L2 , recordamos a denición desta última: |x|h=h|x|2. Na seguinte táboa mostramos a norma euclídea e a norma L2 discreta do erro cometido eh=uh−u sendo u a solución exacta de (1.7) e uh a solución aproximada. O erro relativo coa norma euclídea sería: er,2=|eh|2 |u|2 para os diferentes valores de M e N , obtemos así a seguinte táboa, 4.1. EXEMPLO ANALÍTICO DE COMPROBACIÓN 39 N M h |eh|2|eh|her,2 5 5 0.1667 0.4558 0.0759 0.2026 10 10 0.0909 0.2312 0.0210 0.0560 20 20 0.0476 0.1185 0.0056 0.0150 40 40 0.0244 0.0603 0.0015 0.0039 80 80 0.0123 0.0304 0.0004 0.0010 160 160 0.0062 0.0153 0.0001 0.0003 Táboa 4.1: Erros para o exemplo numérico. Na Figura (4.1) representamos |eh|2 fronte a h obténdose unha recta de pendente 1, o cal indica que o método é converxente de orde 1 coa norma euclídea. Na Figura (4.2) representamos a norma L2 discreta do erro fronte a h en escala lineal. Na Figura (4.3) representamos o mesmo en escala log-log. A pendente da recta na Figura (4.3) é 2, polo que se conrma que o método é converxente de orde 2 na norma L2 discreta, (ver Lui [7], p.32). Figura 4.1: Norma euclídea do erro fronte h . 40 CAPÍTULO 4. EXEMPLOS NUMÉRICOS Figura 4.2: Norma L2 discreta do erro fronte h . Figura 4.3: Norma L2 discreta do erro fronte h en escala log-log. 4.2. EXEMPLOS DE PLACAS 41 4.2. Exemplos de placas A continuación vemos que ocorre cunha placa de silicio ao aplicarlle distintos tipos de forzas: unha forza constante en toda a placa, unha forza arredor do centro da placa e nula no resto, e por último unha forza negativa arredor do centro da placa e positiva no resto. Tomamos unha placa cadrada de 10−3 m de lado. Vemos agora cales son as constantes físicas do material 1 : o módulo de Young E= 1,3×1011 Pa, o coeciente de Poisson µ= 0,25 , o semiespesor ε= 0,5×10−6 m e, por último, tomamos N=M= 40 . Consideramos agora un primeiro exemplo de forza: tomaremos F3=−1000 Pa, constante en toda a placa. Na Figura (4.4) mostramos a solución numérica correspondente. Observamos que a placa exiona ata un desprazamento máximo menor que 10−4 m. Observamos tamén que a condición de empotramento se respeta nos bordes. Figura 4.4: Solución correspondente á forza constante en toda a placa. Consideramos agora unha forza constante arredor do centro da placa, F3=−1000 Pa, e nula no resto (ver Figura 4.5). 1 Pa=Pascal, m=metro. 42 CAPÍTULO 4. EXEMPLOS NUMÉRICOS Figura 4.5: Forza arredor do centro da placa. Na Figura 4.6 vemos a correspondente solución do problema ao aplicarlle a forza anterior. Observamos que a placa exiona ata un desprazamento máximo menor que 4×10−6 m, inferior ao caso anterior. Figura 4.6: Solución uh arredor do centro da placa. Por último consideramos unha forza que toma arredor do centro da placa un valor negativo 4.2. EXEMPLOS DE PLACAS 43 F3=−3000 Pa e un valor positivo no resto F3= 100 Pa (ver Figura 4.7). Na Figura 4.8 móstrase a solución correspondente a esta forza. Notemos que aparece un efecto de alabeado causado pola interacción de forzas de diferente sentido. Figura 4.7: Forza negativa arredor do centro da placa e positiva no resto. Figura 4.8: Solución correspondente á forza negativa arredor do centro da placa e positiva no resto. Apéndice A Programas empregados Os programas que seguen foron modicados a partires dos códigos proporcionados en Joly [6]. Programas para o exemplo de validación %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% An Introduction to Scientific Computing %%%%%%% %%%%%%% I. Danaila, P. Joly, S. M. Kaber & M. Postel %%%%%%% %%%%%%% Springer, 2005 %%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% Modificado por María del Carmen Mosquera Vázquez %%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% Matlab Solution of exercise 1 - project 7 %% ELAS: elastic deformation of a thin plate %% Solution of the plate problem (linear equation) %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% clear all; close all; % % 1) Construction of the linear system % % number of points nx=80;ny=80; 45 Bibliografía [1] Brezis, H., Functional Analysis, Sobolev Spaces and Partial Dierential Equations , 1st ed., Springer, New York, 2011. [2] Ciarlet, P.G. e Destuynder, Ph., A Justication of the Two-Dimensional Linear Plate Model, Part 1: Derivation of the Two-Dimensional Model from the Three-Dimensional Model , ICES REPORT 77-09, The University of Texas at Austin, 1977 [3] Destuyder, Ph., Mathematical analysis of thin plate models , 1st ed., Springer, Berlin, 1996. [4] Evans, L.C., Partial Dierential Equations , 2nd ed, American Mathematical Soc., Berkeley, 2010. [5] Hernandez-Cifre, M. A. e Pastor González, J. A., Un curso de geometría diferencial , 47, Consejo Superior de Investigaciones Cientícas, Madrid, 2010. [6] Danaila, I., Joly, P., Kaber, S.M., Postel, M., An Introduction to Scientic Computing: Twelve Computational Projects Solved with MATLAB , Springer, New York, 2007. [7] Lui, S. H., Numerical Analysis of Partial Dierential Equations , 1st ed., Wiley, Hoboken, New Jersey, 2011. [8] C. Mardare, Méthodes mathématiques en élasticité, https://www.ljll.math.upmc. fr/MathModel/enseignement/polycopies/mardare.pdf . 53