Aproximación de un modelo de cristales líquidos nemáticos con un esquema completamente discreto y penalizado

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Aproximaci´on de un modelo de cristales l´ıquidos nem´aticos con un esquema completamente discreto y penalizado F. Guill´ en-Gonz´ alez1, J.V. Guti´ errez-Santacreu1 1Dpto. E.D.A.N., Universidad de Sevilla, Aptdo. 1160, E-41080 Sevilla. E-mails: [email protected], [email protected]. Palabras clave: cristales l´ıquidos, Navier-Stokes, estabilidad, convergencia, elementos finitos Resumen En esta comunicaci´on proponemos y analizamos un esquema num´erico totalmente discreto con elementos finitos en espacio y diferencias finitas en tiempo para aproximar un modelo de cristales l´ıquidos nem´aticos (de tipo Eriksen-Leslie) por medio de un modelo penalizado de tipo Ginzburg-Landau. Se muestra la estabilidad y convergencia de dicho esquema hacia una soluci´on d´ebil del problema continuo, respecto de los par´ametros de discretizaci´on y del par´ametro de penalizaci´on. 1. Introduci´on Supongamos Ω ⊂IR3un dominio acotado y de frontera ∂Ω poli´edrica tal que el problema de Stokes tenga regularidad W1,∞×L∞en velocidad y presi´on. Denotamos Q= Ω×(0, T) y Σ = Γ×(0, T), donde [0, T] es el intervalo temporal de observaci´on, para T > 0. Las inc´ognitas son: u:Q→IR3el campo de velocidades, p:Q→IR la presi´on del fluido y d:Q→IR3la orientaci´on de las macromol´eculas de cristales l´ıquidos, verificando el siguiente problema en derivadas parciales:            |d|= 1, ∂td+u· ∇d−γ∆d−γ|∇d|2d=0en Q, ∂tu+u· ∇u−ν∆u+∇p+λ∇ · ((∇d)t∇d) = 0en Q, ∇ · u= 0 en Q, u= 0,d=len Σ, u|t=0 =u0,d|t=0 =d0en Ω, (1) donde u0: Ω →IR3yd0: Ω →IR3son las condiciones iniciales, l:∂Ω→IR3es la condici´on de Dirichlet para el vector de orientaci´on d,ν > 0 es una constante dependiente 1 F. Guill´en-Gonz´alez, J. V. Guti´errez-Santacreu de la viscosidad del fluido, λ > 0 es una constante de elasticidad y γ > 0 es una constante del tiempo de relajaci´on. (∇d)tdenota la matriz traspuesta de ∇d. Decir que se supone la condici´on de Dirichlet para dindependiente del tiempo, luego en consecuencia d0¯¯∂Ω=l. En principio los argumentos que desarrollaremos no son v´alidos para el caso dependiente del tiempo. Imponiendo que |d0|= 1 en Qy|l|= 1 sobre Σ, se tiene existencia de soluci´on global en tiempo de (1) con la siguiente regularidad: d∈L∞(0, T;H1(Ω)),u∈L∞(0, T;L2(Ω)) ∩L2(0, T;H1(Ω)). Esta soluci´on se encuentra mediante un proceso asint´otico [5], a partir del modelo penalizado, de tipo Ginzburg-Landau, que se obtiene de (1) relajando la restricci´on |d|= 1 por |d| ≤ 1, y en el sistema para dcambiando el t´ermino m´as no lineales |∇d|2d(dependiente del multiplicador de Lagrange asociado a la restricci´on|d|= 1) por el t´ermino de penalizaci´on fε(d) = ε−2(|d|2−1)d, asociado al par´ametro ε > 0, llegando    |d| ≤ 1, ∂td+u· ∇d+γ(fε(d)−∆d) = 0en Q, ∂tu+u· ∇u−ν∆u+∇p+λ∇ · ((∇d)t∇d) = 0en Q, ∇· u= 0 en Q, (2) junto con las condiciones iniciales y de contorno de (1). Obs´ervese que fε(d) = ∇d(Fε(d)) para cada d∈IR3con Fεla funci´on potencial: Fε(d) = 1 4ε2(|d|2−1)2. Las principales dificultades para dise˜nar esquemas num´ericos estables y convergentes para el modelo penalizado (2) son: aproximar dcon elementos finitos s´olo globalmente continuos aunque la regularidad de des H2y obtener informaci´on de la restricci´on |d| ≤ 1 aunque el esquema num´erico no verifique puntualmente dicha restricci´on. En [3] se introduce un esquema (lineal) con una variable auxiliar para aproximar −∆d, que resulta ser incondicionalmente estable y convergente hacia (2), obteni´endose adem´as estimaciones de error y convergencia de m´etodos iterativos para desacoplar el esquema. Ver tambi´en los trabajos previos [6, 7] para otros esquemas menos eficientes para (2). El objetivo de este trabajo es dise˜nar esquemas num´ericos estables y convergentes hacia (1) a trav´es del modelo penalizado, haciendo tender ε→0 junto con los par´ametros discretos en espacio y tiempo (h, k). Para ello, nos encontramos con las siguientes dificultades principales: 1. Obtener estimaciones de estabilidad independientes de ε, ya que las estimaciones de los esquemas anteriores explotan cuando ε→0. 2. Se pierde la regularidad H2para d, que dificulta primero la compacidad L2para u, que se usa en el paso al l´ımite en el sistema para dy segundo la compacidad L2para ∇dfundamental para pasar al limite en el sistema de momentos. Presentaremos un esquema num´erico lineal aunque completamente acoplado, condicionalmente estable bajo restricciones que involucran los tres par´ametros (h, k, ε) y convergente hacia una soluci´on d´ebil del problema (1). 2 Aproximaci´on de un modelo de cristales l´ıquidos (CEDYA 2007) 2. Esquema num´erico Denotaremos por ³·,·´al producto escalar en L2(Ω). El esquema que proponemos para aproximar las inc´ognitas (velocidad, presi´on y vector de orientaci´on) est´a basado en la siguiente formulaci´on variacional mixta del problema (2) (ver [3] para m´as detalles): ³∂tu,¯ u´+ν³∇u,∇¯ u´+³(u· ∇)u,¯ u´−³(∇d)tw,¯ u´−³p, ∇ · ¯ u´= 0,∀¯ u∈H1 0, ³∂td,¯ w´+³(u· ∇)d,¯ w´+γ³fε(d),¯ w´+γ³w,¯ w´= 0,∀¯ w∈L2, ³∇ · u,¯q´= 0,∀¯q∈L2 0, ³∇d,∇¯ d´−³w,¯ d´= 0,∀¯ d∈H1 0, donde hemos usado previamente la igualdad ∇ · ((∇d)t∇d) = 1 2∇(|∇d|2)+(∇d)t∆d, modificado la presi´on p∼p+1 2|∇d|2e introducido la variable auxiliar w=−∆d. Suponemos por simplicidad, una partici´on uniforme de [0, T] siendo tn=nk, donde k=T/N es el paso de tiempo con N∈N. En cada etapa de tiempo, la velocidad, la presi´on y el vector de orientaci´on (u, p, d) son aproximados en espacios de elementos finitos de funciones globalmente continuas (Xh, Qh,Dh)⊂(H1 0(Ω), L2 0(Ω),H1(Ω)) y el laplaciano del vector de orientaci´on wes aproximado en un espacio de elementos finitos de funciones discontinuas a trozos Wh⊂L2(Ω). Los elementos finitos est´an asociados a una familia regular y quasi-uniforme de triangulaciones {Th}h>0de Ω tales que (Xh, Qh) verifican la condici´on inf-sup discreta y (Wh,Dh) satisfacen: (Xh· ∇)Dh⊂WhyDh⊂Wh. Por ejemplo, la siguiente aproximaci´on verifica las hip´otesis anteriores: P1+bubble/P1para (Xh, Qh), P1+bubble(discontinuos)/P1para (Wh,Dh). Adem´as, suponemos las restricciones (de estabilidad) sobre los par´ametros (k, h, ε): (S) l´ım (h,k,ε)→0 k h2ε6→0 y l´ım (h,k,ε)→0 h ε4→0. Entonces, el algoritmo num´erico que presentamos consiste en: Inicializaci´on: Sea (u0 h,d0 h)∈(Xh,Dh) determinadas aproximaciones de (u0,d0), con d0 h|∂Ω=lhpara lhuna aproximaci´on frontera de l. Etapa n+ 1:Dado (un h,dn h)∈(Xh,Dh) con dn h|∂Ω=lh, encontrar (un+1 h,wn+1 h)∈ Xh×Why (pn+1 h,dn+1 h)∈Qh×Dh(con dn+1 h|∂Ω=lh) tal que: Ãun+1 h−un h k,¯ uh!+³(un h· ∇)un+1 h,¯ uh´+1 2³∇ · un hun+1 h,¯ uh´+ν³∇un+1,∇¯ uh´ −λ³(∇dn h)twn+1 h,¯ uh´+λ³∇ · ¯ uh, Fε(dn h)´−³epn+1 h,∇ · ¯ uh´= 0,∀¯ uh∈Xh, (3) 3 F. Guill´en-Gonz´alez, J. V. Guti´errez-Santacreu ³¯ph,∇ · un+1 h´= 0,∀¯ph∈Qh,(4) Ãdn+1 h−dn h k,¯ wh!+³(un+1 h·∇)dn h,¯ wh´+γ³fε(dn h)+wn+1 h,¯ wh´= 0,∀¯ wh∈Wh,(5) ³∇dn+1 h,∇¯ dh´−³wn+1 h,¯ dh´= 0,∀¯ dh∈D0h:= Dh∩H1 0(Ω),(6) donde epn+1 h=pn+1 h+λFε(dn h) es una presi´on modificada debido al termino estabilizador de tipo potencial λ³∇ · ¯ uh, Fε(dn h)´. Notar que dn+1 h−dn h∈D0hgracias a que el dato de contorno lpara el vector de orientaci´on no depende del tiempo. Ya que el esquema (3)-(6) es un sistema lineal algebraico es f´acil comprobar la existencia y unicidad de soluci´on una vez que las estimaciones a priori sean obtenidas en la Secci´on siguiente. 3. Estabilidad condicional En l´ıneas generales, tomando primero ¯ uh= 2 kun+1 hen (3), ¯ph=pn+1 hen (4), y luego ¯ wh= 2 λ k ¡wn+1 h+Ph(fε(dn h))¢en (5) conjuntamente con ¯ dh=dn+1 h−dn h∈D0hen (6), donde Phes el proyector en L2(Ω) sobre Wh, y usando las igualdades −³(∇dn h)twn+1 h,un+1 h´+³un+1 h· ∇dn h,wn+1 h´= 0, ³un+1 h· ∇dn h,fε(dn h)´+³∇ · un+1 h, Fε(dn h)´= 0, llegamos al siguiente resultado (ver [4] para los detalles). Denotamos por | · | a la norma en L2(Ω). Lema 1 Supongamos que existe una constante Cd>0independiente de (h, k, ε)tal que |un h|2+λ|∇dn h|2≤Cd. Entonces, existen h0>0,k0>0yε0>0tales que para todo h≤h0,k≤k0yε≤ε0 satisfaciendo la hip´otesis (S), la soluci´on correspondiente (un+1 h,dn+1 h,wn+1 h)del problema discreto (3)-(6) verifica la siguiente desigualdad:            ³|un+1 h|2− |un h|2+|un+1 h−un h|2´+νk|∇un+1 h|2 +λ³|∇dn+1 h|2− |∇dn h|2+|∇(dn+1 h−dn h)|2´+λγk|Ph(fε(dn h)) + wn+1 h|2 +2λZΩ (Fε(dn+1 h)−Fε(dn h)) + λ ε2ZΩ³1 4(|dn+1 h|2− |dn h|2)2+|dn+1 h−dn h|2´≤0. (7) A continuaci´on, por un proceso de inducci´on sobre la etapa de tiempo apoy´andonos en el lema anterior establecemos el siguiente resultado: 4 Aproximaci´on de un modelo de cristales l´ıquidos (CEDYA 2007) Teorema 2 Existen h0,k0yε0de modo que para cualquier h≤h0,k≤k0yε≤ε0 satisfaciendo (S), la soluci´on del problema discreto (3)-(6) verifica las estimaciones: i)m´ax 0≤n≤N|un h| ≤ C, ii)k N−1 X n=0 |∇un+1 h|2≤C, iii) N−1 X n=0 |un+1 h−un h|2≤C, iv)m´ax 0≤n≤Nkdn hkH1(Ω) ≤C, v)k N−1 X n=0 |Ph(fε(dn h)) + wn+1 h|2≤Cvi) N−1 X n=0 kdn+1 h−dn hk2 H1(Ω) ≤C, vii)m´ax 0≤n≤NZΩ Fε(dn h)≤Cviii)1 ε2 N−1 X n=0 ZΩ (|dn+1 h|2− |dn h|2)2≤Cix)1 ε2 N−1 X n=0 ZΩ |dn+1 h−dn h|2≤C donde C > 0es una constante independiente de (k, h, ε). 4. Compacidad para uyd Eligiendo como funci´on test ¯ wh=Ph¯ wen (5), con w∈L3(Ω), y usando las estimaciones del Teorema 2 y la estabilidad en L3(Ω) del proyector ortogonal en L2sobre Wh ([2]), obtenemos Lema 3 Bajo las condiciones del Teorema 2, se tiene k N−1 X n=0 ° ° ° ° ° dn+1 h−dn h k° ° ° ° ° 2 L3/2(Ω) ≤C, donde C > 0una constante independiente de (k, h, ε). Como consecuencia de este lema y de las estimaciones de estabilidad del Teorema 2, conseguimos la compacidad de la sucesi´on (dh,k,ε) en Lq(0, T;Lr(Ω)) con 1 ≤r < 6 y 1≤q < ∞, donde dh,k,ε es la funci´on continua y lineal a trozos que vale dh,k,ε(tn) = dn h. Para la compacidad de la velocidad discreta {un+1 h}consideramos el espacio de velocidades de divergencia discreta cero Vh={vh∈Xh:³∇ · vh, qh´= 0,∀qh∈Qh} y consideramos A−1 h:Vh→Vhel inverso del operador de Stokes discreto definido como ³∇A−1 huh,∇vh´=³uh,vh´,∀vh∈Vh. Se puede probar que |∇A−1 huh|ykuhkV0 hson normas equivalentes (V0 hes el dual de Vh). Multiplicando (3) por k, sumando (localmente) para n=m, ..., m−1+ r, seleccionado ¯ uh=A−1 h(um+r h−um h) como funci´on test y finalmente sumando (globalmente) para m= 0, ..., N −r, llegamos a la siguientes estimaciones de una derivada fraccionaria para uh,k,ε (la funci´on constante a trozos que vale uh,k,ε(t) = un+1 hen (tn, tn+1]): Lema 4 Bajo las hip´otesis del Teorema 2, se tiene ZT−δ 0 kuh,k,ε(t+δ)−uh,k,ε(t)k2 V0 hdt ≤C δ1/2∀δ: 0 < δ < T, (8) donde C > 0es independiente de (h, k, ε). 5 F. Guill´en-Gonz´alez, J. V. Guti´errez-Santacreu N´otese que la derivada fraccionaria en tiempo para la velocidad discreta ha sido acotada en la norma V0 hla cual “se mueve” con respecto al par´ametro de espacio h. Luego, no podemos aplicar los resultados de compacidad para espacios Lpen tiempo con valores en espacios de Banach, dados por J. Simon en ([8]). Por tanto, la siguiente idea es encontrar una norma, que no dependa de los par´ametros de discretizaci´on, donde la derivada fraccionaria en tiempo pueda ser acotada. Para ello, consideramos el espacio V={u∈H1 0(Ω) : ∇ · u= 0} y la proyecci´on ortogonal Rh:Vh→Vtal que ³∇(Rhuh−uh),∇v´= 0,∀v∈V. Usando las siguientes propiedades del operador Rh: kRhuhk ≤ kuhk(dependencia continua en H1), |Rhuh−uh| ≤ C h|∇ · uh|(estimaci´on de error en L2), se prueba ([4]) que kRhuhkV0≤C³h|∇ · uh|+kuhkV0 h´. Tomando uh=um+r h−um hy usando (8), se tiene ZT−δ 0 kRhuh,k,ε(t+δ)−Rhuh,k,ε(t)k2 V0dt ≤C δ1/2+C h. Entonces, la convergencia fuerte de sucesi´on (Rhuh,k,ε) en L2(L2) se deduce a partir del siguiente resultado obtenido por P. Az´erad y F. Guill´en-Gonz´alez [1], para la tripleta de espacios V,→H,→V0, donde H={u∈L2(Ω) : ∇ · u= 0,u·n¯¯∂Ω= 0}. Teorema 5 Sea T > 0, y sea X,→B,→Y, con la primera inyecci´on compacta y la segunda continua. Sea {f²}²>0una familia de funciones de Lp(0, T;X), 1 ≤p < ∞tal que (H1) {f²}²>0est´a acotada en Lp(0, T;X), (H2) kf²(t+δ)−f²(t)kLp(0,T;Y)≤ϕ(δ) + ψ(²) con l´ım δ→0ϕ(δ) = 0 y l´ım ²→0ψ(²) = 0. Entonces, la familia {f²}²>0posee un punto de acumulaci´on en Lp(0, T;B) cuando ²→0. Finalmente, gracias a la aproximaci´on (externa) de VhaV, se tiene ([4]) la convergencia fuerte de (uh,k,ε) en L2(L2(Ω)). 5. Convergencia del esquema para d La convergencia para (5) est´a basado en el siguiente resultado: Lema 6 Los dos sistemas siguientes son equivalentes: ∂td+u· ∇d−γ∆d−γ|∇d|2d=0en Q, (9) y |d|= 1, ∂td∧d+u· ∇d∧d−γ∇ · (∇d∧d) = 0en Q. (10) 6 Aproximaci´on de un modelo de cristales l´ıquidos (CEDYA 2007) Que el vector de orientaci´on l´ımite d, es decir, el vector de orientaci´on obtenido como consecuencia del Teorema 2, verifique la restricci´on (10)ase deduce tomando l´ımite en la estimaci´on m´ax nZΩ Fε(dn h) = 1 ε2ZΩ (|dn+1 h|2−1)2≤C. Luego, basta probar (10)b Usando Ph:L2(Ω) →Whel proyector respecto de los producto escalar en L2(Ω) sobre Wh, las hip´otesis sobre los espacios discretos y el hecho que fε(dn h)∧dn h= 0, llegamos a la siguiente versi´on discreta de (10)ba partir de (6): Ãdn+1 h−dn h k∧dn h,¯ dh!+³(un+1 h· ∇)dn h∧dn h,¯ dh´ +γ³wn+1 h∧[dn+1 h+ (dn h−dn+1 h)],¯ dh´+³(Ph(fε(dn h)) −fε(dn h)) ∧dn h,¯ dh´= 0. (11) Por otra parte, queremos obtener la versi´on discreta de la siguiente identidad en continuo −∆d∧d=−∇ · (∇d∧d)− ∇d∧ ∇d=−∇ · (∇d∧d). Ello se obtiene usando P1 h:H1 0(Ω) →D0hel proyector respecto del producto escalar en H1 0(Ω), llegando a ³∇dn+1 h∧dn+1 h,¯ dh´=³wn+1 h∧dn+1 h,¯ dh´+Rn+1 1+Rn+1 2,(12) donde Rn+1 1=³wn+1 h, P1 h[dn+1 h}¯ dh]−[dn+1 h}¯ dh]´, Rn+1 2=³∇e dh,∇[dn+1 h}¯ dh]´, donde }denota el producto vectorial tal que: ³wn+1 h∧dn+1 h,¯ dh´:= ³wn+1 h,dn+1 h}¯ dh´.(13) Luego, usando (12) en (11), multiplicando por ky sumando en (11) para n= 0, . . . , N− 1, conseguimos k N−1 X n=0 (Ãdn+1 h−dn h k∧dn h,¯ dh!+³(un+1 h· ∇)dn h∧dn h,¯ dh´+γ³∇dn+1 h∧dn+1 h,∇¯ dh´) =k N−1 X n=0 (Rn+Rn+1 1+Rn+1 2+Rn+1 3). donde Rn:= ³fε(dn h), Ph(dn h}¯ dh)−dn h}¯ dh´ Rn+1 3:= γ³wn+1 h∧(dn+1 h−dn h),¯ dh´. Teniendo en cuenta las convergencias obtenidas se puede pasar al l´ımite en la formulaci´on d´ebil anterior y mostrar que los t´erminos residuales del segundo miembro van a cero. 7 F. Guill´en-Gonz´alez, J. V. Guti´errez-Santacreu 6. Compacidad para el gradiente de d Definimos zh,k,ε(t) = wn+1 h+Ph(fε(dn h)) y dh,k,ε(t) = dn+1 hsi t∈(tn, tn+1]. Entonces, dh,k,ε(t) es caracterizado por el siguientes problema de optimizaci´on (sin restricciones): Jh,k,ε(dh,k,ε(t)) = m´ın eh∈Dlh Jh,k,ε(eh) (14) donde Jh,k,ε :Dlh→Res definido como Jh,k,ε(eh) = ZΩµ1 2|∇eh|2+Fε(eh)−zh,k,ε(t)·eh¶ y Dlh={dh∈Dh:dh=lhsobre ∂Ω}. Por otro lado, definimos d∗(t) : [0, T]→H1 l(Ω) como la soluci´on del problema de optimizaci´on (con restricciones): J(d∗(t)) = m´ın {e∈Hl,|e|=1|} J(e) (15) donde J(e) = ZΩµ1 2|∇e|2−z(t)·e¶(16) yzes una funci´on limite de zh,k,ε, que existe ya que zh,k,ε est´a acotado en L2(Q) gracias al Teorema 2. Se prueba ([4]) que Jh,k,ε(dh,k,ε) converge a J(d) y que d∗se identifica con el vector de orientaci´on l´ımite dconseguido por las estimaciones de (dh,k,ε). Entonces, tenemos que ∇dh,k,ε converge d´ebil a ∇den L2(L2) y tambi´en converge en norma, de donde deducimos la convergencia fuerte en L2(L2), con la que se puede pasar al l´ımite en el sistema de momentos discreto (3). Agradecimientos Los autores han sido subvencionados por el proyecto MTM2006–07932 Referencias [1] P. Az´ erad, F. Guillen-Gonz´ alez.Mathematical justification of the hidrostatic approximation in the primitive equations of geophysical fluid dynamics. SIAM J. Math. Anal. 33 (2001), no. 4, 847–859. [2] J. J. Douglas, T. Dupont, L. Wahlbin. The stability in Lqof the L2-projection into finite element function spaces. Numer. Math. 23 (1974/75), 193–197. [3] V. Girault, F. Guill´en-Gonz´alez. Mixed formulation, approximation and decoupling algorithm for a nematic liquid crystals model. In preparation. [4] F. Guill´en-Gonz´alez and J.V. Guti´errez-Santacreu. In preparation. [5] F. Guill´en-Gonz´alez and M.A. Rojas-Medar. Global solution of nematic crystals models. C.R.Acad.Sci. Paris, Ser. I 335 (2002) 1085-1090. [6] C. Liu, N.J. Walkington. Mixed methods for the approximation of liquid crystal flows. M2AN Math. Model. Numer. Anal. 36 (2002), no. 2, 205–222. [7] C. Liu, N.J. Walkington. Approximation of liquid crystal flows. SIAM J. Numer. Anal. 37 (2000), no. 3, 725–741. [8] J. Simon. Compact sets in Lp(0, T ;B). Ann. Mat. Pura Appl., s´er. IV, CXLVI (1987), 65–96. 8