Un método paralelo para flujos de partículas rígidas basado en dominios ficticios

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Un m´etodo paralelo para flujos de part´ıculas r´ıgidas basado en dominios ficticios Jordi Blasco1, MaCarmen Calzada2, Mercedes Mar´ ın2 1Dpt. de Matem`atica Aplicada I, Univ. Polit`ecnica de Catalunya, Avgda. Diagonal 647, 08028 Barcelona. E-mail: [email protected]. 2Dpto. de Inform´atica y An´alisis Num´erico, Universidad C´ordoba, Campus de Rabanales, Ed. C2-3, E-14071 C´ordoba. E-mails: [email protected], [email protected]. Palabras clave: Flujo de part´ıculas, Dominios ficticios, Ecuaci´on de Navier-Stokes, Algoritmo paralelo Resumen En esta comunicaci´on presentamos un esquema num´erico paralelo para la simulaci´on del flujo de part´ıculas r´ıgidas en un fluido incompresible viscoso. El m´etodo de dominios ficticios permite resolver las ecuaciones de flujo sobre una malla fija; se emplea aqu´ı, adem´as, una t´ecnica de bajo coste computacional para imponer la codici´on de rigidez sobre las part´ıculas basada en la conservaci´on de los momentos lineal y angular. La paralelizaci´on de los c´alculos se consigue mediante el m´etodo SDI (Simultaneous Directions Implicit). Se presentan tambi´en algunos resultados num´ericos que demuestran el buen comportamiento del m´etodo propuesto. 1 Introducci´on En esta comunicaci´on se presenta un m´etodo num´erico paralelo para la simulaci´on de flujos de part´ıculas r´ıgidas en un fluido incompresible viscoso newtoniano. El movimiento del fluido est´a regido por la ecuaci´on de Navier-Stokes incompresible no estacionaria en el dominio fluido, mientras que el movimiento de las part´ıculas est´a gobernado por las ecuaciones de s´olido r´ıgido. Se han desarrollado recientemente diversos m´etodos para la simulaci´on de flujos de part´ıculas r´ıgidas, algunos de los cuales se emplean en problemas de animaci´on en el campo de la inform´atica gr´afica (ver [3], por ejemplo). Muchos de ellos est´an basados en el m´etodo de dominios ficticios ([7], [9]), que permite el c´alculo del movimiento del fluido en una malla fija sin necesidad de remallar. La idea fundamental de este m´etodo consiste en resolver la ecuaci´on de Navier-Stokes en todo el dominio en cada paso de 1 Jordi Blasco, MaCarmen Calzada, Mercedes Mar´ın tiempo, suponiendo en esa etapa que todo el dominio est´a ocupado por fluido; se impone a continuaci´on la condici´on de s´olido r´ıgido sobre las part´ıculas. Los esquemas resultantes pueden verse como m´etodos de paso fraccionado, consistiendo uno de los pasos en la proyecci´on de un campo de velocidades intermedio sobre el movimiento de s´olido r´ıgido. La proyecci´on de la velocidad intermedia sobre el movimiento de s´olido r´ıgido se realiza habitualmente mediante el m´etodo de multiplicadores de Lagrange distribuidos ([7], [8], [9], [10]) asociados a la condici´on de rigidez sobre las part´ıculas; la idea es an´aloga al papel de la presi´on en problemas de flujo incompresible como multiplicador de Lagrange asociado a la condici´on de incompresibilidad. Por otra parte, en [11], se introduce un m´etodo basado en la conservaci´on de los momentos lineal y angular para calcular dicha proyecci´on que resulta ser m´as r´apido y eficiente que los multiplicadores de Lagrange, reduciendo dicha proyecci´on al c´alculo de dos integrales sobre cada part´ıcula. Emplearemos este ´ultimo m´etodo en nuestro algoritmo. Se han empleado tambi´en diferentes m´etodos de resoluci´on de la ecuaci´on de Navier- Stokes para el c´alculo del movimiento del fluido, la mayor´ıa de ellos basados en esquemas de paso fraccionado secuenciales junto con diferentes aproximaciones espaciales: elementos finitos mixtos ([7], [9]), diferencias finitas ([3]), vol´umenes de control ([11]), etc... Nosotros consideramos un m´etodo de paso fraccionado paralelo estudiado en [1], [2] y [5] que permite separar la incompresibilidad de la convecci´on resolviendo simult´aneamente un problema de Burguers no lineal y un problema de Stokes. Adem´as, tanto la resoluci´on del problema no lineal como la del problema de Stokes se efect´uan mediante esquemas iterativos en los que cada iteraci´on consiste en la resoluci´on de problemas de Dirichlet para la ecuaci´on de Helmholz. El m´etodo SDI (Simultaneous Directions Implicit, ver [2], [4], [5], [6]) permite la resoluci´on de dicho problema en paralelo, calculando simult´aneamente la soluci´on en cada unos de los llamados segmentos integrables horizontales y verticales de manera independiente unos de otros. Este esquema espacial paralelo proporciona un alto nivel de escalabilidad al algoritmo resultante. En el Apartado 2 se plantean las ecuaciones del movimiento tanto del fluido como de las part´ıculas. En el Apartado 3 se describe el m´etodo num´erico empleado para su resoluci´on, mientras que en el Apartado 4 se detalla la paralelizaci´on del mismo. Seguidamente, en el Apartado 5 se presentan algunos resultados num´ericos obtenidos con este algoritmo en diferentes casos test; finalmente, en el Apartado 6 se extraen algunas conclusiones de este trabajo y se plantean posibles l´ıneas futuras de investigaci´on. 2 Ecuaciones de movimiento 2.1 Movimiento del fluido Consideramos un dominio abierto y acotado Ω ⊂IR2con frontera poligonal ∂Ω. Supondremos que Ω est´a ocupado por un fluido incompresible, viscoso, newtoniano y homog´eneo, y por Nppart´ıculas r´ıgidas suspendidas en ´el. Llamaremos Pi(t)⊂Ω al espacio ocupado por la part´ıcula i-´esima (i= 1, . . . , Np) en el instante de tiempo t∈[0, T] (T > 0 es el tiempo final de simulaci´on) y sea P(t) = ∪i=1,...,NpPi(t). El movimiento del fluido est´a descrito por un sistema de ecuaciones en derivadas parciales conocidas como ecuaciones 2 M´etodo paralelo para flujos de part´ıculas de Navier-Stokes: ρf∂u ∂t + (u· ∇)u−µ∆u+∇P =ρfgen (Ω −P(t)), t ∈(0, T) (1) ∇ · u= 0 en (Ω −P(t)), t ∈(0, T ).(2) En ellas, ρf>0 es la densidad del fluido, u(·, t) : (Ω−P(t)) →IR2el campo de velocidades en el instante de tiempo t∈(0, T), P(·, t) : (Ω −P(t)) →IR la presi´on en el instante t∈(0, T), µ > 0 la viscosidad din´amica y g= (0,−g) la aceleraci´on de la gravedad. Llamaremos p=P/ρfa la presi´on cinem´atica y ν=µ/ρfa la viscosidad cinem´atica. Estas ecuaciones se completan con condiciones de contorno apropiadas seg´un sea la naturaleza del problema a resolver. De momento, consideraremos s´olo una condici´on de contorno de tipo Dirichlet homog´enea, es decir, condici´on de no deslizamiento, en la frontera de Ω: u=0en ∂Ω×(0, T).(3) En la interfase fluido-part´ıculas, deben igualarse las velocidades de ambos medios e imponerse la fuerza de tracci´on ejercida por las part´ıculas sobre el fluido: u=uien ∂Pi(t), t ∈(0, T), i = 1, . . . , Np(4) σ·n=tien ∂Pi(t), t ∈(0, T), i = 1, . . . , Np,(5) siendo uila velocidad de las particulas en su superficie, σel tensor de tensiones del fluido, nla normal exterior (al fluido) en la superficie de las part´ıculas y tila tracci´on. Se impone tambi´en una condici´on inicial sobre la velocidad: u(x,0) = u0(x),verificando ∇ · u0= 0 en (Ω −P(0)).(6) 2.2 Movimiento de las part´ıculas El movimiento de las part´ıculas, consideradas s´olidos r´ıgidos homog´eneos, est´a gobernado por las ecuaciones de Newton y Euler. Llamando Xi∈IR2a la posici´on del centro de masas de la part´ıcula i-´esima, θia su posici´on angular respecto a un eje fijo, Ui∈IR2 a la velocidad lineal de su centro de masas y ωia su velocidad angular, y escribiendo ωi= (0,0, ωi), se tiene: dXi dt =Ui, Mi dUi dt =Fh i+Fp i+Fw i+ (ρs,i −ρf)Aig,(7) dθi dt =ωi, Ii dωi dt +ωi×(Iiωi) = Ti,(8) siendo Mila masa de la part´ıcula i,Iisu tensor de inercia, Fp ila fuerza de repulsi´on ejercida por las otras part´ıculas sobre la i,Fw ila fuerza de repulsi´on ejercida por las paredes de Ω, ρs,i la densidad (superficial) de masa, Aiel ´area y Fh iyTila fuerza hidrodin´amica y el impulso angular, respectivamente, ejercidos por el fluido sobre la part´ıcula i: Fh i=−Z∂Pi σ·nd∂Pi,Ti=−Z∂Pi r×(σ·n)d∂Pi. 3 Jordi Blasco, MaCarmen Calzada, Mercedes Mar´ın Aqu´ı, res el vector posici´on desde el centro de masas de la part´ıcula (los vectores de IR2 que aparecen en productos vectoriales se suponen extendidos por 0 a IR3). Supondremos en lo sucesivo que las part´ıculas son circulares, de manera que no es necesario resolver (8) para determinar su posici´on. Llamaremos X,Uyωa los vectores que contienen las posiciones, velocidades y velocidades angulares de todas las part´ıculas. Las fuerzas de repulsi´on entre part´ıculas se calculan de la siguiente manera (ver [8]): Fp i= Np X j=1,j6=i Fp i,j,Fp i,j =   0, di,j > Ri+Rj+δ πR2 iρs,ig pXi−Xj di,j Ri+Rj+δ−di,j δ2 , di,j ≤Ri+Rj+δ, donde di,j es la distancia entre los centros de masas de las part´ıculas iyj,RiyRjsus radios, y δypdos par´ametros algor´ıtmicos. La fuerza de repulsi´on de las paredes se calcula, an´alogamente, de la siguiente manera: Fw i= Nw X j=1 Fw i,j,Fw i,j =     0, d0 i,j >2Ri+δ πR2 iρs,ig wXi−X0 i,j d0 i,j 2Ri+δ−d0 i,j δ2 , d0 i,j ≤2Ri+δ, donde Nwes el n´umero de lados de ∂Ω, d0 i,j es la distancia entre el centro de masas de la part´ıcula iy el de una part´ıcula ficticia de radio Ritangente al lado jyX0 i,j es la posici´on del centro de masas de esta part´ıcula ficticia. 3 Resoluci´on num´erica Se describe a continuaci´on el m´etodo num´erico empleado para la resoluci´on del problema (1)-(8). Se considera, en primer lugar, una integraci´on temporal mediante un m´etodo paralelo de paso fraccionado. As´ı pues, suponemos conocidos X0=X(0), U0=U(0), ω0=ω(0) y u0=u0(en Ω −P(0)); calculamos P0, aproximaci´on de P(0), y ponemos u0=U0 i+ω0 i×ren cada part´ıcula P0 i. Entonces, dado un tama˜no del paso de tiempo ∆t > 0, en cada instante de tiempo tn=n∆t(n= 0, . . . , [T/∆t]−1) y partiendo de Xn, Un,un(en Ω) y pn, calculamos Xn+1,Un+1,ωn+1,un+1 ypn+1 mediante el siguiente algoritmo: 1. Se mueven las particulas: Xn+1 =L(Xn, Un) (ver m´as adelante la definici´on del operador L). Se calcula tambi´en Pn+1, aproximaci´on de P(tn+1). 2. Se avanza el fluido suponiendo que llena todo Ω: ρ ρf∂˜u ∂t + (˜u · ∇)˜u −ν∆˜u +∇p=ρ ρf gen Ω, t ∈(tn, tn+1) ∇ · ˜u = 0 en Ω, t ∈(tn, tn+1) ˜u =0sobre ∂Ω, t ∈(tn, tn+1) siendo ρ=ρf(1 −Hn+1) + ρsHn+1, con Hn+1 = 1 en Pn+1,Hn+1 = 0 en Ω −Pn+1 y ρs(t) = ρs,i en Pi(t). N´otese que se ha dividido la ecuaci´on de Navier-Stokes (1) por la densidad del fluido ρf. Este problema se resuelve mediante el siguiente m´etodo paralelo 4 M´etodo paralelo para flujos de part´ıculas de paso fraccionado (ver [1]): se calculan ˜un+1 ypn+1 resolviendo simult´aneamente el problema no lineal (de Burgers) para un+a ρ ρfun+a−un a∆t+2θ a(un+a· ∇)un+a−νσ∆un+a=ρ ρf g+ν(1 −σ)∆un− ∇pn(9) un+a=0sobre ∂Ω,(10) y el problema lineal (de Stokes generalizado) para un+bypn+b ρ ρf un+b−un b∆t−ν(1 −σ)∆un+b+∇pn+b=ρ ρf g+νσ∆un−ρ ρf 2(1 −θ) b(un· ∇)un (11) ∇ · un+b= 0,(12) un+b=0sobre ∂Ω.(13) siendo a, b > 0 tales que a+b= 2, σ∈(1/2,1] y θ∈[0,1]. Finalmente se coordinan ambas soluciones tomando ˜un+1 = (1/2)(un+a+un+b) y pn+1 =pn+b. 3. Se toma un+1 =˜un+1 en Ω −Pn+1. 4. Se proyecta la velocidad dentro de las part´ıculas sobre la de s´olido r´ıgido. Para ello se invoca la propiedad de conservaci´on del momento lineal y angular, como en [11], tomando, para i= 1, . . . , Np: Un+1 i=ρs,i MiZPn+1 i ˜un+1 dΩ, ωn+1 i=ZPn+1 i rטun+1 dΩZPn+1 i |r|2dΩ, un+1 =Un+1 i+ωn+1 i×ren Pn+1 i. El operador L: el c´alculo de la nueva posici´on de las part´ıculas a partir de su posici´on y velocidad en el paso de tiempo anterior se hace mediante el siguiente esquema iterativo, basado en las leyes del movimiento rectil´ıneo uniformemente acelerado, que parte de Xn+1,0=Xny calcula, para cada k= 1, . . . , K: ∗ X n+1,k =Xn+1,k−1+Un∆t K, Xn+1,k =∗ X n+1,k +1 2F(∗ X n+1,k )+F(Xn+1,k−1)(∆t)2 2M K2, tomando, finalmente, Xn+1 =Xn+1,K =: L(Xn, Un). Aqu´ı, F=Fp+Fwes la suma de las fuerzas de repulsi´on entre part´ıculas y de las paredes. 4 Algoritmo paralelo SDI Tanto la resoluci´on del problema no lineal (9)-(10) como la del problema de Stokes generalizado (11)-(12)-(13) se efect´uan mediante esquemas iterativos en los que cada iteraci´on consiste en la resoluci´on de problemas de Dirichlet para la ecuaci´on de Helmholz como el siguiente: hallar vtal que αv −β∆v=fen Ω, v = ¯vsobre ∂Ω.(14) 5 Jordi Blasco, MaCarmen Calzada, Mercedes Mar´ın Figura 1: Velocidad vertical de la part´ıcula. Figura 2: Ordenada del centro de la part´ıcula. El algoritmo SDI (Simultaneous Directions Implicit, ver [2], [4], [5], [6]) permite la resoluci´on de (14) en paralelo mediante un esquema iterativo que resuelve el problema simult´aneamente en cada direcci´on espacial: partiendo de una aproximaci´on inicial v0, en cada iteraci´on jse calculan dos nuevas aproximaciones vj+1,1yvj+1,2haciendo: (I+τ(α 2I−β∂2 ∂x2 1 ))vj+1,1= (I−τ(α 2I−β∂2 ∂x2 2 ))vj+τf en Ω, vj+1,1=v0sobre ∂Ω,(15) y an´alogamente para vj+1,2, siendo τel par´ametro de evoluci´on. Si se discretiza el dominio Ω mediante un malla cartesiana, la resoluci´on de (15) puede realizarse de manera independiente en cada uno de los llamados segmentos integrables, tanto horizontales como verticales, lo cual proporciona un alto nivel de paralelizaci´on al algoritmo resultante: eventualmente, podr´ıan usarse tantos procesadores en paralelo como segmentos integrables halla en total en el mallado del dominio y resolver el problema de manera independiente en cada segmento. La resoluci´on de (15) se lleva a cabo mediante diferencias finitas centradas, lo cual conduce a un sistema de ecuaciones lineales tridiagonal en cada segmento. Finalmente, se coordinan ambas soluciones tomando vj+1 = (γ/2) (vj+1,1+vj+1,2) + (1−γ)vj, donde γ∈(0,1) es el par´ametro de coordinaci´on. 5 Resultados num´ericos Se presentan a continuaci´on algunos resultados num´ericos obtenidos con el algoritmo descrito en un caso test. 5.1 Sedimentaci´on de una part´ıcula circular Resolvimos el problema de la sedimentaci´on de una part´ıcula circular de radio R= 0.125 en una cavidad rectangular de dimensiones 2 ×6, cuyo centro se sit´ua inicialmente en el punto (1,4) con velocidad y velocidad angular iniciales nulas. La viscosidad del fluido es 6 M´etodo paralelo para flujos de part´ıculas Figura 3: Campo de velocidades y posici´on de la part´ıcula en t= 0.3, 0.5, 0.7y0.9. ν= 0.1, su densidad ρ= 1 y la densidad de la part´ıcula es ρs= 1.25. Este problema fue ampliamente estudiado, entre otros art´ıculos, en [12]. En la resoluci´on de este problema usamos un mallado uniforme de tama˜no h= 1/100 y un paso de tiempo ∆t= 0.001, calculando la soluci´on hasta T= 1. La ejecuci´on fue secuencial. En las Figuras 1 y 2 se representa la evoluci´on temporal de la velocidad vertical de la part´ıcula y de la ordenada de su centro de masas. La part´ıcula sedimenta en t= 0.7810. La velocidad terminal obtenida en este caso (maximo temporal de kuk) fue de 5.522, con un 0.6% de error respecto al valor hallado en [12] que fue de 5.488. En la Figura 3 se muestra el campo de velocidades del fluido y la posici´on de la part´ıcula en diferentes instantes de tiempo; puede observarse que en este caso no se rompe la simetr´ıa del problema. 6 Conclusiones Se ha desarrollado un algoritmo basado en el m´etodo de dominios ficticios para el c´alculo del movimiento de part´ıculas r´ıgidas en un fluido incompresible viscoso. La imposici´on del movimiento de s´olido r´ıgido en cada part´ıcula se lleva a cabo mediante un m´etodo de bajo coste computacional basado en la conservaci´on del momento lineal y angular. Se utiliza el m´etodo SDI para la resoluci´on de las ecuaciones del fluido el cual permite la paralelizaci´on del c´alculo de las nuevas velocidades en cada uno de los diferentes segmentos integrables. Los primeros resultados obtenidos con este algoritmo en un caso test sencillo han sido satisfactorios; se prevee una pr´oxima implementaci´on paralela del mismo que permitir´a el estudio de su escalabilidad. 7 Jordi Blasco, MaCarmen Calzada, Mercedes Mar´ın Agradecimientos Este trabajo ha sido financiado por el Proyecto MCyT BFM2003-06446-C02-02. Agradecemos a Isidoro Albarreal, Jos´e Luis Cruz, Enrique Fern´andez-Cara y Jos´e Rom´an Galo su inestimable ayuda en la realizaci´on de este trabajo. Referencias [1] I.I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fern´andez-Cara, J.R. Galo, M. Mar´ın. Convergence analysis and error estimates for a parallel algorithm for solving the Navier-Stokes equations. Numerische Mathematik, 93 (2002), 201-221. [2] I.I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fern´andez-Cara, J.R. Galo, M. Mar´ın. Time and space parallelization of the Navier-Stokes equations. Journal of Computational and Applied Mathematics, 24 (2005), 115-130. [3] M. Carlson, P. Mucha, G. Turk. Rigid fluid: animating the interplay between rigid bodies and fluid. Conference Proceedings of ACM SIGGRAPH 2004. [4] J.R. Galo, I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fern´andez-Cara, M. Mar´ın. Simultaneous directions parallel methods for elliptic and parabolic systems. C.R. Acad. Sci. Paris, Ser. I, 339 (2004), 145-150. [5] J.R. Galo, I. Albarreal, M.C. Calzada, J.L. Cruz, E. Fern´andez-Cara, M. Mar´ın. A Simultaneous directions parallel algorithm for the Navier-Stokes equations. C.R. Acad. Sci. Paris, Ser. I, 339 (2004), 235-241. [6] J.R. Galo, M.C. Calzada, J.L. Cruz, M. Mar´ın, I. Albarreal, E. Fern´andez-Cara. The smoothing effect of a simultaneous directions parallel method as applied to Poisson problems. Numerical Methods for Partial Differential Equations, 22 (2006), 414–434. [7] R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph. A distributed Lagrange multiplier/fictitious domain method for particulate flows. International Journal of Multiphase Flow, 25 (1999), 755-794. [8] R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph, J. Periaux. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow. Journal of Computational Physics, 169 (2001), 363-426. [9] N.A. Patankar, P. Singh, D.D. Joseph, R. Glowinski, T.-W. Pan. A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flow. International Journal of Multiphase Flow, 26 (2000), 1509-1524. [10] N. Sharma, Y. Chen, N.A. Patankar. A distributed Lagrange multiplier based computational method for the simulation of particulate-Stokes flow. Computer Methods in Applied Mechanics and Engineering, 194 (2005), 4716-4730. [11] N. Sharma, N.A. Patankar. A fast computation technique for the direct numerical simulation of rigid particulate flows. Journal of Computational Physics, 205 (2005), 439-457. [12] D. Wan, S. Turek. An efficient multigrid-FEM method for the simulation of solid-liquid two phase flows. Journal of Computational and Applied Mathematics, 203 (2007), 561-580. 8