Application to the tube chock problem

En el libro An Introducctión to Scientific Computing Twelve Computational Projects Solved with MATLAB  presentan doce proyectos, que si bien estan resueltos, son un buen ejercicio eso días en los que no hay deberes. Entre mis intereses están los métodos numéricos para resolver las ecuaciones que gobiernan la dinámica de los fluidos (CFD-Dinámica de fluidos Computacional). Los proyectos 10 y 12 tratan estas ecuaciones. Hay algunas cuestiones importantes en la solución numérica de estas ecuaciones, que tienen que ver con la matemática de las ecuaciones, soluciones discontinuas, condiciones iniciales discontinuas; ver la matemática necesaria para entender estas ecuaciones es una de mis expectativas del curso de EDP 2 que voy a llevar este ciclo. 

No tengo MATLAB por el momento, así que me las arregle, use SAGE y lo programe en Python; estoy obteniendo los resultados.

class Tubechock:
    def __init__(self):
        self.gamma=1.4
        self.cfl=0.95
        self.g1=line([])
        self.g2=line([])
        self.g3=line([])
        self.g4=line([])
        self.M=81
        self.x0=0.5
        self.rhoL=8.0
        self.rhoR=1.0
        self.pL=10.0/self.gamma
        self.pR=1.0/self.gamma
        self.UL=0
        self.UR=0
        self.rhoUL=self.rhoL*self.UL
        self.rhoUR=self.rhoR*self.UR
        self.EL=self.pL/(self.gamma-1)
        self.ER=self.pR/(self.gamma-1)
        self.dx=1/(self.M-1)
        self.dt=0
        self.x=[i*self.dx for i in range(self.M)]
        self.t=0
        self.w=[[0 for j in range(3)]for i in range(self.M)]
        self.F=[[0 for j in range(3)]for i in range(self.M)]
        self.FF=[[0 for j in range(3)]for i in range(self.M-1)]
        for i in range(self.M):
            if (self.x[i]<=self.x0):
                self.w[i]=[self.rhoL,self.rhoUL,self.EL]
            if (self.x[i]>self.x0):
                self.w[i]=[self.rhoR,self.rhoUR,self.ER]
    
    
    def plot_d(self):        
        xvsd=[[self.x[i],self.w[i][0]] for i in range(self.M)]
        self.g1=line(xvsd,rgbcolor=(3/4,1/2,5/8))
        self.g1.axes_labels(['$x$',r'$\rho$'])
        self.g1.fontsize(15)
        
    
    
    def plot_v(self):
        xvsv=[[self.x[i],self.w[i][1]/self.w[i][0]] for i in range(self.M)]        
        self.g2=line(xvsv,rgbcolor=(1/4,1/4,5/8))
        self.g2.axes_labels(['$x$','$v$'])
        self.g2.fontsize(15)
    
    
    def plot_p(self):
        xvsp=[]
        for i in range(self.M):
            k=(self.w[i][2]-0.5*self.w[i][1]*self.w[i][1]/self.w[i][0])*(self.gamma-1)
            xvsp.append([self.x[i],k])
        self.g3=line(xvsp,rgbcolor=(6/8,3/4,3/8))
        self.g3.axes_labels(['$x$','$p$'])
        self.g3.fontsize(15)
 
 
    def plot_T(self):
        xvsT=[]
        for i in range(self.M):
            k=(self.w[i][2]-0.5*self.w[i][1]*self.w[i][1]/self.w[i][0])*(self.gamma-1)
            xvsT.append([self.x[i],k/self.w[i][0]])
        self.g4=line(xvsT,rgbcolor=(3/4,3/5,1/3))
        self.g4.axes_labels(['$x$','$T$'])
        self.g4.fontsize(15)
        
    
    def graf(self):
        self.plot_d()
        self.plot_v()
        self.plot_p()
        self.plot_T()
        g=graphics_array([[self.g1,self.g2],[self.g3,self.g4]])
        g.show(frame=True, axes=True, figsize=[12,8])
        
    
    
    def Dt(self):
        uloc=[(self.w[i][1]/self.w[i][0]) for i in range(self.M)]
        press=[(self.gamma-1)*(self.w[i][2]-0.5*self.w[i][1]*uloc[i]) for i in range(self.M)]
        k=self.gamma*press[i]/self.w[i][0]        
        a=[((k)^(0.5)) for i in range(self.M)]
        max=[(abs(uloc[i])+a[i]) for i in range(self.M)]
        max.sort()
        self.dt=0.95*self.dx/max[self.M-1]
        
        
    
    def Fw(self):
        for i in range(self.M):            
            uloc  = self.w[i][1]/self.w[i][0]
            rhou2 = self.w[i][1]*uloc
            press = (0.4)*(self.w[i][2]-0.5*rhou2)
            self.F[i]=[self.w[i][1],rhou2+press,(self.w[i][2]+press)*uloc]
            
    def FFw(self):
        for i in range(self.M-1):
            pas=[0.5*(self.w[i][j]+self.w[i+1][j])-0.5*(self.dt/self.dx)*(self.F[i+1][j]-self.F[i][j]) for j in range(3)]
            uloc  = pas[1]/pas[0]
            rhou2 = pas[1]*uloc
            press = (0.4)*(pas[2]-0.5*rhou2)
            self.FF[i]=[pas[1],rhou2+press,(pas[2]+press)*uloc]
                
    
    
    def sol(self,tf=0.2):
        while (self.t<tf):                    
            self.Dt()
            self.Fw()
            self.FFw()
            for i in range(1,self.M-1):
                self.w[i][0]=self.w[i][0]-(self.dt/self.dx)*(self.FF[i][0]-self.FF[i-1][0])
                self.w[i][1]=self.w[i][1]-(self.dt/self.dx)*(self.FF[i][1]-self.FF[i-1][1])
                self.w[i][2]=self.w[i][2]-(self.dt/self.dx)*(self.FF[i][2]-self.FF[i-1][2])
            self.t=self.t+self.dt            
        
 
