Análise e resolución numérica de modelos de fluxo de persoas

Traballo Fin de Grao Análise e resolución numérica de modelos de uxo de persoas Alfonso Salgado Pais 2020/2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Análise e resolución numérica de modelos de uxo de persoas Alfonso Salgado Pais Xullo 2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Matemática Aplicada Título: Análise e resolución de modelos de uxo de persoas Breve descrición do contido A análise do uxo peonil vense estudando desde os anos 60 e nos últimos anos o movemento de multitudes presenta moitas aplicacións desde o deseño de espazos ata as evacuacións, e na actualidade a denición da distancia social. Un dos obxectivos do TFG é a análise de algúns modelos continuos de uxo de persoas nos que se realizan hipóteses nas que os peóns son axentes activos e intelixentes e teñen un comportamento anisótropo, é dicir, non están igualmente afectados por estímulos procedentes de todas as direccións do espazo. Entre os modelos continuos analizados están os que se corresponden a leis de conservación nos que a variable conservada é o uxo de persoas, proporase a resolución numérica de un modelo coa metodoloxía de volumes nitos desenvolvendo un código especíco a partir dun inicial facilitado polas titoras en MATLAB. iii iv Agradecementos Gustaríame agradecer aquí a grande axuda recibida por parte de Elena Vázquez Cendón, a titora deste proxecto, por todo o tempo adicado a resolver as miñas dúbidas, pola súa proximidade e tranquilidade e pola paixón que sente polas matemáticas. Igualmente, debo darlle as grazas tamén a Saray Busto Ulloa, pola imprescindible axuda cos códigos en Matlab, por axudarme a entender ben o método de volumes nitos e por permitirme volver a Italia. v Índice xeral Resumen ix Introducción xi 1. Un modelo do uxo de persoas 1 1.1. O modelo e as súas variables . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2. Conservación do número de peóns . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3. A velocidade coma función da densidade . . . . . . . . . . . . . . . . . . . . 4 1.4. Aecuaciónimaxe................................. 5 1.5. Omodelodesegundaorde............................ 7 2. Nocións básicas e leis de conservación 9 2.1. Denición e exemplos de leis de conservación . . . . . . . . . . . . . . . . . 9 2.2. O problema de Cauchy e o concepto de solución clásica . . . . . . . . . . . . 15 2.3. Soluciónsfebles.................................. 21 3. O método de volumes nitos 25 3.1. Descricióndométodo .............................. 25 3.2. Converxencia ................................... 29 3.3. Esquemasnuméricos............................... 32 3.4. Programación en Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4. Resultados numéricos 39 A. Código en Matlab 55 Bibliografía 65 vii 2 CAPÍTULO 1. UN MODELO DO FLUXO DE PERSOAS cisións acorde co seu obxectivo, chegar a un punto determinado, e son capaces de avaliar o seu entorno para decidir o camiño óptimo ata dito obxectivo. Comportamento anisotrópico: Os peóns non están afectados por igual polos estímulos externos procedentes de direccións distintas, e vense maiormente condicionados por aqueles que se producen dentro do seu campo visual ou que ocorren na dirección e sentido no que se moven. Comportamento non local: Os peóns reaccionan dun xeito non local ás condicións externas tanto en espazo, posto que reciben estímulos da rexión sensorial, coma en tempo, xa que precisan dun tempo de percepción para captar o estímulo e dun tempo de reacción para decidir como afrontar o mesmo. Estas tres hipóteses, que semellan dunha clara intuición lóxica, refírense a situacións normais nas que os peóns adoitan realizar as súas actividades cunha certa planicación e dun xeito natural. Por outra banda, en situacións de pánico coma unha evacuación, un sismo ou un incendio aparece un comportamento competitivo guiado polo pánico: É habitual que aparezan os nervios e as persoas traten de moverse o máis rápido posible. A interacción entre os peóns é de natureza física, facendo desaparecer a distancia mínima habitual que normalmente manteñen os individuos entre eles ou coas paredes e os obstáculos. As presións elevadas que se rexistran nos atascos de persoas (superiores nalgúns casos aos 4500N/m2 ) poden ser sucientes para derribar unha parede. Diminúe o uxo nas saídas do espazo debido ás aglomeracións ou en ocasións ás persoas feridas ou que se atopan no chan. Xurde o efecto rabaño, que fai que as persoas deixen de pensar por si mesmas e se limiten a seguir o comportamento da maioría e das que teñen máis próximas. Normalmente isto provoca que se centre a atención nas saídas principais, deixando de lado as saídas de emerxencia. Aínda non se conta con modelos que reproduzan con exactitude os comportamentos das persoas, ben porque os individuos poden actuar de xeito compulsivo, poden descoñecer a estratexia óptima debido por exemplo a que se distraen co entorno que os rodea, e nalmente porque calquera modelo que busque reproducir o comportamento humano terá 1.1. O MODELO E AS SÚAS VARIABLES 3 un pequeno grado de imperfección. Estas desviacións do comportamento óptimo trátanse coma utuacións da conduta. Porén, de cara a unha maior simplicidade, o noso estudo basearase case na súa totalidade en situacións normais. 1.1. O modelo e as súas variables De cara a estudar o modelo do uxo de persoas debemos ter en conta, ademais do tempo, t (en segundos), tres variables: a velocidade, v (en metros por segundo); o uxo, q (en peóns por segundo), que representa a cantidade de peóns por unidade de tempo; e a densidade, ρ (en peóns por metro cadrado), que indica a cantidade de peóns por unidade de espazo. Ademais, como veremos posteriormente, ditas variables están intimamente relacionadas entre si. Medir individualmente as magnitudes de cada un dos individuos pode resultar complicado e ademais carece de interese práctico, polo que imos adoptar un enfoque macroscópico e ter en conta a xa comentada hipótese do continuo, que claricaremos cun exemplo: Exemplo 1.1. Un transeúnte situado nun punto x decide contar o número de peóns que pasan polo seu carón. Supoñamos que pasan 4 peóns nos instantes t , t+ 2 , t+ 4.5 e t+ 7 , de feito que a cantidade de peóns por segundo que pasaron polo seu carón no intervalo [t, t + 7] é 4 7 . Así, se lle preguntasen polo uxo no instante t+ 3 , o transeúnte contestaría 4 7 , a pesar de que nese instante non pasou ningún peón polo seu carón. A velocidade e a densidade calcularíanse de xeito similar, tendo en conta o acontecido en intervalos de tempo pequenos, mais dun tamaño suciente para non ter conta de cada peón individualmente. Asumindo isto, podemos denir as nosas variables como campos escalares, que adoptarán un valor único en cada punto x e en cada instante t : v(x, t) , q(x, t) e ρ(x, t) . 1.2. Conservación do número de peóns A idea básica do noso modelo, así como dos métodos numéricos que imos empregar, é a conservación do número de peóns, é dicir, que se partimos dun dominio pechado, as persoas no interior do mesmo serán unha constante. Polo tanto, podemos asumir que o cambio no número de persoas por unidade de espazo dependerá unicamente do uxo entrante e saínte na fronteira do noso dominio. Para explicar esta idea dun xeito máis matemático, basearémonos no traballo sobre o uxo de tráco [7], posto que as ideas neste punto son practicamente as mesmas se falamos de coches coma se o facemos de peóns. 4 CAPÍTULO 1. UN MODELO DO FLUXO DE PERSOAS Así, desenvolveremos o razoamento matemático nunha dimensión espacial e sen ter en conta inicialmente a regularidade das funcións escalares q e ρ , impoñendo logo certas restricións sobre elas para chegar a leis de conservación máis precisas. Matematicamente, o feito de que o número de peóns se conserva pode expresarse mediante a chamada lei de conservación en forma integral: d dt Zb a ρ(x, t)dx =q(a, t)−q(b, t), (1.1) isto é, a variación da cantidade total de peóns coincide coa variación do uxo entre a entrada, q(a, t) , e a saída, q(b, t) . Posto que non puxemos restrición sobre a función uxo q , supoñamos agora que esta é derivable respecto da variable espacial x . Tense así q(a, t)−q(b, t) = −Zb a ∂q ∂x(x, t)dx. (1.2) Se ademais a función densidade ρ é derivable respecto do tempo t, pódese introducir o operador diferencial na integral, obtendo d dt Zb a ρ(x, t)dx =Zb a ∂ρ ∂t (x, t)dx. (1.3) Substituíndo (1.2) e (1.3) en (1.1) chegamos a Zb a ∂ρ ∂t (x, t)dx =−Zb a ∂q ∂x(x, t)dx ⇒Zb a∂ρ ∂t (x, t) + ∂q ∂x(x, t)dx = 0. (1.4) Agora ben, posto que [a, b] é un intervalo arbitrario, verifícase ∂ρ ∂t (x, t) + ∂q ∂x(x, t) = 0. (1.5) Esta última expresión é o que se coñece coma lei de conservación en forma diferencial, válida para todo x∈[a, b] , e como o seu nome nos indica, para chegar a ela é preciso supoñer que as funcións q e ρ son diferenciables. No capítulo 2 veremos que estas condicións non sempre se cumpren en xeral, pero para a dinámica de peóns si que podemos supoñer que q e ρ son funcións de clase un a pedazos, e polo tanto poderemos establecer o modelo de comportamento xeral do uxo de peóns grazas á ecuación (1.5) unida á unha lei que nos especique o que acontece no conxunto de puntos onde se perde a diferenciabilidade. 1.3. A velocidade coma función da densidade Resulta bastante intuitivo entender a relación entre a velocidade dos peóns e a densidade do lugar onde se moven. Supoñamos que nos estamos movendo por unha praza baleira, 1.4. A ECUACIÓN IMAXE 5 de xeito que podemos adaptar a nosa velocidade atendendo exclusivamente a preferencias individuais, e a medida que aumenta a densidade de peóns na praza teremos que ir reducindo a nosa velocidade por mor de evitar choques, ata chegar a unha situación de densidade máxima na cal non nos poderemos apenas mover, algo así coma a densidade de saturación da praza, e que na vida real se alcanza en eventos moi multitudinarios coma por exemplo o 25 de xullo ás 00:00 na praza do Obradoiro en Santiago de Compostela. Así, a intuición física dinos que o módulo da velocidade é función da densidade, V(ρ) , que neste traballo deniremos baseándonos en [13]: V: [0, ρmax]−→ R+ ρ7−→ V(ρ) = vmaxe−α(ρ/ρmax)2, (1.6) onde ρmax é a densidade de saturación comentada anteriormente, vmax é a velocidade libre que un peón pode acadar, supoñendo a situación inicial de ausencia doutros peóns que imposibiliten o seu percorrido, e α é unha constante positiva. Deste xeito, tense que o módulo da velocidade é decrecente con respecto á densidade, e tomará valores entre o seu mínimo vmaxe−α , cando se alcance a densidade de saturación, e o seu máximo vmax na situación ideal na que a densidade sexa nula. 1.4. A ecuación imaxe Con vistas a seguir expoñendo a forma que vai adoptar o noso modelo, cómpre explicar a idea da ecuación imaxe, a cal, baseándonos no feito en [13], vainos permitir denir o campo de vectores que representa a dirección desexada. Partimos dun dominio conexo bidimensional Ω⊂R2 que representa a área obxecto do noso estudo, a cal conta ademais cunha ou varias saídas polas que poden entrar e saír peóns, aumentando ou diminuíndo o uxo q . A fronteira do dominio, ∂Ω=Γ0∪Γp , consiste na unión de ditas saídas polas que entran e saen os peóns, Γ0 , así como do conxunto dos demais límites do espazo en cuestión, Γp , é dicir, os lugares onde remata Ω polos que non se produce intercambio de uxo entre Ω e o exterior. Nese dominio bidimensional, o modelo macroscópico derívase da ecuación de conservación de masa ∂tρ+ div(q)=0, t > 0,(x, y)∈Ω, (1.7) onde ρ=ρ(x, y, t) representa a densidade e q=q(x, y, t)∈R2 o uxo. Se a esta ecuación de conservación de masa lle engadimos unha relación pechada entre velocidade, uxo e densidade, xa teríamos o noso modelo de primeira orde, o cal soamente habería que complementar cunha ecuación de momento adicional (que supón ter en conta a aceleración) 6 CAPÍTULO 1. UN MODELO DO FLUXO DE PERSOAS para chegar ao modelo de segunda orde. Ademais, a partir de agora escribiremos x para referirnos a (x, y) . No que ás relacións se rere, o campo de velocidades ven dado por v=V(ρ)µ , sendo V(ρ) o denido en (1.6) e µ o campo de vectores da dirección desexada polo peón, froito do proceso de avaliación do entorno en busca da ruta óptima para alcanzar o seu obxectivo, Γ0, e que caracteriza o comportamento intelixente que comentamos na primeira sección. Para denir µ dun xeito adecuado, seguimos a caracterización adoptada por Hughes en [5] deníndoo coma a dirección unitaria oposta ao gradiente dun potencial escalar φ µ=−∇φ ||∇φ|| (1.8) Pola súa parte, a función φ é a solución da chamada eikonal equation ou ecuación imaxe, e que ven dada por (||∇φ|| =c(t, x ) para t > 0, x ∈Ω φ= 0 para x ∈Γ0. (1.9) A función c≡c(t, x )∈[1,+∞) é o coeciente de custo, posto que fai referencia ao traballo que debe facer o peón nun determinado instante de tempo t para pasar por un punto concreto x ∈Ω . A súa escolla caracteriza a solución da ecuación imaxe φ , por exemplo se tomamos c= 1 a solución de (1.9) será a distancia entre o punto x e o obxectivo Γ0 , é dicir, φ( x ) = dist( x ,Γ0) . Esta elección precisa facer unha pequena restrición sobre Ω , o cal debe ser un dominio convexo, xa que do contrario, o camiño máis curto unindo x e Γ0 pode non estar completamente contido en Ω , de feito que a aplicación distancia non estaría ben denida. Outras posibilidades á hora de escoller a función custo c buscan expresar dun xeito máis claro a intención dos peóns de evitar as zonas máis transitadas, tomando a función custo c≡c(ρ(t, x )) coma unha función crecente con respecto da densidade. Este é o caso do modelo de Hughes, no cal se ten que c≡c(t, x )=1/V (ρ(t, x )) . A maior crítica que se lle pode facer á esa elección para c é que esta supón que os peóns son capaces de analizar de antemán a rexión por onde camiñan, obtendo un coñecemento absoluto de todo o dominio Ω en calquera instante de tempo, o cal resulta imposible, especialmente se consideramos un dominio o sucientemente extenso. Coa intención de superar este suposto, en [15] introdúcese o concepto de efecto de memoria, segundo o cal os peóns seguen o camiño máis curto ata o seu destino grazas á información que a súa memoria recopilou de anteriores ocasións no mesmo lugar ou en zonas e situacións con densidades similares. Isto leva á elección dun campo de direccións 1.5. O MODELO DE SEGUNDA ORDE 7 µ algo máis complexo: µ=−∇φ−ω∇[V(ρ)−1+g(ρ)] || − ∇φ−ω∇[V(ρ)−1+g(ρ)]||, (1.10) onde φ resolve a ecuación (1.9), ω > 0 é o coeciente de repulsión e g(ρ) unha función que busca modelar do xeito máis adecuado posible a incomodidade do peón en situacións de densidade elevada. Normalmente na expresión de g adoita aparecer un coeciente β positivo que dependerá da persoa en cuestión e representa a sensibilidade individual fronte a estes aspectos. Sen embargo, nos problemas test que resolveremos numericamente imos supoñer que a dirección unitaria µ é un dato do problema. 1.5. O modelo de segunda orde Unha vez denido o campo de velocidades, é preciso establecer tamén a relación entre o uxo, a velocidade e densidade. A elección que faremos neste traballo será o produto de ambas variables, de xeito que q(ρ, x) = ρv , e tendo en conta (1.6) e (1.8), así coma a elección de c(t, x) , chegamos ao seguinte: q(ρ, x) = ρv =ρV (ρ)µ=−ρV (ρ)∇φ ||∇φ|| =−ρV 2(ρ)∇φ=−ρv2 maxe−2α(ρ ρmax )2∇φ. (1.11) Ademais, tendo en conta esta relación, podemos reescribir a ecuación de conservación de masa (1.7) como ∂tρ+div(ρv)=0, t > 0, x ∈Ω, (1.12) e a resolución da mesma é completamente factible se ven acompañada de condicións iniciais e condicións de contorno. Os modelos de segunda orde, empregados por moitos autores para problemas similares a este, ou para a resolución de problemas con características comúns, como no problema do uxo de tráco ou o das aguas someras, foron readaptados a modelos de uxo de persoas por Jiang et al [6]. Así, a ecuación de conservación de masa (1.12) compleméntase cunha ecuación de momento, na cal aparece un produto tensorial, ademais dun termo difusivo no segundo membro que representa unha forza externa actuando sobre os peóns, e que supón incluír os efectos do campo vectorial da aceleración. No noso caso, dita ecuación de momento será a seguinte: ∂t(ρv) + div(ρv ⊗v) + ∇P(ρ) = 1 τρ(V(ρ)µ−v), (1.13) con τ un parámetro de relaxación que permite modelar a forma en que os peóns adaptan a súa velocidade v á velocidade desexada V(ρ)µ . O termo do segundo membro representa 8 CAPÍTULO 1. UN MODELO DO FLUXO DE PERSOAS unha forza repulsiva de xeito similar a como acontece coa lei de enerxía na dinámica de gases, e denotarémolo por A debido ás relacións obvias deste termo coa aceleración dos peóns. Pola súa parte, P(ρ) = p0ργ é unha función dependente da densidade onde p0 e γ son coecientes de repulsión maiores que 0 e 1, respectivamente. Cómpre ademais observar que se trata dunha ecuación vectorial, con dúas compoñentes, posto que o produto tensorial supón unha matriz 2×2 e o último sumando do segundo membro representa o gradiente da composición da función P coa función ρ , e esta última depende de x e de y . Capítulo 2 Nocións básicas e leis de conservación Nesta sección introduciremos as nocións principais acerca dos sistemas de ecuacións hiperbólicos no cal se enmarca o modelo de uxo de persoas. Conceptos importantes serán o de solución clásica, solución feble e problema de Cauchy, que servirán de base teórica para os métodos numéricos que imos empregar. Igualmente, desenvolveremos o concepto de lei de conservación introducido no primeiro capítulo, o cal resulta de vital importancia nos métodos numéricos, e as súas distintas versións dependendo das restricións baixo as que teñamos as funcións. Case todos os conceptos aquí explicados foron extraídos do libro de Vázquez Cendón, [14]. 2.1. Denición e exemplos de leis de conservación Temos xa adquirida a intuición física de que nunha dimensión unha magnitude se conserva cando non se crea nin se destrúe, é dicir, cando os cambios se deben unicamente ás entradas e saídas de dita magnitude nos extremos do intervalo a considerar. Se queremos xeralizar isto a n variables, a intuición será parecida, posto que os cambios dependerán do uxo entrante e saínte na fronteira do dominio. Mostramos así a denición formal de sistema de leis de conservación en varias variables dada por E. Godlewski y P. A. Raviart en [2]. Denición 2.1. Sexa Ω un subconxunto de Rp e sexan fj , con 1≤j≤d , d funcións regulares de Ω en Rp : fj: Ω −→ Rp,1≤j≤d. 9 10 CAPÍTULO 2. NOCIÓNS BÁSICAS E LEIS DE CONSERVACIÓN A forma xeral xeral dun sistema de leis de conservación é ∂w ∂t + d X j=1 ∂ ∂xj (fj(w)) = 0, x = (x1, ... , xd)∈Rd, t > 0, (2.1) onde w= (w1, ... , wp) , chamado vector das variables conservativas, é a función vectorial w:Rd×[0,∞]−→ Ω (x1, ... , xd, t)7−→ (w1, ... , wp). (2.2) O conxunto Ω denomínase conxunto dos estados, mentres que as funcións fj= (f1j, ... , fpj) reciben o nome de funcións uxo. Deste xeito, temos que a ecuación de conservación de masa (1.12) entra claramente na denición de lei de conservación, simplemente identicando o operador lineal diverxencia, é dicir, div(f(w)) = d X j=1 ∂ ∂xj fj(w) . O sistema (2.1) representa formalmente a idea de conservación das cantidades w1, ... , wp . Para velo dun xeito máis claro e similar á forma integral introducida en (1.1), consideremos un dominio arbitrario Ω∈Rd e integremos a ecuación (2.1) sobre Ω . Teremos así: d dt ZΩ w(x, t)dx + d X j=1 Z∂Ω fj(w)·~ηjds = 0, onde ∂Ω representa a fronteira de Ω e ~ηj= (η1, ... ηd) é o vector normal unitario a ∂Ω que apunta cara ó exterior do dominio. Así, a variación de RΩw(x, t)dx co tempo coincide co uxo que entra a través da súa fronteira, −Pd j=1 R∂Ωfj(w)·~ηjds , onde o signo negativo ven de considerar coma positivo o sentido que apunta cara o exterior de Ω , en concordancia coa elección de ~η . Facendo as restricións necesarias sobre o vector das variables conservativas w e sobre as funcións uxo fj , así como tendo en conta que se poden eliminar as integrais posto que o dominio Ω foi escollido de forma arbitraria, chegamos á lei de conservación en forma diferencial, ou simplemente lei de conservación, (2.1). Introducimos agora algúns exemplos de leis de conservación: Exemplo 2.2. A ecuación do transporte é un dos exemplos máis sinxelos que podemos atopar (por exemplo, en Toro [12]), posto que Ω se considera un conxunto de R e só temos unha ecuación de uxo, f(w) = λw , con λ∈R unha constante que pode ser interpretada como a velocidade de propagación. Así, p=d= 1 e a ecuación ven dada por ∂w ∂t (x, t) + λ∂w ∂x (x, t) = 0. (2.3) 2.1. DEFINICIÓN E EXEMPLOS DE LEIS DE CONSERVACIÓN 11 Exemplo 2.3. O modelo do uxo de tráco, tamén con d=p= 1 , considera w=ρ a densidade e o uxo coma unha función non linear da variable conservativa, f(ρ) = ρvm(1 −ρ ρm) , sendo vm e ρm a velocidade e densidade máximas. Exemplo 2.4. As ecuacións das augas someras, ou shallow water equations , na literatura inglesa, modelan o comportamento da auga en ríos, canles e zonas pouco profundas en xeral, e constitúen un dos exemplos máis importantes de leis de conservación. Se traballamos nunha dimensión e supoñemos a profundidade constante (este modelo introdúcese en Stoker [11]) podemos caracterizalo como segue: d= 1, p = 2, w = (h, q)T, f(w) = q, q2 h+gh2 2T , onde h representa a altura da auga para cada punto x e cada instante t , g a gravidade, a variable conservativa q=hu é o caudal por unidade de lonxitude e u a velocidade promediada en vertical. A ecuación das primeiras compoñentes dos vectores correspóndese coa lei de conservación de masa: ∂h ∂t +∂q ∂x = 0, mentres que a das segundas compoñentes, coa ecuación de conservación do momento: ∂q ∂t +∂ ∂x q2 h+gh2 2= 0. Obsérvense as similitudes entre estas ecuacións e as ecuacións de conservación de masa e de conservación de momento do noso modelo, (1.12) e (1.13), a pesar de que aquí non aparecen termos fonte. Estendendo este problema a dúas dimensións chegamos ás ecuacións das augas someras bidimensionais: d= 2, p = 3, w = (h, q1, q2)T, f1(w) = q1,q2 1 h+gh2 2,q1q2 hT , f2(w) = q2,q1q2 h,q2 2 h+gh2 2T , mantendo a notación empregada para o caso unidimensional. Exemplo 2.5. A ecuación de Burgers, ∂w ∂t +w∂w ∂x = 0, foi introducida por Burgers coma o modelo diferencial máis simple para un uído, e presenta certas peculiaridades dun sistema non linear nunha variable que nos servirán posteriormente para introducir algunhas ideas relacionadas co concepto de solución e solución feble. 18 CAPÍTULO 2. NOCIÓNS BÁSICAS E LEIS DE CONSERVACIÓN Despexando e substituíndo adecuadamente, chegamos a ∂w ∂t =−F(ξ)w0 0(ξ) 1 + F0(ξ)t, (2.11) ∂w ∂x =w0 0(ξ) 1 + F0(ξ)t. (2.12) É inmediato ver que estas expresións verican a lei de conservación hiperbólica: ∂w ∂t +F(ξ)∂w ∂x = 0. A importancia da ecuación do transporte, así coma do método das características neste exemplo, radica principalmente no feito de que moitos problemas asociados a sistemas de leis de conservación hiperbólicos poden reducirse a este caso, ou a un número nito dos mesmos. Os sistemas de leis de conservación lineais son un exemplo disto. Denición 2.11. Un sistema de leis de conservación hiperbólico dise que é lineal se é da forma ∂w ∂t +A∂w ∂x = 0, (2.13) sendo A unha matriz con coecientes constantes de orde p×p e diagonalizable, de tal xeito que o vector de uxos é f(w) = Aw . Coñecendo a solución da ecuación do transporte, pódese calcular a solución deste sistema lineal sen máis que considerar a matriz diagonal B formada polos autovalores λk , k= 1, ... p, da matriz A e introducindo as variables características u= (u1, ... , up)T obtidas a partir da matriz de cambio de base P formada polos autovectores de A , w=Pu . Tendo en conta esta relación entre w e u , e que B=P−1AP chegamos a que o sistema (2.13) é equivalente a ∂u ∂t +B∂u ∂x = 0. Como ademais B é unha matriz diagonal, obtemos n ecuacións do transporte da forma ∂uk ∂t +λk ∂uk ∂x = 0, k = 1, ... , p das cales coñecemos a solución exacta: uk(x, t) = u0 k(x−λkt), k = 1, ... , p. Desfacendo nalmente o cambio de variable chegamos á solución do sistema lineal en cada instante t e cada punto x , w(x, t) = Pu(x, t) . 2.2. O PROBLEMA DE CAUCHY E O CONCEPTO DE SOLUCIÓN CLÁSICA 19 Co exposto anteriormente na aplicación do método das características sabemos que a solución e constante ao longo das mesmas, pero aínda non sabemos se esa solución está denida en calquera punto do plano xt . Veremos a continuación un exemplo, sacado de [7] que mostrará que a resposta a esa cuestión é negativa. Exemplo 2.12. Consideremos a ecuación de Burgers introducida anteriormente, ∂w ∂t +w∂w ∂x = 0, coa condición inicial dada por w0(x) =        1, se x < 0, 1−x, se 0<x<1, 0, se x > 1, (2.14) O diagrama das características para este problema pódese observar na gura 2.2. -2 -1,5 -1 -0,5 0 0,5 1 1,5 0,4 0,8 1,2 1,6 2 t=1 t Figura 2.2: Diagrama das características para a ecuación de Burgers coa condición inicial (2.14). Á vista da gura xa observamos que as características se cortan entre sí. Isto signica que w toma tantos valores diferentes no mesmo punto (x, t) coma características se cortan nel, e obviamente isto non ten signicado físico, polo que vai supoñer un problema. De feito, obsérvase coma unhas características son máis rápidas que as outras, é dicir, a pendente da solución aumenta a medida que pasa o tempo. Este exemplo mostra o que non sempre é posible denir a solución para calquera tempo t , e pon de manifesto un problema de vital importancia: a rotura da solución pasado un certo período de tempo. Este fenómeno é o que se coñece coma formación de ondas de 20 CAPÍTULO 2. NOCIÓNS BÁSICAS E LEIS DE CONSERVACIÓN choque, e o que realmente sucede non é que a solución tome distintos valores nun mesmo punto, pois isto carecería de sentido físico, senón que os valores elevados do gradiente motivan a aparición dunha discontinuidade que se mantén no tempo e que se pode desprazar ao longo do espazo. Fisicamente, o que impide que a solución se rompa é a inuencia de certos factores que poden desprezarse para valores pequenos do gradiente, polo que adoitan omitirse á hora de modelizar o problema coma un sistema de leis de conservación. No caso da dinámica de gases, estes factores poden ser a viscosidade ou a condución de calor, mentres que se buscamos situacións máis ligadas ao noso modelo de uxo de persoas, como á formación de atascos ou de situacións de grandes densidades de peóns en lugares públicos, estes factores poderían ser, entre outros, a capacidade dos peóns de adiantarse aos acontecementos ou o tempo de reacción. O momento no que a solución deixa de estar denida coñécese co nome de tempo de rotura, e o seu cálculo pode resultar de interese analítico. Unha condición necesaria para que w tome varios valores á vez nun mesmo punto é que a súa derivada respecto de x, (2.12), deixe de estar denida, o cal vai ocorrer cando o denominador de dita expresión se anule: 1 + F0(ξ)t= 0. Deste xeito, se unha característica pasa por (ξ, 0) , a derivada parcial da solución respecto de x deixa de estar denida en t=−1 F0(ξ). Así, o tempo crítico ou de rotura alcanzarase na característica que minimice dita expresión, tB= m´ın ξ∈R −1 F0(ξ). Exemplo 2.13. Para deixar isto algo máis claro, imos amosar o tempo de rotura na ecuación de Burgers, substituíndo a condición inicial (2.14) por unha máis regular, w0(x) = 2 1 + ex . Temos que F0(ξ) = −2eξ (1 + eξ)2. O máximo desta función atópase en ξ= 0 , porén tB= m´ın ξ∈R −1 F0(ξ)=−1 F0(0) = 2, de xeito que será no instante t= 2 no que as rectas características comenzarán a cortarse. Este problema suxire que vai ser necesario atopar outro tipo de solucións para o problema de Cauchy de modo que o tempo de rotura non supoña un impedimento. Vai ser o caso das solucións febles, que pasamos a denir na sección seguinte. 2.3. SOLUCIÓNS FEBLES 21 2.3. Solucións febles Acabamos de ver na sección anterior que as solucións clásicas de leis de conservación poden deixar de estar denidas a partir dun certo tempo. Co obxectivo de enmarcar estas situacións, entre as que se atopan as do modelo de uxo de persoas, na nosa teoría, é preciso introducir unha nova formulación para o que entendemos por solución do problema de Cauchy que poidamos idealizar cunha discontinuidade. Nesta formulación equivalente para a solución non poden aparecer as derivadas parciais da solución, posto que se esta non vai ser continua, tampouco poderá ser derivable. Denición 2.14. Sexa f unha función uxo dada por unha lei de conservación integral, e consideremos o espazo das funcións test, que é o seguinte conxunto: V=C1 0(Rd×[0,∞)) = nϕ:Rd×[0,∞)/ϕ ∈C1 e ten soporte compacto } Dise que unha función w:Rd×[0,∞) medible e localmente limitada é unha solución feble do problema de Cauchy asociado á lei de conservación integral se satisfai a ecuación ZRd w( x ,0) ·ϕ( x ,0)d x +Z∞ 0ZRd w∂ϕ ∂t + d X j=1 fj(w)∂ϕ ∂xj d x dt = 0 ∀ϕ∈V (2.15) así coma a condición inicial (2.8), é dicir, w( x ,0) = w0( x )∀ x ∈Rd, onde w0 é unha función coñecida. É importante remarcar o feito de que toda solución feble vai asociada a unha lei de conservación integral. Ademais, para a ecuación (2.15), non temos unicidade de solución. Introducimos agora esta proposición que formaliza a relación entre solución clásica e solución feble, algo que semella bastante intuitivo a priori, coa súa demostración con d = 1, posto que resulta máis comprensible nunha dimensión e estendela logo a varias dimensións resulta trivial. Proposición 2.15. Toda solución clásica do problema de Cauchy (2.9) é tamén unha solución feble. Demostración. Supoñamos que w é unha solución clásica do problema de Cauchy, ou o que é o mesmo, unha función C1(R×[0,∞)) vericando as condicións (2.1) e (2.8). 22 CAPÍTULO 2. NOCIÓNS BÁSICAS E LEIS DE CONSERVACIÓN Se multiplicamos a lei de conservación diferencial (2.1) por unha función test ϕ∈V e integramos respecto de x e de t en R×[0,∞) , chegamos ó seguinte: Z∞ 0ZR∂w ∂t +∂f(w) ∂x ϕdxdt = 0. Como ϕ é unha función con soporte compacto, existe un δ > 0 de forma que sop(ϕ)⊂ U= (−δ, δ)×[0, δ) . Deste xeito, a función ϕ anúlase fóra de dito conxunto U, polo que podemos reescribir a integral anterior coma ZU∂w ∂t +∂f(w) ∂x ϕdxdt = 0. Podemos empregar a regra da cadea sobre as funcións w , f e ϕ , posto que as tres son funcións de clase C1 , polo que escribimos ∂w ∂t ϕ=∂ ∂t(wϕ)−w∂ϕ ∂t , ∂f(w) ∂x ϕ=∂ ∂x(f(w)ϕ)−f(w)∂ϕ ∂x , obtendo deste xeito ZU∂ ∂t(wϕ) + ∂ ∂x(f(w)ϕ)dxdt −ZUw∂ϕ ∂t +f(w)∂ϕ ∂x dxdt = 0. (2.16) Aplicamos agora o teorema de Green ó primeiro sumando, converténdose este en Z∂U (f(w)ϕ, wϕ)·~ndγ, onde ~n denota ó vector unitario normal a ∂U , que é a fronteira de U . Nótese agora que U∩∂U = (−δ, δ)×{0} , e coma sop(ϕ)⊂U , concluímos que ϕ vale cero sobre ∂U −(−δ, δ)× 0, polo que podemos simplicar aínda máis a expresión anterior: Z(−δ,δ)×0 (f(w)ϕ, wϕ)·~ndγ. O valor do vector normal unitario nese anaco de fronteira é ~n = (0,−1) , e tendo iso en conta chegamos nalmente a ZU∂ ∂t(wϕ) + ∂ ∂x(f(w)ϕ)dxdt =−Zδ −δ w(x, 0)ϕ(x, 0)dx. Substituíndo nalmente esa expresión en (2.16) chegamos a −Zδ −δ w(x, 0)ϕ(x, 0)dx −ZUw∂ϕ ∂t +f(w)∂ϕ ∂x dxdt = 0, 2.3. SOLUCIÓNS FEBLES 23 e ampliamos agora estas integrais a todo espazo, razoando de xeito inverso a coma xemos antes, xa que ϕ vale cero fóra de U : −Z∞ −∞ w(x, 0)ϕ(x, 0)dx −ZR×[0,∞)w∂ϕ ∂t +f(w)∂ϕ ∂x dxdt = 0. Esta última expresión proba que w é unha solución feble do problema de Cauchy. Denición 2.16. Un caso particular do problema de Cauchy é o problema de Riemann, aquel no que a condición inicial ven dada por unha función w0 , continua en todo R agás nun punto arbitrario, x0 : w0(x) = (wl, se x≤x0 wr, se x>x0, (2.17) onde wl , e wr son constantes arbitrarias que reciben o nome de estados a esquerda e a dereita de x0 , respectivamente. Esta condición inicial discontinua vai ser de gran importancia para o modelo de uxo de persoas, posto que será así coma supoñamos a densidade inicial, e no noso código estaremos resolvendo problemas de Riemann. En particular, supoñeremos a densidade inicial do seguinte modo: ρ0=(ρl, se x≤x0 ρr, se x>x0, (2.18) sendo x0 o valor intermedio do intervalo do dominio no eixo das X. Se temos que ρl> ρr o estado da esquerda propágase a unha velocidade menor que o estado da dereita e estaremos ante o que se coñece coma unha onda de rarefacción. No caso contrario, é dicir se ρl< ρr , será o estado da esquerda o que se propague a unha maior velocidade, e diremos que se forma unha onda de choque. Mostramos isto gracamente nas guras 2.3 e 2.4. Ademais, no noso caso as demais compoñentes do vector de variables conservativas serán inicialmente 0 , supoñendo que os peóns inician o seu movemento dende o repouso. Polo tanto teremos que wl= (ρl,0,0)T e wr= (ρr,0,0)T , de xeito que o valor inicial de w será un ou outro dependendo de se x é maior ou menor que x0 . 24 CAPÍTULO 2. NOCIÓNS BÁSICAS E LEIS DE CONSERVACIÓN -2,5 -2 -1,5 -1 -0,5 0 0,5 1 1,5 2 2,5 0,5 1 1,5 2 2,5 3 fan t Figura 2.3: Diagrama das características para unha onda de enrarecemento. -2 -1,5 -1 -0,5 0 0,5 1 1,5 2 2,5 3 x=st Figura 2.4: Diagrama das características para unha onda de choque. Capítulo 3 O método de volumes nitos Nesta sección imos explicar o funcionamento dos métodos numéricos que nos permitan aproximar a solución das ecuación do modelo hiperbólico de uxo de persoas, así coma veri- car a súa converxencia baixo certas condicións e repasar brevemente a súa implementación en Matlab. 3.1. Descrición do método A idea principal do método de volumes nitos consiste en discretizar o dominio computacional Ω en pequenas celas, Ci e aproximar a integral da variable conservativa sobre cada unha delas no instante t=tn , w(·, tn) mediante a constante wn i : wn i≃1 A(Ci)ZCi w(x, tn)dA, (3.1) para despois empregar ese valor aproximado no seguinte paso de tempo e calcular con eles o uxo na fronteira das celas, que nos permitirá nalmente actualizar o valor da integral de ω . Explicaremos primeiro o método nunha dimensión espacial por simplicidade, estendendo a continuación a explicación a dúas dimensións empregando un mallado estruturado. Asemade teremos en conta que debemos denir dous uxos numéricos en lugar de un. No caso unidimensional, o dominio Ω consiste simplemente nun intervalo [a, b] con a, b ∈R , e as celas quedan determinadas por unha malla de puntos {xi}M i=1 . Cada cela Ci defínese como Ci= (xi−1 2, xi+1 2), con xi−1 2=xi−xi−1 2 e xi+1 2=xi+1−xi 2 . Neste traballo, a pesar de que o métodos tamén é aplicable para mallas arbitrarias, limitarémonos a mallas uniformes, é dicir, particións do intervalo [a, b] nas que os puntos 25 26 CAPÍTULO 3. O MÉTODO DE VOLUMES FINITOS {xi}M i=1 son equidistantes e a distancia que separa dous puntos consecutivos denotámola por ∆x . Porén, a expresión das celas pode simplicarse como segue: Ci=xi−∆x 2, xi+∆x 2. Partimos da lei de conservación (2.1) con d= 1 , ∂w ∂t +∂f(w) ∂x = 0, e comezamos integrando esta expresión sobre cada cela Ci : ZCi ∂w ∂t dt =−ZCi ∂f(w) ∂x dx, o cal é equivalente a d dt Zxi+1 2 xi−1 2 w(x, t)dx =f(w(xi−1 2, t)) −f(w(xi+1 2, t)). (3.2) Introducimos agora o tempo, e consideramos dous instantes tn e tn+1 distintos e tales que tn< tn+1 , integrando a ecuación anterior entre eles. Deixamos no primeiro membro só o primeiro sumando, xa que estamos en busca dun método explícito: Zxi+1 2 xi−1 2 w(x, tn+1)dx =Zxi+1 2 xi−1 2 w(x, tn)dx −Ztn+1 tn f(w(xi+1 2, t))dt −Ztn+1 tn f(w(xi−1 2, t))dt. (3.3) Se recordamos a expresión (3.1) aparece no denominador o termo A(Ci) que neste caso será ∆x , polo que imos dividir por ∆x a expresión (3.3): 1 ∆xZxi+1 2 xi−1 2 w(x, tn+1)dx =1 ∆xZxi+1 2 xi−1 2 w(x, tn)dx −1 ∆xZtn+1 tn f(w(xi+1 2, t))dt −Ztn+1 tn f(w(xi−1 2, t))dt. (3.4) Agora necesitamos unha aproximación dos valores medios do uxo nos puntos xi±1 2 durante o intervalo de tempo [tn, tn+1] : fn i±1 2 ≃1 ∆tZtn+1 tn f(w(xi±1 2, t))dt. (3.5) 3.1. DESCRICIÓN DO MÉTODO 27 Se temos en conta esta expresión, así coma a aproximación de w(·, tn) dada por (3.1), a ecuación (3.3) transfórmase en wn+1 i=wn i−∆t ∆xfi+1 2−fi−1 2. Nos problemas hiperbólicos que imos resolver, dispoñemos sempre dunha condición inicial, polo que podemos asegurar que no momento de calcular wn+1 i coñeceremos o valor de wn i , xa que basta con retroceder no algoritmo ata chegar a w0 i . Por outra banda, pode non parecer claro a priori a forma de obter as aproximacións fn i±1 2 . É aquí onde entra a idea principal sobre a que se sustenta o método de volumes nitos: a información propágase a unha velocidade nita se traballamos con problemas hiperbólicos (por ser os autovalores reais), polo que podemos obter a aproximación do uxo na fronteira entre dúas celas baseándonos en aproximacións dos valores promedio de w en cada unha delas. Se temos isto en conta, as aproximacións fn i±1 2 faranse mediante a chamada función uxo numérico: fi+1 2= Φ(wn i, wn i+1), fi−1 2= Φ(wn i−1, wn i). A elección desta función uxo numérico será o que caracterice cada un dos diferentes esquemas do método de volumes nitos. En resumo, mediante esta notación chegamos á que se coñece como forma conservativa do método: wn+1 i=wn i−∆t ∆xΦ(wn i, wn i+1)−Φ(wn i−1, wn i), i = 2,··· , M −1. (3.6) Teremos ademais M−2 condicións iniciais, posto que os valores de w0 i son coñecidos para i= 2,··· , M −1 . Sen embargo, imos ter que implementar algún tipo de condición de contorno para os extremos da malla, C1 e CM e describimos a continuación as alternativas máis habituais para facer isto: Supoñendo que a condición inicial é constante no tempo nas proximidades dos extremos, pódese tomar wn+1 1=wn 1 e wn+1 M=wn M . Isto pode non ser adecuado nalgúns casos xa que o método deixará de ser consistente se chega algunha perturbación ós extremos. Outra posibilidade habitual consiste en dar condicións de contorno periódicas. Para isto, suponse que o problema é cíclico, de tal xeito que as fronteiras x1−1 2 e xM+1 2 pasan a ser a mesma e se calcula o uxo nesa fronteira mediante os valores de wn M e wn 1 . 34 CAPÍTULO 3. O MÉTODO DE VOLUMES FINITOS mentres que a derivada espacial de f aproxímase de forma centrada: ∂f(w) ∂x ≃f(wn i+1)−f(wn i−1) 2∆x. Aínda que a primeira vista pode non apreciarse que se trata dun método conservativo, facendo algunhas manipulacións alxébricas chegamos á seguinte función uxo numérico: Φ(u, v) = f(u) + f(v) 2−∆x 2∆t(v−u). Trátase dun método consistente, pois verica a condición (3.10), e condicionalmente estable para o problema lineal. Q-esquemas Para poder caracterizar de xeito preciso estes esquemas, cómpre primeiro denir algúns conceptos, introducidos por Harten, Lax e van Leer en [3]. Denición 3.2. Sexa x∈R , chamamos parte positiva e parte negativa de x ás seguintes cantidades: x+=x+|x| 2, x−=x− |x| 2. Denición 3.3. Un esquema conservativo dise descentrado se ∀u, v estados próximos, a función uxo numérico Φ verica as seguintes propiedades: a) Existe outro estado, u∗ , próximo tamén a u e a v de xeito que: Φ(u, v) = f(u∗) + f0 +(u∗)(u−u∗) + f0 −(u∗)(v−u∗) + o(|u−u∗|+|v−u∗|), onde f0 + e f0 − son as partes positiva e negativa da derivada da función uxo. b) Se f0(u) e f0(v) son ambas positivas, entón Φ(u, v) = f(u) . c) Se f0(u) e f0(v) son ambas negativas, entón Φ(u, v) = f(v) . Unha vez explicados estes conceptos, podemos englobar os Q-esquemas coma un subconxunto dos esquemas descentrados nos que o uxo numérico ven determinado por unha matriz Q : Φ(u, v) = f(u) + f(v) 2−1 2|Q(u, v)|(v−u). (3.13) Dita matriz Q , que depende continuamente de u e v , será o que diferencie os distintos Q-esquemas. A matriz |Q| obtense a partir de Q tendo os mesmos autovectores e os seus 3.3. ESQUEMAS NUMÉRICOS 35 autovalores son os valores absolutos dos autovalores de Q . Algúns exemplos de Q-esquemas son o de van Leer, que considera Q coma a matriz xacobiana do uxo avaliada no punto medio de ambos estados, Q(u, v) = ∂f ∂w u+v 2, ou o Q-esquema de Roe, no cal Q ven denida como segue: f(v)−f(u) = Q(u, v)(v−u). Se analizamos a función uxo (3.13) vemos que ten dous sumandos. O primeiro sumando, f(u) + f(v) 2 é a parte centrada do uxo, mentres que o termo −1 2|Q(u, v)|(v−u) representa a parte descentrada e recibe o nome de viscosidade numérica. Un problema que pode suceder cos Q-esquemas é que se algún dos autovalores da matriz Q é nulo, entón o termo da viscosidade numérica anúlase e pérdese a converxencia do método. Así, para asegurar a existencia de viscosidade numérica, redefínese o valor absoluto dos autovalores na parte correspondente do uxo numérico. Explicamos a continuación algún exemplo desta situación. Fluxo de Rusanov O modelo de Rusanov coincide con aplicar a regularización local de Lax-Friedrichs aos Q-esquemas no caso escalar, de xeito que a corrección do valor absoluto dos autovalores da matriz xacobiana do uxo en cada fronteira das celas se leve a cabo tomando o máximo dos valores dos mesmos a cada lado da cela. No caso de sistemas de leis de conservación, o uxo de Rusanov emprega dito máximo no termo correspondente á viscosidade numérica en lugar da matriz |Q| . Deste xeito, a función uxo numérico para o método ven dada por: Φ(u, v) = f(u) + f(v) 2−1 2|λ(u, v)|(v−u), (3.14) sendo |λ(u, v)|= m´ax k=1,···n{|λk(u)|,|λk(v)|}, para unha matriz Q cadrada de orde n . Regularización parabólica de Harten A regularización parabólica de Harten tamén ten como obxectivo garantir a existencia dunha determinada viscosidade numérica. Para iso, redene o valor absoluto da matriz Q do modo seguinte: 36 CAPÍTULO 3. O MÉTODO DE VOLUMES FINITOS |λ|ε=|λ|+1 2(1 + sig(ε− |λ|)) λ2+ε2 2ε− |λ|, (3.15) A idea é introducir un parámetro moi pequeno, ε de xeito que se |λ|< ε , entón |λ| substitúese por λ2+ε2 2ε e o valor de |λ|ε é ε 2 . Un inconvinte desta regularización é que o termo ε ten que ser axustado para cada problema. Existen outros uxos numéricos que se empregan no método de volumes nitos, coma por exemplo as técnicas de descomposición do uxo, que non desenvolveremos aquí. Sen embargo, os métodos aquí desenvoltos serán sucientes para o noso obxectivo e para unha primeira aproximación á idea xeral do método de volumes nitos. Pasamos polo tanto a programar en Matlab un modelo xeral do método de volumes nitos, empregando o uxo numérico de Rusanov. 3.4. Programación en Matlab Para desenvolver o modelo de uxos de persoas e atopar unha solución do mesmo faremos uso do método de Rusanov. O código de volumes nitos desenvolto permitiría implementar outros uxos numéricos dos mencionados na memoria ou dos existentes na literatura. Ademais, o código do método de Rusanov en dúas dimensións atópase no Apéndice A. Argumentos de entrada Para comezar debemos introducir os datos que caracterizan o noso problema. Denotamos por argumentos de entrada a todos estes valores, como os extremos dos intervalo de denición, xL e xR no eixo horizontal e yL e yR no eixo vertical, así coma o tempo no que comenzamos, o tempo nal ó que queremos chegar e os coecientes propios do noso modelo, γ , vmax e p0 , denidos no capítulo 2. Ademais, tamén indicamos as condicións iniciais das variables conservativas. Nos diferentes test considerados resolvemos problemas de Riemann, nos que caracterizamos a densidade a ambos lados dunha curva na que se dene a discontinuidade e asumimos que o dato de entrada do uxo é nulo. Denición das mallas A continuación especicamos o número de nodos que teremos no eixo horizontal, IMAX, e no eixo vertical, JMAX, obtendo os valores de ∆x e ∆y , denindo deste xeito a malla estruturada regular. Tamén introducimos o número de Courant para calcular no código o valor de ∆t a partir dos valores de ∆x e ∆y e o máximo dos autovalores en todos os nodos. 3.4. PROGRAMACIÓN EN MATLAB 37 Inicialización das variables conservativas nas matrices Q e Qn A notación empregada no código para as variables conservativas é Q . Non entra en condición coa detallada na memoria por empregar o uxo numérico de Rusanov que non precisa a matriz |Q| . Inicializamos as variables en dous tempos asignándolle o valor cero e polo tanto introducindo dúas matrices de ceros de dimensión 3×IMAX ×JMAX : Q , que representa a solución no instante tn , e Qn , que representa a solución no instante tn+1 , posto que debido á metodoloxía do método de volumes nitos serán necesarias as dúas. Ademais, tamén introduciremos a densidade inicial na primeira columna da matriz Q , dependendo do test que desexemos realizar. As outras dúas columnas serán 0 , posto que, como xa foi mencionado, supoñeremos que o movemento comeza dende o repouso, de xeito que a velocidade inicial vai ser nula, e polo tanto as dúas compoñentes do uxo tamén. Paso temporal Feito isto, iniciaremos o bucle temporal e calcularemos ∆t en cada iteración tendo en conta o número de Courant, ∆x , ∆y e os autovalores das matriz Q en cada instante de tempo, que calcularemos mediante unha función Matlab. Deste xeito, teremos que dt =CFL amaxx dx +amaxy dy , onde CFL denota o número de Courant e amaxx e amaxy os valores absolutos dos autovalores máximos da matriz Q na dirección horizontal e vertical, respectivamente. Ademais, asegurarémonos de que se alcanza o tempo nal cunha condición if . Método de volumes nitos Dentro do bucle temporal iniciado no paso anterior, iniciaremos bucles espaciais no eixo X e no eixo Y para o cálculo dos uxos en ambas direccións, a partir dunha función Matlab. Calcularemos ademais o termo fonte en cada iteración, tamén mediante unha función externa, e nalmente actualizamos a solución envorcando Qn en Q e o paso temporal. Cálculo de erros No caso de que se coñeza a solución exacta ou que exista unha solución de referencia pódense calcular os erros absoluto e relativo cometidos en cada paso de tempo, imprimindo por pantalla os seus máximos en norma L2 . Representación gráca Finalmente, represéntase a solución aproximada obtida en tres grácas diferentes, unha para cada compoñente da matriz de variables conservativas. Isto consistira nunha animación en dúas dimensións, para poder visualizar como evoluciona a solución en calquera instante de tempo e en calquera punto (x, y) do dominio de denición. 38 CAPÍTULO 3. O MÉTODO DE VOLUMES FINITOS Capítulo 4 Resultados numéricos Neste último capítulo amosaremos os resultados obtidos nalgunhas simulacións feitas en Matlab, empleando o uxo numérico de Rusanov. Para iso, faremos distintos tests, onde algúns coecientes permanecerán sempre xos, e outros variarán para intentar simular situacións diferentes. Baseándonos nas simulacións feitas en [13] tomaremos os seguintes valores dos parámetros, que serán comúns a todos os tests que fagamos: p0= 0.005, γ = 2, vmax = 2m/s, τ = 0.61, α = 7.5. Ademais, supoñeremos inicialmente que a densidade máxima é ρmax = 10 persoas por metro cadrado, e nun test en concreto variaremos este valor para simular unhas condicións nas que os peóns deben manter unha distancia de seguridade entre elas. Traballaremos cun dominio de 50 ×50 , en metros, polo que será un cadrado, e en xeral consideraremos 100 nodos en ambas direccións, aínda que isto tamén variará dependendo do test. Tamén pararemos as simulacións en t= 100 s. Por outra banda, o número de Courant que tomaremos será 0.7 , tendo en conta que para que haxa estabilidade este debe ser menor que 1 . Non escollemos un número máis preto de 1 , coma 0.9 ou 0.8 porque a condición (3.11) ven dada en situacións xerais, pero neste caso a existencia do termo fonte pode afectar á condición de estabilidade. Para as condicións iniciais, supoñeremos que a velocidade inicial é nula, como se non houbese movemento no instante t0 , e daremos as densidades, ρl e ρr , en dúas seccións distintas e complementarias do dominio. Deste xeito, pretendemos simular unha situación na que os peóns se atopan nunha sección do dominio, e que non se lles permite o paso á outra ata o instante t0 , coma se houbese unha valla ou cinta de separación. Finalmente, as condicións de contorno serán de tipo Neumann, supoñendo por exemplo que o uxo a esquerda nas celas C1,j é o uxo na propia fronteira, aplicando a consistencia 39 40 CAPÍTULO 4. RESULTADOS NUMÉRICOS dos uxos numéricos cos correspondentes uxos físicos. Faremos isto nos catro bordos de Ω . Test 1: Movemento horizontal Neste primeiro caso, dividiremos o dominio pola sección x= 25 m, e a dirección de movemento dos peóns será completamente horizontal: µ= (1,0) . Ademais, partimos dunha densidade inicial dada por ρl= 4 , ρr= 0 , polo que estamos a supoñer que non hai peóns na parte dereita do dominio, e que os peóns da esquerda se desprazan cara alí. Mostramos gracamente a densidade inicial na gura 4.1. Figura 4.1: Densidade inicial nos tests horizontais supoñendo ρl= 4 e ρr= 0 . Os resultados obtidos, que vemos na gura 4.2 conrman con certa claridade a intuición que tiñamos antes de facela simulación: a densidade vaise distribuíndo cara a parte dereita do dominio, coma se os peóns que se atopaban na parte esquerda se movesen cara alá. Ademais, os peóns poden saír pola parte dereita unha vez chegan ó bordo dereito de Ω , polo que pasado un tempo a densidade total dentro do dominio se reduce. Co uxo horizontal, q1 , vese claramente coma o entorno de maior movemento é a sección x= 25 , é dicir, o lugar onde estaría posta inicialmente a cinta que prohibe o paso. Ademais, é maior conforme nos movemos cara x= 50 que cara x= 0 . O uxo vertical q2 é constantemente igual a 0 , xa que o movemento é só horizontal. A velocidade horizontal é maior canto máis nos acerquemos a x= 50 , especialmente para t preto de 0 , xa que os peóns que pasaron á parte dereita non teñen restricións á hora de moverse, mentres que os da parte esquerda deben esperar a que os peóns que teñen diante deixen espazo. 41 Fixando o valor de y=y0 e un instante de tempo determinado, podemos ver mellor como varían as solucións dependendo do punto x : Figura 4.2: Densidade ( ρ ), en azul; uxo horizontal ( q1 ), en vermello; e velocidade horizontal ( v1 ), en amarelo; nos instantes t= 25 s, t= 50 s e t= 75 s con y=y0 xo. Pódese extraer bastante información se analizamos as tres grácas da gura 4.2 en profundidade, pero a modo de resumo diremos que a densidade se move claramente cara a dereita conforme pasa o tempo, mentres que o uxo estabilízase bastante en 3 peóns por segundo e a velocidade parece ser sempre crecente con respecto a x , aínda que isto deixa de suceder a partir de t= 75 . Con todo, cando nos acercamos ao tempo t= 100 aparecen pequenas utuacións tanto no uxo coma na velocidade, tal vez debidas á escala empregada nos eixos, coma se pode apreciar nas grácas da gura 4.3: 42 CAPÍTULO 4. RESULTADOS NUMÉRICOS Figura 4.3: Densidade, uxo horizontal e velocidade horizontal no instante t= 100 . Para mostrar agora como poden variar os resultados cunha malla máis ou menos re- nada, repetiremos o este test dúas veces en idénticas condicións, salvo por unha pequena variación: na primeira imos considerar unha malla con 50 nodos en cada dirección, mentres que na segunda consideraremos 200 nodos en cada dirección. 43 Figura 4.4: Densidade ( ρ ), en azul; uxo horizontal ( q1 ), en vermello; e velocidade horizontal ( v1 ), en amarelo; nos instantes t= 25 s, t= 50 s e t= 75 s con y=y0 xo e 50 nodos en cada eixo. 50 CAPÍTULO 4. RESULTADOS NUMÉRICOS Figura 4.11: Densidade, uxo vertical e velocidade horizontal para o test diagonal no instante t= 25 . Mostramos a continuación na gura 4.12 a velocidade vertical, para dar unha evidencia gráca de que é a mesma que a velocidade horizontal. Figura 4.12: Velocidade vertical para o test diagonal no instante t= 25 . 51 Variando parámetros coma o vector de dirección ou as condicións iniciais poden facerse outros tests que simulan situacións diferentes. Unha opción que contemplamos pero que nalmente non incluímos aquí é un test no que a condición inicial veña dada de xeito radial, coma se a discontinuidade se atopase na circunferencia de radio 10. 52 CAPÍTULO 4. RESULTADOS NUMÉRICOS Conclusión Antes de nalizar co traballo, cómpre resumir de xeito breve o feito nel: Comezamos expoñendo o modelo de uxo de persoas unidimensional e bidimensional, explicando como variar algún dos parámetros podía permitir adaptarse mellor ó comportamento humano, tan variable e difícil de prever. Seguimos explicando a teoría matemática na que se atopaba dito modelo, e desenvolvemos tamén a metodoloxía que nos ía permitir facer as simulacións pertinentes. Finalmente, puidemos atopar certas solucións para o problema en situacións sinxelas nas que a dirección de movemento é un dato do problema. Recordemos ademais que a aplicación deste modelo debe limitarse a situacións normais, posto que a modelización do comportamento humano en situacións de pánico é dunha grande complexidade. Tamén é preciso destacar que o presente traballo serviu para afondar no estudo das ecuacións en derivadas parciais hiperbólicas, así coma para desenvolver o método de volumes nitos, metodoloxía con grandes aplicacións en diversos ámbitos. Finalmente, aínda que aquí non o chegamos a desenvolver, a idea da ecuación imaxe e a súa resolución pode ser de grande interese e contar con moitas aplicacións, polo que animamos a quen o desexe que continúe este Traballo de Fin de Grado desenvolvendo ese concepto e resolvendo dita ecuación en busca da dirección da velocidade en cada paso temporal. 53 54 CAPÍTULO 4. RESULTADOS NUMÉRICOS Apéndice A Código en Matlab Mostramos aquí o código Matlab do problema de uxo de persoas empregando o uxo de Rusanov. % Resolutor para o modelo de fluxo de persoas en 2D clear all close all clc global p0 global gamma global vmax global eps % Definimos o problema físico p0 = 0.005; % presión gamma = 2; % expoñente adiabático vmax = 2; % velocidade desexada eps = 10^-6; tempo = 0; tf = 100; rho_l = 4; rho_r = 0; % densidades a esquerda e dereita % Definimos o dominio computacional, que vai ser regular 55 56 APÉNDICE A. CÓDIGO EN MATLAB xL = 0; xR = 50; % bordos a esquerda e dereita yL = 0; yR = 50; % bordos arriba e abaixo IMAX = 100; JMAX = 100; % número de celas en ambas direccións dx = (xR-xL)/IMAX; dy = (yR-yL)/JMAX; % incrementos espaciais x = linspace(xL+dx/2,xR+dx/2,IMAX); y = linspace(yL+dy/2,yR+dy/2,JMAX); CFL = 0.7; % Condición de Courant % Damos a densidade inicial (rho_l ou rho_r) e a velocidade inicial (0,0) Q = zeros(3,IMAX,JMAX); Qn = zeros(3,IMAX,JMAX); for i = 1:IMAX for j = 1:JMAX %if x(i) < (xR-xL)/2 % test 1 %if y(j) < (yR-yL)/2 % test 2 %if x(i) + y(j) < 50 % test 3 Q(:,i,j) = [rho_l; 0; 0]; else Q(:,i,j) = [rho_r; 0; 0]; end end end % Representacion gráfica da densidade inicial surf(x,y,squeeze(Q(1,:,:))') axis([xL xR yL yR 0 5]) xlabel ('X(m)') ylabel ('Y(m)') zlabel ('densidade(persoas/m^2)') title (sprintf('Densidade inicial')) % Definimos os vectores para almacenar as velocidades v_x = zeros(1,IMAX,JMAX); v_y = zeros(1,IMAX,JMAX); % Método de volumes finitos NMAX = 10000; 57 for n = 1:NMAX % iniciamos o bucle temporal amax_x = 0; amax_y = 0; for i = 1:IMAX for j = 1:JMAX amax_x = max(amax_x,max(abs(lambdax(Q(:,i,j))))); amax_y = max(amax_y,max(abs(lambday(Q(:,i,j))))); end end dt = CFL/(amax_x/dx + amax_y/dy); % Chegada ó tempo final if (tempo + dt > tf) dt = tf - tempo; end % Criterio de parada if (tempo >= tf) break end for i = 1:IMAX % bucle no eixo X for j = 1:JMAX % bucle no eixo Y % fluxo na dirección X if (i==1) % fronteira da esquerda QL = Q(:,i,j); fp = Rusanovx(Q(:,i,j),Q(:,i+1,j)); fm = Rusanovx(QL,Q(:,i,j)); elseif (i==IMAX) % fronteira da dereita QR = Q(:,i,j); fp = Rusanovx(Q(:,i,j),QR); fm = Rusanovx(Q(:,i-1,j),Q(:,i,j)); else % dentro do dominio fp = Rusanovx(Q(:,i,j),Q(:,i+1,j)); fm = Rusanovx(Q(:,i-1,j),Q(:,i,j)); end %fluxo na dirección Y if (j==1) % fronteira de abaixo QL = Q(:,i,j); 58 APÉNDICE A. CÓDIGO EN MATLAB gp = Rusanovy(Q(:,i,j),Q(:,i,j+1)); gm = Rusanovy(QL,Q(:,i,j)); elseif (j==JMAX) % fronteira de arriba QR = Q(:,i,j); gp = Rusanovx(Q(:,i,j),QR); gm = Rusanovx(Q(:,i,j-1),Q(:,i,j)); else % dentro do dominio gp = Rusanovx(Q(:,i,j),Q(:,i,j+1)); gm = Rusanovx(Q(:,i,j-1),Q(:,i,j)); end % Definimos o termo fonte s = S(Q(:,i,j)); % Solución no instante n+1 Qn(:,i,j) = Q(:,i,j)- dt/dx *(fp-fm) - dt/dy *(gp-gm) + dt*s; end end % Actualizamos o tempo e a solución Q = Qn; tempo = tempo + dt; % Representacion gráfica figure(1) surf(x,y,squeeze(Q(1,:,:))') % densidade axis([xL xR yL yR 0 5]) title(sprintf('Densidade no tempo t = %f',tempo)) drawnow figure(2) surf(x,y,squeeze(Q(2,:,:))') % fluxo horizontal axis([xL xR yL yR]) title(sprintf('Fluxo horizontal no tempo = %f',tempo)) drawnow figure(3) 59 surf(x,y,squeeze(Q(3,:,:))') % fluxo vertical axis([xL xR yL yR]) title(sprintf('Fluxo vertical no tempo = %f',tempo)) drawnow % Introducimos unha sentencia condicional para poder representar tamén % as velocidades. if (Q(1,:,:) <= eps) v_x(1,:,:) = 0; v_y(1,:,:) = 0; else v_x(1,:,:) = Q(2,:,:)./Q(1,:,:); v_y(1,:,:) = Q(3,:,:)./Q(1,:,:); end figure (4) surf (x,y,squeeze(v_x)') % velocidade horizontal axis([xL xR yL yR]) title(sprintf('Velocidade horizontal no tempo =%f',tempo)) drawnow figure(5) surf (x,y,squeeze(v_y)') % velocidade vertical axis([xL xR yL yR]) title(sprintf('Velocidade vertical no tempo =%f',tempo)) drawnow end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fluxo numérico na dirección x function fluxo = Rusanovx(Qm,Qp) s_max = max(max(abs(lambdax(Qm)),abs(lambdax(Qp)))); fluxo = 0.5*(f(Qm)+f(Qp))- 0.5*s_max*(Qp-Qm); 66 BIBLIOGRAFÍA [12] E. F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction. Springer, 1997. [13] Monika Twarogowska, Paola Goatin, Régis Duvigneau. Numerical study of macroscopic pedestrian ow models.[Research Report] RR-8340,INRIA. hal-00849587, 2013. [14] M.Elena Vázquez Cendón. Introducción al método de volúmenes nitos . Manuais universitarios(10), 2008. [15] Y. Xia, S. C. Wong, C. Shu, Dynamic continuum pedestrian ow model with memory eect. Physics Revision E 79: 066113, 2009.