Análisis matemático y simulación numérica con métodos puramente Lagrangianos y semi-Lagrangianos de problemas de la mecánica de medios continuos
Abstract
El objetivo final de la tesis es desarrollar herramientas propias e independientes, suficientemente robustas y fiables, que contribuyan a la simulación numérica de modelos de la mecánica de medios continuos. En particular, se considerarán problemas de fluidos e interacción fluido-estructura. La primera parte de la tesis tiene como finalidad la resolución numérica de las ecuaciones de Navier-Stokes (3D) y la implementación de los algoritmos mediante programas de ordenador. Dado el carácter evolutivo del problema que se quiere analizar, la importancia de la convección en el mismo y la inclusión de fronteras libres que provocan deformaciones en el dominio computacional , se propone desarrollar métodos de características puramente Lagrangianos y semi-Lagrangianos en combinación con métodos de elementos finitos, para resolver numéricamente el problema. Puesto que uno de los objetivos de este proyecto es el estudio de la evolución de un fluido para tiempos grandes, es necesario incorporar métodos de características de orden dos, pero manteniendo la estabilidad. Con ello se pretende obtener métodos de orden dos en tiempo y en espacio y, por tanto, la posibilidad de considerar pasos de tiempo mayores con objeto de acortar el tiempo de cálculo necesario para la simulación.
Full text
TESIS DE DOCTORADO ANÁLISIS MATEMÁTICO Y SIMULACIÓN NUMÉRICA CON MÉTODOS PURAMENTE LAGRANGIANOS Y SEMI-LAGRANGIANOS DE PROBLEMAS DE LA MECÁNICA DE MEDIOS CONTINUOS Pedro Fontán Muiños ESCUELA DE DOCTORADO INTERNACIONAL DE LA UNIVERSIDAD DE SANTIAGO DE COMPOSTELA PROGRAMA DE DOCTORADO EN MÉTODOS MATEMÁTICOS Y SIMULACIÓN NUMÉRICA EN INGENIERÍA Y CIENCIAS APLICADAS SANTIAGO DE COMPOSTELA 2021
DECLARACIÓN DO AUTOR/A DA TESE D./Dna. Pedro Fontán Muiños Título da tese: ANÁLISIS MATEMÁTICO Y SIMULACIÓN NUMÉRICA CON MÉTODOS PURAMENTE LAGRANGIANOS Y SEMI-LAGRANGIANOS DE PROBLEMAS DE LA MECÁNICA DE MEDIOS CONTINUOS Presento a miña tese, seguindo o procedemento axeitado ao Regulamento, e declaro que: 1) A tese abarca os resultados da elaboración do meu traballo. 2) De ser o caso, na tese faise referencia ás colaboracións que tivo este traballo. 3) Confirmo que a tese non incorre en ningún tipo de plaxio doutros autores nin de traballos presentados por min para a obtención doutros títulos. 4) A tese é a versión definitiva presentada para a súa defensa e coincide a versión impresa coa presentada en formato electrónico E comprométome a presentar o Compromiso Documental de Supervisión no caso de que o orixinal non estea na Escola. En Santiago, 02 de August de 2021. Sinatura electrónica
AUTORIZACIÓN DO DIRECTOR/TITOR DA TESE D./Dna. Marta Benítez García En condición de: Co-directora Título da tese: ANÁLISIS MATEMÁTICO Y SIMULACIÓN NUMÉRICA CON MÉTODOS PURAMENTE LAGRANGIANOS Y SEMI-LAGRANGIANOS DE PROBLEMAS DE LA MECÁNICA DE MEDIOS CONTINUOS INFORMA: Que a presente tese, correspóndese co traballo realizado por D/Dna Pedro Fontán Muiños, baixo a miña dirección/titorización, e autorizo a súa presentación, considerando que reúne os requisitos esixidos no Regulamento de Estudos de Doutoramento da USC, e que como director/titor desta non incorre nas causas de abstención establecidas na Lei 40/2015. En A Coruña a 3 de Agosto de 2021 Sinatura electrónica
AUTORIZACIÓN DO DIRECTOR/TITOR DA TESE D./Dna. Alfredo Bermúdez de Castro López-Varela En condición de: Co-director Título da tese: ANÁLISIS MATEMÁTICO Y SIMULACIÓN NUMÉRICA CON MÉTODOS PURAMENTE LAGRANGIANOS Y SEMI-LAGRANGIANOS DE PROBLEMAS DE LA MECÁNICA DE MEDIOS CONTINUOS INFORMA: Que a presente tese, correspóndese co traballo realizado por D/Dna Pedro Fontán Muiños, baixo a miña dirección/titorización, e autorizo a súa presentación, considerando que reúne os requisitos esixidos no Regulamento de Estudos de Doutoramento da USC, e que como director/titor desta non incorre nas causas de abstención establecidas na Lei 40/2015. En Santiago de Compostela a 3 de agosto de 2021 Sinatura electrónica
Agradecimientos En primer lugar quisiera agradecer a mis directores de tesis, a la profesora Marta Benítez García y al profesor Alfredo Bermúdez de Castro López-Varela , su tiempo, dedicación, trabajo y su infinita paciencia. Ha sido un privilegio poder trabajar bajo su tutela. También me gustaría agradecer a la profesora Pilar Salgado y Branca García su ayuda y colaboración en el planteamiento y resolución de algunos de los problemas industriales contenidos en este trabajo. Por último me gustaría mencionar a los profesores, colegas y compañeros que han participado en mi etapa formativa y profesional hasta la fecha. Las enseñanzas y experiencias acumuladas en este periodo son los elementos habilitantes que han permitido afrontar este trabajo. Muchas gracias a todos. v
vi
Índice general 1. Introducción 1 1.1. Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2. Estructura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 I Métodos no Eulerianos para problemas de frontera libre 5 Introducción 7 2. Fundamentos de la mecánica de los medios continuos 9 3. Problema de Cauchy escalar 15 3.1. Formulaciones Eulerianas y no Eulerianas del problema . . . . . . . . . . . . . . . . . . . 16 3.1.1. Formulación Euleriana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.2. Formulación general no Euleriana . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1.3. Formulación Lagrangiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2. Resolución numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.1. Discretización temporal: métodos de características . . . . . . . . . . . . . . . . . . 19 4. Problema de Cauchy vectorial 21 4.1. Formulaciones Eulerianas y no Eulerianas del problema . . . . . . . . . . . . . . . . . . . . 21 4.1.1. Formulación Euleriana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1.2. Formulación general no Euleriana . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1.3. Formulación Lagrangiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2. Resolución numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5. Remallado y reinicialización 27 5.1. Remallado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 vii
4Capítulo 1. Introducción problemas de transporte escalar y vectorial, describiendo las estrategias necesarias para su transformación desde una configuración Euleriana a otra no Euleriana general. Para cada uno de los problemas se indicará un caso particular de configuración no Euleriana, la Lagrangiana, por ser una configuración singular de aplicación general. Por último en el Capítulo 5 se abordará el problema del remallado y la reinicialización de las variables del problema en las situaciones en las que la malla del dominio computacional sufre deformaciones que producen elementos degenerados. En la segunda parte de la tesis se aplicarán las técnicas desarrolladas en la primera parte a dos problemas particulares de la mecánica de medios continuos. El primer problema abordado es el de las ecuaciones de Navier-Stokes para un fluido Newtoniano incompresible (Capítulo 6). En este capítulo se plantean los problemas que resultan de la transformación del problema original definido en una configuración Euleriana, a configuraciones no Eulerianas. Este problema es un caso particular de problema de Cauchy vectorial tratado en el Capítulo 4. Tras su transformación a una configuración no Euleriana, se proponen diferentes estrategias de discretización temporal empleando las fórmulas de Newmark. Según la elección de los valores de los parámetros que las definen se pueden obtener métodos de diferente orden. Los diferentes métodos considerados y los programas informáticos desarrollados son evaluados y validados utilizando ejemplos test numéricos habituales en dos dimensiones. Los códigos computacionales se desarrollan a partir de bibliotecas de propósito general. El contenido de este capítulo reproduce los resultados publicados en [25] y [26]. En el Capítulo 7 se plantea un problema derivado de la industria de conformado de piezas metálicas mediante la técnica de electrorrecalcado. Este problema involucra la resolución de diferentes fenómenos físicos. Requiere de tres modelos: mecánico para un sólido visco-elasto-plástico, electromagnético y térmico. Los tres problemas presentan un acople fuerte entre ellos: la mayor parte de los parámetros que definen el comportamiento del material son dependientes de la temperatura, el término fuente de la ecuación de transferencia de calor proviene del efecto Joule originado por las corrientes eléctricas que recorren la pieza y los modelos implicados dependen de la geometría de la pieza, que varía en el tiempo por la deformación que sufre la pieza. El problema completo se aborda en este capítulo incluyendo algunas simplificaciones. En concreto, el problema electromagnético se transforma en un problema eléctrico considerando que la corriente que se emplea es corriente continua. De este modo este problema pasa a ser un caso particular del problema de Cauchy escalar considerado en el Capítulo 3. El problema de transporte de calor también es un caso particular de este problema general de Cauchy escalar. En el caso del problema mecánico, éste es un caso particular del problema general de Cauchy vectorial considerado en el Capítulo 4 con una ley constitutiva para un sólido visco-elasto-plástico. Se obtiene una descripción del problema completo en una configuración Lagrangiana y se plantea su resolución de forma acoplada. Para la discretización temporal, en este caso, se opta por emplear integradores de tipo Runge-Kutta implícitos con mecanismos para el uso de pasos de tiempo variables. Esta característica es determinante para la resolución del problema, puesto que durante la simulación completa del proceso existen periodos en los que las variables propias del problema presentan cambios muy abruptos y se requieren tamaños de pasos temporales muy pequeños. Sin embargo, en otros periodos, los cambios son suaves y se pueden usar tamaños de pasos de tiempo más grandes, reduciendo de este modo el tiempo de cálculo necesario para la simulación.
Parte I Métodos no Eulerianos para problemas de frontera libre 5
Introducción En esta parte se describen las estrategias para la transformación de problemas generales, definidos en una configuración Euleriana, a una configuración de referencia general no Euleriana. Los problemas generales abordados serán dos problemas de transporte, uno escalar y otro vectorial, definidos dentro del marco de la teoría de medios continuos. Para obtener estas transformaciones se introducen en el Capítulo 2 los conceptos de la mecánica de medios continuos que posteriormente se emplearán en el tratamiento de los problemas generales. En los Capítulos 3 y 4 se abordarán los problemas de Cauchy escalar y vectorial, respectivamente. Finalmente, en el Capítulo 5, se describirá el problema de remallado y reinicialización que surge del empleo de métodos definidos en configuraciones no Eulerianas. Se propondrán estrategias para abordarlos. Estas estrategias permiten afrontar los problemas de frontera móvil de una manera natural. Los elementos desarrollados en esta parte serán empleados en los problemas planteados en la segunda parte de la tesis. 7
8
Capítulo 2 Fundamentos de la mecánica de los medios continuos El objetivo de este capítulo es recordar los elementos fundamentales de la mecánica de los medios continuos. Para ello se han seguido casi siempre las notaciones del libro de Gurtin [27]. Sea Ωun dominio acotado en Rd(d= 2,3) con frontera Lipschitziana Γdividida en dos partes: Γ=ΓD∪ΓR, con ΓD∩ΓR=∅. Sea X:¯ Ω×R→Rdun movimiento en el sentido de Gurtin [27], es decir X∈C3(¯ Ω×R)y para cada t∈R,X(·, t)es una función inyectiva tal que det F(p, t)>0,∀(p, t)∈¯ Ω×R,(2.1) donde F(·, t)es el gradiente del movimiento X(·, t). Nótese que por la compacidad de ¯ Ωy la continuidad de X(·, t)se deduce que X(¯ Ω, t)es una región cerrada para todo t. En general, esta condición no se verifica y debe ser añadida a las anteriores. Para cualquier A ⊂ ¯ Ω, se define At:= X(A, t). En la práctica se considera un intervalo temporal acotado, t∈[t0, tf], siendo t0,tfdos números no negativos. Se asume que X(p, t0) = p∀p∈¯ Ω, es más, en muchos casos se asumirá que el cuerpo se encuentra en reposo hasta el tiempo inicial (X(p, t) = p∀t≤ t0∀p∈¯ Ω) y en esta situación la velocidad y la aceleración son nulas para t < t0. La trayectoria del movimiento se define como T:= {(x, t) : x∈¯ Ωt, t ∈[t0, tf]}.(2.2) Para cada t, se denota por P(·, t) : ¯ Ωt→¯ Ωa la inversa de la aplicación X(·, t), es decir, X(P(x, t), t) = x,P(X(p, t), t) = p,∀(x, t)∈ T ,∀(p, t)∈¯ Ω×[t0, tf].(2.3) Esta aplicación P:T → ¯ Ω, se denomina aplicación de referencia del movimiento X. Para un (x, t)∈ T dado, P(x, t)es la posición inicial del punto que ocupa la posición xen el instante t. Dado un instante t0≤τ≤tf, se denota por Xτ:¯ Ωτ×[t0, tf]→Rdal movimiento relativo a la configuración en el instante τ: Xτ(y, t) := X(P(y, τ), t),∀(y, t)∈¯ Ωτ×[t0, tf].(2.4) Por tanto, Xτ(y, t)es la posición en el instante tdel punto que ocupa la posición yen el instate τ. Nótese que se denotan por plos puntos en ¯ Ω, por ylos puntos en ¯ Ωτy por xlos puntos en ¯ Ωt. 9
10 Capítulo 2. Fundamentos de la mecánica de los medios continuos Se pueden definir los campos asociados con la configuración ¯ Ωτ, por ejemplo Fτ:= grad Xτ,en ¯ Ωτ×[t0, tf].(2.5) A las funciones con dominio Tse les llama campos espaciales oEulerianos mientras que las funciones definidas en ¯ Ω×[t0, tf]se les llama campos materiales oLagrangianos. Un campo espacial Ψtambien se puede definir en ¯ Ωτ×[t0, tf]: Ψτ(y, t) := Ψ(Xτ(y, t), t),∀(y, t)∈¯ Ωτ×[t0, tf].(2.6) Para τ=t0,Ψt0es la descripción material oLagrangiana de Ψque se denota también como Ψm, es decir, Ψm(p, t) := Ψ(X(p, t), t),∀(p, t)∈¯ Ω×[t0, tf].(2.7) Estas funciones se representan en la Figura 2.1. p y x X(·, t) X(·, τ) P(·, τ) P(·, t) Xτ(·, t) ¯ Ω ¯ Ωτ¯ Ωt R RdLin Ψ(·, t) Ψτ(y, t) := Ψ(Xτ(y, t), t)Ψm(p, t) := Ψ(X(p, t), t) Figura 2.1: Transformaciones asociadas a las distintas configuraciones. A continuación se introducen algunas definiciones y notaciones sobre campos espaciales y materiales que se emplearán a lo largo del documento. •Sea Φun campo material escalar, vectorial o tensorial regular. Se denota por ˙ Φa la derivada parcial con respecto al segundo argumento (t). Si Φes un campo material escalar o vectorial regular, ∇Φ denota el gradiente con respecto al primer argumento (p). El símbolo Div denota el operador divergencia para campos materiales regulares que toman valores en Rd, es decir, Div Φ = tr (∇Φ). •Sea Ψun campo espacial escalar, vectorial o tensorial regular. Por Ψ0se indica la derivada parcial con respecto al segundo argumento (t), es decir, Ψ0(x, t) := ∂Ψ ∂t (x, t). Si Ψes un campo espacial escalar o vectorial regular, se denota por grad Ψ el gradiente con respecto al primer argumento (x). El símbolo div denota el operador divergencia para campos espaciales regulares que toman valores en Rd, es decir, div Ψ = tr ( grad Ψ).
11 •Para campos escalares regulares definidos en ¯ Ωτ×[t0, tf]se considera la notación empleada para campos espaciales, concretamente, grad Ψτ=∂Ψτ ∂y,div Ψτ= tr ( grad Ψτ),Ψ0 τ:= ∂Ψτ ∂t .(2.8) •En algunas situaciones donde están involucrados distintos operadores previamente introducidos, se puede especificar la variable de diferenciación con objeto de clarificar las expresiones empleadas, por ejemplo, ∇pΦ,grad xΨ,grad yΨτ,Div pΦ,div xΨ,div yΨτ. De las notaciones anteriores, empleando la regla de la cadena se pueden obtener las siguientes igualdades: grad yΨτ= ( grad xΨ)τFτ⇒( div xΨ)τ= tr ( grad yΨτF−1 τ) = grad yΨτF−1 τ·I=F−T τ·grad yΨτ, (2.9) para cualquier campo vectorial espacial regular Ψque toma valores en Rd. •Si Ψes un campo espacial escalar, vectorial o tensorial regular, ˙ Ψdenota la derivada material de Ψcon respecto al tiempo, es decir: ˙ Ψ(x, t) = ∂ ∂t(Ψ(X(p, t), t))|p=P(x,t),∀(x, t)∈ T .(2.10) En particular, si Ψes un campo escalar espacial regular se puede usar la regla de la cadena para obtener ˙ Ψ=Ψ0+ grad Ψ ·v,en T.(2.11) Análogamente, si Ψes un campo vectorial espacial regular se verifica ˙ Ψ=Ψ0+ ( grad Ψ)v,en T.(2.12) Se definen uel desplazamiento material, vla velocidad espacial y ala aceleración espacial, de la manera siguiente, u(p, t) = X(p, t)−p,v(x, t) = ˙ X(p, t)|p=P(x,t),a(x, t) = ¨ X(p, t)|p=P(x,t),(2.13) para (p, t)∈¯ Ω×[t0, tf]y(x, t)∈ T . Con las notaciones anteriores se obtienen las siguientes igualdades vτ=X0 τ,aτ=v0 τ=X00 τ,en ¯ Ωτ×[t0, tf],(2.14) vm=˙ X=˙ u,am=˙ vm=¨ X=¨ u,en ¯ Ω×[t0, tf].(2.15) Además, uτdenota el campo de desplazamientos relativo a la configuración en el instante τ: uτ(y, t) := Xτ(y, t)−y,∀(y, t)∈¯ Ωτ×[t0, tf].(2.16) Nótese que se verifican las siguientes igualdades vτ=u0 τ,aτ=v0 τ=u00 τ.(2.17) Se indicarán ahora algunas igualdades y lemas que se emplearán en capítulos posteriores para obtener formulaciones de los problemas considerados en distintas configuraciones.
12 Capítulo 2. Fundamentos de la mecánica de los medios continuos •Se tiene la siguiente relación entre el vector unitario normal a la frontera del dominio deforado en el instante t,Γt, y dirigido hacia el exterior de este, que será denotado por n(·, t), y el vector unitario normal a la frontera del dominio deformado en el instante τ,Γτ, y dirigido hacia el exterior de este, que será denotado por mτ: n(x, t) = F−T τ(y, t)mτ(y) |F−T τ(y, t)mτ(y)|y=Xt(x,τ) ,(2.18) para x∈Γtyt∈[t0, tf]. Nótese que por (2.4) se tiene que Xt(x, τ) = X(P(x, t), τ). En particular, tomando τ=t0en (2.18), se tiene: n(x, t) = F−T(p, t)mt0(p) |F−T(p, t)mt0(p)|p=P(x,t) ,(2.19) para x∈Γtyt∈[t0, tf]. •Se introducen los siguientes lemas: Lema 2.0.1. Dado un campo escalar espacial suficientemente regular ψ, se verifica, ZΩt div x[A(x, t) grad xψ(x, t)] dx =ZΩτ div ydet Fτ(y, t)F−1 τ(y, t)Aτ(y, t)F−T τ(y, t) grad yψτ(y, t)dy,(2.20) siendo Aun campo tensorial espacial con suficiente regularidad. En particular, si τ=t0, ZΩt div x[A(x, t) grad xψ(x, t)] dx =ZΩ Div det F(p, t)F−1(p, t)Am(p, t)F−T(p, t)∇pψm(p, t)dp.(2.21) Demostración. En primer lugar, aplicando el teorema de la divergencia, se tiene ZΩt div x[A(x, t) grad xψ(x, t)] dx=ZΓt A(x, t) grad xψ(x, t)·n(x, t)dAx. Aplicando ahora el cambio de variable x=Xτ(y, t)y usando la regla de la cadena y el teorema de la divergencia se deduce el resultado. En efecto, ZΓt A(x, t) grad xψ(x, t)·n(x, t)dAx =ZΓτ det Fτ(y, t)Aτ(y, t) grad xψ(Xτ(y, t), t)·F−T τ(y, t)mτ(y)dAy =ZΓτ det Fτ(y, t)F−1 τ(y, t)Aτ(y, t)F−T τ(y, t) grad yψτ(y, t)·mτ(y)dAy =ZΩτ div ydet Fτ(y, t)F−1 τ(y, t)Aτ(y, t)F−T τ(y, t) grad yψτ(y, t)dy.
13 Lema 2.0.2. Dado un campo escalar espacial suficientemente regular ψ, se verifica, div x[A(x, t) grad xψ(x, t)] = div ydet Fτ(y, t)F−1 τ(y, t)Aτ(y, t)F−T τ(y, t) grad yψτ(y, t)1 det Fτ(y, t)y=Xt(x,τ)=X(P(x,t),τ) , (2.22) para x∈Ωtyt∈[t0, tf], siendo Aun campo tensorial espacial con suficiente regularidad. En particular, si τ=t0, div x[A(x, t) grad xψ(x, t)] = Div det F(p, t)F−1(p, t)Am(p, t)F−T(p, t)∇ψm(p, t)1 det F(p, t)p=P(x,t) .(2.23) Demostración. Procediendo de forma análoga a la demostración del resultado anterior, se deduce ZB(x0,λ) div x[A(x, t) grad xψ(x, t)] dx =ZXt(B(x0,λ),τ) div ydet Fτ(y, t)F−1 τ(y, t)Aτ(y, t)F−T τ(y, t) grad yψτ(y, t)dy, para x0∈Ωt,t∈[t0, tf]yλ > 0. Recordemos que por definición Xt(B(x0, λ), τ) = X(P(B(x0, λ), t), τ). Aplicando el cambio de variable y=Xt(x, τ) = X(P(x, t), τ), resulta ZB(x0,λ) div x[A(x, t) grad xψ(x, t)] dx =ZB(x0,λ) div ydet Fτ(Xt(x, τ), t)F−1 τ(Xt(x, τ), t)A(x, t)F−T τ(Xt(x, τ), t) grad yψτ(Xt(x, τ), t)] 1 det Fτ(Xt(x, τ), t)dx, para x0∈Ωt,t∈[t0, tf]yλ > 0. Por último, utilizando el teorema de localización (ver [27] página 5) se obtiene el resultado. Lema 2.0.3. Dado un campo tensorial espacial suficientemente regular H, se verifica, ZΩt div H(x, t)dx=ZΩτ div ydet Fτ(y, t)Hτ(y, t)F−T τ(y, t)dy,(2.24) div H(x, t) = div ydet Fτ(y, t)Hτ(y, t)F−T τ(y, t)1 det Fτ(y, t)y=Xt(x,τ)=X(P(x,t),τ) ,(2.25) para x∈Ωtyt∈[t0, tf]. Demostración. El procedimiento para obtener estas igualdades es análogo al empleado para deducir (2.20) y (2.22).
20 Capítulo 3. Problema de Cauchy escalar instante en un dominio computacional conocido y, dependiendo del esquema de integración temporal empleado para resolver el problema (3.12), utilizar estrategias de alto orden para mejorar la precisión y permitir un aumento del tamaño del paso de tiempo. Conocido el movimiento Xτse puede obtener el gradiente de deformaciones a través de la igualdad Fτ= grad yXτo planteando el problema de valor inicial que verifica Fτ: F0 τ(y, t) = L(Xτ(y, t), t)Fτ(y, t) = Lτ(y, t)Fτ(y, t), Fτ(y, τ) = I,(3.13) donde L(x, t)es el campo tensorial espacial grad v(x, t). Para la resolución numérica de los problemas de valor inicial (3.10), (3.12) y (3.13) se pueden utilizar distintos métodos como colocación, multipaso lineales o Runge-Kutta. En [24] se emplea un método de Runge-Kutta explícito de orden dos para obtener el movimiento y el gradiente de deformaciones.
Capítulo 4 Problema de Cauchy vectorial Tal vez el ejemplo más importante de problema de Cauchy vectorial es el problema general de movimiento de un medio continuo. Para caracterizar un medio concreto será necesario relacionar el tensor de tensiones con ciertas magnitudes cinemáticas mediante leyes constitutivas. Las ecuaciones de Navier-Stokes incompresibles y de Navier-Lamé son casos particulares para fluidos Newtonianos incompresibles y sólidos elásticos lineales, respectivamente. En [32] se utilizan estos modelos para tratar problemas de interacción fluido-estructura. Partiendo del problema de Cauchy vectorial para el movimiento en un dominio con diferentes tipos de materiales, particularizándolo en cada dominio mediante la ley constitutiva correspondiente y planteando los subproblemas en una configuración no Euleriana con el campo de desplazamientos como incógnita, es posible definir un problema acoplado de una manera natural, que maneja en todos los subdominios el mismo tipo de ecuaciones. Al igual que ocurría en el caso escalar, el tratamiento del término convectivo y las fronteras libres presentan dificultades cuando se emplea una configuración Euleriana del problema vectorial. Sin embargo, si se plantean en configuraciones no Eulerianas, estos elementos se abordan de manera natural, facilitando su tratamiento. En este capítulo se desarrollarán los cálculos para la transformación del problema de Cauchy para la ecuación general del movimiento de un medio continuo, desde configuraciones Eulerianas a configuraciones no Eulerianas, sin concretar ninguna ley constitutiva. Más adelante se utilizarán los desarrollos de este capítulo para plantear problemas concretos que involucran el movimiento acoplado con otros fenómenos. En cada caso se definirán las leyes constitutivas necesarias para el correcto tratamiento de los fenómenos a modelar. 4.1. Formulaciones Eulerianas y no Eulerianas del problema En esta sección se introducirá la forma general Euleriana de la ecuación del movimiento, sin tener en cuenta leyes constitutivas que definan el tensor de tensiones, que será denotado por T. Posteriormente se emplearán diferentes equivalencias para transformar el problema Euleriano en una familia de problemas no Eulerianos. Para estos problemas definidos en configuraciones no Eulerianas, se definirán también los 21
22 Capítulo 4. Problema de Cauchy vectorial problemas variacionales asociados. El caso puramente Lagrangiano se tratará como un caso particular de la formulación no Euleriana general. 4.1.1. Formulación Euleriana Se considera el siguiente problema de valor inicial-frontera: Problema 4.1.1. Problema fuerte vectorial en configuración Euleriana. Encontrar dos funciones v:T → RdyT:T → Sym tales que ρ(x, t)˙ v(x, t) = div T(x, t) + b(x, t),(4.1a) para x∈Ωtyt∈[t0, tf], sujetas a las condiciones de contorno v(x, t) = vD(x, t),∀x∈ΓD t,(4.1b) T(x, t)n(x, t) = h(x, t),∀x∈ΓR t,(4.1c) para t∈[t0, tf], y a la condición inicial v(x, t0) = v0(x)∀x∈¯ Ωt0.(4.1d) En estas ecuaciones T:T −→ Sym representa el tensor de tensiones de Cauchy siendo Sym el espacio de los tensores simétricos, ρ:T −→ R+es la densidad de masa, b:T −→ V es el término fuente (una densidad volumétrica de fuerza), h(·, t) : ΓR t→Rdes la densidad superficial de fuerza en la frontera ΓR t donde se da una condición de contorno de tipo Neumann, vD(·, t):ΓD t→Rdes la condición Dirichlet impuesta sobre la velocidad y v0:¯ Ωt0→Rdes la velocidad inicial. Para completar la definición del problema es necesario definir el tensor de tensiones de Cauchy T mediante una ley constitutiva. Estas leyes constitutivas pueden añadir nuevas ecuaciones e incógnitas que serán incorporadas al problema original. En los Capítulos 6 y 7 se introducirán las leyes constitutivas consideradas para varios tipos de materiales: sólido elástico e isótropo, fluido Newtoniano y sólido viscoelasto-plástico. 4.1.2. Formulación general no Euleriana A continuación se desarrollarán las expresiones necesarias para escribir el Problema 4.1.1 en una configuración no Euleriana Ωτ×[t0, tf], con τ∈[t0, tf]. A partir de la definición de derivada material (2.10), usando la regla de la cadena y teniendo en cuenta las relaciones entre la velocidad y el movimiento dadas en (2.13) y (2.14), se obtiene la siguiente igualdad: ˙ v(x, t) = ∂vτ ∂t (y, t)|y=Xt(x,τ)≡v0 τ(y, t)|y=Xt(x,τ),∀(x, t)∈ T .(4.2) Se evalúa la ecuación (4.1a) en el punto x=Xτ(y, t)y se utiliza la ecuación (4.2) para obtener
4.1. Formulaciones Eulerianas y no Eulerianas del problema 23 ρ(Xτ(y, t), t)v0 τ(y, t) = div xT(Xτ(y, t), t) + b(Xτ(y, t), t),(4.3) para (y, t)∈Ωτ×[t0, tf]. Utilizando ahora la igualdad (2.25), con H=T, se deduce la ecuación del movimiento en formulación no Euleriana: ρτ(y, t)v0 τ(y, t) = div ydet Fτ(y, t)Tτ(y, t)F−T τ(y, t)1 det Fτ(y, t)+bτ(y, t),(4.4) para (y, t)∈Ωτ×[t0, tf]. Se realiza el cambio de variable x=Xτ(y, t)en (4.1b)-(4.1d) para conseguir las condiciones de contorno e iniciales en la configuración de referencia intermedia, vτ= (vD)τen ΓD τ×[t0, tf],(4.5) TτF−T τmτ=|F−T τmτ|hτen ΓR τ×[t0, tf],(4.6) vτ(·, t0) = v0(P(·, τ)) en ¯ Ωτ,(4.7) donde se ha utilizado (2.18). Entonces, a partir de estos resultados se deduce la siguiente formulación no Euleriana del Problema 4.1.1. Problema 4.1.2. Problema fuerte vectorial en configuración no Euleriana, Ωτ×[t0, tf]. Encontrar dos funciones vτ:¯ Ωτ×[t0, tf]→RdyTτ:¯ Ωτ×[t0, tf]→Sym tales que ρτv0 τ= div ydet FτTτF−T τ1 det Fτ +bτ.(4.8a) en Ωτ×[t0, tf], sujetas a las condiciones de contorno vτ= (vD)τ,en ΓD τ×[t0, tf],(4.8b) TτF−T τmτ=|F−T τmτ|hτ,en ΓR τ×[t0, tf],(4.8c) y a la condición inicial vτ(·, t0) = v0(P(·, τ)) en ¯ Ωτ.(4.8d) Para obtener una formulación débil del Problema 4.1.2 basta con multiplicar (4.8a) por det Fτy por una función test z∈H1 ΓD τ(Ωτ), integrar en Ωτy aplicar la fórmula de Green adecuada junto con la condición de contorno (4.8c). De esta forma se deduce, ZΩτ det Fτρτv0 τ·zdy=ZΓR τ det Fτ|F−T τmτ|hτ·zdAy−ZΩτ det FτTτF−T τ·grad yzdy +ZΩτ det Fτbτ·zdy,(4.9) para todo z∈H1 ΓD τ(Ωτ).
24 Capítulo 4. Problema de Cauchy vectorial 4.1.3. Formulación Lagrangiana A partir del Problema 4.1.2, particularizándolo para el instante τ=t0, se obtiene la formulación fuerte del Problema 4.1.1 en la configuración de referencia Ω=Ωt0. En concreto, Problema 4.1.3. Problema fuerte vectorial en configuración Lagrangiana. Encontrar dos funciones vm:¯ Ω×[t0, tf]→RdyTm:¯ Ω×[t0, tf]→Sym tales que ρ0˙ vm= Div det FTmF−T+ det Fbm,(4.10a) en Ω×[t0, tf], sujetas a las condiciones de contorno vm= (vD)m,en ΓD×[t0, tf],(4.10b) TmF−Tmt0=|F−Tmt0|hm,en ΓR×[t0, tf],(4.10c) y a la condición inicial vm(·, t0) = v0en ¯ Ω.(4.10d) Nótese que se ha utilizado el principio de conservación de la masa en la forma siguiente (véase [27]) ρm(p, t) det F(p, t) = ρ0(p),(4.11) donde ρ0es la densidad en la configuración de referencia. Tomando τ=t0en (4.9) y utilizando el principio de conservación de la masa (4.11) se obtiene la siguiente formulación débil del Problema 4.1.3: ZΩ ρ0˙ vm·zdp=ZΓR det F|F−Tmt0|hm·zdAp−ZΩ det FTmF−T· ∇zdp+ZΩ det Fbm·zdp,(4.12) para todo z∈H1 ΓD(Ω). 4.2. Resolución numérica Al contrario de lo que ocurría en el problema general escalar planteado en la Sección 3, ahora la velocidad no es un dato sino que es una incógnita del problema. De modo que el cálculo de Xτ, necesario para obtener Fτy realizar las transformaciones de los campos espaciales, no es independiente de la resolución del Problema 4.1.2.
4.2. Resolución numérica 25 En concreto, las relaciones entre vτ,Xτyuτdadas en (2.14) y (2.17) permiten escribir el Problema 4.1.2 de la siguiente manera: ρτv0 τ= div ydet FτTτF−T τ1 det Fτ +bτ,en Ωτ×[t0, tf],(4.13a) u0 τ=vτ,en ¯ Ωτ×[t0, tf],(4.13b) Xτ(y, t) = y+uτ(y, t),∀(y, t)∈¯ Ωτ×[t0, tf],(4.13c) Fτ= grad yXτ,en ¯ Ωτ×[t0, tf],(4.13d) vτ= (vD)τen ΓD τ×[t0, tf],(4.13e) TτF−T τmτ=|F−T τmτ|hτen ΓR τ×[t0, tf],(4.13f) uτ(·, τ) = 0,en ¯ Ωτ,(4.13g) vτ(·, t0) = v0(P(·, τ)),en ¯ Ωτ.(4.13h) Nótese que se puede plantear este problema tomando como incógnita el desplazamiento en lugar de la velocidad. En efecto, ρτu00 τ= div ydet FτTτF−T τ1 det Fτ +bτ,en Ωτ×[t0, tf],(4.14a) Xτ(y, t) = y+uτ(y, t),∀(y, t)∈¯ Ωτ×[t0, tf],(4.14b) Fτ= grad yXτ,en ¯ Ωτ×[t0, tf],(4.14c) vτ= (vD)τen ΓD τ×[t0, tf],(4.14d) TτF−T τmτ=|F−T τmτ|hτen ΓR τ×[t0, tf],(4.14e) uτ(·, τ) = 0,en ¯ Ωτ,(4.14f) vτ(·, t0) = v0(P(·, τ)),en ¯ Ωτ.(4.14g) De esta manera, tras la discretización espacial, se obtiene un problema de valor inicial de segundo orden. Estas formulaciones fuertes admiten, respectivamente, las siguientes formulaciones variacionales que son equivalentes a la dada en (4.9): ZΩτ det Fτρτv0 τ·zdy= ZΓR τ det Fτ|F−T τmτ|hτ·zdAy−ZΩτ det FτTτF−T τ·grad yzdy+ZΩτ det Fτbτ·zdy,∀z∈H1 ΓD τ(Ωτ), (4.15a) ZΩτ u0 τ·˜ udy=ZΩτ vτ·˜ udy,∀˜ u∈L2(Ωτ)(4.15b) y ZΩτ det Fτρτu00 τ·zdy= ZΓR τ det Fτ|F−T τmτ|hτ·zdAy−ZΩτ det FτTτF−T τ·grad yzdy+ZΩτ det Fτbτ·zdy,∀z∈H1 ΓD τ(Ωτ), (4.16a)
26 Capítulo 4. Problema de Cauchy vectorial con las correspondientes condiciones de contorno e iniciales, y definiciones dadas en (4.13e), (4.13g), (4.13h), (4.13c) y (4.13d). Nótese que, cuando la incógnita principal es el campo de desplazamientos en lugar de la velocidad, será necesario obtener las condiciones iniciales y de contorno para dicho campo a partir de las dadas para la velocidad. Tras considerar una discretización espacial de estas formulaciones débiles mediante el método de los elementos finitos, se obtiene un problema implícito de valores iniciales del siguiente tipo: f1(t, y,y0,y00) = f2(t, y)⇒f(t, y,y0,y00) = f1(t, y,y0,y00)−f2(t, y) = 0. Adicionalmente, las leyes constitutivas introducirán nuevas ecuaciones en forma de restricciones. Las leyes constitutivas que se consideran en este trabajo admiten la siguiente escritura: G(Tτ,Fτ,grad vτ,wτ) = 0,(4.17) con sus correspondientes incógnitas wτ. En ese caso se obtendrá un problema semidiscretizado en espacio de valores iniciales, definido por un sistema de ecuaciones diferenciales y algebraicas (DAE por sus siglas en inglés, Differential Algebraic Equations).
Capítulo 5 Remallado y reinicialización Los métodos no Eulerianos que se han descrito en las secciones precedentes permiten resolver el problema planteado en configuraciones espaciales diferentes de la del instante que se pretende resolver. En estos métodos habitualmente se fija una configuración de referencia y se mantiene hasta un determinado instante o por un periodo constante de tiempo. En todo caso, cuando se considera el momento adecuado para cambiar la configuración de referencia y reinicializar el problema, si ese instante no coincide con el instante final de la simulación, será necesario desarrollar métodos para trasladar la información de una configuración a otra con objeto de poder continuar con la simulación. En el caso de los métodos puramente Lagrangianos la decisión de cambiar de configuración de referencia está determinada por el valor del determinante del tensor de gradiente de deformaciones Fτ. Este tensor que está definido en (2.5) debe cumplir con la condición (2.1) para que el movimiento Xτsea válido. Durante la simulación puede ocurrir que en un instante determinado se viole la condición (2.1) y entonces será necesario definir una nueva configuración de referencia previa a ese instante. Pueden existir otros criterios más estrictos que lleven a tomar esta decisión. El cambio de configuración de referencia implica un remallado con complejidad variable, dependiendo de las características del problema que se está resolviendo. En la Sección 5.1 se incidirá en este tema. Cuando se procede a un cambio de configuración de referencia es necesario convertir los campos asociados al problema, a la nueva configuración. Esta acción se denominará reinicialización de los campos y se detallará en la Sección 5.2. 5.1. Remallado En esta sección se introduce el proceso de remallado distinguiendo entre el caso en el que el dominio del problema sea fijo o variable con el tiempo. En el segundo caso, la nueva configuración de referencia es habitualmente desconocida y se aproxima utilizando una aproximación del desplazamiento. Los campos aproximados se denotarán con el subíndice h. 5.1.1. Frontera fija Cuando la frontera del dominio es fija se puede reutilizar la malla original en cada nuevo reinicio del método pues todas las configuraciones coinciden: Ωt= Ωs∀s, t ∈[t0, tf]). 27
28 Capítulo 5. Remallado y reinicialización A pesar de esta circunstancia es posible que en algunos casos esté justificado el remallado del dominio: •Malla adaptativa, refinada en las zonas con cambios bruscos de los campos que son incógnita. •Reutilizar los vértices originales desplazados para evitar los errores de interpolación que se añaden en el proceso de reinicialización de los campos. Si se reutilizan los vértices desplazados, es conveniente llevar a cabo un cambio de conectividades para intentar mejorar la calidad de la malla. Esta mejora de la calidad de la malla a través de un nuevo esquema de conectividades, no siempre es posible. 5.1.2. Frontera móvil Si el problema involucra alguna frontera libre, el cambio de configuración de referencia implica una deformación de la frontera del dominio original. Para obtener dicha frontera deformada se aplica una aproximación del movimiento Xτ,h sobre la frontera original Γτ. Por simplicidad, se asume que Ωτ⊂R2y que tiene una frontera poligonal. Se puede definir una familia de particiones de ¯ Ωτdenotada por Tτ,h consistente en un conjunto de triángulos {Ki}ital que [ K∈Tτ,h K=¯ Ωτ,(5.1) |K| ≤ h, ∀K∈Tτ,h,(5.2) donde |K|representa el diámetro del triángulo K. Además, una arista es compartida, como máximo, por 2 triángulos. Los puntos de la triangulación se denominan vértices. Sean dos instantes temporales τi, τj∈[t0, tf]con τi< τj. Para obtener una aproximación del dominio deformado Ωτj=Xτi(Ωτi, τj), siendo conocidos el dominio Ωτio una aproximación de este y una aproximación de uτi(·, τj), se llevan a cabo las siguientes tareas: •Se construye una nueva triangulación de ¯ Ωτi, denotada por Tr τi,h, cuyos vértices son los nudos que están en las aristas asociados al campo de desplazamientos uτi,h(·, τj), que es una función continua en ¯ Ωτiy polinómica sobre cada elemento de la triangulación original Tτi,h. •Se aplica a cada vértice de Tr τi,h el movimiento correspondiente al desplazamiento aproximado uτi,h(·, τj)y se obtiene la malla deformada ˜ Tr τj,h y una aproximación del dominio ¯ Ωτjcomo se ilustra en la Figura 5.1. En concreto, si denotamos por {yh,τi l}ly{˜ yh,τj l}llos vértices de Tr τi,h y˜ Tr τj,h, respectivamente, se tiene ˜ yh,τj l=yh,τi l+uτi,h(yh,τi l, τj),∀l, (5.3) ¯ Ωτj∼¯ ˜ Ωτj:= [ ˜ K∈˜ Tr τj,h ˜ K, (5.4) donde los vértices de cada elemento de ˜ Tr τj,h son vértices desplazados de un elemento de Tr τi,h. •Utilizando como frontera aproximada de Ωτjla frontera de ˜ Ωτj, se construye una nueva triangulación de Ωτj. La triangulación ˜ Tr τj,h también es una triangulación de ¯ Ωτjpero, en general, la calidad de esta malla no es buena. Este procedimiento se ilustra en la Figura 5.2. En el caso tridimensional, el procedimiento es análogo.
5.1. Remallado 29 Figura 5.1: Malla original superpuesta a la malla deformada que se obtiene aplicando a cada vértice de la malla original el movimiento correspondiente al desplazamiento uti,h(·, tj).
36 constitutiva para un sólido visco-elasto-plástico. Se obtiene una descripción del problema completo en una configuración Lagrangiana y se plantea su resolución de forma acoplada. Por la naturaleza del problema, se considera un caso de simetría cilíndrica como el considerado en el ejemplo numérico, permitiendo su resolución en una configuración 2D. Para la discretización temporal, en este caso, se opta por emplear integradores de tipo Runge-Kutta implícitos con mecanismos para el uso de pasos de tiempo variables. Esta característica es determinante para la resolución del problema, puesto que durante la simulación completa del proceso existen periodos en los que los cambios de las variables propias del problema son muy abruptos y requieren tamaños de pasos temporales muy pequeños. En otras situaciones los cambios son suaves pudiendo usar tamaños de pasos de tiempo más grandes, reduciendo el número total de pasos de tiempo.
Capítulo 6 Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles* *Los resultados de este capítulo ya han sido publicados como: Benítez M.a, Bermúdez A.b,c, Fontán P.d(2018) A Second-Order Linear Newmark Method for Lagrangian Navier-Stokes Equations. In: Doubova A., González-Burgos M., Guillén-González F., Marín Beltrán M. (eds) Recent Advances in PDEs: Analysis, Numerics and Control. SEMA SIMAI Springer Series, vol 17. Springer, Cham. https://doi.org/10.1007/978-3-319-97613-6_3 Benítez, M.a, Bermúdez, A.b,c and Fontán, P.dNon-Eulerian Newmark Methods: A Powerful Tool for Free-Boundary Continuum Mechanics Problems. J Sci Comput 83, 44 (2020). https://doi.org/10.1007/s10915-020-01207-y aDepartamento de Matemáticas, Universidade da Coruña, Centro de Investigación CITIC, Elviña s/n, 15071 A Coruña, Spain bDepartamento de Matemática Aplicada and Instituto de Matemáticas, Universidade de Santiago de Compostela, c/ Lope Gómez de Marzoa s/n, 15782 Santiago de Compostela, Spain cInstituto Tecnológico de Matemática Industrial (ITMATI), Rúa Constantino Candeira s/n, 15782 Santiago de Compostela, Spain dREPSOL Technology Lab, Autovía de Extremadura s/n, 28935 Móstoles, Madrid, Spain El objetivo principal de este capítulo es desarrollar una familia de métodos precisos y estables para la resolución numérica de problemas con frontera libre en el ámbito de la mecánica de fluidos. Los métodos y resultados de este capítulo han sido publicados en [25, 26]. La familia de métodos indicada hace uso de formulaciones no Eulerianas y de métodos de Newmark para la discretización temporal. La resolución de problemas de la mecánica de fluidos incompresibles con fronteras móviles planteados en configuración Euleriana presenta dos principales dificultades: el tratamiento del término convectivo para altos números de Reynolds [33] y la modelización y seguimiento de la frontera móvil. En el caso del término convectivo, es posible emplear diferentes técnicas para conseguir una formulación estable [33]. Una de ellas es emplear el método de las características [6] de tal manera que el término convectivo se incluye dentro de la derivada material con respecto al tiempo, como se indica en los Capítulos 37
38 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. 3 y 4. La descripción Lagrangiana, típicamente empleada en mecánica de sólidos, proporciona un tratamiento natural de problemas con dominios dependientes del tiempo tal y como son los problemas de frontera libre e interacción fluido-estructura. Inicialmente, el problema principal al emplear los métodos Lagrangianos en problemas de mecánica de fluidos residía en la dificultad del proceso de remallado necesario en estos métodos cuando se enfrentan a grandes deformaciones (ver [34, 35]). Como consecuencia del desarrollo de algoritmos de mallado eficientes que permiten estrategias adaptativas automáticas, los métodos puramente Lagrangianos han adquirido un importante interés en el ámbito de la mecánica de fluidos (ver por ejemplo [36, 37]). En este capítulo se detallarán algunas técnicas para la resolución numérica de problemas en mecánica de fluidos con fronteras móviles mediante el uso de formulaciones no Eulerianas como las introducidas en los capítulos anteriores. En concreto, se considerará un fluido viscoso, Newtoniano e incompresible en un dominio dependiente del tiempo que puede presentar grandes deformaciones. Los cambios topológicos en las interfaces, como pueden ser roturas de la interfaz, no se contemplarán en esta tesis. Este tipo de problemas aparecen de forma recurrente en el ámbito de la ingeniería, por ejemplo en problemas de frontera libre y problemas de interacción fluido-estructura. En la primera parte del capítulo se recordarán las ecuaciones de Navier-Stokes en configuración Euleriana y, a continuación, se establecerá su formulación en configuración no Euleriana. Posteriormente se plantearán las discretizaciones temporales y espaciales mediante el método de Newmark y elementos finitos, respectivamente. Por último, se mostrarán resultados numéricos obtenidos mediante los distintos esquemas propuestos. 6.1. Formulaciones del problema Las ecuaciones de Navier-Stokes pertenecen a la familia de ecuaciones vectoriales introducidas en el Capítulo 4. A partir de estas ecuaciones y utilizando la ley constitutiva de los fluidos Newtonianos incompresibles se consigue caracterizar la ley del movimiento para este tipo de fluidos. A continuación se indicará la ley constitutiva aplicada. Para un desarrollo más profundo se puede consultar [38]. 6.1.1. Leyes constitutivas para un fluido Newtoniano incompresible Para fluidos Newtonianos incompresibles el tensor de tensiones de Cauchy es de la siguiente forma: T=−πI+µgrad v+ grad vT,en T,(6.1) donde π:T → Res la presión y µ:T → R+la viscosidad dinámica. Además, dado que el fluido es incompresible, es necesario añadir una restricción a la velocidad: div v= 0,en T.(6.2)
6.1. Formulaciones del problema 39 Recordemos cómo aparece la condición (6.2). Un cuerpo Ωes incompresible si solo admite movimientos isocóricos, es decir, movimientos que conservan el volumen. Dada una región material Pdenotamos por Pt=X(P, t), su posición en el instante t. El movimiento Xes isocórico si d dtVol(Pt) = 0,(6.3) donde Vol(Pt) = RPtdxdenota el volumen ocupado por los puntos materiales de Pen el instante t. El teorema del transporte de Reynolds, d dt ZPt φ(x, t)dx=ZPt (˙ φ(x, t) + φ(x, t) div v(x, t))dx,(6.4) permite evaluar la derivada temporal de Vol(Pt), simplemente considerando φ= 1 y tomando v(x, t)la velocidad espacial asociada al movimiento X: d dtVol(Pt) = d dt ZPt dx=ZPt div v(x, t)dx= 0.(6.5) Como la igualdad anterior es válida para cualquier parte Pde B, entonces, por el teorema de localización, se deduce que un movimiento será isocórico si y solo si div v= 0. Además, teniendo en cuenta que se verifica (ver [27] página 27) ( det F)·= det Ftr ˙ FF−1, y que ˙ F=LmFcon L= grad v, se tienen las siguientes equivalencias Xes isocórico ⇔div v= 0 ⇔( det F)·= 0.(6.6) Si un movimiento es isocórico, de la ecuación de conservación de la masa y de la ecuación (6.4) para φ=ρ, se deduce d dt ZPt ρ(x, t)dx=ZPt ( ˙ρ(x, t) + ρ(x, t) div v(x, t))dx=ZPt ˙ρ(x, t)dx= 0,(6.7) y, de nuevo por el teorema de localización, se deduce que ˙ρ(x, t)=0lo que significa que los puntos materiales conservan su densidad a lo largo del movimiento. 6.1.2. Formulación Euleriana A partir del Problema 4.1.1, utilizando la ley constitutiva del tensor de Cauchy (6.1) y la restricción asociada con los movimientos isocóricos se define el siguiente problema de valores iniciales y de contorno: Problema 6.1.1. Problema fuerte de Navier-Stokes en configuración Euleriana. Encontrar dos funciones v:T → Rdyπ:T → Rtales que ρ(x, t)˙ v(x, t) + grad π(x, t)−div {µ(x, t)( grad v(x, t) + ( grad v(x, t))T}=b(x, t),(6.8a) div v(x, t) = 0,(6.8b)
40 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. para x∈Ωtyt∈[t0, tf], sujetas a las siguientes condiciones de contorno, v(x, t) = vD(x, t),∀x∈ΓD t,(6.8c) −π(x, t)I+µ(x, t)( grad v(x, t) + ( grad v(x, t))T)n(x, t) = h(x, t),∀x∈ΓR t,(6.8d) para t∈[t0, tf], y a la condición inicial v(x, t0) = v0(x),∀x∈¯ Ωt0.(6.8e) 6.1.3. Formulación general no Euleriana Para poder escribir el problema en una configuración no Euleriana Ωτcon τ∈[t0, tf], se parte del problema de Cauchy vectorial (4.8a)-(4.8d) y de la expresión de la ley constitutiva en una configuración no Euleriana. En concreto, utilizando la primera igualdad de (2.9) con Ψ = v, se obtiene la formulación no Euleriana en Ωτ×[t0, tf]de la ley constitutiva dada en (6.1): Tτ=−πτI+µτgrad yvτF−1 τ+F−T τ( grad yvτ)T.(6.9) Utilizando las igualdades (2.25) y (2.9) se tiene, div xT(x, t) = −1 det Fτ(y, t)div yπτ(y, t) det Fτ(y, t)F−T τ(y, t)(6.10) +1 det Fτ(y, t)div yµτ(y, t)grad yvτ(y, t)F−1 τ(y, t) + F−T τ(y, t)( grad yvτ(y, t))Tdet Fτ(y, t)F−T τ(y, t), con y=Xt(x, τ) = X(P(x, t), τ),x∈Ωt,y∈Ωτ, donde X(p, t)es el movimiento asociado a la velocidad v(x, t)y verifica (2.13). De esta manera la ecuación (6.8a) se escribe en configuración no Euleriana Ωτ× [t0, tf]como sigue ρτ(y, t)v0 τ(y, t) + 1 det Fτ(y, t)div yπτ(y, t) det Fτ(y, t)F−T τ(y, t) −1 det Fτ(y, t)div yµτ(y, t)grad yvτ(y, t)F−1 τ(y, t) + F−T τ(y, t)( grad yvτ(y, t))Tdet Fτ(y, t)F−T τ(y, t)= bτ(y, t).(6.11) En el caso de la restricción asociada con la incompresibilidad (6.8b), utilizando (2.9) con Ψ = vse obtiene su forma no Euleriana, F−T τ(y, t)·grad yvτ(y, t) = 0.(6.12)
6.1. Formulaciones del problema 41 Nótese que en este trabajo se asume que X(p, t0) = p∀p∈¯ Ω, y por tanto det F(p, t0) = 1 ∀p∈¯ Ω. Teniendo en cuenta esta condición inicial para det Fy la segunda equivalencia de (6.6), se deduce div v= 0 ⇔det F= 1 ⇔det Fτ= 1.(6.13) Con el objetivo de desarrollar procedimientos que son válidos en general, estas formulaciones Lagrangiana y no Euleriana de la restricción de incompresibilidad no se tienen en consideración en las formulaciones Lagrangiana y no Euleriana del problema. Para obtener las condiciones de contorno e inicial, (6.8c)-(6.8e), en descripción no Euleriana se aplica (4.8b)-(4.8d) para el caso particular dado en (6.9), deduciéndose vτ(y, t)=(vD)τ(y, t),∀(y, t)∈ΓD τ×[t0, tf],(6.14) −πτ(y, t)I+µτ(y, t)grad yvτ(y, t)F−1 τ(y, t) + F−T τ(y, t)( grad yvτ(y, t))TF−T τ(y, t)mτ(y) =|F−T τ(y, t)mτ(y)|hτ(y, t),∀(y, t)∈ΓR τ×[t0, tf],(6.15) vτ(y, t0) = v0(P(y, τ)),∀y∈¯ Ωτ,(6.16) siendo mτel vector normal unitario a Γτy dirigido hacia el exterior del dominio. Finalmente el problema en configuración no Euleriana se puede escribir de la siguiente manera: Problema 6.1.2. Problema fuerte de Navier-Stokes en configuración no Euleriana, Ωτ×[t0, tf]. Encontrar dos funciones vτ:¯ Ωτ×[t0, tf]→Rdyπτ:¯ Ωτ×[t0, tf]→Rtales que ρτv0 τ+1 det Fτ div yπτdet FτF−T τ −1 det Fτ div yµτgrad yvτF−1 τ+F−T τ( grad yvτ)Tdet FτF−T τ=bτ,en Ωτ×[t0, tf],(6.17a) F−T τ·grad yvτ= 0,en Ωτ×[t0, tf],(6.17b) sujetas a las siguientes condiciones de contorno vτ= (vD)τ,en ΓD τ×[t0, tf],(6.17c) −πτI+µτgrad yvτF−1 τ+F−T τ( grad yvτ)TF−T τmτ=|F−T τmτ|hτ,en ΓR τ×[t0, tf],(6.17d) y a la condición inicial vτ(·, t0) = v0(P(·, τ)),en ¯ Ωτ.(6.17e) Para obtener la formulación variacional de este problema se procede de la siguiente manera: •En primer lugar, se considera la ecuación (6.17a), se multiplica por det Fτy por la función test z∈H1 ΓD τ(Ωτ)y se integra en Ωτ:
42 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. ZΩτ det Fτρτv0 τ·zdy +ZΩτ div yπτdet FτF−T τ−µτgrad yvτF−1 τ+F−T τ( grad yvτ)Tdet FτF−T τ·zdy =ZΩτ det Fτbτ·zdy. •Se aplica la fórmula de Green adecuada y se utiliza la condición de contorno (6.17d) en la segunda integral de arriba obteniéndose, ZΩτ div yπτdet FτF−T τ−µτgrad yvτF−1 τ+F−T τ( grad yvτ)Tdet FτF−T τ·zdy =ZΓN τ det FτπτI−µτgrad yvτF−1 τ+F−T τ( grad yvτ)TF−T τmτ·zdAy −ZΩτ πτdet FτF−T τ·grad yzdy+ZΩτ µτgrad yvτF−1 τ+F−T τ( grad yvτ)Tdet FτF−T τ·grad yzdy =−ZΓR τ det Fτ|F−T τmτ|hτ·zdAy −ZΩτ πτdet FτF−T τ·grad yzdy+ZΩτ µτgrad yvτF−1 τ+F−T τ( grad yvτ)Tdet FτF−T τ·grad yzdy. •Por último, se multiplica la ecuación (6.17b) por det Fτy por la función test q∈L2(Ωτ)y se integra en Ωτ: ZΩτ det FτF−T τ·grad τvτqdy= 0. Por tanto, se obtiene el siguiente problema variacional en configuración no Euleriana: Problema 6.1.3. Problema débil de Navier-Stokes en configuración no Euleriana, Ωτ×[t0, tf]. Encontrar dos funciones vτ: [t0, tf]→H1(Ωτ)yπτ: [t0, tf]→L2(Ωτ)tales que ZΩτ det Fτρτv0 τ·zdy−ZΩτ πτdet FτF−T τ·grad yzdy +ZΩτ µτgrad yvτF−1 τ+F−T τ( grad yvτ)Tdet FτF−T τ·grad yzdy =ZΩτ det Fτbτ·zdy+ZΓR τ det Fτ|F−T τmτ|hτ·zdAy,∀z∈H1 ΓD τ(Ωτ),(6.18a) ZΩτ det FτF−T τ·grad yvτqdy= 0,∀q∈L2(Ωτ),(6.18b) sujetas a la condición de contorno de Dirichlet (6.17c) y a la condición inicial (6.17e).
6.2. Resolución numérica 43 6.1.4. Formulación Lagrangiana La formulación fuerte Lagrangiana o material del Problema 6.1.1 se obtiene particularizando el Problema 6.1.2 para el caso τ=t0, de modo que el dominio de referencia es Ωt0= Ω, y el cambio de variable aplicado x=X(p, t). Problema 6.1.4. Problema fuerte de Navier-Stokes en configuración Lagrangiana. Encontrar dos funciones vm:¯ Ω×[t0, tf]→Rdyπm:¯ Ω×[t0, tf]→Rtales que ρm˙ vm+1 det FDiv πmdet FF−T(6.19a) −1 det FDiv µm∇vmF−1+F−T(∇vm)Tdet FF−T=bm,en Ω×[t0, tf], F−T· ∇vm= 0,en Ω×[t0, tf],(6.19b) sujetas a las siguientes condiciones de contorno, vm= (vD)m,en ΓD×[t0, tf],(6.19c) −πmI+µm∇vmF−1+F−T(∇vm)TF−Tmt0=|F−Tmt0|hm,en ΓR×[t0, tf],(6.19d) y a la condición inicial vm(·, t0) = v0en ¯ Ω.(6.19e) La formulación débil del Problema 6.1.4 se obtiene particularizando, para el caso τ=t0, el procedimiento descrito para definir el Problema 6.1.3: Problema 6.1.5. Problema débil de Navier-Stokes en configuración Lagrangiana. Encontrar dos funciones vm: [t0, tf]→H1(Ω) yπm: [t0, tf]→L2(Ω) tales que ZΩ det Fρm˙ vm·zdp−ZΩ πmdet FF−T· ∇zdp +ZΩ µm∇vmF−1+F−T(∇vm)Tdet FF−T· ∇zdp =ZΩ det Fbm·zdp+ZΓR det F|F−Tmt0|hm·zdAp,∀z∈H1 ΓD(Ω),(6.20a) ZΩ det FF−T· ∇vmqdp= 0,∀q∈L2(Ω),(6.20b) sujetas a la condición de contorno de Dirichlet (6.19c) y a la condición inicial (6.19e). 6.2. Resolución numérica En esta sección se aplicarán distintos esquemas de discretización temporal descritos en el Anexo A.4 a la resolución numérica del Problema 6.1.3.
44 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. 6.2.1. Discretización temporal La notación utilizada en esta parte es la considerada a lo largo del documento. Recordemos que N es el número de pasos de tiempo, ∆t=tf−t0 Nes el tamaño del paso de tiempo y tn=t0+n∆t(para n= 0, . . . , N) son los instantes de la discretización temporal. Para cualquier campo espacial o material ϕ, se introduce la notación, ϕl:= ϕ(·, tl), yϕl ∆tdenotará una aproximación de ϕlobtenida mediante una semi-discretización temporal. Se considera la semi-discretización temporal del problema obtenida mediante la aplicación del método de Newmark [39] y que ha sido empleada en [25] para el Problema 6.1.5 y en [26] para el Problema 6.1.3. Este método pertenece a la familia de métodos de Runge-Kutta-Nyström y ha sido definido en el Anexo A.4 mediante la Tabla de Butcher (A.65). A partir de la definición del método de Runge-Kutta-Nyström (A.58) y de la tabla de Butcher que define el método de Newmark (A.65), se pueden definir las siguientes relaciones para un campo ψsolución de la ecuación diferencial de segundo orden, ∂2ψ ∂t2=ft, ψ,∂ψ ∂t : ∂ψ ∂t n ∆t =∂ψ ∂t n−1 ∆t + ∆t(1 −γ)kn−1 1+γkn−1 2,(6.21) ψn ∆t=ψn−1 ∆t+ ∆t∂ψ ∂t n−1 ∆t + ∆t21−2β 2kn−1 1+βkn−1 2,(6.22) donde kn−1 1ykn−1 2se definen de la forma kn−1 1=ftn−1,ψn−1 ∆t,∂ψ ∂t n−1 ∆t=∂2ψ ∂t2n−1 ∆t ,(6.23) kn−1 2=ftn,ψn−1 ∆t+ ∆t∂ψ ∂t n−1 ∆t + ∆t21−2β 2kn−1 1+kn−1 2, ∂ψ ∂t n−1 ∆t + ∆t(1 −γ)kn−1 1+γkn−1 2.(6.24) Utilizando (6.21) y (6.22) en (6.24) resulta, kn−1 2=ftn,ψn ∆t,∂ψ ∂t n ∆t=∂2ψ ∂t2n ∆t .(6.25) Finalmente, gracias a las igualdades previas se llega a las siguientes expresiones: ψn ∆t=ψn−1 ∆t+ ∆t∂ψ ∂t n−1 ∆t + ∆t2"1−2β 2∂2ψ ∂t2n−1 ∆t +β∂2ψ ∂t2n ∆t#,(6.26) ∂ψ ∂t n ∆t =∂ψ ∂t n−1 ∆t + ∆t"(1 −γ)∂2ψ ∂t2n−1 ∆t +γ∂2ψ ∂t2n ∆t#.(6.27)
6.2. Resolución numérica 45 Estas expresiones también pueden obtenerse mediante el siguiente desarrollo de Taylor cuyas fórmulas se emplearán de manera explícita en la discretización de los problemas planteados: ψn=ψn−1+ ∆t∂ψ ∂t n−1 + ∆t2 β∂2ψ ∂t2n +1 2−β∂2ψ ∂t2n−1! +O(∆t3)1 6−β+O(∆t4),(6.28) ∂ψ ∂t n =∂ψ ∂t n−1 + ∆t γ∂2ψ ∂t2n + (1 −γ)∂2ψ ∂t2n−1! +O(∆t2)1 2−γ+O(∆t3).(6.29) Si se aplican estas fórmulas a ψ=uτse obtienen las siguientes expresiones para el desplazamiento y la velocidad: un τ=un−1 τ+ ∆tvn−1 τ+ ∆t2βan τ+1 2−βan−1 τ+O(∆t3)1 6−β+O(∆t4),(6.30) vn τ=vn−1 τ+ ∆tγan τ+ (1 −γ)an−1 τ+O(∆t2)1 2−γ+O(∆t3).(6.31) A partir de (6.31), si γ6= 0 se puede obtener la siguiente expresión para la aceleración aτ: an τ=1 ∆tγ (vn τ−vn−1 τ)−1 γ−1an−1 τ+O(∆t)1 2−γ+O(∆t2).(6.32) A continuación se puede plantear la semidiscretización temporal del Problema 6.1.3 mediante las fórmulas (6.30) y (6.31), y obtener el esquema general de Newmark en coordenadas no Eulerianas en el que se contempla la posibilidad de que τpueda variar con n, (por simplicidad, se tomará siempre τ0=t0): Problema 6.2.1. Esquema general de Newmark para el problema de Navier-Stokes en configuraciones no Eulerianas, {Ωτn× {tn}}N n=0 con τn< tn∀n≥1. Encontrar cuatro sucesiones de funciones, {an τn,∆t}N n=0,{vn τn,∆t}N n=0,{un τn,∆t}N n=0 y{πn τn,∆t}N n=1,tales
52 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. Encontrar cuatro sucesiones de funciones {an τn,∆t,h}N n=0 ∈QN n=0 Xp1 h(Ωτn),{vn τn,∆t,h}N n=0 ∈QN n=0 Xp1 h(Ωτn), {un τn,∆t,h}N n=0 ∈QN n=0 Xp1 h(Ωτn)y{πn τn,∆t,h}N n=1 ∈QN n=1 Vp2 h(Ωτn)tales que, ZΩτn det Fn τn,∆t,hρn◦Xn τn,∆t,han τn,∆t,h ·zhdy−ZΩτn πn τn,∆t,h det Fn τn,∆t,h(Fn τn,∆t,h)−T·grad yzhdy +ZΩτn µn◦Xn τn,∆t,h grad yvn τn,∆t,h(Fn τn,∆t,h)−1+ (Fn τn,∆t,h)−T( grad yvn τn,∆t,h)T. . . . . . det Fn τn,∆t,h(Fn τn,∆t,h)−T·grad yzhdy =ZΩτn det Fn τn,∆t,hbn◦Xn τn,∆t,h ·zhdy+ZΓR τn det Fn τn,∆t,h|(Fn τn,∆t,h)−Tmτn|hn◦Xn τn,∆t,h ·zhdAy, ∀zh∈Xp1 0h(Ωτn),(6.48a) ZΩτn det Fn τn,∆t,h(Fn τn,∆t,h)−T·grad yvn τn,∆t,hqhdy= 0 ∀qh∈Vp2 h(Ωτn),(6.48b) ZΩτn un τn,∆t,h ˜ uhdy=ZΩτnun−1 τn,∆t,h + ∆tvn−1 τn,∆t,h + ∆t2βan τn,∆t,h +1 2−βan−1 τn,∆t,h˜ uhdy ∀˜ uh∈Xp1 h(Ωτn),(6.48c) ZΩτn vn τn,∆t,h ˜ vhdy=ZΩτnhvn−1 τn,∆t,h + ∆tγan τn,∆t,h + (1 −γ)an−1 τn,∆t,hi˜ vhdy∀˜ vh∈Xp1 h(Ωτn),(6.48d) para 1≤n≤N, con las condiciones de contorno de Dirichlet y las condiciones iniciales, vn τn,∆t,h =vn D◦Xn τn,∆t,h,para todo nudo en ΓD τn,(6.48e) u0 τ0,∆t,h =0,en ¯ Ωτ0=¯ Ω,(6.48f) v0 τ0,∆t,h =γt0,h(v0),en ¯ Ωτ0=¯ Ω,(6.48g) a0 τ0,∆t,h =γt0,h(a0),en ¯ Ωτ0=¯ Ω,(6.48h) y donde las aproximaciones del movimiento Xn τn,∆t,h y del gradiente de deformaciones Fn τn,∆t,h se definen de la forma, Xn τn,∆t,h(y) := y+un τn,∆t,h(y),y∈Ωτn,(6.48i) Fn τn,∆t,h|K:= I+ grad yun τn,∆t,h|K, K ∈ Tτn,h,(6.48j) para 1≤n≤N. En el caso del Problema 6.2.3, su discretización completa es la siguiente: Problema 6.2.10. Esquema de Newmark-Galerkin lineal para el problema de Navier-Stokes en configuraciones no Eulerianas, {Ωτn× {tn}}N n=0 con τn< tn∀n≥1. Encontrar dos sucesiones de funciones {vn τn,∆t,h}N n=0 ∈QN n=0 Xp1 h(Ωτn)y{πn τn,∆t,h}N n=1 ∈QN n=1 Vp2 h(Ωτn) tales que
6.2. Resolución numérica 53 ZΩτn det Fn τn,∆t,hρn◦Xn τn,∆t,h 1 ∆tγ (vn τn,∆t,h −vn−1 τn,∆t,h)−1 γ−1an−1 τn,∆t,h·zhdy −ZΩτn πn τn,∆t,h det Fn τn,∆t,h(Fn τn,∆t,h)−T·grad yzhdy +ZΩτn µn◦Xn τn,∆t,h grad yvn τn,∆t,h(Fn τn,∆t, h)−1+ (Fn τn,∆t,h)−T( grad yvn τn,∆t,h)T. . . . . . det Fn τn,∆t,h(Fn τn,∆t,h)−T·grad yzhdy =ZΩτn det Fn τn,∆t,hbn◦Xn τn,∆t,h ·zhdy +ZΓR τn det Fn τn,∆t,h|(Fn τn,∆t,h)−Tmτn|hn◦Xn τn,∆t,h ·zhdAy,∀zh∈Xp1 0h(Ωτn),(6.49a) ZΩτn det Fn τn,∆t,h(Fn τn,∆t,h)−T·grad τn,∆t,hvn τn,∆t,hqhdy= 0,∀qh∈Vp2 h(Ωτn),(6.49b) para 1≤n≤N, con las condiciones de contorno de Dirichlet y las condiciones iniciales, vn τn,∆t,h =vn D◦Xn τn,∆t,h,para todo nudo en ΓD τn,(6.49c) u0 τ0,∆t,h =0,en ¯ Ωτ0=¯ Ω,(6.49d) v0 τ0,∆t,h =γt0,h(v0),en ¯ Ωτ0=¯ Ω,(6.49e) a0 τ0,∆t,h =γt0,h(a0),en ¯ Ωτ0=¯ Ω,(6.49f) y donde las aproximaciones del desplazamiento un τn,∆t,h, del movimiento Xn τn,∆t,h y del gradiente de deformaciones Fn τn,∆t,h se definen de la forma, un τn,∆t,h =un−1 τn,∆t,h + ∆tvn−1 τn,∆t,h +∆t2 2an−1 τn,∆t,h en Ωτn,(6.49g) Xn τn,∆t,h(y) = y+un τn,∆t,h(y),y∈Ωτn,(6.49h) Fn τn,∆t,h|K=I+ grad yun τn,∆t,h|K, K ∈ Tτn,h,(6.49i) para 1≤n≤N. La actualización de la aceleración para el instante tnse obtiene mediante la fórmula (6.34j): an τn,∆t,h =1 ∆tγ (vn τn,∆t,h −vn−1 τn,∆t,h)−1 γ−1an−1 τn,∆ten Ωτn,(6.49j) siendo 1≤n≤N. En el caso del Problema 6.2.6, su discretización completa se puede obtener modificando en el Problema 6.2.10 la actualización para un τn,∆t,h:
54 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. Problema 6.2.11. Esquema de Newmark-Galerkin lineal modificado para el problema de Navier-Stokes en configuraciones no Eulerianas, {Ωτn× {tn}}N n=0 con τn< tn∀n≥1. Encontrar dos sucesiones de funciones {vn τn,∆t,h}N n=0 ∈QN n=0 Xp1 h(Ωτn)y{πn τn,∆t,h}N n=1 ∈QN n=1 Vp2 h(Ωτn) que satisfagan las ecuaciones, condiciones de contorno y condiciones iniciales (6.49a)-(6.49f), donde las aproximaciones del movimiento y del gradiente de deformaciones se definen de la forma, Xn τn,∆t,h(y) := y+un−1 τn,∆t,h(y)+∆tvn−1 τn,∆t,h(y) + ∆t2 2an−1 τn,∆t,h(y),y∈Ωτn,(6.50a) Fn τn,∆t,h|K:= grad yXn τn,∆t,h|K, K ∈ Tτn,h.(6.50b) El valor de an τn,∆t,h se calcula como en (6.49j) y el de un τn,∆t,h se actualiza con la fórmula (6.30) para un βcualquiera: un τn,∆t,h =un−1 τn,∆t,h + ∆tvn−1 τn,∆t,h + ∆t2βan τn,∆t,h +1 2−βan−1 τn,∆t,h,en Ωτn,(6.50c) para 1≤n≤N. Para la resolución del Problema 6.2.9 en el caso de elegir τn=tn−1∀n(reinicialización en cada paso de tiempo), se siguen los pasos definidos en el Algoritmo 6.1. En el caso de elegir τn=t0∀n, método Lagrangiano, se siguen los pasos definidos en el Algoritmo 6.2. Algoritmo 6.1 Resolución del Problema 6.2.9 para {τn}N n=1 ={tn−1}N n=1 Dadas las condiciones iniciales u0 τ0,∆t,h =0,v0 τ0,∆t,h =γt0,h(v0)ya0 τ0,∆t,h =γt0,h(a0)(obtenidas aplicando las ecuaciones (6.48f)-(6.48h) donde se tiene en cuenta que τ0=t0), unos valores para las constantes del método de Newmark βyγy teniendo en cuenta que un−1 τn,∆t,h =un−1 tn−1,∆t,h =0, para cada instante tn, n∈[1, N]se realizan los siguientes pasos: 1. Resolver el sistema de ecuaciones (6.48a)-(6.48d), sujeto a la condición de contorno (6.48e) y a las igualdades (6.48i)-(6.48j), para obtener un τn,∆t,h,vn τn,∆t,h,an τn,∆t,h yπn τn,∆t,h. 2. Actualizar el dominio computacional ¯ Ωτn+1 =¯ Ωtnaplicando el procedimiento descrito en la Sección 5.1 donde se parte de la configuración en el instante τn=tn−1y para lo que se requiere el desplazamiento aproximado calculado en el punto anterior un τn,∆t,h =un tn−1,∆t,h. 3. Reinicializar las variables vn τn+1,∆t,h =vn ∆t,h yan τn+1,∆t,h =an ∆t,h aplicando el procedimiento descrito en la Sección 5.2 para lo que se requieren los campos aproximados un τn,∆t,h =un tn−1,∆t,h,vn τn,∆t,h = vn tn−1,∆t,h yan τn,∆t,h =an tn−1,∆t,h. Por último, tener en cuenta que un τn+1,∆t,h =un tn,∆t,h =0.
6.2. Resolución numérica 55 Algoritmo 6.2 Resolución del Problema 6.2.9 para τn=t0∀n Dadas las condiciones iniciales u0 ∆t,h =0,v0 m,∆t,h =γt0,h(v0)ya0 m,∆t,h =γt0,h(a0), unos valores para las constantes del método de Newmark βyγ, y un valor de tolerancia εF∈[0,1], para cada instante tn, n∈[1, N]se realizan los siguientes pasos: 1. Resolver el sistema de ecuaciones (6.48a)-(6.48d), sujeto a la condición de contorno (6.48e) y a las igualdades (6.48i)-(6.48j), para obtener un ∆t,h,vn m,∆t,h,an m,∆t,h yπn m,∆t,h. 2. Si minK∈Tt0,h det (Fn ∆t,h|K)≤εF, interrumpir el bucle en pasos de tiempo para reinicializar el problema a la configuración de referencia ¯ Ωtn0eligiendo n0=n: •Actualizar el dominio computacional ¯ Ωtn0aplicando el procedimiento descrito en la Sección 5.1 donde se parte de la configuración inicial ¯ Ωy para lo que se requiere el desplazamiento aproximado calculado en el punto anterior un ∆t,h. •Reinicializar las variables vn0 tn0,∆t,h =vn0 ∆t,h yan0 tn0,∆t,h =an0 ∆t,h aplicando el procedimiento descrito en la Sección 5.2 para lo que se requieren los campos aproximados un0 ∆t,h,vn0 m,∆t,h y an0 m,∆t,h (calculados en el punto anterior). Por último, tener en cuenta que un0 tn0,∆t,h =0y retomar el bucle en pasos de tiempo tomando τn=tn0∀n > n0. Para la resolución del Problema 6.2.10 en el caso de elegir τn=tn−1∀n(reinicialización en cada paso de tiempo) se siguen los pasos definidos en el Algoritmo 6.3. En el caso de elegir τn=t0∀n, método Lagrangiano, se siguen los pasos descritos en el Algoritmo 6.4. Algoritmo 6.3 Resolución del Problema 6.2.10 para {τn}N n=1 ={tn−1}N n=1 Dadas las condiciones iniciales u0 τ0,∆t,h =0,v0 τ0,∆t,h =γt0,h(v0)ya0 τ0,∆t,h =γt0,h(a0)(donde se tiene en cuenta que τ0=t0), un valor para las constante del método de Newmark γy teniendo en cuenta que un−1 τn,∆t,h =un−1 tn−1,∆t,h =0, para cada instante tn,n∈[1, N]se realizan los siguientes pasos: 1. Obtener la aproximación un τn,∆t,h mediante la fórmula explícita (6.49g). 2. Resolver el sistema de ecuaciones (6.49a)-(6.49b), sujeto a la condición de contorno (6.49c) y a las igualdades (6.49h)-(6.49i), para obtener vn τn,∆t,h yπn τn,∆t,h. 3. Obtener an τn,∆t,h mediante la fórmula (6.49j). 4. Actualizar el dominio computacional ¯ Ωτn+1 =¯ Ωtnaplicando el procedimiento descrito en la Sección 5.1 donde se parte de la configuración en el instante τn=tn−1y para lo que se requiere el desplazamiento aproximado calculado en el primer punto un τn,∆t,h =un tn−1,∆t,h. 5. Reinicializar las variables vn τn+1,∆t,h =vn ∆t,h yan τn+1,∆t,h =an ∆t,h aplicando el procedimiento descrito en la Sección 5.2 para lo que se requieren los campos aproximados un τn,∆t,h =un tn−1,∆t,h,vn τn,∆t,h = vn tn−1,∆t,h yan τn,∆t,h =an tn−1,∆t,h. Por último, tener en cuenta que un τn+1,∆t,h =un tn,∆t,h =0.
56 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. Algoritmo 6.4 Resolución del Problema 6.2.10 para τn=t0∀n Dadas las condiciones iniciales u0 ∆t,h =0,v0 m,∆t,h =γt0,h(v0)ya0 m,∆t,h =γt0,h(a0), un valor para la constante del método de Newmark γ, y un valor de tolerancia εF∈[0,1], para cada instante tn,n∈[1, N] se realizan los siguientes pasos: 1. Obtener la aproximación un ∆t,h mediante la fórmula explícita (6.49g). 2. Resolver el sistema de ecuaciones (6.49a)-(6.49b), sujeto a la condición de contorno (6.49c) y a las igualdades (6.49h)-(6.49i), para obtener vn m,∆t,h yπn m,∆t,h. 3. Obtener an m,∆t,h mediante la fórmula (6.49j). 4. Si minK∈Tt0,h det (Fn ∆t,h|K)≤εF, interrumpir el bucle en pasos de tiempo para reinicializar el problema a la configuración de referencia ¯ Ωtn0eligiendo n0=n: •Actualizar el dominio computacional ¯ Ωtn0aplicando el procedimiento descrito en la Sección 5.1 donde se parte de la configuración inicial ¯ Ωy para lo que se requiere el desplazamiento aproximado calculado en el primer punto un ∆t,h. •Reinicializar las variables vn0 tn0,∆t,h =vn0 ∆t,h yan0 tn0,∆t,h =an0 ∆t,h aplicando el procedimiento descrito en la Sección 5.2 para lo que se requieren los campos aproximados un0 ∆t,h,vn0 m,∆t,h y an0 m,∆t,h (calculados en los pasos 1 y 2). Por último, tener en cuenta que un0 tn0,∆t,h =0y retomar el bucle en pasos de tiempo tomando τn=tn0∀n > n0. Para la resolución del Problema 6.2.11, en el caso de τn=tn−1∀n(reinicialización en cada paso de tiempo), se emplea el Algoritmo 6.5 que es una modificación del Algoritmo 6.3 obtenida considerando una nueva actualización de un τn,∆t,h. Para resolver la descripción Lagrangiana de este problema, es decir, el que se obtiene tomando τn=t0∀nen el Problema 6.2.11, se emplea el Algoritmo 6.6. Algoritmo 6.5 Resolución del Problema 6.2.11 para {τn}N n=1 ={tn−1}N n=1 Dadas las condiciones iniciales u0 τ0,∆t,h =0,v0 τ0,∆t,h =γt0,h(v0)ya0 τ0,∆t,h =γt0,h(a0)(donde se tiene en cuenta que τ0=t0), unos valores para las constantes del método de Newmark βyγy teniendo en cuenta que un−1 τn,∆t,h =un−1 tn−1,∆t,h =0, para cada instante tn,n∈[1, N]se realizan los siguientes pasos: 1. Obtener las aproximaciones Xn τn,∆t,h yFn τn,∆t,h mediante las fórmulas explícitas (6.50a)-(6.50b). 2. Resolver el sistema de ecuaciones (6.49a)-(6.49b), sujeto a las condición de contorno (6.49c), para obtener vn τn,∆t,h yπn τn,∆t,h. 3. Obtener an τn,∆t,h mediante la fórmula (6.49j). 4. Actualizar el desplazamiento un τn,∆t,h mediante la fórmula (6.50c), empleando el parámetro β. 5. Actualizar el dominio computacional ¯ Ωτn+1 =¯ Ωtnaplicando el procedimiento descrito en la Sección 5.1 donde se parte de la configuración en el instante τn=tn−1y para lo que se requiere el desplazamiento aproximado calculado en el punto anterior un τn,∆t,h =un tn−1,∆t,h. 6. Reinicializar las variables vn τn+1,∆t,h =vn ∆t,h yan τn+1,∆t,h =an ∆t,h aplicando el procedimiento descrito en la Sección 5.2 para lo que se requieren los campos aproximados un τn,∆t,h =un tn−1,∆t,h,vn τn,∆t,h = vn tn−1,∆t,h yan τn,∆t,h =an tn−1,∆t,h. Por último, tener en cuenta que un τn+1,∆t,h =un tn,∆t,h =0.
6.2. Resolución numérica 57 Algoritmo 6.6 Resolución del Problema 6.2.11 para τn=t0∀n Dadas las condiciones iniciales u0 ∆t,h =0,v0 m,∆t,h =γt0,h(v0)ya0 m,∆t,h =γt0,h(a0), unos valores para las constantes del método de Newmark βyγ, y un valor de tolerancia εF∈[0,1], para cada instante tn, n∈[1, N]se realizan los siguientes pasos: 1. Obtener las aproximaciones Xn ∆t,h yFn ∆t,h mediante las fórmulas explícitas (6.50a)-(6.50b). 2. Resolver el sistema de ecuaciones (6.49a)-(6.49b), sujeto a la condición de contorno (6.49c), para obtener vn m,∆t,h yπn m,∆t,h. 3. Obtener an m,∆t,h mediante la fórmula (6.49j). 4. Actualizar el desplazamiento un ∆t,h mediante la fórmula (6.50c), empleando el parámetro β. 5. Si minK∈Tt0,h det (Fn ∆t,h|K)≤εF, interrumpir el bucle en pasos de tiempo para reinicializar el problema a la configuración de referencia ¯ Ωtn0eligiendo n0=n: •Actualizar el dominio computacional ¯ Ωtn0aplicando el procedimiento descrito en la Sección 5.1 donde se parte de la configuración inicial ¯ Ωy para lo que se requiere el desplazamiento aproximado calculado en punto anterior un ∆t,h. •Reinicializar las variables vn0 tn0,∆t,h =vn0 ∆t,h yan0 tn0,∆t,h =an0 ∆t,h aplicando el procedimiento descrito en la Sección 5.2 para lo que se requieren los campos aproximados un0 ∆t,h,vn0 m,∆t,h y an0 m,∆t,h (calculados en los pasos 4, 2 y 3, respectivamente). Por último, tener en cuenta que un0 tn0,∆t,h =0y retomar el bucle en pasos de tiempo tomando τn=tn0∀n > n0. Observación 6.2.12. La resolución del problema de Darcy (6.44) requiere emplear una combinación inf-sup estable de los espacios de elementos finitos utilizados para representar las incógnitas a0yπ0(ver [40]), como por ejemplo la combinación de elementos de Raviart-Thomas para a0y de tipo L2para π0. Otra alternativa es añadir términos de estabilización a la formulación variacional (6.44) que permitan emplear cualquier combinación de espacios de elementos finitos (ver [41]). Esta estrategia permite elegir un espacio de elementos finitos para a0que facilite la posterior asignación de los grados de libertad del espacio a an τn,∆t,h. La formulación variacional estabilizada empleada para resolver el problema (6.44) es la descrita en [41], es decir, ZΩ ρ0a0·wdp−ZΩ π0Div wdp−ZΩ k1∇π0+ρ0a0·ρ0wdp+ZΩ k2Div a0·Div wdp= ZΩ Div µ0∇v0+ (∇v0)T·wdp−ZΓRµ0∇v0+ (∇v0)Tmt0−h0mt0·wdAp+ZΩ b0·wdp −ZΩ k1b0+ Div µ0∇v0+ (∇v0)T·ρ0wdp+ZΩ k2∇v0·(∇v0)TDiv wdp,∀w∈HΓD( div ,Ω), (6.51) ZΩ Div a0qdp+ZΩ k1(∇π0+ρ0a0)· ∇qdp= ZΩ ∇v0·(∇v0)Tqdp+ZΩ k1b0+ Div µ0∇v0+ (∇v0)T· ∇qdp,∀q∈L2(Ω),(6.52) siendo k1yk2dos factores de estabilización que se fijan a los valores k1=1 2ρ0yk2= 1.
58 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. 6.3. Resultados numéricos En esta sección se evaluarán las características de los métodos numéricos descritos en la sección anterior mediante test numéricos habituales en dos dimensiones. En todos los test, para la discretización espacial se empleará una combinación estable de espacios de elementos finitos continuos sobre mallas triangulares; en concreto se emplearán los elementos MINI [42] (lineal a trozos continuo más la burbuja para la velocidad/aceleración/desplazamiento, y lineal a trozos continuo para la presión) . Para la implementación de estos métodos se emplea como base la plataforma para la resolución de ecuaciones en derivadas parciales FEniCS [15], que permite una representación directa de las formulaciones variacionales mediante el lenguaje UFL (Unified Form Language) [43]. Con esta representación se pueden obtener de manera analítica las derivadas parciales de la formulación variacional con respecto a las incógnitas del problema, habilitando la resolución de problemas no lineales mediante representaciones explícitas de las matrices Jacobianas. Para la resolución de los sistemas lineales y no lineales se emplea, de forma explícita, el paquete PETSc [44]. Las mallas triangulares iniciales se han generado con Gmsh [21] y para el remallado y la reinicialización del método se ha programado una biblioteca en Fortran y Python, que hace uso de la interfaz meshpy [22] para la generación de mallas triangulares. 6.3.1. Problema con solución analítica En este ejemplo se pretende evaluar los órdenes de convergencia de los distintos esquemas propuestos, mediante un problema cuyo dominio espacial es Ω = [0,1]×[0,1],t0= 0 ytf= 1. Se utiliza una viscosidad dinámica µ= 0.1y una densidad ρ= 1. Se ajustan las funciones b,g,vDyv0de tal manera que la solución exacta del problema sea, u(p1, p2, t)=[t3et, t3etcos(p1−1)],(p1, p2, t)∈¯ Ω×[t0, tf],(6.53) v(x1, x2, t) = [(3t2+t3)et,(3t2+t3)etcos(x1−t3et−1)],(x1, x2, t)∈ T ,(6.54) a(x1, x2, t) = [(6t+ 6t2+t3)et,(6t+ 6t2+t3)etcos(x1−t3et−1)],(x1, x2, t)∈ T ,(6.55) π(x1, x2, t) = etsin(x1−t3et−1),(x1, x2, t)∈ T .(6.56) Este problema se ha resuelto empleando las estrategias no Eulerianas definidas en los Algoritmos 6.2, 6.4 y 6.6, para distintas combinaciones de los parámetros (γ, β). En este ejemplo se calculan los errores para el desplazamiento, la velocidad, la aceleración y la presión en la norma l∞(L2 h(Ω)) que se define como || ˆ w||l∞(L2 h(Ω)) := m´axN n=1 ||wn||L2 h(Ω), siendo L2 h(Ω) la aproximación de la norma L2(Ω) mediante fórmulas de cuadratura sobre los vértices de la malla. En las Figuras 6.1, 6.2, 6.3 y 6.4 se muestran los errores asociados con el desplazamiento, la velocidad, la aceleración y la presión frente al número de pasos de tiempo (en escala logarítmica) utilizando una malla espacial fija de 201 ×201 vértices. La elección de los parámetros de Newmark (γ, β) = (1/2,1/6) es la que proporciona los mejores órdenes de convergencia, acorde con los órdenes que se obtienen de los desarrollos de Taylor (6.30) y (6.31).
6.3. Resultados numéricos 59 N: número de pasos de tiempo 101102103104 Error 10-12 10-10 10-8 10-6 10-4 10-2 100 102 Curva de error para el desplazamiento l∞(Lh 2(Ω)) Algoritmo 6.6 (γ,β)=(1,1/6) curva de orden 1 Algoritmo 6.4 γ=1/2 Algoritmo 6.2 (γ,β)=(1/2,1/4) curva de orden 2 Algoritmo 6.6 (γ,β)=(1/2,1/6) curva de orden 3 Figura 6.1: Problema con solución analítica: Orden de convergencia temporal para el desplazamiento considerando una malla espacial de 201 ×201 vértices. N: número de pasos de tiempo 101102103104 Error 10-12 10-10 10-8 10-6 10-4 10-2 100 102 Curva de error para la velocidad l∞(Lh 2(Ω)) Algoritmo 6.6 (γ,β)=(1,1/6) curva de orden 1 Algoritmo 6.4 γ=1/2 Algoritmo 6.2 (γ,β)=(1/2,1/4) curva de orden 2 Algoritmo 6.6 (γ,β)=(1/2,1/6) curva de orden 3 Figura 6.2: Problema con solución analítica: Orden de convergencia temporal para la velocidad considerando una malla espacial de 201 ×201 vértices.
60 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. N: número de pasos de tiempo 101102103104 Error 10-10 10-8 10-6 10-4 10-2 100 102 Curva de error para la aceleración l∞(Lh 2(Ω)) Algoritmo 6.6 (γ,β)=(1,1/6) curva de orden 1 Algoritmo 6.4 γ=1/2 Algoritmo 6.2 (γ,β)=(1/2,1/4) Algoritmo 6.6 (γ,β)=(1/2,1/6) curva de orden 2 Figura 6.3: Problema con solución analítica: Orden de convergencia temporal para la aceleración considerando una malla espacial de 201 ×201 vértices. N: número de pasos de tiempo 101102103 Error 10-10 10-8 10-6 10-4 10-2 100 102 Curva de error para la presión l∞(Lh 2(Ω)) Algoritmo 6.6 (γ,β)=(1,1/6) curva de orden 1 Algoritmo 6.4 γ=1/2 Algoritmo 6.2 (γ,β)=(1/2,1/4) Algoritmo 6.6 (γ,β)=(1/2,1/6) curva de orden 2 Figura 6.4: Problema con solución analítica: Orden de convergencia temporal para la presión considerando una malla espacial de 201 ×201 vértices.
6.3. Resultados numéricos 61 En las Figuras 6.5, 6.6, 6.7 y 6.8 se muestran los errores asociados con el desplazamiento, la velocidad, la aceleración y la presión frente al número de grados de libertad en cada dirección espacial principal para un paso de tiempo fijo ∆t=1 2·104. El orden de convergencia espacial observado es igual a dos para el desplazamiento, la velocidad y la aceleración e igual a uno para la presión. Nx=Ny: número de grados de libertad en cada dirección espacial (h=1/Nx) 101102 Error 10-7 10-6 10-5 10-4 10-3 10-2 Curva de error para el desplazamiento l∞(Lh 2(Ω)) con el Algoritmo 6.6 (γ,β)=(1,1/6) (γ,β)=(1/2,0) (γ,β)=(1/2,1/4) (γ,β)=(1/2,1/6) curva de orden 2 Figura 6.5: Problema con solución analítica: Orden de convergencia espacial para el desplazamiento considerando un paso temporal fijo ∆t=1 2·104.
68 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. Figura 6.12: Problema de cavidad con dos paredes móviles: Líneas de corriente correspondientes a Re = 2000 para diferentes condiciones iniciales. 6.3.3. Modos de oleaje en un recinto confinado Este ejemplo se ha elegido para mostrar las capacidades del esquema (incondicionalmente estable y con precisión óptima) en problemas de chapoteo (sloshing) de gran amplitud. El análisis de los modos de oleaje en un recinto confinado (sloshing waves) ha sido estudiado en varios artículos. Por ejemplo, en [52, 53] se trata este problema empleando una descripción Lagrangian-Euleriana arbitraria (ALE). En [8] se emplea un esquema lineal Lagrangiano. Estos esquemas lineales limitan en gran medida el tamaño del paso de tiempo que se puede emplear en problemas con grandes amplitudes de onda. Por esta razón se ha utilizado el esquema no lineal definido por el Algoritmo 6.2, seleccionando la combinación A-estable definida por los parámetros (γ, β) = (1/2,1/4). El dominio empleado es un tanque rectangular de 0.8metros de ancho y la profundidad inicial del líquido es de 0.3metros. Se considera un fluido incompresible de densidad ρ= 1000kg/m3y una viscosidad µ= 0.001kg/(m·s). Inicialmente el fluido se encuentra en reposo y es excitado mediante una aceleración
6.3. Resultados numéricos 69 horizontal sinusoidal y la aceleración vertical de la gravedad: b(x, t) = ρ(A·g·sin(ωt),−g), donde Aes una constante arbitraria que gobierna la amplitud de la excitación, ges la aceleración de la gravedad y ωes la frecuencia de la excitación. En este ejemplo, A= 0.01,g= 9.8m/s2yω= 5.642rad/s. Se considera una condición de contorno de deslizamiento en las fronteras verticales y en la frontera inferior horizontal, y una condición de Neumann nula en la frontera horizontal superior (frontera libre). El problema se resuelve utilizando el Algoritmo 6.2 con la combinación de parámetros (γ, β) = (1/2,1/4). Para la discretización espacial se emplea una malla de 20×20 vértices. En la Figura 6.13 se muestra, para diferentes pasos de tiempo, el desplazamiento vertical de los vértices situados en las esquinas superiores a lo largo del tiempo. En estas simulaciones no ha sido necesario aplicar la estrategia de reinicialización y remallado. Los resultados obtenidos son semejantes a los publicados en [8, 52, 53]. Cabe destacar que, para pasos de tiempo grandes, el método propuesto es estable aunque la precisión en instantes temporales alejados del instante inicial no es buena.
70 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. Figura 6.13: Problema de modos de oleaje: altura de la ola adimensionalizada en las paredes izquierda (gráfica superior) y derecha (gráfica inferior) frente al tiempo para una malla espacial de 20 ×20 vértices.
6.3. Resultados numéricos 71 Para poder simular periodos de tiempo más largos es necesario aplicar la estrategia de remallado y reinicialización varias veces, tal como indica el Algoritmo 6.2. En las Figuras 6.14 y 6.15 se muestran los resultados obtenidos con esta estrategia. En particular, en la Figura 6.15, se muestra una captura de una configuración del dominio junto con los vectores que representan la velocidad sobre las isoregiones del campo de presiones. Figura 6.14: Problema de modos de oleaje: altura de la ola adimensionalizada en las paredes izquierda y derecha frente al tiempo para ∆t= 0.02 y una malla espacial de 20 ×20 vértices. 6.3.4. Rotura de presa En este ejemplo se considera el colapso de una columna de agua en dos dimensiones. Este problema se emplea a menudo como caso test para problemas con frontera libre, pues las condiciones de contorno y la configuración inicial son muy simples (ver Figura 6.16). Se impone una condición de deslizamiento en la pared vertical izquierda y en la horizontal inferior y condición Neumann nula (frontera libre) en las fronteras horizontal superior y vertical derecha. El ancho de la columna de agua es L= 3.5cm y la altura H= 7cm. La aceleración de la gravedad es g= 980cm/s2, la densidad ρ= 1g/cm3y la viscosidad ν= 0.01cm2/s.
72 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. Figura 6.15: Problema de modos de oleaje: Solución calculada en el instante t= 16.24 s. Las flechas indican el campo de velocidades. v2= 0 v1= 0 v1=v2= 0 t= 0 Figura 6.16: Problema de la rotura de presa: condiciones iniciales y de contorno.
6.3. Resultados numéricos 73 El problema se resuelve empleando el Algoritmo 6.4 y los parámetros (γ, β) = (1/2,0) sin reinicialización con una malla espacial de 20 ×20 vértices y un paso de tiempo ∆t= 0.01s. En la Figura 6.17 se representa la posición horizontal del nodo asociado con la esquina inferior derecha de la configuración inicial frente al tiempo. Estos resultados se comparan con los resultados numéricos obtenidos en [35, 54-56] y con los resultados experimentales recogidos en [3, 57]. En estos resultados se distinguen dos grupos de curvas, cada uno de ellos asociado con uno de los resultados experimentales [3, 57]. Los resultados obtenidos con el Algoritmo 6.4 se asemejan a los datos experimentales de [3]. tiempo t (2g/L)1/2 0 0.5 1 1.5 2 2.5 3 3.5 4 x(t)/L 1 1.5 2 2.5 3 3.5 4 4.5 Algoritmo 6.4 γ=1/2 W.A. Wall et al. [50] E. Walhorn et al. [49] P. Hansbo [48] B. Ramaswamy et al. [28] C.W. Hirt et. al [51] (exp.) J.C. Martin et. al [52] (exp.) Figura 6.17: Problema de la rotura de presa: Desplazamiento horizontal del nodo inferior derecho. En la Figura 6.18 se representan tres configuraciones deformadas del dominio y el campo de presiones asociados con tres instantes temporales. En este ejemplo también se ha comprobado numéricamente que el esquema conserva el volumen con una convergencia de segundo orden en tiempo. En la Figura 6.19 se representa la norma l∞del error del área frente al número de pasos de tiempo para una malla espacial fija. 6.3.5. Flujo 2D alrededor de un cilindro En este ejemplo se considera el flujo alrededor de un cilindro infinito. El problema se modela en una configuración 2D donde el cilindro es un círculo y el dominio computacional es un rectángulo en el que en su centro se ha eliminado el círculo. Este ejemplo se utiliza habitualmente para validar métodos para la resolución de los problemas en mecánica de fluidos. El campo de velocidades que se obtiene es simétrico para valores moderados del número de Reynolds. A medida que el número de Reynolds empieza a aumentar, típicamente en el entorno de Re = 90 el flujo empieza a separarse detrás del cilindro provocando el desprendimiento de vórtices. Se considera un fluido incompresible, una longitud del dominio computacional de 20 metros y una altura de 0,5 metros. El radio del cilindro es de 0,1 metros y se sitúa
74 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. en el centro del rectángulo. Se imponen condiciones Dirichlet de velocidad vertical nula y horizontal constante igual a la unidad en las fronteras horizontales del dominio, y condiciones periódicas en las fronteras verticales. Se fija el valor de la densidad a ρ= 1kg/m3y variando el valor de la viscosidad se consigue el cambio del número de Reynolds para cada ensayo. Para la resolución de este problema se emplea el Algoritmo 6.4 reinicializando cuando la deformación provoca una degeneración en el dominio computacional deformado. Para valores de viscosidad µ= 1.0yµ= 0.01 se obtienen soluciones estables simétricas (ver Figura 6.20). Para un valor de µ= 0.001 los desprendimientos de vórtices a los dos lados del cilindro son notorios y la solución ya no es estacionaria. En la Figura 6.21 se muestra la solución del problema para un valor de µ= 0.001 y cuatro instantes diferentes.
6.3. Resultados numéricos 75 Figura 6.18: Problema de la rotura de presa: Campo de presiones y dominio deformado calculado en los instantes t= 0.01,t= 0.04 yt= 0.09 .
76 Capítulo 6. Mecánica de fluidos. Ecuaciones de Navier-Stokes incompresibles. N: número de pasos de tiempo 102103 Error 10-5 10-4 10-3 10-2 10-1 100 Algoritmo 6.4 γ=1/2 curva de orden 2 curva de orden 1 Figura 6.19: Problema de la rotura de presa: error del área medido con la norma l∞frente al número de pasos de tiempo.
6.3. Resultados numéricos 77 Figura 6.20: Problema del flujo alrededor de un cilindro: líneas de corriente para valores de viscosidad de µ= 1.0(superior) y µ= 0.01 (inferior) obtenidas con una malla espacial de 38910 vértices y ∆t= 0.001.
84 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado Problema 7.1.5. Problema eléctrico fuerte de corriente continua en configuración no Euleriana, Ωτ× [t0, tf]. Encontrar una función Vτ:¯ Ωτ×[t0, tf]→Rtal que div ydet FτF−1 τˇσ(Θτ)F−T τgrad yVτ= 0,(7.13a) para (y, t)∈Ωτ×[t0, tf], sujeta a las condiciones de contorno, Vτ= (VD)τ,en ΓD τ×[t0, tf],(7.13b) −ˇσ(Θτ)F−T τgrad yVτ·F−T τmτ=|F−T τmτ|gτ,en ΓR τ×[t0, tf].(7.13c) Para la configuración Lagrangiana (caso τ=t0) se obtiene el siguiente problema: Problema 7.1.6. Problema eléctrico fuerte de corriente continua en configuración Lagrangiana. Encontrar una función Vm:¯ Ω×[t0, tf]→Rtal que Div det FF−1ˇσ(Θm)F−T∇Vm= 0,(7.14a) en Ω×[t0, tf], sujeta a las condiciones de contorno Vm= (VD)men ΓD×[t0, tf],(7.14b) −ˇσ(Θm)F−T∇Vm·F−Tmt0=|F−Tmt0|gmen ΓR×[t0, tf].(7.14c) Para este problema se considera la correspondiente formulación débil que se obtiene de (3.9): ZΩ det FF−1ˇσ(Θm)F−T∇Vm· ∇ ˜ V dp+ZΓR |F−Tmt0|gmdet F˜ V dAp= 0,∀˜ V∈H1 ΓD(Ω).(7.15) Se considera ahora el caso de simetría cilíndrica. En el sistema de coordenadas cilíndricas (véase A.3) el potencial es una función de la forma Vm(p, t) = Vm(rcos θ, r sin θ, p3, t) = ˆ Vm(r, p3, t),(7.16) y la formulación débil (7.15), considerando funciones test de esta misma forma, se reescribe como sigue: Zˆ Ω det ˆ ˆ F1 + 1 rˆu1ˆ ˆ F−1ˇσ(ˆ Θm)ˆ ˆ F−Tˆ ∇ˆ Vm·ˆ ∇˜ V rdrdp3 +Zˆ ΓR |ˆ ˆ F−Tˆ mt0|ˆgmdet ˆ ˆ F1 + 1 rˆu1˜ V rdl(r, p3)=0,∀˜ V∈H1 ˆ ΓD(ˆ Ω).(7.17)
7.1. Modelización del problema 85 Observación 7.1.7. El término fuente en el problema térmico Q(x, t)es debido al efecto Joule: Q(x, t) = J(x, t)·E(x, t),(7.18) y los términos de esta expresión pueden darse en función del potencial eléctrico Q(x, t) = ˇσ◦Θ(x, t) grad V(x, t)·grad V(x, t).(7.19) Utilizando la regla de la cadena se obtiene la misma expresión del término fuente en una configuración no Euleriana Qτ(y, t) = ˇσ◦Θτ(y, t) grad yVτ(y, t)F−1 τ(y, t)·grad yVτ(y, t)F−1 τ(y, t),(7.20) y en el caso de una configuración Lagrangiana: Qm(p, t) = ˇσ◦Θm(p, t)∇Vm(p, t)F−1(p, t)· ∇Vm(p, t)F−1(p, t),(7.21) y si se considera ahora la simetría cilíndrica (véase A.3) ˆ Qm(r, p3, t) = ˇσ◦ˆ Θm(r, p3, t)ˆ ∇ˆ Vm(r, p3, t)ˆ ˆ F−1(r, p3, t)·ˆ ∇ˆ Vm(r, p3, t)ˆ ˆ F−1(r, p3, t).(7.22) 7.1.3. Modelo mecánico Se plantea a continuación la parte mecánica del problema de electrorrecalcado. Para ello, partiendo del Problema 4.1.1 del movimiento, se considerará una ley constitutiva viscoplástica bajo la hipótesis de pequeñas deformaciones y elasticidad lineal isótropa. A continuación, se describirán los distintos elementos que componen esta ley constitutiva de una forma simplificada. Para una descripción con mayor detalle puede consultarse, por ejemplo, [58]. Cinemática En esta sección nos remitimos a los conceptos de la mecánica de los medios continuos introducidos en el Capítulo 2. Se considera un movimiento Xy se denota su gradiente por F. El determinante del gradiente de deformaciones relaciona el volumen de una región infinitesimal en la configuración deformada con el volumen en la configuración de referencia. Cuando det F= 1 la deformación es isocórica. Por otra parte, Fadmite una descomposición polar: F=RU =VR;V=RURT,(7.23)
86 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado donde Res el tensor ortogonal local de rotación y donde los tensores simétricos definidos positivos Uy Vson, respectivamente, los tensores derecho e izquierdo de estiramiento. Se definen los tensores de deformación de Cauchy-Green derecho e izquierdo de la forma C=FTF=U2,(7.24) B=FFT=V2.(7.25) Otro tensor importante es el de Green-Lagrange (también llamado de Green-Saint Venant) G=1 2(C−I). La medida de la deformación local de un cuerpo establece la tensión. En un movimiento, las deformaciones provocan cambios en las distancias entre dos puntos. Las traslaciones y rotaciones, al no alterar las distancias, no producen tensiones. Dejando por un momento el rigor matemático damos una explicación sencilla sobre el significado de estos tensores: como la distancia infinitesimal entre dos puntos en la configuración deformada viene dada por dx =Fdp, se puede cuantificar el estiramiento, de la siguiente manera: kdxk2=Fdp ·Fdp =FTFdp ·dp =U2kdpk2=Ckdpk2.(7.26) Por otra parte, como F=I+∇use tiene C=I+∇u+ (∇u)T+ (∇u)T∇u,(7.27) B=I+∇u+ (∇u)T+∇u(∇u)T,(7.28) G=1 2(∇u+ (∇u)T+ (∇u)T∇u).(7.29) Si las deformaciones son infinitesimales se pueden despreciar los términos de segundo orden, (∇u)T∇u≈ ∇u(∇u)T≈0, y se obtienen la siguientes aproximaciones: C≈B≈I+∇u+ (∇u)T,(7.30) G≈E=1 2∇u+ (∇u)T.(7.31) siendo Eel llamado tensor de deformaciones infinitesimales. En este caso, la medida de la deformación adquiere la forma kdxk2= (I+ 2E)kdpk2+o(∇u). El vector de tensiones de Cauchy s(x, t, n)es un campo espacial dependiente de la normal na la superficie considerada en el punto x. Puede definirse a través del tensor de tensiones de Cauchy: s(x, t, n) = T(x, t)n. El tensor de tensiones de Cauchy es un campo tensorial simétrico que se define en la configuración Euleriana (es decir, en la deformada) a partir del cual pueden definirse otros tensores de tensiones en configuración Lagrangiana:
7.1. Modelización del problema 87 •Primer tensor de Piola-Kirchhoff: S= det FTmF−T. •Segundo tensor de Piola-Kirchhoff: P= det FF−1TmF−T. •Tensor de tensiones de Kirchhoff: τ= det FTm. Nótese que T,Pyτson tensores simétricos, mientras que el primer tensor de Piola-Kirchoff Sno lo es en general. El campo espacial Ldefinido como L= grad v,(7.32) se conoce como el gradiente de la velocidad. Su parte simétrica D=1 2(L+LT), se denomina tensor de la tasa de deformación (strain rate tensor). Elasticidad lineal A nivel fenomenológico un sólido elástico se caracteriza por una serie de comportamientos claramente definidos. Un sólido elástico se deforma de manera reversible, es decir, al retirar las cargas el sólido recupera su forma original. La deformación del sólido depende únicamente del esfuerzo al que está sometido, y no de la velocidad de la carga ni de la historia. En el caso de elasticidad lineal, la tensión es una función lineal de la deformación. Este comportamiento se puede observar en un ensayo de tensión uniaxial representado en la Figura 7.2, donde una pieza cilíndrica de sección inicial A0y longitud inicial L0es sometida a una fuerza axial Pque provoca una deformación tal que la longitud final es Ly la sección final es A. Después, se deja de aplicar la fuerza axial y la pieza recupera su forma inicial. Los ejes representan s=P A0y e=L L0−1. s e Figura 7.2: Ensayo de tensión uniaxial para un material elástico lineal. El tensor de tensiones de Cauchy Tpara un cuerpo elástico admite la siguiente expresión: T(x, t) = ˇ T(F(p, t),p)|p=P(x,t),∀(x, t)∈ T ,(7.33)
88 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado y en consecuencia S(p, t) = ˇ S(F(p, t),p),P(p, t) = ˇ P(F(p, t),p),∀(p, t)∈¯ Ω×[t0, tf].(7.34) Dado p∈¯ Ω, se introduce el tensor de elasticidad en el punto material pcomo la aplicación lineal C(p) : Lin −→ Lin, definida por C(p) = DFˇ S(I,p).(7.35) Lema 7.1.1. En ausencia de tensiones residuales, es decir, si ˇ T(I,p) = 0se tiene C(p) = DFˇ T(I,p) = DFˇ P(I,p).(7.36) Demostración. La demostración de la primera igualdad se puede consultar en [27] (pág. 194) y de forma análoga se demuestra la segunda igualdad. Concretamente, en primer lugar derivamos la siguiente ecuación con respecto a F: ˇ S(F,p) = Fˇ P(F,p),(7.37) y obtenemos DFˇ S(F,p)[H] = Hˇ P(F,p) + FDFˇ P(F,p)[H],∀H∈Lin,(7.38) donde hemos utilizado la regla del producto. Finalmente el resultado se deduce evaluando esta expresión en F=Iy teniendo en cuenta que las tensiones residuales son nulas. 2 Nótese que las ecuaciones constitutivas para cuerpos elásticos dadas en (7.33)-(7.34) se pueden linealizar despreciando los términos o(∇u)en los siguientes desarrollos de Taylor: Tm(p, t) = ˇ T(I+∇u(p, t),p) = ˇ T(I,p) + C(p)[∇u(p, t)] + o(∇u(p, t)) = C(p)[E(p, t)] + o(∇u(p, t)), (7.39) S(p, t) = ˇ S(I+∇u(p, t),p) = ˇ S(I,p) + C(p)[∇u(p, t)] + o(∇u(p, t)) = C(p)[E(p, t)] + o(∇u(p, t)), (7.40) P(p, t) = ˇ P(I+∇u(p, t),p) = ˇ P(I,p) + C(p)[∇u(p, t)] + o(∇u(p, t)) = C(p)[E(p, t)] + o(∇u(p, t)), (7.41) para (p, t)¯ Ω×[t0, tf], donde se ha asumido que las tensiones residuales son nulas. En las ecuaciones anteriores se ha utilizado el Lema 7.1.1 y que C(p)[W] = 0,∀W∈Skw donde Skw representa el espacio de los tensores antisimétricos (ver [27], pág. 195). Entonces, en ausencia de tensiones residuales y asumiendo que los términos o(∇u)pueden ser despreciados, se obtienen las siguientes aproximaciones: Tm(p, t)≃S(p, t)≃P(p, t)≃ C(p)[E(p, t)],∀(p, t)∈¯ Ω×[t0, tf].(7.42) Para el caso de un material isótropo el tensor de elasticidad tiene la siguiente expresión (ver [27], pág. 196)
7.1. Modelización del problema 89 C(p)[H] = 2µ(p)H+λ(p)( tr H)I,(7.43) para H∈Sym, y donde µyλson campos materiales con valores escalares que se denominan coeficientes de Lamé del material. Entonces, en las situaciones en las que los términos o(∇u)pueden ser despreciados y el material es elástico e isótropo, podemos expresar el tensor Sen función de E, o viceversa, como sigue: S= 2µE+λ( tr E)I,(7.44) E=1 2µS−λ 2µ+dλ( tr S)I,(7.45) y además se verifica (recordemos que des la dimensión espacial) tr S= (2µ+dλ) tr E.(7.46) Nótese que esta igualdad se escribe en función del módulo de compresibilidad K=1 d(2µ+dλ)de la siguiente manera: tr S=d·K( tr E).(7.47) En algunas ocasiones es conveniente descomponer el tensor Sen partes esférica y desviatoria e introducir una nueva variable πasociada a la parte esférica que se puede interpretar como la presión hidrostática a la que está sometida el cuerpo. En concreto, S=Sd+1 d( tr S)I=Sd−πI,(7.48) siendo π=−1 d( tr S)y donde el superíndice ddenota la parte desviatoria del tensor correspondiente. Empleando estas notaciones, es inmediato comprobar que se verifica la siguiente equivalencia (nótese que tr E= Div u): S= 2µE+λ( tr E)I⇐⇒ S= 2µEd−πI, Div u=−π K.(7.49) Elastoplasticidad para pequeñas deformaciones Los materiales elastoplásticos son aquellos que tras ser sometidos a determinadas cargas presentan deformaciones permanentes (plásticas) cuando se retiran todas las cargas aplicadas. Estos materiales, entre otros los metales, presentan una serie de características que pueden ser identificadas en un ensayo de tensión uniaxial (ver [38]). Típicamente, estos ensayos uniaxiales sobre metales dúctiles presentan una curva tensión-deformación como la que se aprecia en la Figura 7.3, en la cual se representa el ensayo en el que una pieza cilíndrica de sección inicial A0y longitud inicial L0es sometida a una fuerza P. Los
90 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado ejes representan s=P A0ye=L L0−1. En su deformación, el material transita por los puntos O,AyB; en ese punto se retira la fuerza, retornando la pieza a la situación C, obteniéndose como resultado una deformación permanente ep. Si se repite el ensayo, el punto de partida es Cy si en el ensayo no se supera el punto B, al retirar la fuerza la pieza retorna al punto C. Este ensayo muestra que existe un dominio elástico definido por la tensión de fluencia, que si no se supera, las deformaciones son reversibles. Si se supera la tensión de fluencia se alcanza un régimen plástico de deformaciones permanentes que provocan, a su vez, una evolución de la tensión de fluencia que se denomina endurecimiento. En el ejemplo las tensiones de fluencia están indicadas por los puntos AyB. s e O A B C epee Figura 7.3: Ensayo de tensión uniaxial para material elástoplástico. El fenómeno descrito en este ensayo puede modelizarse matemáticamente para obtener un conjunto de ecuaciones para el caso uniaxial. En el caso de pequeñas deformaciones se considera una descomposición aditiva de la deformacion infinitesimal uniaxial ε=εe+εp,(7.50) donde εerepresenta la parte elástica de la deformación infinitesimal uniaxial y εpla parte plástica o irreversible. Considerando el modelo de elasticidad lineal para la parte elástica, la tensión uniaxial σ puede expresarse como σ=Eεe,(7.51) para un parámetro Ellamado módulo de elasticidad de Young. El comportamiento del material vendrá determinado por el régimen en que se encuentra en cada momento, ya sea elástico o plástico y la característica que define la frontera entre los dos es la tensión de fluencia σy. Se define ahora un criterio de fluencia en base a una función de fluencia que permite determinar el régimen del material. Así, la función de fluencia Φ(σ, σy) = |σ| − σy,(7.52) permite caracterizar el dominio elástico como
7.1. Modelización del problema 91 E={σ|Φ(σ, σy)<0}.(7.53) Además, la tensión máxima admisible es la tensión de fluencia de manera que siempre se verifica la desigualdad Φ(σ, σy)≤0.(7.54) En el dominio elástico solo se puede producir deformación elástica y, por tanto, ˙εp= 0 si Φ(σ, σy)<0, 6= 0 si Φ(σ, σy) = 0,(7.55) de modo que se puede considerar la siguiente ecuación para definir la evolución temporal de la deformación plástica (regla del flujo plástico) ˙εp= ˙γsign(σ),(7.56) donde ˙γse denomina multiplicador plástico y toma valores no negativos. Además, satisface la condición de complementariedad Φ(σ, σy)˙γ= 0, asociada con la condición (7.55). El último elemento necesario para definir el comportamiento elastoplástico es una ley de endurecimiento de modo que la tensión de fluencia σyno sea una constante sino una función de la deformación plástica acumulada ¯εp(t) = Zt 0 |˙εp(s)|ds. El desarrollo de este modelo junto con la generalización a dos y tres dimensiones pueden consultarse en [58, Cap. 6]). Elasto-visco-plasticidad para pequeñas deformaciones El comportamiento elasto-visco-plástico se puede observar en materiales como metales sometidos a altas temperaturas. Se caracteriza por un comportamiento similar a la elastoplasticidad pero dependiente de la velocidad de la deformación. En un ensayo uniaxial con velocidad de deformación constante, para diferentes velocidades de deformación se pueden observar las curvas tensión-deformación de la Figura 7.4. Cada uno de los ensayos se realiza con una nueva pieza y una velocidad de deformación diferente. Puede observarse que cada curva del ensayo tiene un comportamiento equivalente a la curva obtenida para materiales elastoplásticos. La modelización de este comportamiento es similar a la de la elastoplasticidad, compartiendo los mismo elementos e introduciendo una nueva ecuación para definir la evolución temporal de la deformación plástica. En este caso, la expresión toma la siguiente forma
92 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado s e ˙e1 ˙e2 ˙e3 Figura 7.4: Ensayo uniaxial con velocidad de deformación constante para un material elasto-visco-plástico. ˙εp= ˙γ(σ, σy)sign (σ),(7.57) y el multiplicador plástico ˙γpuede ser definido de varias formas. A modo de ejemplo, se reproduce a continuación el modelo empleado en [58] ˙γ(σ, σy) = 1 ν"|σ| σy1 −1#Φ(σ, σy)≥0, 0 Φ(σ, σy)<0, (7.58) donde las constantes del material son ν, parámetro relacionado con las viscosidad, y el parámetro adimensional . Estos parámetros pueden depender de la temperatura. Este modelo, al igual que en el caso de la elastoplasticidad, puede generalizarse al caso multidimensional, cuya definición puede consultarse en [58, Cap. 11]. Elasto-visco-plasticidad: modelo aditivo de Anand para pequeñas deformaciones A continuación se introduce un modelo elasto-visco-plástico que se obtiene asumiendo pequeñas deformaciones, considerando una descomposición aditiva del tensor de deformaciones en parte elástica y plástica, y utilizando el modelo de Anand para describir el comportamiento plástico del material [59]. Nótese que en estas circunstancias, el tensor de tensiones no corresponde a la totalidad del tensor de deformaciones sino solo a una parte elástica. En concreto, se tiene E=Ee+Ep, donde el primer sumando representa el tensor de las deformaciones elásticas y el segundo el tensor de las deformaciones plásticas, ambos en su versión infinitesimal. Se considera la siguiente ley constitutiva: S(p, t) = C(p) (Ee(p, t)) = C(p) (E(p, t)−Ep(p, t)) ,∀(p, t)∈¯ Ω×[t0, tf],(7.59)
7.1. Modelización del problema 93 siendo C(p)el tensor de elasticidad dado en (7.43) y donde la evolución temporal del tensor de deformaciones plásticas ˙ Epse modela como sigue: ˙ Ep(p, t) = ˙p(p, t)3 2 Sd(p, t) q(Sd(p, t)),(7.60) ˙p(p, t) = Aexp −Q R·Θm(p, t)sinh ξq(Sd(p, t)) s(p, t)1 m ,(7.61) ˙s(p, t) = h01−s(p, t) s∗(p, t) a sign1−s(p, t) s∗(p, t)˙p(p, t),(7.62) s∗(p, t) = es˙p(p, t) Aexp Q R·Θm(p, t)n ,(7.63) para (p, t)∈¯ Ω×[t0, tf], siendo q:Lin −→ R+la aplicación definida por q(H) = r3 2H·H,∀H∈Lin.(7.64) De este modo, para modelar la evolución temporal del tensor de deformaciones plásticas, se añade una única variable interna s(p, t). Nótese que, a partir de (7.60) y asumiendo que el tensor de deformaciones plásticas en el instante inicial tiene traza nula, se deduce que dicho tensor tienen traza nula en todo instante de tiempo y para todo punto material (ya que tr Sd= 0). Utilizando esta propiedad y que el tensor de elasticidad tiene la forma dada en (7.43), se deduce S= 2µ(E−Ep) + λDiv uI ⇐⇒ S= 2µEd−Ep−πI, Div u=−π K,(7.65) Sd= 2µ(Ed−Ep).(7.66) Incorporando esta ley constitutiva a la ecuación del movimiento en configuración Lagrangiana definida en el Problema 4.1.3 se obtienen las siguientes ecuaciones ρ0(p)¨ u(p, t)−Div (S(p, t)) = det F(p, t)bm(p, t),(7.67) S(p, t) = 2µ(p)1 2∇u(p, t) + ∇uT(p, t)−Ep(p, t)+λ(p) Div u(p, t)I,(7.68) o equivalentemente, ρ0(p)¨ u(p, t)−Div (S(p, t)) = det F(p, t)bm(p, t),(7.69) S(p, t) = 2µ(p)1 2∇u(p, t) + ∇uT(p, t)−Ep(p, t)−1−2µ(p) d·K(p)π(p, t)I,(7.70) Div u(p, t) = −π(p, t) K(p),(7.71)
100 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado una nueva ecuación al problema que relaciona la derivada temporal del desplazamiento con la velocidad, para que el problema solo incluya derivadas primeras temporales (ver Anexo A.4). 7.2.1. Discretización espacial En esta sección se aborda la semidiscretización espacial de los problemas (7.8), (7.17) y (7.110) mediante el método de los elementos finitos. Sea ˆ Ωun dominio acotado en R2con frontera poligonal a trozos. Se considera una familia de particiones de ¯ ˆ Ωdenotada por Thconsistente en elementos simpliciales (triángulos) K, de diámetro ≤h. Se asume que estas particiones son compatibles con la división de la frontera de ˆ Ωen ˆ ΓDyˆ ΓR. Esta división de las fronteras es diferente para cada uno de los subproblemas, así Tˆ ΓD,Eˆ ΓDyMˆ ΓDcorresponden a las fronteras Dirichlet de los subproblemas térmico, eléctrico y mecánico respectivamente y Tˆ ΓR,Eˆ ΓR yMˆ ΓRa las fronteras Robin. Se propone una discretización espacial mediante el empleo de espacios de elementos finitos Xp1 h(ˆ Ω) para el desplazamiento ˆ uy la velocidad ˆ v,Vp2 h(ˆ Ω) para la temperatura ˆ Θ, potencial eléctrico ˆ Vy la resistencia a la deformación ˆsyXp3 h(ˆ Ω) para el tensor de deformaciones plásticas ˆ Ep, donde los números enteros positivos p1,p2yp3representan el “orden de aproximación” en el siguiente sentido: Existen tres operadores de interpolación γh:C0(¯ ˆ Ω) →Xp1 h(ˆ Ω),ζh:C0(¯ ˆ Ω) → Vp2 h(ˆ Ω) y ξh:C0(¯ ˆ Ω) → Xp3 h(ˆ Ω) que satisfacen las siguientes propiedades: ||γh(ϕ)−ϕ||s,ˆ Ω≤Q1hr−s||ϕ||r,ˆ Ω,∀ϕ∈C0(¯ ˆ Ω) ∩Hr(ˆ Ω),0≤r≤p1+ 1,(7.111) ||ζh(ϕ)−ϕ||s,ˆ Ω≤Q2hr−s||ϕ||r,ˆ Ω,∀ϕ∈C0(¯ ˆ Ω) ∩Hr(ˆ Ω),0≤r≤p2+ 1,(7.112) ||ξh(M)−M||s,ˆ Ω≤Q3hr−s||M||r,ˆ Ω,∀M∈C0(¯ ˆ Ω) ∩Hr(ˆ Ω),0≤r≤p3+ 1,(7.113) donde s= 0,1,Q1,Q2yQ3son tres constantes positivas independientes de h,H0(ˆ Ω) := L2(ˆ Ω),H0(ˆ Ω) := L2(ˆ Ω),H0(ˆ Ω) := L2(ˆ Ω) y|| · ||r,ˆ Ωes la norma usual en Hr(ˆ Ω),Hr(ˆ Ω) yHr(ˆ Ω). Se definen los siguientes espacios: Xp1 0h(ˆ Ω) = nˆ zh∈Xp1 h(ˆ Ω) : ˆ zh=0en Mˆ ΓDo, iVp2 0h(ˆ Ω) = nˆzh∈ Vp2 h(ˆ Ω) : ˆzh= 0 en iˆ ΓDo, ∀i∈ {T, E}. Para obtener una discretización espacial de los Problemas 7.8, 7.17 y 7.110 se emplean los espacios: Xp1 h(ˆ Ω) para aproximar el espacio funcional asociado al desplazamiento ˆ uy a la velocidad ˆ v,Vp2 h(ˆ Ω) para el asociado a la temperatura ˆ Θ, el potencial eléctrico ˆ Vy la resistencia a la deformación ˆsyXp3 h(ˆ Ω) para el asociado al tensor de deformaciones plásticas ˆ Ep. De esta manera la discretización completa del problema acoplado se define como:
7.2. Resolución numérica 101 Zˆ Ω ˆρ0ˇc(ˆ Θm,h)˙ ˆ Θm,h ˜ Θhrdrdp3+Zˆ Ω det ˆ ˆ Fh1 + 1 rˆu1hˆ ˆ F−1 hˇ k(ˆ Θm,h)ˆ ˆ F−T hˆ ∇ˆ Θm,h ·ˆ ∇˜ Θhrdrdp3 +ZTˆ ΓR|ˆ ˆ F−T hˆ mt0|hαmˆ Θm,h −ˆ T∞,m+mσ0ˆ Θ4 m−ˆ T4 ∞,midet ˆ ˆ Fh1 + 1 rˆu1h˜ Θhrdˆ l =Zˆ Ω ˆ Qmdet ˆ ˆ Fh1 + 1 rˆu1h˜ Θhrdrdp3,∀˜ Θh∈TV1 0h(ˆ Ω),(7.114a) Zˆ Ω det ˆ ˆ Fh1 + 1 rˆu1hˆ ˆ F−1 hˇσ(ˆ Θm,h)ˆ ˆ F−T hˆ ∇ˆ Vm,h ·ˆ ∇˜ Vhrdrdp3 +ZEˆ ΓR|ˆ ˆ F−T hˆ mt0|ˆgmdet ˆ ˆ Fh1 + 1 rˆu1h˜ Vhrdl(r, p3) = 0,∀˜ Vh∈EV1 0h(ˆ Ω),(7.114b) Zˆ Ω ˆρ0¨ ˆ uh·ˆ zhdrdp3+Zˆ Ω ˆ ˆ Sh·ˆ ∇ˆ zhdrdp3+Zˆ Ω 1 r(ˆ S11h−ˆ Sθθh)ˆz1h+ˆ S21hˆz2hdrdp3 =Zˆ Ω ˆ bm·ˆ zhdet ˆ ˆ Fh1 + 1 rˆu1hdrdp3 +ZMˆ ΓR|ˆ ˆ F−T hˆ mt0|det ˆ ˆ Fh1 + 1 rˆu1hˆ hm·ˆ zhdˆ l, ∀ˆ zh∈X1 0h(ˆ Ω),(7.114c) Zˆ Ω ˙ ˆ Ep h·ˆ Hhdrdp3=3 2Zˆ Ω ˆ ˙p h ˆ Sd h q(ˆ Sd h)·ˆ Hhdrdp3,∀ˆ Hh∈ X0 h(ˆ Ω),(7.114d) Zˆ Ω ˙ ˆshˆwhdrdp3=Zˆ Ω(h01−ˆsh ˆ s∗ h a sign 1−ˆsh ˆ s∗ h!)ˆ ˙p hˆwhdrdp3,∀ˆwh∈ V0 h(ˆ Ω),(7.114e) siendo ˆ ˆ Sh= 2ˆµ1 2ˆ ∇ˆ uh+ˆ ∇ˆ uT h−ˆ ˆ Ep h+ˆ λˆ Div ˆ uh+1 rˆu1hI, ˆ Sθθh = 2ˆµ1 rˆu1h−ˆ Ep θθh+ˆ λˆ Div ˆ uh+1 rˆu1h,(7.114f) ˆ ˙p h=Aexp −Q Rˆ Θm,h !"sinh ξq(ˆ Sd h) ˆsh!#1 m ,(7.114g) ˆ s∗ h=es"ˆ ˙p h Aexp Q Rˆ Θm,h !#n .(7.114h) Nótese que en esta semidiscretización ya se ha sustituido el contacto unilateral por una condición de Dirichlet de no penetración. 7.2.2. Discretización temporal El problema semidiscretizado espacialmente mediante un método de elementos finitos constituye un sistema de ecuaciones algebraico-diferenciales (DAE) que se puede resolver con integradores apropiados para este tipo de sistemas.
102 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado Resolución numérica con integradores Runge-Kutta implícitos. El problema 7.114 se puede representar mediante el siguiente sistema de DAEs: ˜ M(ˆ yd,ˆ ya)˙ ˆ yd=˜ f(ˆ yd,ˆ ya, t),(7.115a) 0 = ˜ g(ˆ yd,ˆ ya),(7.115b) ˆ yd(·, t) = ¯ yd(·, t),en ˆ ΓD,(7.115c) ˆ ya(·, t) = ¯ ya(·, t),en ˆ ΓD,(7.115d) ˆ yd(·, t0) = y0 d,en ¯ ˆ Ω,(7.115e) siendo ˜ Muna matriz no singular. Si se introduce la notación ˆ y=ˆ yd ˆ ya,¯ y=¯ yd ¯ ya,f=˜ f ˜ g,M=˜ M0 0 0 ,(7.116) el problema se puede reescribir de una forma más compacta: M(ˆ y)˙ ˆ y=f(ˆ y, t),(7.117a) ˆ y(·, t) = ¯ y(·, t),en ˆ ΓD,(7.117b) ˆ yd(·, t0) = y0 d,en ¯ ˆ Ω.(7.117c) Para la resolución de estos problemas acoplados se empleará un método de Runge-Kutta como los descritos en el Anexo A.4, en concreto se utilizarán los métodos implícitos definidos en la herramienta PETSc. La biblioteca PETSc [20, 44, 62] incluye un módulo de integración temporal denominado time stepping ode solver [16] que permite la resolución numérica de ODEs y DAEs mediante distintos esquemas numéricos (Runge-Kutta, Rosembrock-Wanner, BDF). Para la resolución del problema tipo (7.117) se emplea la familia de esquemas denominada ARKIMEX, constituida por métodos de Runge-Kutta aditivos donde es posible separar el problema en dos partes, una tratada de manera implícita y otra explícita. En el caso de DAEs se emplea únicamente la estrategia implícita que conlleva el empleo de métodos DIRK. En concreto, para el problema (7.117), para cada instante tnconocido ˆ yn−1≈ˆ y(·, tn−1), se resuelve el siguiente sistema de ecuaciones no lineales para cada etapa i∈[1, s]del método: ˇ f(Yn,i,˙ Yn,i, tn−1+ci∆t) = M(Yn,i)˙ Yn,i −f(Yn,i, tn−1+ci∆t) = 0,(7.118a) Z=ˆ yn−1+ ∆t i−1 X j=1 aij ˙ Yn,j,(7.118b) ˙ Yn,i =1 ∆t(Yn,i −Z),(7.118c) cuya matriz jacobiana puede obtenerse empleando la regla de la cadena:
7.3. Resultados numéricos 103 ∂ ∂Yn,i (ˇ f(Yn,i,˙ Yn,i(Yn,i),·)) = ∂ˇ f ∂˙ Yn,i (Yn,i,˙ Yn,i(Yn,i),·)∂˙ Yn,i ∂Yn,i +∂ˇ f ∂Yn,i (Yn,i,˙ Yn,i(Yn,i),·),(7.119) donde ∂˙ Yn,i ∂Yn,i es una función constante, con todas las componentes iguales, que depende del propio método de discretización temporal. Este parámetro lo proporciona el módulo de PETSc. De esta manera, para resolver el problema será necesario definir dos funciones: •ˇ f(y,˙ y, t), •∂ˇ f ∂y(y,˙ y, t, b), donde bes el escalar asociado a ∂˙ Yn,i ∂Yn,i y proporcionado por el método. En el módulo se incluyen estrategias de variación del tamaño del paso de tiempo basadas en estimaciones del error local (ver Anexo A.4). 7.3. Resultados numéricos En este capítulo se va a considerar un problema de conformado de piezas metálicas mediante el proceso de electrorrecalcado. Para este caso se asumirán ciertas simplificaciones que faciliten el tratamiento del problema. Se empleará la formulación variacional semidiscretizada en espacio (7.114). Se emplearán elementos L1 1⊕B3(lineal a trozos continuo más la burbuja) para cada componente del desplazamiento y de la velocidad, L1 0(constante por elemento) para cada componente del tensor de deformaciones plásticas y L1 1(lineal a trozos continuo) para los campos de resistencia a la deformación, temperatura y potencial eléctrico. Para la programación del problema semidiscreto en espacio se emplea la plataforma FEniCS [15] que permite la resolución de ecuaciones en derivadas parciales mediante elementos finitos. Utilizando el lenguaje UFL (Unified Form Language) [43] se implementa la formulación variacional, de tal manera que se emplea una función vectorial no lineal equivalente a ˇ f(y,˙ y, t)definida en la Subsección 7.2.2. Esta representación también permite calcular las derivadas asociadas a la resolución de la integración temporal (∂ˇ ∂y). Para la integración temporal se ha empleado un esquema de orden 3, el indicado con la etiqueta 3 dentro de los métodos ARKIMEX definidos en PETSc [44]. Este es un método DIRK de orden 3 con un esquema encajado de orden 2 para el empleo del paso de tiempo variable. En el Apéndice A.4 se detallan estos procedimientos. Se considera la barra cilíndrica metálica de 1.33 metros de longitud y 0.02875 metros de radio. Al emplear la simetría cilíndrica el dominio computacional se reduce a un rectángulo de 1.33 metros de alto y0.02875 metros de ancho. La malla empleada es una malla de triángulos generada con el programa Gmsh [21] compuesta por 709 triángulos.En este dominio se resolverán de manera acoplada los tres subproblemas indicados (térmico, eléctrico y mecánico). Se consideran las fronteras indicadas en la Figura 7.5 donde el vértice inferior izquierdo se sitúa en el origen de coordenadas. La frontera 2 comienza en un vértice situado en las coordenadas p3= 0.16,r= 0.02875 y tiene una longitud de 0.08 metros. En cada subproblema se aplicarán las condiciones de contorno adecuadas que se indicarán a continuación.
104 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado Figura 7.5: Fronteras consideradas en el ejemplo de proceso de electrorrecalcado. 7.3.1. Datos del problema térmico Se consideran las siguientes expresiones para las propiedades térmicas dependientes de la temperatura, donde la función que proporciona el calor específico ˇc[J kgK ]se ha obtenido de [63] y la de la conductividad térmica ˇ k[W mK ]de [64], ˇc(ˆ Θh) =3.15ˆ Θh−1.68 ·1030.00366ˆ Θh−12+ 1.55 ·1030.00366ˆ Θh−13 −684 0.00366ˆ Θh−14+ 152 0.00366ˆ Θh−15−16.50.00366ˆ Θh−16 + 0.699 0.00366ˆ Θh−17−506 + 892e−199(0.001ˆ Θh−1)2 π0.5,
7.3. Resultados numéricos 105 ˇ k(ˆ Θh) =0.0464ˆ Θh−10.20.00366ˆ Θh−12 + 2.25 0.00366ˆ Θh−13−0.155 0.00366ˆ Θh−14+ 21.3. En las fronteras 1, 2 y 3 se impone la siguiente condición de contorno Robin: −ˇ k(ˆ Θm)ˆ ˆ F−Tˆ ∇ˆ Θm·ˆ ˆ F−Tˆ mt0=|ˆ ˆ F−Tˆ mt0|αˆ Θm−ˆ T∞,m+σ0ˆ Θ4 m−ˆ T4 ∞,m, donde se considera una temperatura ambiente T∞,m = 303 [K], un coeficiente de transferencia térmica α= 20 [W K m2], una emisividad = 0.9y la constante de Stefan-Boltzmann σ0= 5.669 ·10−8[W m2K4]. En la frontera 5 se considera la siguiente condición Robin de contacto con resistencia entre la barra y el yunque: −ˇ k(ˆ Θm)ˆ ˆ F−Tˆ ∇ˆ Θm·ˆ ˆ F−Tˆ mt0=|ˆ ˆ F−Tˆ mt0|5·103ˆ Θm−573. En la frontera 4 se considera una condición de Neumann nula. Como condición inicial para la temperatura se utiliza Θ0= 293 K. 7.3.2. Datos del problema eléctrico de corriente continua Se considera la siguiente expresión para la función que proporciona la conductividad eléctrica ˇσ[1 Ohm m] del material a partir de la temperatura: ˇσ(ˆ Θh) = 1 1.08 ·10−9ˆ Θh−3.23 ·10−80.00366ˆ Θh−12−9.44 ·10−8 . En cuanto a las condiciones de contorno, se considera una condición Dirichlet ˆ V= 0 en la frontera 5. En la frontera 2 se considera una corriente entrante dependiente del tiempo en Amperios: I(t) = 50.62 ·103,si t∈[0,19), 35 ·103,si t∈[19,50), 20 ·103,si t∈[50,60), 20 ·103,si t∈[60,70), 30 ·103,si t∈[70,80), 25 ·103,si t≥80, y como consecuencia, la componente normal de la densidad de corriente será: g(x, t) = −1 2πR0l2 I(t),
106 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado donde R0es el radio del cilindro (0.02875 metros) y l2la longitud de la frontera 2 (0.08 metros). En el resto de fronteras se considera condición de Neumann nula. 7.3.3. Datos del problema mecánico Para el modelo mecánico se consideran los siguientes parámetros del material: •Se utiliza una densidad homogénea ρ(x, t) = 7799 [ kg m3]. •Para los parámetros asociados con la elasticidad lineal se utilizan las Tablas 7.1 y 7.2 que identifican, respectivamente, el valor del módulo de Young Ey del coeficiente de Poisson νcon la temperatura, dado que estos parámetros son dependientes de esta magnitud. En la literatura es habitual encontrar estas tablas para estos parámetros de elasticidad. Se utilizará una interpolación lineal para construir los valores de los parámetros en temperaturas no incluidas en las tablas. Para obtener los coeficientes de Lamé utilizados en la formulación del problema se emplean las siguientes igualdades: µ=E 2(1 + ν), λ=νE (1 + ν)(1 −2ν). •Los parámetros asociados al modelo de elasto-visco-plasticidad se indican en la Tabla 7.3 y han sido obtenidos de [59]. En cuanto a las condiciones de contorno, se imponen condiciones de Dirichlet nulas en la componente normal del desplazamiento en las fronteras 4 y 5. La condición en la frontera 4 viene determinada por la simetría cilíndrica y la condición en la frontera 5 es una forma simplificada de modelar el contacto sin rozamiento. Sobre la frontera 3 se aplica una fuerza normal dependiente del tiempo: F3(t) = 0,si t∈[0,55), t−55 5·105,si t∈[55,60), 105,si t≥60, de tal manera que la función asociada con la condición de contorno Neumann en la frontera 3 será ˆ hm="0 −F3 πR2 0#, donde R0es el radio de la pieza cilíndrica (0.02875 metros). En las fronteras 1 y 2 se aplican las condicones Neumann nulas. Se emplean las siguientes condiciones iniciales: ˆ u(·,0) = 0en ˆ Ω, ˆ Ep(·,0) = 0en ˆ Ω, ˆs(·,0) = 43 ·106en ˆ Ω.
7.3. Resultados numéricos 107 Temperatura Θ [oC]Módulo de Young E[P a] −200 2.17 ·1011 −100 2.17 ·1011 −50 2.15 ·1011 −25 2.14 ·1011 0 2.14 ·1011 20 2.12 ·1011 50 2.10 ·1011 100 2.07 ·1011 150 2.03 ·1011 200 1.99 ·1011 250 1.96 ·1011 300 1.92 ·1011 350 1.88 ·1011 400 1.84 ·1011 450 1.79 ·1011 500 1.75 ·1011 550 1.70 ·1011 600 1.64 ·1011 650 1.55 ·1011 1000 43 ·109 1250 20 ·109 1400 15 ·109 1500 15 ·109 Tabla 7.1: Módulo de Young a diferentes temperaturas para el acero (valores extraídos de [59]).
108 Capítulo 7. Aplicación a un problema de multifísica: proceso industrial de electro recalcado Temperatura Θ [oC]coeficiente de Poisson ν −200 0.282 −100 0.282 −50 0.283 −25 0.284 0 0.284 20 0.285 50 0.286 100 0.287 150 0.288 200 0.289 250 0.291 300 0.293 350 0.295 400 0.297 450 0.299 500 0.302 550 0.306 600 0.311 650 0.318 700 0.318 Tabla 7.2: Coeficiente de Poisson a diferentes temperaturas para el acero. A108 Q R[K] 32514 ξ1.15 h0[Pa] 13.29 ·108 ˜s[Pa] 1.476 ·108 m0.147 a1 n0.06869 Tabla 7.3: Parámetros del modelo elasto-visco-plástico para el acero extraídos de [59].
7.3. Resultados numéricos 109 7.3.4. Resultados En la Figura 7.6 se muestra el dominio computacional original y una vista de detalle, centrada en la zona de interés, que se empleará para representar los diferentes campos. También se indica el nodo que se empleará como punto de medida, pues es el nodo que sufre mayor deformación en la dirección radial. Figura 7.6: Dominio computacional y vista en detalle. En la Figura 7.7 se muestra la evolución temporal del campo de temperaturas sobre el dominio deformado. Cada una de las imágenes se corresponde con un instante temporal (0, 19, 50, 70, 80, 84.42 segundos) ordenados de izquierda a derecha y de arriba a abajo. En la Figura 7.8 se muestra la evolución temporal de la magnitud del tensor de deformaciones plásticas Epen los instantes temporales 83, 84, 84.2, 84.3, 84.4 y 84.42. En la Figura 7.9 se puede observar el cambio de la magnitud del efecto Joule a causa del aumento de la sección transversal del cilindro en la zona de máxima deformación radial. En los instantes representados en la figura la corriente aplicada es constante pero la densidad de corriente varía al modificarse el área de la sección. Este efecto se puede apreciar también en la Figura 7.10 donde se representa la evolución temporal del efecto Joule y el desplazamiento en la dirección rmedidos en el nodo de referencia. Claramente se observa un descenso del valor del efecto Joule cuando se incrementa el desplazamiento. En la parte inferior del gráfico se incluye la evolución temporal de la fuerza y corriente aplicadas como referencia. Los cambios en el efecto Joule en el tramo comprendido entre 0 y 19 segundos son debidos al cambio de temperatura y su efecto en la conductividad eléctrica. Además, se observa un cambio brusco del desplazamiento a partir del instante 80, que se acelera en el entorno del segundo 84. Debido a estos cambios tan repentinos de las magnitudes resulta conveniente el empleo de un método de integración temporal con estrategias robustas de selección de paso de tiempo como el considerado en este trabajo. En este problema concreto, en los instantes iniciales es posible emplear pasos de tiempo del tamaño de 5 segundos, pero en las etapas finales los pasos han de ser del orden de 10−3segundos.
116 Capítulo 8. Conclusiones
Capítulo 9 Trabajo futuro Se identifican posibles líneas de desarrollo de los métodos y aplicaciones introducidos en este trabajo: •Desarrollo y aplicación de nuevas estrategias de reconstrucción de las fronteras deformadas que mejoren el proceso de remallado y reinicialización (ver, por ejemplo, [65]). Actualmente una de las estrategias de remallado que se aplica consiste en desplazar los vértices de la frontera utilizando el campo de deformaciones para obtener la frontera deformada y el dominio computacional deformado, y generar una malla de dicho dominio. En este proceso de deformación de las fronteras es posible que los vértices desplazados tiendan a presentar una distribución no homogénea, pudiendo incluso colapsar. •Aplicación de métodos de Runge-Kutta de alto orden y estrategias adaptativas de elección del paso de tiempo al problema de Navier-Stokes considerado en el Capítulo 6. Cabe destacar que los métodos de discretización temporal de Newmark aplicados al problema de Navier-Stokes se han mostrado muy eficientes, permitiendo desarrollar métodos lineales de alto orden. Sin embargo, teniendo en cuenta el desempeño observado de los métodos Runge-Kutta encajados (ver Anexo A.4) en la resolución del problema de electrorrecalcado (Capítulo 7), se plantea la opción de aplicar esta metodología al problema de Navier-Stokes y realizar análisis comparativos entre los diferentes métodos. Además, en algunos problemas es posible que el empleo de un paso de tiempo fijo no sea una estrategia factible de resolución, en términos de coste computacional. •Implementación informática en el problema de electrorrecalcado del contacto unilateral sin rozamiento con un obstáculo rígido a través de técnicas del tipo semismooth Newton. Nótese que este problema ha sido rigurosamente planteado pero el problema de contacto no ha sido incluido en el algoritmo de resolución para faciliar su resolución. •Simulación del proceso de conformado de piezas metálicas considerado en el Capítulo 7 pero en régimen de corriente alterna. En este caso el problema eléctrico 7.1.5 se sustituye por un problema electromagnético cuya solución se emplea para la obtención de la fuente de calor en el modelo térmico. Nótese que dicha fuente de calor es debida al efecto Joule y su expresión depende del tipo de corriente suministrada al sistema (continua o alterna). 117
118 Capítulo 9. Trabajo futuro
Apéndice A Apéndice 119
120 Apéndice A. Apéndice A.1. Espacios funcionales Sea Aun dominio acotado en Rd. Se denotarán los espacios de Lebesgue por Lr(A)y los espacios de Sobolev por Wm,r(A), siendo r= 1,2, ..., ∞ym∈N. En concreto, Lr(A) := ϕ:A −→ Rmedible,ZA |ϕ|r<∞,1≤r < ∞,(A.1) L∞(A) := {ϕ:A −→ Rmedible,|ϕ| ≤ Kcasi por doquier en A} ,(A.2) Wm,r(A) := nϕ:A −→ Rmedible, ∂βϕ∈Lr(A),β∈Nd,|β| ≤ m, 1≤r≤ ∞o,(A.3) siendo |β|:= |(β1, β2, . . . , βd)|=β1+β2+· · · +βdy∂βϕla derivada parcial ∂βϕ:= ∂|β|ϕ ∂xβ1 1∂xβ2 2· · · ∂xβd d . Para r= 2, se tienen los espacios de Hilbert Hm(A) = Wm,2(A). Además, se denota por H1 ΓP(A)el subespacio cerrado de H1(A)definido por H1 ΓP(A) := ϕ∈H1(A), ϕ|ΓP= 0,(A.4) siendo ΓPuna parte de la frontera de Acon medida no nula. Cuando las aplicaciones tengan como conjunto de llegada el espacio vectorial Rd, escribiremos los espacios en negrita, es decir, Lr(A),Wm,r(A),Hm(A) yH1 ΓP(A). Análogamente, cuando las aplicaciones tengan como conjunto de llegada el espacio Lin, escribiremos, Lr(A),Wm,r(A)yHm(A). A.2. Espacios de elementos finitos En esta sección se fijará la notación y se definirán algunos conceptos relativos a la aproximación de espacios de funciones mediante espacios de elementos finitos. Este apartado corresponde a un resumen de [40, Cap. 2], adaptando partes de la notación. Para un estudio más riguroso se puede consultar esa referencia. En esta primera parte de esta sección se concretan algunos espacios introducidos en la sección anterior y se complementan algunos aspectos. Primeramente, en el espacio de funciones L2(A)de cuadrado integrable en el sentido de Lebesgue: L2(A) := v:A → Rmedible,ZA v2dx < ∞,(A.5) se define la norma kvkA:= ZA v2dx1 2 ,(A.6)
A.2. Espacios de elementos finitos 121 que deriva del producto escalar (u, v)A:= ZA u(x)v(x)dx. (A.7) Análogamente, en el espacio de Sobolev Hm(A) := nv:A → Rmedible, ∂βv∈L2(A),β∈Nd,|β| ≤ mo,(A.8) se considera la norma kvkm,A:= X |β|≤mZA |∂βv|2dx 1 2 ,(A.9) asociada con el producto escalar (u, v)m,A:= X |β|≤mZA ∂βu ∂βvdx. (A.10) Los espacios anteriores son espacios de Hilbert con los productos escalares definidos. Para v∈H1(A), si ∂Aes suficientemente regular (por ejemplo Lipschitziana), es posible definir la traza γv =v|∂A. Entonces se puede definir el subespacio cerrado de H1(A), H1 0(A) := v∈H1(A), v|∂A= 0.(A.11) y, análogamente, H1 ΓP(A) := v∈H1(A), v|ΓP= 0,(A.12) siendo ΓP⊂∂A. Sea Th(A) = {KR}m R=1 una partición de Aen una familia de subdominios: A=Sm R=1 KR.KRse denomina elemento y en este trabajo será un triángulo (para d= 2) o un tetraedro (para d= 3). Se denotan por eilas aristas/caras de cada elemento y por eij =∂Ki∩∂Kj,(A.13) la interfaz entre los elementos KiyKj. Dada una partición del dominio Aen elementos poligonales o poliédricos, una aproximación conforme de H1(A)es un espacio de funciones continuas con un número finito de grados de libertad (DOF por sus siglas en inglés, Degrees Of Freedom). Para obtener un subespacio de dimensión finita habitualmente
122 Apéndice A. Apéndice se emplean funciones polinómicas a trozos sobre los elementos de la partición. La continuidad en Ase garantiza con una correcta elección de la posición de los grados de libertad. Se denota por Pk(K)el espacio de los polinomios de grado ≤ksobre el elemento K. Además, se introducen los siguientes espacios: Rk(∂K) := φ∈L2(∂K), φ|ei∈Pk(ei)∀ei,(A.14) Tk(∂K) := φ∈Rk(∂K)∩C0(∂K).(A.15) Elemento de referencia. Sea ˆ K⊂Rd. Se denota por ∂ˆ Ksu frontera. Sea FKuna transformación regular (al menos C1) de Rden Rdtal que K=FK(ˆ K). Cuando FKes afín, se tiene FK(ˆx) = xK+BKˆx. Si ˆves una función definida en ˆ K, entonces v:= ˆv◦F−1 K,(A.16) es una función definida sobre K. El elemento ˆ Kse denomina elemento de referencia. Un elemento se dice afín si FKes afín. Elemento finito. Un elemento finito se define mediante 3 características: •La geometría: el elemento de referencia ˆ Ky el cambio de variable FK(ˆx). •El espacio de polinomios definidos sobre ˆ K:ˆ P. •Un conjunto de grados de libertad:ˆ Σ. Concretamente, un conjunto de formas lineales {li}1≤i≤dim ˆ P y tales que li:ˆ P→R. Para elementos finitos Lagrangianos, el conjunto de grados de libertad ˆ Σse define mediante un conjunto {ˆai}1≤i≤dim ˆ Pde puntos de ˆ K: li(ˆp) = ˆp(ˆai),1≤i≤dim ˆ P. (A.17) A los puntos ˆaise les denomina nodos. Estos elementos finitos permiten aproximar los espacios de Sobolev de orden uno. Coordenadas baricéntricas. Sea Kun elemento simplicial (es decir, un triángulo para d=2o un tetraedro para d=3). Sean xi(1≤i≤nv), sus vértices. Existe un conjunto de funciones lineales sobre K, λK i(x),1≤i≤nv, llamadas coordenadas baricéntricas tales que λK i(xj) = δij,X i λK i(x) = 1.(A.18)
A.2. Espacios de elementos finitos 123 Las coordenadas baricéntricas λK ison las funciones de base del elemento afín P1(K). Para el elemento finito de grado dos, las funciones de base asociadas con los vértices xison λK i(2λK i−1), y para los nodos situados en el punto medio de la arista que une xicon xj, 4λK iλK j. Se pueden definir ahora aproximaciones del espacio H1(A). Sea Sk(K)un subespacio de Pk(K). Para la partición Th(A)se define L1(Sk,Th) := v∈H1(A), v|K∈Sk(K),(A.19) se puede construir una notación más compacta si no existe ambigüedad L1 k=L(Pk,Th).(A.20) Así, por ejemplo, se representarán por L1 1las aproximaciones del espacio H1(A)mediante polinomios de grado menor o igual que 1. Función burbuja. Las funciones burbuja, definidas sobre el elemento K, son aquellas que se anulan en ∂K y su función base se puede definir mediante las coordenadas baricéntricas: bK d= (d+ 1)d+1 d+1 Y i=1 λK i,en Rd.(A.21) Dado un subespacio Sk(K)⊂H1 0(K)se puede usar la notación ya definida para denotar el espacio asociado a la función burbuja: B(Sk) = L1(Sk,Th),(A.22) y, cuando no exista ambigüedad se puede emplear la siguiente notación: Bk=B(Pk∩H1 0).(A.23)
124 Apéndice A. Apéndice A.3. Sistema de coordenadas cilíndricas Se considera un caso con simetría cilíndrica. En particular, el recinto inicial que ocupa el cuerpo, Ω, se obtiene al revolucionar alrededor de un eje la región del plano limitada por una curva, ˆ Γ, y dicho eje que la corta en dos puntos. Dicha región se denotará por ˆ Ω. Suponiendo que el eje de revolución es el correspondiente a la tercera coordenada cartesiana, p3, se tiene Ω = {(r, θ, p3) : θ∈[0,2π),(r, p3)∈ˆ Ω}, y Γ := ∂Ω = {(r, θ, p3) : θ∈[0,2π),(r, p3)∈ˆ Γ}. El conjunto ˆ Ωse denomina sección meridional de Ω. Nótese que la frontera de ˆ Ωtiene dos partes bien diferenciadas: ∂ˆ Ω = ˆ Λ∪ˆ Γ, siendo ˆ Λ := {(r, p3)∈∂ˆ Ω : r= 0}. Además, las fronteras ΓD,ΓRyΓCtambién se obtienen al revolucionar unas curvas contenidas en ˆ Γ,ˆ ΓD,ˆ ΓRyˆ ΓCalrededor del eje, respectivamente. En adelante, denotamos por er(θ),eθ(θ)yep3los vectores de la base física del sistema de coordenadas cilíndrico, es decir, er(θ) = cos θe1+ sin θe2,eθ(θ) = −sin θe1+ cos θe2,ep3=e3, siendo eiel i-ésimo vector de la base canónica de R3. Entonces, se verifica mt0=mrer+mp3ep3,ˆ mt0=mrˆ e1+mp3ˆ e2, siendo ˆ mt0el vector normal unitario exterior a ˆ Γyˆ eiel i-ésimo vector de la base canónica de R2. Además, si φ(escalar), ϑ(vectorial) y Ψ(tensorial) son campos materiales con simetría de revolución, se definen ˆ φ(r, p3, t) := φ(rcos θ, r sin θ, p3, t),(A.24) ˆ ϑ1(r, p3, t)er(θ) + ˆ ϑ2(r, p3, t)ep3:= ϑ(rcos θ, r sin θ, p3, t),(A.25) ˆ ϑ(r, p3, t) := ˆ ϑ1(r, p3, t)ˆ e1+ˆ ϑ2(r, p3, t)ˆ e2, ˆ Ψ(r, θ, p3, t) = ˆ Ψ11(r, p3, t)er(θ)⊗er(θ) + ˆ Ψ12(r, p3, t)er(θ)⊗ep3 +ˆ Ψ21(r, p3, t)ep3⊗er(θ) + ˆ Ψ22(r, p3, t)ep3⊗ep3 +ˆ Ψθθ(r, p3, t)eθ(θ)⊗eθ(θ) := Ψ(rcos θ, r sin θ, p3, t),(A.26) ˆ ˆ Ψ(r, p3, t) := 2 X i,j=1 ˆ Ψij(r, p3, t)ˆ ei⊗ˆ ej.(A.27)
A.3. Sistema de coordenadas cilíndricas 125 Por tanto, la matriz de coordenadas del tensor ˆ Ψ(r, θ, p3, t)es hˆ Ψ(r, p3, t)i= ˆ Ψ11(r, p3, t) 0 ˆ Ψ12(r, p3, t) 0ˆ Ψθθ(r, p3, t) 0 ˆ Ψ21(r, p3, t) 0 ˆ Ψ22(r, p3, t) .(A.28) Conviene recordar que si wes un campo vectorial y Wes un campo tensorial se tiene ∇w=∂wr ∂r er⊗er+1 r ∂wr ∂θ −1 rwθer⊗eθ+∂wr ∂p3 er⊗ep3 +∂wθ ∂r eθ⊗er+1 r ∂wθ ∂θ +1 rwreθ⊗eθ+∂wθ ∂p3 eθ⊗ep3 +∂wp3 ∂r ep3⊗er+1 r ∂wp3 ∂θ ep3⊗eθ+∂wp3 ∂p3 ep3⊗ep3,(A.29) Div w=1 r ∂ ∂r(rwr) + 1 r ∂wθ ∂θ +∂wp3 ∂p3 ,(A.30) Div W=1 r∂ ∂r(rWrr) + ∂Wrθ ∂θ +r∂Wrp3 ∂p3 −Wθθer(A.31) +∂Wθr ∂r +1 r ∂Wθθ ∂θ +∂Wθp3 ∂p3 +1 rWrθ +1 rWθreθ(A.32) +1 r∂ ∂r(rWp3r) + ∂Wp3θ ∂θ +r∂Wp3p3 ∂p3ep3.(A.33) En particular, si ϑyΨson campos de la forma dada en (A.25) y (A.27), respectivamente, ∇ϑ=∂ˆ ϑ1 ∂r er⊗er+∂ˆ ϑ1 ∂p3 er⊗ep3+1 rˆ ϑ1eθ⊗eθ+∂ˆ ϑ2 ∂r ep3⊗er+∂ˆ ϑ2 ∂p3 ep3⊗ep3 Div ϑ= tr (∇ϑ) = ˆ Div ˆ ϑ+1 rˆ ϑ1,(A.34) Div Ψ=1 r"∂ ∂r(rˆ Ψ11) + r∂ˆ Ψ12 ∂p3 −ˆ Ψθθ#er+1 r"∂ ∂r(rˆ Ψ21) + r∂ˆ Ψ22 ∂p3#ep3(A.35) =(ˆ Div ˆ ˆ Ψ)1+1 rˆ Ψ11 −1 rˆ Ψθθer+(ˆ Div ˆ ˆ Ψ)2+1 rˆ Ψ21ep3,(A.36) y, por tanto, si ues de la forma dada en (A.25), se tiene F=I+∇u=1 + ∂ˆu1 ∂r er⊗er+∂ˆu1 ∂p3 er⊗ep3+1 + 1 rˆu1eθ⊗eθ+∂ˆu2 ∂r ep3⊗er +1 + ∂ˆu2 ∂p3ep3⊗ep3, det F=1 + ∂ˆu1 ∂r 1 + 1 rˆu11 + ∂ˆu2 ∂p3−∂ˆu1 ∂p31 + 1 rˆu1∂ˆu2 ∂r =1 + 1 rˆu1det (I+ˆ ∇ˆ u) = 1 + 1 rˆu1det ˆ ˆ F.
132 Apéndice A. Apéndice
Bibliografía 1. Tezduyar, T. E. Interface-tracking and interface-capturing techniques for finite element computation of moving boundaries and interfaces. Computer Methods in Applied Mechanics and Engineering 195. Incompressible CFD, 2983-3000. issn: 0045-7825 (2006). 2. Harlow, F. H. y Welch, J. E. Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface. Physics of Fluids 8, 2182 (1965). 3. Hirt, C. W. y Nichols, B. D. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of computational physics 39, 201-225 (1981). 4. Osher, S. y Sethian, J. A. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics 79, 12-49 (nov. de 1988). 5. Hou, G., Wang, J. y Layton, A. Numerical Methods for Fluid-Structure Interaction — A Review. Communications in Computational Physics 12, 337-377 (ago. de 2012). 6. Ewing, R. y Wang, H. A summary of numerical methods for time-dependent advection-dominated partial differential equations. J. Comput. Appl. Math. 128, 423-445 (2001). 7. Allievi, A. y Bermejo, R. Finite element modified method of characteristics for the Navier–Stokes equations. International Journal for Numerical Methods in Fluids 32, 439-463 (2000). 8. Benítez, M. y Bermúdez, A. Pure Lagrangian and semi-Lagrangian finite element methods for the numerical solution of Navier–Stokes equations. Applied Numerical Mathematics 95, 62-81 (2015). 9. Douglas, J. y Russell, T. Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. SIAM J. Numer. Anal. 19, 871-885 (1982). 10. Pironneau, O. On the Transport-Diffusion Algorithm and Its Applications to the Navier-Stokes Equations. Numer. Math. 38, 309-332 (1982). 11. Renard, Y. y Poulios, K. GetFEM: Automated FE modeling of multiphysics problems based on a generic weak form language. ACM Transactions on Mathematical Software (TOMS) 47, 1-31 (2020). 12. Hecht, F. New development in FreeFem++. J. Numer. Math. 20, 251-265. issn: 1570-2820 (2012). 13. Kirk, B. S., Peterson, J. W., Stogner, R. H. y Carey, G. F. libMesh: A C++ Library for Parallel Adaptive Mesh Refinement/Coarsening Simulations. Engineering with Computers 22. https://doi. org/10.1007/s00366-006-0049-3, 237-254 (2006). 14. Rathgeber, F. y col. Firedrake: automating the finite element method by composing abstractions. ACM Transactions on Mathematical Software (TOMS) 43, 1-27 (2016). 15. Alnæs, M. S. y col. The FEniCS Project Version 1.5. Archive of Numerical Software 3. doi:10. 11588/ans.2015.100.20553 (2015). 133
134 BIBLIOGRAFÍA 16. Abhyankar, S. y col. Petsc/ts: A modern scalable ode/dae solver library. arXiv preprint arXiv:1806.01437 (2018). 17. Hindmarsh, A. C. y col. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software (TOMS) 31, 363-396 (2005). 18. Guennebaud, G., Jacob, B. y col. Eigen v3 http://eigen.tuxfamily.org. 2010. 19. Davis, T. A. Algorithm 832: UMFPACK V4. 3—an unsymmetric-pattern multifrontal method. ACM Transactions on Mathematical Software (TOMS) 30, 196-199 (2004). 20. Balay, S. y col. PETSc Web page https://www.mcs.anl.gov/petsc. 2019. https://www.mcs.anl. gov/petsc. 21. Geuzaine, C. y Remacle, J.-F. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. International journal for numerical methods in engineering 79, 1309-1331 (2009). 22. Klöckner, A. meshpy: 2D/3D simplicial mesh generator interface for Python (Triangle, TetGen, gmsh) https://github.com/inducer/meshpy. 2020. 23. Dapogny, C., Dobrzynski, C. y Frey, P. Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems. Journal of Computational Physics 262, 358-378. issn: 0021-9991 (2014). 24. Benitez, M. y Bermudez, A. Pure Lagrangian and semi-Lagrangian finite element methods for the numerical solution of convection-diffusion problems. International Journal of Numerical Analysis & Modeling 11 (2014). 25. Benítez, M., Bermúdez, A. y Fontán, P. en Recent Advances in PDEs: Analysis, Numerics and Control 33-59 (Springer, 2018). 26. Benítez, M., Bermúdez, A. y Fontán, P. Non-Eulerian Newmark Methods: A Powerful Tool for Free- Boundary Continuum Mechanics Problems. Journal of Scientific Computing 83 (2020). 27. Gurtin, M. An Introduction to Continuum Mechanics. Mathematics in Science and Engineering 158 (Academic Press, San Diego, 1981). 28. Süli, E. Convergence and nonlinear stability of the Lagrange-Galerkin method for the Navier-Stokes equations. Numer. Math. 53, 459-483 (1988). 29. Ferretti, R. Equivalence of semi-Lagrangian and Lagrange-Galerkin schemes under constant advection speed. Journal of Computational Mathematics, 461-473 (2010). 30. Staniforth, A. y Côté, J. Semi-Lagrangian integration schemes for atmospheric models—A review. Monthly weather review 119, 2206-2223 (1991). 31. Benítez, M. y Bermúdez, A. Numerical analysis of a second order pure Lagrange–Galerkin method for convection-diffusion problems. Part I: Time discretization. SIAM Journal on Numerical Analysis 50, 858-882 (2012). 32. Benítez, M. y Bermúdez, A. Second-Order Pure Lagrange–Galerkin Methods for Fluid-Structure Interaction Problems. SIAM Journal on Scientific Computing 37, B744-B777 (2015). 33. Oñate, E. y Manzan, M. Stabilization techniques for finite element analysis of convection-diffusion problems. Developments in Heat Transfer 7, 71-118 (2000). 34. Hirt, C., Cook, J. y Butler, T. A Lagrangian method for calculating the dynamics of an incompressible fluid with free surface. Journal of Computational Physics 5, 103-124 (1970). 35. Ramaswamy, B. y Kawahara, M. Lagrangian finite element analysis applied to viscous free surface fluid flow. International Journal for Numerical Methods in Fluids 7, 953-984 (1987).
BIBLIOGRAFÍA 135 36. Idelsohn, S. R., Oñate, E. y Pin, F. D. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. International journal for numerical methods in engineering 61, 964-989 (2004). 37. Radovitzky, R. y Ortiz, M. Lagrangian finite element analysis of newtonian fluid flows. International Journal for Numerical Methods in Engineering 43, 607-619 (1998). 38. Gurtin, M. E., Fried, E. y Anand, L. The mechanics and thermodynamics of continua (Cambridge University Press, 2010). 39. Newmark, N. M. A method of computation for structural dynamics. Journal of the engineering mechanics division 85, 67-94 (1959). 40. Boffi, D., Brezzi, F., Fortin, M. y col. Mixed finite element methods and applications (Springer, 2013). 41. Barrios, T. P., Cascón, J. M. y González, M. A posteriori error analysis of an augmented mixed finite element method for Darcy flow. Computer Methods in Applied Mechanics and Engineering 283, 909-922 (2015). 42. Arnold, D. N., Brezzi, F. y Fortin, M. A stable finite element for the Stokes equations. Calcolo 21, 337-344 (1984). 43. Alnæs, M. S., Logg, A., Ølgaard, K. B., Rognes, M. E. y Wells, G. N. Unified form language: A domain-specific language for weak formulations of partial differential equations. ACM Transactions on Mathematical Software (TOMS) 40, 1-37 (2014). 44. Balay, S. y col. PETSc Users Manual inf. téc. ANL-95/11 - Revision 3.14 (Argonne National Laboratory, 2020). https://www.mcs.anl.gov/petsc. 45. Erturk, E., Corke, T. C. y Gökçöl, C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. International journal for Numerical Methods in fluids 48, 747-774 (2005). 46. Ghia, U., Ghia, K. N. y Shin, C. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of computational physics 48, 387-411 (1982). 47. Fortin, A., Jardak, M., Gervais, J. y Pierre, R. Localization of Hopf bifurcations in fluid flow problems. International Journal for Numerical Methods in Fluids 24, 1185-1210 (1997). 48. Gervais, J., Lemelin, D. y Pierre, R. Some experiments with stability analysis of discrete incompressible flows in the lid-driven cavity. International Journal for Numerical Methods in Fluids 24, 477-492 (1997). 49. Kuhlmann, H., Wanschura, M. y Rath, H. Flow in two-sided lid-driven cavities: non-uniqueness, instabilities, and cellular structures. Journal of Fluid Mechanics 336, 267-299 (1997). 50. Perumal, D. A. y Dass, A. K. Multiplicity of steady solutions in two-dimensional lid-driven cavity flows by lattice Boltzmann method. Computers & Mathematics with Applications 61, 3711-3721 (2011). 51. Wahba, E. Multiplicity of states for two-sided and four-sided lid driven cavity flows. Computers & Fluids 38, 247-253 (2009). 52. Huerta, A. y Liu, W. K. Viscous flow with large free surface motion. Computer methods in applied mechanics and engineering 69, 277-324 (1988). 53. Souli, M. y Zolesio, J. Arbitrary Lagrangian–Eulerian and free surface methods in fluid mechanics. Computer methods in applied mechanics and engineering 191, 451-466 (2001).
136 BIBLIOGRAFÍA 54. Hansbo, P. The characteristic streamline diffusion method for the time-dependent incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering 99, 171-186 (1992). 55. Walhorn, E., Kölke, A., Hübner, B. y Dinkler, D. Fluid–structure coupling within a monolithic model involving free surface flows. Computers & structures 83, 2100-2111 (2005). 56. Wall, W. A., Genkinger, S. y Ramm, E. A strong coupling partitioned approach for fluid–structure interaction with free surfaces. Computers & Fluids 36, 169-183 (2007). 57. Martin, J. C. y col. Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 244, 312-324 (1952). 58. De Souza Neto, E. A., Peric, D. y Owen, D. R. Computational methods for plasticity: theory and applications (John Wiley & Sons, 2011). 59. Koric, S. y Thomas, B. G. Thermo-mechanical models of steel solidification based on two elastic visco-plastic constitutive laws. journal of materials processing technology 197, 408-418 (2008). 60. Ciarlet, P. G. y Nečas, J. Unilateral problems in nonlinear, three-dimensional elasticity. Archive for rational mechanics and analysis 87, 319-338 (1985). 61. Yang, H. y Cai, X.-C. Parallel two-grid semismooth Newton-Krylov-Schwarz method for nonlinear complementarity problems. Journal of Scientific Computing 47, 258-280 (2011). 62. Balay, S., Gropp, W. D., McInnes, L. C. y Smith, B. F. Efficient Management of Parallelism in Object Oriented Numerical Software Libraries en Modern Software Tools in Scientific Computing (eds. Arge, E., Bruaset, A. M. y Langtangen, H. P.) (Birkhäuser Press, 1997), 163-202. 63. Serendisnky, F. Peformance analysis and optimization of the plate-rolling process en Mathematical Process Models in Iron and Steelmaking (Amsterdam, 1973). 64. Krzyzanowski, M. y Beynon, J. Finite element model of steel oxide failure during tensile testing under hot rolling conditions. Materials science and technology 15, 1191-1198 (1999). 65. Marchandise, E., Remacle, J.-F. y Geuzaine, C. Optimal parametrizations for surface remeshing. Engineering with Computers 30, 383-402 (2014). 66. Hairer, E., Nørsett, S. y Wanner, G. Solving Ordinary Differential Equations I: Nonstiff Problems isbn: 9783540566700 (Springer Berlin Heidelberg, 2008). 67. Hairer, E. y Wanner, G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems isbn: 9783642052217 (Springer Berlin Heidelberg, 2010).
BIBLIOGRAFÍA 137 Los resultados del Capítulo 6 ya han sido publicados como: Benítez M.a, Bermúdez A.b,c, Fontán P.d(2018) A Second-Order Linear Newmark Method for Lagrangian Navier-Stokes Equations. In: Doubova A., González-Burgos M., Guillén-González F., Marín Beltrán M. (eds) Recent Advances in PDEs: Analysis, Numerics and Control. SEMA SIMAI Springer Series, vol 17. Springer, Cham. https://doi.org/10.1007/978-3-319-97613-6_3 Contribución específica en la publicación Formulaciones del problema, desarrollo e implementación del código para la resolución numérica de los problemas, ensayos numéricos. Índices de calidad Esta publicación presenta un índice CiteScore de 1.0, índice SNIP de 0.777 e índice SJR de 0.191, métricas calculadas por Scopus para el año 2020. Se sitúa en el cuartil 4 (Q4) en la categoría de Applied Mathematics en el año 2020 calculado por Scimago: https://www.scimagojr.com/journalsearch.php?q=21100834919&tip=sid&clean=0 Autorización de la revista La publicación SEMA SIMAI Springer Series perteneciente a la editorial Springer Nature, permite la reutilización del artículo por parte del autor como parte de su tesis: https://www.springer.com/gp/rights-permissions/obtaining-permissions/882 Benítez, M.a, Bermúdez, A.b,c and Fontán, P.dNon-Eulerian Newmark Methods: A Powerful Tool for Free-Boundary Continuum Mechanics Problems. J Sci Comput 83, 44 (2020). https://doi.org/10.1007/s10915-020-01207-y Contribución específica en la publicación Formulaciones del problema, desarrollo e implementación del código para la resolución numérica de los problemas, ensayos numéricos. Índices de calidad Esta publicación presenta un índice CiteScore de 4.1, índice SNIP de 1.488 e índice SJR de 1.53, métricas calculadas por Scopus para el año 2020. Se sitúa en el cuartil 1 (Q1) en la categoría de Applied Mathematics en el año 2020 calculado por Scimago: https://www.scimagojr.com/journalsearch.php?q=23490&tip=sid&clean=0 Autorización de la revista La publicación Journal of Scientific Computing perteneciente a la editorial Springer Nature, permite la reutilización del artículo por parte del autor como parte de su tesis: https://www.springer.com/gp/rights-permissions/obtaining-permissions/882 aDepartamento de Matemáticas, Universidade da Coruña, Centro de Investigación CITIC, Elviña s/n, 15071 A Coruña, Spain bDepartamento de Matemática Aplicada and Instituto de Matemáticas, Universidade de Santiago de Compostela, c/ Lope Gómez de Marzoa s/n, 15782 Santiago de Compostela, Spain cInstituto Tecnológico de Matemática Industrial (ITMATI), Rúa Constantino Candeira s/n, 15782 Santiago de Compostela, Spain dREPSOL Technology Lab, Autovía de Extremadura s/n, 28935 Móstoles, Madrid, Spain