SIMULACION EN MATLAB DE UN PROBLEMA

...

(2.8)[pic 22]

Viscosidad adimensional

(2.9)[pic 23]

Como ya he comentado he elegido un valor de tiempo característico que hace que se me adimensionalice la viscosidad.

ECUACIÓN DE TRANSPORTE ADIMENSIONALIZADA

(2.10)[pic 24]

Donde la variación de la función de corriente adimensional con respecto a las coordenadas espaciales adimensionales pueden ponerse de la forma:

(2.11)[pic 25]

(2.12)[pic 26]

Para ver estas relaciones anteriores debemos recordar las expresiones 1.5 y 1.6

ECUACIÓN DE POISSON ADIMENSIONALIZADA

(2.13)[pic 27]

- MÉTODOS NUMÉRICOS

Para la resolución de las ecuaciones adimensionales 2.10 y 2.13 voy a implementar dos métodos numéricos distintos para cada una. Ello es debido a que son dos tipos de ecuaciones diferentes una es hiperbólica () y otra es elíptica (). Además las derivadas parciales que hay en una con respecto al tiempo, no aparecen en la ecuación de Poisson. Por lo que los métodos elegidos para cada ecuación son los siguientes.

MÉTODO DE DUFORT-FRANKEL PARA LA ECUACIÓN DE TRANSPORTE

En este método la derivadas de primer orden, tanto temporales como espaciales, se discretizan de forma centrada.

(3.1)[pic 28]

(3.2)[pic 29]

(3.3)[pic 30]

(3.4)[pic 31]

(3.5)[pic 32]

Las derivadas de segundo orden espaciales se realizan también centradas, pero para el nodo central se escoge la media entre los instantes anterior y posterior. Con ello se elimina la inestabilidad del término difusivo.

(3.6)[pic 33]

(3.7)[pic 34]

MÉTODO DE LIEBMAN PARA LA ECUACIÓN DE POISSON

Para el caso de que los nodos espaciales sean de distinto tamaño, se puede aplicar en resumen la ecuación siguiente para calcular la función de corriente en el instante n+1.

(3.8)[pic 35]

Donde es:[pic 36]

(3.9)[pic 37]

- CONDICIONES DE CONTORNO Y CONDICIONES INICIALES

Dependiendo del tipo de superficie se aplicaran unas determinadas condiciones (ver ilustración 1).

PARED SUPERIOR

Se trata de una pared sólida móvil con velocidad en la dirección x, U.

RESTO DE PAREDES

Se trata de un conjunto de paredes sólidas fijas.

Para ambos tipos de condiciones, al ser paredes sólidas, se les pueden aplicar las mismas ecuaciones discretizadas. En lo que se diferenciarían sería en el valor de la velocidad de la pared.

En las paredes la función de corriente es constante. Por simplicidad para todas las paredes elegimos que el valor de la función de corriente sea nulo.

(4.1)[pic 38]

Mediante el método de Woods, que aproxima la vorticidad en la pared mediante la expansión en serie de Taylor de la función de corriente para la pared. Esta aproximación de segundo orden se obtiene sustituyendo la derivada de la vorticidad frente a la coordenada requerida y sustituyendo la derivada tercera obtenida de la función de corriente.

(4.2)[pic 39]

Donde el subíndice p indica el nodo de la pared y la distancia sobre la normal a la pared del nodo p al nodo p1. El signo – es para las paredes situadas en la parte superior de un dominio y el signo + para las inferiores (dentro de una coordenada espacial específica).[pic 40][pic 41]

CONDICIONES INICIALES

Inicialmente no hay ningún movimiento en el fluido, este está en reposo. Por lo que tanto la vorticidad como la función de corriente serán nulas en este instante.

- IMPLEMENTACIÓN EN MATLAB

A continuación se expondrá el código con los comentarios oportunos para su explicación:

%Para ejecutar la función le indico a esta el número de divisiones en %x y en y, además de el paso de tiempo y el número de iteraciones de %tiempo

function ejercicio3(Nx,Ny,Dt,N)

%PARÁMETRoS

U=1;

H=1;

L=1;

%visc=1;

Re=100;

Dy=1/Ny;

Dx=1/Nx;

beta=Dx/Dy;

evalx=floor(Nx/1.05);

evaly=floor(Ny/1.05);

convx=floor(Nx/1.5);

convy=floor(Ny/1.5);

%creación de arrays

x=0:Dx:L;

y=0:Dy:H;

tauf=N*Dt;

tau=0:Dt:tauf;

%creación de matrices de vorticidad y función de corriente adimensionales y

%condiciones iniciales

vort=zeros(length(x),length(y),length(tau));

fc=zeros(length(x),length(y),length(tau));

for i=3:length(tau) %bucle en tiempo

%inicializo la condición del while en cada tiempo

cond=true;