mysol=Tubechock()
mysol.sol()
mysol.graf()

Esquema numérico de Lax-Wendroff y Mac-Cormack

El método de Euler o de Runge Kutta no son apropiados para calcular soluciones discontinuas(y es que las soluciones no son siempre tan bonitas) desde que estos generan oscilaciones sin significado físico o desbalances de cantidades invariantes(masa por ejemplo).

Por ejemplo, para el problema "the shock tube", para la que hay una solución discontinua (velocidad, presión, densidad)

Velocidad

Densidad

tenemos que resolver las ecuaciones de Euler. Entonces usamos los esquemas de Lax-Wendroff y Mac-Cormack.

Contaminantes

Transporte de Contaminantes

Existen principalmente dos tipos de procesos físicos por los cuales los elementos químicos son transportados mediante fluidos en el medio ambiente: advección y difusión.

ADVECCIÓN.-

Este proceso se debe al movimiento del fluido, ya sea aire o agua (convección, término similar, normalmente se asocia al movimiento vertical de advección debido a diferencias de densidad). Luego un elemento químico presente en el aire o en el agua será pasivamente llevado por este movimiento advectivo de masas.

El movimiento advectivo es descrito matemáticamente por la dirección y la magnitud de su velocidad, dado que a pesar de la ocurrencia de dispersión, el centro de masa del elemento químico que es transportado por advección, se mueve a la velocidad promedio del fluido, siempre y cuando no se produzca adsorción y retardo .

DIFUSIÓN.-(transporte difusivo o “Fickiano”)

En este segundo tipo de proceso, el elemento químico se mueve desde un lugar donde su concentración es relativamente alta hacia otro donde es menor, por efecto de un movimiento aleatorio de las moléculas (difusión molecular), a un movimiento aleatorio del aire o agua que acarrea al elemento químico (difusión turbulenta) o por una combinación de ambos.

BALANCE DE MASA.-

Esta ecuación, que describe el principio de conservación de masas en un volumen infinitesimal, establece que la tasa de cambio de almacenamiento del compuesto o elemento químico en cualquier punto del espacio, dC/dt, es igual a la razón de ingreso y salida del químico por medios físicos, más la razón de producción interna (fuentes menos sumideros). Los ingresos y salidas ocurren por medios físicos, advección y difusión (transporte “Fickiano”), y se expresan en términos de la velocidad del fluido (v), el coeficiente de difusión/dispersión (D), y el gradiente de concentración del químico en el fluido (dC/dx). Las entradas y salidas asociadas con fuentes o sumideros internos se denotan por r. Así, en una dimensión, la ecuación para un punto fijo es:

donde
C : concentración del compuesto o elemento químico
v  : velocidad del fluido
D : coeficiente de difusión dispersión

Esta ecuación nos dice que la concentración del elemento o compuesto puede cambiar (dC/dt) si existe una concentración diferente en alguna parte del fluido y esta es transportada al punto de interés(v*dC/dx), o debido al transporte “Fickiano” si es que existe una variación espacial en la concentración del fluido ((d/dx)(D*(dC/dx)). Cambios en la concentración también pueden ocurrir como consecuencia de reacciones biológicas o químicas, que remueven o introducen el compuesto o elemento de interés (r).

Para una situación tridimensional, la ecuación de advección-dispersión-reacción puede ser sucintamente representada usando notación vectorial, donde ∇ es el operador gradiente, de la siguiente forma, asumiendo que D es igual en todas las direcciones:

Fuente:

http://www.aulados.net/Temas_ambientales/Contaminantes_aguas_subterraneas/Transporte_contaminantes.pdfhttp://www.aulados.net/Temas_ambientales/Contaminantes_aguas_subterraneas/Transporte_contaminantes.pdf