Análisis numérico de un problema de contacto viscoelástico en piezoelectricidad

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) An´alisis num´erico de un problema de contacto viscoel´astico en piezoelectricidad M. Barboteu1, J.R. Fern´ andez2, Y. Ouafik1 1Laboratoire de Math´ematiques et Physique pour les Syst`emes (MEPS), Universit´e de Perpignan, Bˆatiment B3, case courrier 12, 52 Avenue Paul Alduy, 66860 Perpignan, Francia. E-mails: [email protected], [email protected]. 2Dpto. de Matem´atica Aplicada, Universidade de Santiago de Compostela, Facultade de Matem´aticas, Campus Sur s/n, 15782 Santiago de Compostela. E-mail: [email protected]. Palabras clave: Piezoelectricidad, Viscoelasticidad, Respuesta normal, Estimaciones del error, Simulaciones num´ericas. Resumen En este trabajo consideramos un problema de contacto sin rozamiento entre un cuerpo electro-viscoel´astico y un obst´aculo deformable. La ley constitutiva electroviscoel´astica se utiliza para modelar el material y, para el contacto, la cl´asica condici´on de contacto de tipo respuesta normal. La formulaci´on variacional de este problema se escribe como un problema acoplado de una ecuaci´on variacional no lineal, de primer orden en tiempo, para el campo de desplazamientos, y una ecuaci´on variacional lineal para el potencial el´ectrico. Entonces, presentamos aproximaciones discretas basadas en el m´etodo de los elementos finitos para aproximar la variable espacial y un esquema de Euler para discretizar la derivada temporal. Demostramos un teorema general de estimaci´on del error del que, bajo condiciones de regularidad adicionales, deducimos la convergencia lineal de estas aproximaciones. Finalmente, presentamos algunos ejemplos bidimensionales que demuestran la eficacia del algoritmo. 1. Introducci´on En este trabajo estudiamos, desde el punto de vista num´erico, un problema de contacto cuasiest´atico sin rozamiento entre un cuerpo viscoel´astico y un obst´aculo deformable, incluyendo los efectos piezoel´ectricos. La piezoelectricidad es la capacidad de ciertos cristales, como el cuarzo (tambi´en algunos materiales cer´amicos (BaTiO3, KNbO3, LiNbO3, etc) o incluso los huesos humanos) para producir corriente el´ectrica cuando se encuentran sometidos a esfuerzos internos. A 1 M. Barboteu, J.R. Fern´andez, Y. Ouafik nivel nanosc´opico, el fen´omeno piezoel´ectrico aparece por una distribuci´on de carga no uniforme en las c´elulas. Cuando el cristal se deforma, las centros de carga positiva y negativa se desplazan debido a la diferencia de potencial, de tal forma que el cristal permanece el´ectricamente neutro. Sin embargo, la diferencia en el desplazamiento de la carga produce una polarizaci´on el´ectrica dentro del cristal. El efecto piezoel´ectrico se caracteriza por el acoplamiento entre las propiedades mec´anicas y el´ectricas del material: se ha observado que la aparici´on de cargas el´ectricas en algunos cristales era debida a la acci´on de fuerzas vol´umicas y/o superficiales y, rec´ıprocamente, la acci´on del potencial el´ectrico provocaba deformaciones o tensiones dentro del cuerpo. Este tipo de materiales aparece muchas veces en la industria como interruptores en electroac´ustica o equipos de medici´on (un ejemplo son los conocidos “thermistor”). En los ´ultimos cincuenta a˜nos se han desarrollado diferentes modelos para describir esta interacci´on entre los campos el´ectricos y mec´anicos (v´eanse, por ejemplo, [1, 5, 7, 10, 11] y las referencias all´ı citadas). Adem´as, recientemente se han estudiado algunos problemas de contacto que inclu´ıan materiales electro-el´asticos ([2, 9]) o electro-viscoel´asticos ([8]). 2. Problema mec´anico y formulaci´on variacional Denotemos por Sdel espacio de tensores de segundo orden sim´etricos en Rdy por “·” yk · k los productos interiores y las normas eucl´ıdeas en ambos espacios. Sea Ω ⊂Rd,d= 1,2,3, un dominio ocupado por un cuerpo electro-viscoel´astico con frontera Γ = ∂Ω. Denotamos por νel vector unitario normal y exterior a Γ y suponemos que esta frontera se descompone en tres partes medibles y disjuntas ΓD, ΓFy ΓC, por un lado, y en dos partes medibles y disjuntas ΓAy ΓB, por el otro lado, tales que med (ΓD)> 0, med (ΓA)>0, y ΓC⊆ΓB. Finalmente, sea [0, T], T > 0, el intervalo temporal donde estudiaremos la evoluci´on del cuerpo (v´ease la Figura 1). f Foundation Γ ΓΓΓ Γ Ω A B f0 C D F qF g ν u = 0 F Figura 1: Un cuerpo electro-viscoel´astico en contacto con un obst´aculo deformable. Denotamos por uel campo de desplazamientos, σel tensor de tensiones, ε(u) = (εij(u))d i,j=1 el tensor de deformaciones linealizado dado por εij(u) = 1 2µ∂ui ∂xj +∂uj ∂xi¶, 2 An´alisis num´erico de un problema de contacto viscoel´astico en piezoelectricidad yϕel potencial el´ectrico. El cuerpo se supone electro-viscoel´astico y satisfaciendo la siguiente ley constitutiva ([3, 8]): σ=Aε(˙ u) + Bε(u)− E∗E(ϕ), donde AyBson los tensores de cuarto orden de viscosidad y elasticidad, respectivamente, E(ϕ) = (Ei(ϕ))d i=1 representa el campo el´ectrico definido por Ei(ϕ) = −∂ϕ ∂xi , i = 1,...,d, yE∗= (e∗ ijk)d i,j,k=1 denota el traspuesto del tensor de tercer orden piezoel´ectrico E= (eijk)d i,j,k=1. Siguiendo [1] se ha utilizado la siguiente ley constitutiva para el potencial el´ectrico: D=Eε(u) + βE(ϕ), donde Des el campo de desplazamientos el´ectricos y βel tensor de permitividad el´ectrica. Sean f0la densidad de fuerzas vol´umicas que act´uan en Ω y q0la densidad de volumen de las cargas el´ectricas libres. En la parte de la frontera ΓDsuponemos que el cuerpo est´a fijado (y as´ı el desplazamiento se anula all´ı), mientras que en ΓFact´ua una densidad de fuerzas superficiales denotada por fF. En la parte ΓCel cuerpo puede entrar en contacto con un obst´aculo deformable. Siguiendo [6], hemos utilizado la siguiente condici´on de contacto (conocida como condici´on de respuesta normal): −σν=p(uν−g) en ΓC×(0, T), donde σν=σν ·νes la tensi´on normal, uν=u·νdenota el desplazamiento normal, g representa la distancia entre el cuerpo y el obst´aculo medida en la direcci´on del vector normal exterior y pes una funci´on constitutiva dada. Finalmente, suponemos que el contacto es sin rozamiento y, por tanto, στ=σ−σνν=0. Para las condiciones frontera del potencial el´ectrico, suponemos que Ω est´a sometido a un potencial el´ectrico ϕAprescrito en ΓAy a una densidad de fuerzas superficiales qF en ΓB. El problema mec´anico se escribe, entonces, de la forma siguiente: Problema P.Encontrar el campo de desplazamientos u: Ω ×(0, T)→Rd, el tensor de tensiones σ: Ω ×(0, T)→Sd, el potencial el´ectrico ϕ: Ω ×(0, T)→Ry el campo de desplazamientos el´ectricos D: Ω ×(0, T)→Rdtales que: σ=Aε(˙ u) + Bε(u)− E∗E(ϕ)en Ω×(0, T), D=Eε(u) + βE(ϕ)en Ω×(0, T), Div σ+f0=0en Ω×(0, T), div D=q0en Ω×(0, T), u=0en ΓD×(0, T), σν =fFen ΓF×(0, T), στ=0,−σν=p(uν−g)en ΓC×(0, T), ϕ=ϕAen ΓA×(0, T), D·ν=qFen ΓB×(0, T), u(0) = u0en Ω, 3 M. Barboteu, J.R. Fern´andez, Y. Ouafik donde u0representa una condici´on inicial para el campo de desplazamientos y, en la ecuaci´on de equilibrio para el campo de desplazamientos, hemos suprimido los efectos de inercia porque consideramos el problema cuasiest´atico. Para poder escribir la formulaci´on variacional del Problema P, definimos los siguientes espacios variacionales V,QyW, y el conjunto convexo WA: V={v∈[H1(Ω)]d;v=0en ΓD}, Q={τ= (τij)d i,j=1 ∈[L2(Ω)]d×d;τij =τji, i, j = 1, . . . , d}, W={ψ∈H1(Ω) ; ψ= 0 en ΓA}, WA={ψ∈H1(Ω) ; ψ=ϕAen ΓA}, y denotamos por H= [L2(Ω)]d. Suponemos que los tensores de viscosidad y permitividad, Ayβ, son sim´etricos, acotados y definidos positivo, que los tensores de elasticidad y piezoelectricidad, ByE, son sim´etricos y acotados, que la funci´on de respuesta normal pes de Lipschitz, positiva y mon´otona, y, finalmente, que las fuerzas vol´umicas y superficiales son continuas en tiempo y que la funci´on distancia ges positiva. Usando el Teorema de Riesz, definimos las siguientes formas lineales: (f(t),w)V=ZΩ f0(t)·wdx+ZΓF fF(t)·wdΓ,∀w∈V, (q(t), ψ)W=ZΩ q0(t)ψ dx+ZΓB qF(t)ψ dΓ,∀ψ∈W. Sea j:V×V→Rel funcional de respuesta normal dado por: j(u,v) = ZΓC p(uν−g)vνdΓ,∀u,v∈V, donde, para todo v∈V, denotamos por vν=v·ν. Aplicando la f´ormula de Green y las condiciones frontera enunciadas en el Problema P, obtenemos la siguiente formulaci´on variacional. Problema PV.Encontrar un campo de desplazamientos u: [0, T]→Vy un potencial el´ectrico ϕ: [0, T]→WAtales que u(0) = u0y para todo t∈[0, T ]: (Aε(˙ u(t)),ε(w))Q+ (Bε(u(t)),ε(w))Q+ (E∇ϕ(t),ε(w))Q+j(u(t),w) = (f(t),w)V,∀w∈V, (β∇ϕ(t),∇ψ)H−(Eε(u(t)),∇ψ)H= (q(t), ψ)W,∀ψ∈W. Bajo las hip´otesis anteriores, utilizando resultados cl´asicos de ecuaciones variacionales no lineales parab´olicas y ecuaciones en derivadas parciales, en [8] se demostr´o que el Problema PV tiene una ´unica soluci´on con la regularidad: u∈C1([0, T]; V), ϕ ∈C([0, T]; WA). 3. Aproximaciones discretas: estimaciones del error Para simplificar, suponemos que ϕA= 0 (y entonces WA=W). La discretizaci´on del Problema PV se realiza en dos pasos. En primer lugar, consideramos dos espacios de 4 An´alisis num´erico de un problema de contacto viscoel´astico en piezoelectricidad dimensi´on finita Vh⊂VyWh⊂Waproximando los espacios variacionales, respectivamente. h > 0 denota el par´ametro de discretizaci´on espacial. En segundo lugar discretizamos las derivadas temporales. Definimos una partici´on uniforme del intervalo [0, T], denotada por 0 = t0< t1< . . . < tN=T, y sea kel paso de tiempo, k=T/N. Para una funci´on continua, denotamos por fn=f(tn). Utilizando un esquema de Euler, obtenemos el siguiente problema discreto. Problema PVhk.Encontrar un campo de desplazamientos discreto uhk ={uhk n}N n=0 ⊂ Vhy un potencial el´ectrico discreto ϕhk ={ϕhk n}N n=0 ⊂Whtales que uhk 0=uh 0y para todo n= 1,...,N: (Aε((uhk n−uhk n−1)/k),ε(wh))Q+ (Bε(uhk n),ε(wh))Q+ (E∗∇ϕhk n,ε(wh))Q +j(uhk n,wh) = (fn,wh)V,∀wh∈Vh, (β∇ϕhk n,∇ψh)H−(Eε(uhk n),∇ψh)H= (qn, ψh)W,∀ψh∈Wh, donde uh 0es una aproximaci´on adecuada de la condici´on inicial u0. Utilizando resultados cl´asicos de ecuaciones variacionales y teniendo en cuenta las hip´otesis anteriores se obtiene f´acilmente que el Problema P V hk tiene una ´unica soluci´on. El inter´es de esta secci´on se centra en obtener estimaciones de los errores num´ericos kun−uhk nkVykϕn−ϕhk nkW. Se tiene el siguiente resultado general de estimaci´on del error. Teorema 1. Bajo las hip´otesis enunciadas en la secci´on anterior, sean (u, ϕ)y(uhk, ϕhk) las soluciones de los problemas PV yPV hk, respectivamente. Entonces se tiene la siguiente estimaci´on del error para todos wh={wh j}N j=1 ⊂Vhyψh={ψh j}N j=1 ⊂Wh: m´ax 1≤n≤N{kun−uhk nk2 V+kϕn−ϕhk nk2 W} ≤ c³m´ax 1≤n≤Nkϕn−ψh nk2 W + N X j=1 khk˙ uj−(uj−uj−1)/kk2 V+kuj−wh jk2 Vi+ m´ax 1≤n≤Nkun−wh nk2 V +ku0−uh 0k2 V+1 k N−1 X j=1 kuj−wh j−(uj+1 −wh j+1)k2 V´. La estimaci´on precedente es la base para el an´alisis del orden de convergencia del algoritmo. Por ejemplo, supongamos que Ω es un dominio poligonal y denotemos por Th una triangulaci´on tipo elementos finitos regular de Ω compatible con la partici´on de la frontera Γ = ∂Ω en ΓD, ΓFy ΓCpor un lado, y en ΓAy ΓB, por el otro. Definimos los espacios de elementos finitos VhyWhcomo: Vh={wh∈[C(Ω)]d;wh |T r ∈[P1(Tr)]dTr ∈ T h,wh=0en ΓD}, Wh={ψh∈C(Ω) ; ψh |T r ∈P1(Tr)Tr ∈ T h, ψh= 0 en ΓA}, y construimos la aproximaci´on discreta de la condici´on inicial u0en la forma uh 0= Πhu0, donde Πh= (πh)d i=1 : [C(Ω)]d→Vh, y πh:C(Ω) →Bhes el operador de interpolaci´on por elementos finitos usual. En estas condiciones, podemos demostrar que bajo las condiciones de regularidad adicionales: u∈H2(0, T;V)∩C([0, T]; [H2(Ω)]d), ϕ ∈C([0, T]; H2(Ω)), 5 M. Barboteu, J.R. Fern´andez, Y. Ouafik el problema PV hk es linealmente convergente; es decir, existe c > 0, independiente de hy k, tal que: m´ax 1≤n≤N{kun−uhk nkV+kϕn−ϕhk nkW} ≤ c(h+k). Observaci´on 1. La resoluci´on num´erica del Problema PV hk se realiza en dos pasos. En primer lugar se utiliza una penalizaci´on para modelar la ley de respuesta normal y, en segundo lugar, se emplea el m´etodo de Newton para resolver el problema no lineal, acoplando el algoritmo con un esquema predictor-corrector. 4. Resultados num´ericos Como un ejemplo bidimensional del Problema P, consideramos el cuerpo Ω con frontera Γ que puede entrar en contacto con un s´olido deformable (v´ease la situaci´on f´ısica descrita en el lado izquierdo de la Figura 2). Para describir la geometr´ıa, fijemos los puntos P1= (0,1), P2= (1,0), P3= (3,1), P4= (1,5,1), P5= (1,1,5) y P6= (1,4). Definimos las fronteras ΓB= ΓC= [P1, P2], ΓD= [P4, P5] y ΓF= Γ\(ΓC∪ΓD); {X2, X3}denota la base ortonormal can´onica. Adem´as, utilizamos las siguientes constantes de viscosidad y elasticidad: Elasticidad (GPa) viscosidad (GPa ·s) b22 b23 b33 b44 a22 a23 a33 a44 210 105 211 42.5 21 10.5 21.1 4.25 Tabla 1: Coeficientes de viscosidad y elasticidad del material. Piezoelectricidad (C/m2) Permitividad (C2/Nm2) e32 e33 e24 β22/²0β33/²0 -0.61 1.14 -0.59 -8.3 -8.8 Tabla 2: Coeficientes piezoel´ectricos del material. Supongamos que el cuerpo est´a fijado en ΓD. En las simulaciones hemos empleado los siguientes datos: fF=0N/m2,f0=0N/m3, q0= 0 C/m3, g = 10−4m, T = 10 s, u0=0m, ϕA=   20 Ven [P3, P4], 0Ven [P5, P6]. Este problema nos permitir´a mostrar que la aparici´on de deformaciones y tensiones en el cuerpo es debida a la acci´on del campo el´ectrico. El intervalo de tiempo se discretiza con una partici´on uniforme (k= 0,001) y el cuerpo se discretiza con el mallado tipo elementos finitos que se puede observar en la Figura 2 (derecha). Como podemos ver en la Figura 3, la acci´on de la diferencia del potencial el´ectrico entre [P3, P4] y [P5, P6] induce una deformaci´on en el cuerpo. Esto hace que este entre 6 An´alisis num´erico de un problema de contacto viscoel´astico en piezoelectricidad 6 x x 2 3 Ω Foundation g ϕ Γ Γ ν ϕ = 0 D C PP P P P P 1 2 34 5 Figura 2: Problema f´ısico (izquierda) y mallado tipo elementos finitos (derecha). Figura 3: Mallado tipo elementos finitos, en los instantes inicial y final, y fuerzas de contacto. en contacto con el obst´aculo en ΓC. Adem´as, observamos que algunos nodos de contacto est´an en deslizamiento y las fuerzas de contacto se orientan en la normal saliente a ΓC (ya que no hay rozamiento). Adem´as, la Figura 4 muestra la distribuci´on del campo de desplazamientos el´ectricos y el potencial el´ectrico (izquierda) y la norma von Mises de las tensiones en la configuraci´on deformada del cuerpo. Podemos apreciar que existe una cierta correspondencia entre la distribuci´on de estos desplazamientos y la norma von Mises de las tensiones. Esto es debido a que los mayores valores del campo el´ectrico se localizan en las zonas donde las tensiones son mayores. 7 M. Barboteu, J.R. Fern´andez, Y. Ouafik 22.00 19.74 18.05 16.36 14.67 12.97 11.28 9.590 7.897 6.205 4.513 2.821 1.128 0.0000E+00 32.53 29.59 27.46 25.34 23.21 21.09 18.96 16.84 14.71 12.59 10.46 8.336 6.211 4.086 1.962 0.0000E+00 Figura 4: Campo de desplazamientos el´ectricos (flechas) y potencial el´ectrico (isovalores) en la configuraci´on deformada (izquierda) y norma von Mises del campo de tensiones (derecha). Agradecimientos El trabajo de J.R. Fern´andez fue parcialmente financiado por el Ministerio de Educaci´on y Ciencia (Proyecto MTM2006-13981). Referencias [1] R.C. Batra y J.S. Yang. Saint-Venant’s principle in linear piezoelectricity. J. Elasticity, 38 (1995), 209–218. [2] P. Bisegna, F. Lebon y F. Maceri. The unilateral frictional contact of a piezoelectric body with a rigid support. Contact mechanics (Praia da Consola¸c˜ao, 2001), 347–354, Solid Mech. Appl., 103, Kluwer Acad. Publ., Dordrecht, 2002. [3] G. Duvaut y J.L. Lions. Inequalities in mechanics and physics. Springer Verlag, Berlin, 1976. [4] S. H¨ueber, A. Matei y B.I. Wohlmuth. A mixed variational formulation and an optimal a priori error estimate for a frictional contact problem in elasto-piezoelectricity. Bull. Math. Soc. Sci. Math. Roumanie, 48(96) (2005), 209–232. [5] T. Ideka, Fundamentals of piezoelectricity, Oxford University Press, Oxford, 1990. [6] A. Klarbring, A. Mikeli´c y M. Shillor. Frictional contact problems with normal compliance. Internat. J. Engrg. Sci., 26 (1988), 811–832. [7] R.D. Mindlin. Elasticity, piezoelasticity and crystal lattice dynamics. J. Elasticity, 4(1972), 217–280. [8] M. Sofonea y E.-H. Essoufi. Quasistatic frictional contact of a viscoelastic piezoelectric body. Adv. Math. Sci. Appl., 14(1) (2004), 25–40. [9] M. Sofonea y Y. Ouafik. A piezoelectric contact problem with normal compliance. Appl. Math. 32 (2005), 425–442. [10] R.A. Toupin. The elastic dielectrics. J. Rational Mech. Anal., 5(1956), 849–915. [11] N. Turb´e y G.A. Maugin. On the linear piezoelectricity of composite materials. Math. Methods Appl. Sci., 14 (6) (1991), 403–412. 8