Análisis numérico de un modelo de remodelación ósea

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 modelo de remodelaci´on ´osea J.R. Fern´ andez, R. Mart´ ınez, J.M. Via˜ no Departamento de Matem´atica Aplicada, Universidade de Santiago de Compostela, Facultade de Matem´aticas, Campus Sur s/n, 15782 Santiago de Compostela. E-mails: [email protected], [email protected], [email protected]. Palabras clave: Remodelaci´on ´osea, aproximaciones discretas, estimaciones del error, soluciones num´ericas. Resumen En este trabajo se estudia, desde el punto de vista num´erico, un modelo de remodelaci´on ´osea. Este modelo se caracteriza por una ecuaci´on variacional el´ıptica para el campo de desplazamientos y una ecuaci´on diferencial ordinaria de primer orden en tiempo para describir el proceso fisiol´ogico de remodelado ´oseo. Utilizando el m´etodo de elementos finitos para aproximar la variable espacial y un esquema de Euler para discretizar las derivadas temporales, obtenemos aproximaciones discretas de este problema variacional y probamos un resultado de estimaci´on del error. Bajo condiciones de regularidad adecuadas, deducimos la convergencia lineal del algoritmo respecto de los par´ametros de discretizaci´on. Finalmente, presentamos algunos resultados num´ericos, en un ejemplo bidimensional, para mostrar la validez del algoritmo. 1. Formulaci´on mec´anica y variacional Denotemos por ·el producto interior en Rdy su correspondiente norma por |·|. Sea Sdel espacio de los tensores sim´etricos de segundo orden en Rd, o equivalentemente, el espacio de las matrices sim´etricas de orden d, y sea : su producto interior y |·|su norma. Sea Ω ⊂Rd,d= 1,2,3, un dominio abierto y acotado y sea Γ = ∂Ω su frontera, que suponemos Lipschitz continua y dividida en dos partes disjuntas ΓDy ΓN. El cuerpo est´a sometido a una fuerza vol´umica de densidad f, fijado en ΓDy sometido a una fuerza de tracci´on superficial de densidad gactuando en ΓN. Finalmente, denotemos por νel vector normal exterior unitario (v´ease la Figura 1). Denotamos por x= (xi)d i=1 un punto gen´erico de Ω y por t∈[0, T], T > 0, un instante cualquiera en el intervalo de tiempo [0, T]. Sea u=u(x, t) el campo de desplazamientos en el punto xy en el instante t,ε(u) = (εij(u))d i,j=1 el tensor de deformaciones linealizado dado por εij(u) = 1 2µ∂ui ∂xj +∂uj ∂xi¶. 1 J.R. Fern´andez, R. Mart´ınez, J.M. Via˜no Ω Γ D g Γ N ν Γ N Γ N f Figura 1: Problema de remodelaci´on ´osea. Denotamos por e=e(x, t) la medida del cambio en la fracci´on vol´umica del hueso desde la configuraci´on de referencia y que denominaremos funci´on de remodelaci´on ´osea. Denotamos por σ=σ(x, t) el tensor de tensiones en el punto xy en el instante t. El modelo de remodelaci´on ´osea que consideramos aqu´ı fue postulado en [2, 5]. Se trata de un modelo no lineal con una ley constitutiva que se escribe como σ= (ξ0+e)C(e)ε(u), donde ξ0representa la fracci´on vol´umica de referencia y C(e) = (Cijkl(e))d i,j,k,l=1 los coeficientes de elasticidad que dependen de e. N´otese que si ξ0= 1 y e= 0, la ley constitutiva corresponde a la ley de Hooke cl´asica. La evoluci´on de la funci´on de remodelaci´on ´osea se obtiene de la siguiente ecuaci´on diferencial ordinaria donde ˙edenota la derivada temporal: ˙e=a(e) + A(e) : ε(u), donde a(e) es una funci´on constitutiva y A(e) = (Aij(e))d i,j=1 incluye los coeficientes de remodelaci´on ´osea. Observaci´on 1. Notemos que las funciones C(e),a(e)yA(e)caracterizan las propiedades del material y los datos experimentales determinan su forma. En algunos art´ıculos, se emplean aproximaciones polin´omicas (v´ease, por ejemplo, [6]): Cijkl(e) = 1 ξ0+e(ξ0C0 ijkl +C1 ijkle), a(e) = a0+a1e+a2e2, Aij(e) = A0 ij +A1 ije, donde C0 ijkl,C1 ijkl,a0,a1,a2,A0 ij and A1 ij son constantes que dependen de las propiedades del material. Definimos el siguiente operador de truncamiento ΦL:R→[−L, L] por ΦL(r) =      −Lsi −L≤r, rsi −L≤r≤L, Lsi r≥L. Finalmente, suponemos que el proceso es cuasiest´atico y, en consecuencia, los efectos de inercia son ignorados. Adem´as, denotamos por e0la funci´on de remodelaci´on ´osea inicial. 2 An´alisis num´erico de un modelo de remodelaci´on ´osea El problema mec´anico, obtenido de las leyes de la mec´anica de medios continuos en el caso de peque˜nas deformaciones, es el siguiente (v´ease [4, 6]): Problema P. Encontrar el campo de desplazamientos u: Ω ×(0, T)→Rd, el tensor de tensiones σ: Ω ×(0, T)→Sdy la funci´on de remodelaci´on ´osea e: Ω ×(0, T)→Rtal que e(0) = e0y para todo t∈(0, T): σ(t) = (ξ0+e(t))C(e(t))ε(u(t)) en Ω×(0, T), ˙e(t) = a(e(t)) + A(e(t)) : ε(u(t)) en Ω×(0, T), −Div σ(t) = γ(ξ0+ ΦL(e(t)))f(t)en Ω×(0, T), u(t) = 0en ΓD×(0, T), σ(t)ν=g(t)en ΓN×(0, T). donde γ > 0 es la densidad del material que ocupa ¯ Ω (se considera constante para simplificar). Intentemos ahora obtener una formulaci´on variacional del problema P. Denotemos por Y=L2(Ω) y H= [L2(Ω)]dy definimos los siguientes espacios variacionales: V={w∈[H1(Ω)]d;w=0en ΓD}, Q={τ= (τij)d i,j=1 ∈[L2(Ω)]d×d;τij =τji,1≤i, j ≤d}. Suponemos las siguientes hip´otesis sobre los datos: los coeficientes de elasticidad Cijkl(e) son continuamente diferenciables con respecto aey satisfacen las siguientes propiedades: (a) Cijkl(e) = Cjikl(e) = Cklij (e) para i, j, k, l = 1, . . . , d. (b) Existe mC>0 tal que (ξ0+e)C(e)τ:τ≥mC|τ|2.(1) La funci´on constitutiva a(e) y los coeficientes de remodelaci´on ´osea Aij(e) son continuamente diferenciables con respecto a e. La fracci´on vol´umica de referencia ξ0satisface las siguientes condiciones para alg´un ξm 0, ξM 0<1 : ξ0∈C1(Ω),0< ξm 0≤ξ0(x)≤ξM 0<1 para todo x∈Ω.(2) Las fuerzas de densidad tienen la siguiente regularidad: f∈C([0, T]; [C(Ω)]d),g∈C([0, T]; [C(ΓN)]d),(3) y el valor inicial de la funci´on de remodelaci´on e0verifica que e0∈C1(Ω).(4) 3 J.R. Fern´andez, R. Mart´ınez, J.M. Via˜no Dado e∈L∞(Ω), definimos la siguiente forma bilineal c(e;·,·) : V×V→R: c(e;u,v) = ZΩ (ξ0+e)C(e)ε(u) : ε(v)dx∀u,v∈V, y la forma lineal L(e;·) : V→Rdada por: L(e;v) = ZΩ γ(ξ0+ ΦL(e))f·vdx+ZΓN g·vdΓ∀v∈V. Mediante un procedimiento habitual, obtenemos la siguiente formulaci´on variacional del problema mec´anico P. Problema PV. Encontrar el campo de desplazamientos u: [0, T]→Vy la funci´on de remodelaci´on ´osea e: [0, T]→L∞(Ω) tales que: c(e(t); u(t),v) = L(e(t); v)∀v∈V, c.p.t. t ∈(0, T), ˙e(t) = a(e(t)) + A(e(t)) : ε(u(t)), en D0(0, T;L2(Ω)), e(0) = e0. En [7] se demuestra el siguiente resultado de existencia y unicidad de soluci´on del problema PV. Teorema 1. Supongamos que se cumplen las hip´otesis (1)-(4) Entonces, existe una ´unica soluci´on del problema PV con la siguiente regularidad: u∈C1([0, T]; [H2(Ω)]d), e ∈C1([0, T]; C(Ω)). 2. Aproximaciones num´ericas: estimaciones del error Introducimos ahora un algoritmo basado en el m´etodo de elementos finitos para aproximar las soluciones del problema PV y obtener una estimaci´on del error a partir de ellas. La discretizaci´on del problema PV se hace como sigue. Primero, consideramos dos espacios de dimensi´on finita Vh⊂VyBh⊂L∞(Ω) ⊂Y, aproximaci´on de los espacios V yL∞(Ω), respectivamente. Observaci´on 2. En las simulaciones num´ericas que se presentan en la siguiente secci´on, VhyBhson espacios de funciones continuas afines a trozos y funciones constantes a trozos, respectivamente; es decir: Vh={wh∈[C(Ω)]d;wh |T r ∈[P1(Tr)]d, Tr ∈ T h,wh=0en ΓD},(5) Bh={ξh∈L∞(Ω) ; ξh |T r ∈P0(Tr), Tr ∈ T h},(6) donde Ωes un dominio poligonal, Thdenota una triangulaci´on tipo elementos finitos de Ω, y Pq(Tr),q= 0,1, representa el espacio de los polinomios de grado menor o igual que qen Tr. 4 An´alisis num´erico de un modelo de remodelaci´on ´osea Para discretizar las derivadas en tiempo, utilizamos una partici´on uniforme de [0, T], denotada por 0 = t0< t1< . . . < tN=T, y sea kel tama˜no del paso de tiempo, k=T/N. Adem´as, para una funci´on continua f(t) sea fn=f(tn). En esta secci´on, cdenota una constante positiva que depende de los datos del problema, pero es independiente de los par´ametros de discretizaci´on hyk. Utilizando el esquema de Euler expl´ıcito, la discretizaci´on del problema variacional PV resulta de la forma siguiente. Problema PVhk.Encontrar el campo de desplazamientos uhk ={uhk n}N n=0 ⊂Vhy una funci´on de remodelaci´on ´osea discreta ehk ={ehk n}N n=0 ⊂Bhtales que ehk 0=eh 0, uhk 0=uh 0son adecuadas aproximaciones de e0yu0respectivamente y para n= 1, . . . , N: uhk n∈V, c(ehk n;uhk n,vh) = L(ehk n;vh)∀vh∈Vh, ehk n−ehk n−1 k=a(ehk n−1) + A(ehk n−1)) : ε(uhk n−1). De las propiedades (1) es inmediato obtener la existencia y unicidad de la soluci´on discreta como enunciamos a continuaci´on. Teorema 2. Supongamos que se verifican las hip´otesis (1)-(4). Entonces, el problema PVhk tiene una ´unica soluci´on (uhk, ehk)⊂Vh×Bh. El objetivo de esta secci´on es obtener estimaciones del error para kun−uhk nkVy ken−ehk nkY. Tenemos el siguiente resultado para la estimaci´on del error principal. Teorema 3. Supongamos que se cumplen las hip´otesis (1)-(4). Sean (u, e)y(uhk, ehk) las respectivas soluciones de los problemas PV y PVhk. Bajo la condici´on de regularidad u∈C([0, T]; [W1,∞(Ω)]d),tenemos, para todo {vh n}N n=0 ⊂Vh: m´ax 0≤n≤N{ken−ehk nk2 Y+kun−uhk nk2 V} ≤ ke0−eh 0k2 Y+c³k N X j=1 hk˙ej−(ej−ej−1)/kk2 Y +kuj−uj−1k2 Vi+k2+ m´ax 1≤n≤Nkun−vh nk2 V+ku0−uhk 0k2 V´.(7) Las estimaciones del error (7) son b´asicas para el an´alisis de la convergencia del algoritmo. A continuaci´on presentamos un ejemplo. Sea Ω un dominio poli´edrico y denotemos por Thuna triangulaci´on de Ω compatible con la partici´on de la frontera Γ = ∂Ω en ΓDy ΓN. Sean VhyBhdefinidos por (5) y (6), respectivamente, y supongamos que la condici´on inicial discreta eh 0se obtiene por eh 0=πhe0,donde πh:C(Ω) →Bhes el operador de interpolaci´on de elementos finitos (v´ease, e.g., [1]) y uh 0=πhu0donde Πh= (πh i)d i=1 : [C(Ω)]d→Vh. Supongamos la siguiente condici´on de regularidad en la soluci´on continua: e∈C([0, T]; H2(Ω)) ∩H2(0, T;Y).(8) Del Teorema 3 se tiene que u∈C([0, T]; [H2(Ω)]d) (por tanto, u0=u(0) ∈[H2(Ω)]d) El siguiente resultado se sigue de la estimaci´on (7). Corolario 1. Bajo las hip´otesis del Teorema 3 y la condici´on de regularidad (8), el esquema discretizado basado en el m´etodo de elementos finitos descrito anteriormente es 5 J.R. Fern´andez, R. Mart´ınez, J.M. Via˜no linealmente convergente, es decir, existe una constante positiva c, independiente de hyk, tal que m´ax 0≤n≤Nnkun−uhk nkV+ken−ehk nkYo≤c(h+k). 3. Resultados num´ericos En esta secci´on describimos brevemente el esquema num´erico que hemos implementado, y presentamos un ejemplo bidimensional para mostrar su comportamiento. 3.1. Esquema num´erico Para aproximar los espacios VyL∞(Ω) utilizamos los espacios de elementos finitos VhyBhdefinidos por (5) y (6), respectivamente. En primer lugar, se˜nalemos que, en la pr´actica, uhk 0se obtiene resolviendo el siguiente problema variacional: uhk 0∈Vh, c(eh 0;uhk 0,vh) = L(eh 0;vh),∀vh∈Vh. Se trata de un problema de elasticidad lineal discreto equivalente a un sistema lineal que se resuelve utilizando el m´etodo de Cholesky. Para n∈ {1, . . . , N}se tiene que uhk n−1yehk n−1son conocidos. La funci´on de remodelaci´on ´osea discreta ehk nse calcula expl´ıcitamente como: ehk n=ehk n−1+ka(ehk n−1) + kA(ehk n−1) : ε(uhk n−1). Llevando esto a la ecuaci´on (7), el campo de desplazamientos discreto uhk nse obtiene resolviendo la siguiente ecuaci´on variacional (problema de elasticidad lineal): uhk n∈Vhc(ehk n;uhk n,vh) = L(ehk n;vh),∀vh∈Vh. De nuevo, esto nos lleva a un sistema lineal que resolvemos por el m´etodo de Cholesky. El esquema num´erico fue implementado en lenguaje MATLAB en un PC con 3.2Ghz. Una ejecuci´on de un ejemplo 1D necesita alrededor de 3.5 segundos de tiempo CPU y una ejecuci´on de un ejemplo 2D alrededor de 10 minutos. 3.2. Resultados num´ericos en un problema bidimensional Como ejemplo bidimensional, consideramos el dominio Ω = (0,1,2)×(0,6) que est´a siendo sometido a una fuerza de compresi´on creciente de forma lineal actuando en la frontera horizontal superior mientras que la frontera horizontal inferior permanece fija. No hay fuerzas de volumen en Ω (v´ease la Figura 2). Se han utilizado los siguientes datos en las simulaciones num´ericas de este ejemplo: T= 80 d´ıas, k = 0,01,C(e) = 1 ξ0+e(C0+C1e), a(e) = a0+a1e+a2e2, A(e) = A0+A1e, ξ0= 0,892, γ = 1740 Kg/m3,f=0N/m3, a0=−1296 ×10−4(100 d´ıas)−1, a1=−1296 ×10−2(100 d´ıas)−1, a2= 216 ×10−2(100 d´ıas)−1, 6 An´alisis num´erico de un modelo de remodelaci´on ´osea Ω ΓN g ΓD Figura 2: Descripci´on f´ısica. donde los tensores de cuarto orden C0= (C0 ijkl)2 i,j,k,l=1 yC1= (C1 ijkl)2 i,j,k,l=1 y el tensor de segundo orden A0= (A0 ij)2 i,j=1 yA1= (A1 ij)2 i,j=1 tienen las siguientes componentes: C0 1111 = 25,69 (GPa)−1, C0 2211 = 11,67 (GPa)−1, C0 2222 = 25,69 (GPa)−1, C0 1211 =C0 1222 = 0 (GPa)−1, C0 1212 = 7 (GPa)−1, C1 1111 = 252,08 (GPa)−1, C1 2211 = 114,58 (GPa)−1, C1 2222 = 252,08 (GPa)−1, C0 1211 =C0 1222 = 0 (GPa)−1, C1 1212 = 68,75 (GPa)−1, A0 11 = 216 (100d´ıas)−1, A0 22 =−216 (100d´ıas)−1, A0 12 =A0 21 = 0, A1 11 = 216 (100d´ıas)−1, A1 22 = 216 (100d´ıas)−1, A0 12 =A0 21 = 0. La fuerza aplicada tiene la forma: g(x, y, t) = −28 (0, x)MPa si y= 6 y tomamos como funci´on de remodelaci´on ´osea inicial e0= 0. Utilizando como par´ametro de discretizaci´on temporal k= 0,01, en las Figuras 3 y 4 se muestran los desplazamientos (multiplicados por 20) y la funci´on de remodelaci´on ´osea en el instante final. Como podemos observar, el desplazamiento decrece a causa de la remodelaci´on ´osea. Adem´as, esta funci´on toma valores positivos donde el cuerpo est´a sometido a una compresi´on y valores negativos donde se produce una extensi´on. Agradecimientos Este trabajo fue parcialmente financiado por el Ministerio de Educaci´on y Ciencia (Proyecto MTM2006-13981). Referencias [1] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, en Handbook of Numerical Analysis, Volumen II, Parte 1, Eds. P.G. Ciarlet and J.L. Lions, North Holland (1991), 17–352. [2] S.C. Cowin y D.H. Hegedus, Bone remodeling I: theory of adaptive elasticity, J. Elasticity 6 (3) (1976) 313–326. [3] S.C. Cowin y R.R. Nachlinger,Bone remodeling III: uniqueness and stability in adaptive elasticity theory, J. Elasticity 8 (3) (1978) 285–295. 7 J.R. Fern´andez, R. Mart´ınez, J.M. Via˜no Figura 3: Configuraci´on inicial y desplazamientos (x20) en el instante inicial (izquierda) y t=80 (derecha). −0.01 −0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 Figura 4: Funci´on de remodelaci´on ´osea en el instante final (t=80). [4] M.A. Fern´andez, Resoluci´on num´erica de modelos de elasticidad adaptativa en formaci´on de huesos, Tesina de licenciatura, Universidad de Santiago, 1999. [5] I.N. Figueiredo, Approximation of bone remodeling models, J. Math. Pures Appl. 84 (2005) 1794– 1812. [6] D.H. Hegedus y S.C. Cowin, Bone remodeling II: small strain adaptive elasticity, J. Elasticity 6(4) (1976) 337–352. [7] J. Monnier y L. Trabucho, Existence and uniqueness of a solution to an adaptive elasticity model, Math. Mech. Solids 3 (1998) 217–228. 8