Equação de Blasius Camada Limite MATLAB

Prévia do material em texto

%PUC-Rio Pontificia Universidade Católica do Rio de Janeiro
%Departamento de Mecânica
%Mestrado em Engenharia Mecânica
%Mecânica dos Fluidos II
%Aluno: Juliano Leodonio Xavier Lima
%Problema 11
Clc %limpa memória/ limpa tela
%-----------------------------------------------
%Equação de Blasius : 2f"+ff'=0
%-----------------------------------------------
%Condicao de contorno
%1. f'(0) = 0
%2. f'(inf) = 0
%3. f(0) = 0
%Algumas sugestoes para a solução inicial do metodo
%4. f"(0) = 0.332
%-----------------------------------------------
%find: f,f" and f'''
%-----------------------------------------------
%valor da funcao em forma de matriz
%eta = x
%f = y0,f'=y1,f"=y2,f'''= -(1/2)*y0*y2
%-----------------------------------------------
func1 = inline('y1','x','y0','y1','y2');
%-----------------------------------------------
func2 = inline('y2','x','y0','y1','y2');
%-----------------------------------------------
func3 = inline('-0.5*y0*y2','x','y0','y1','y2');
%-----------------------------------------------
%Implantacao da linguagem com input estabelecido usualmente (tirar o coments)
%input: x = 0 , y0 = 0 , y1= 0
% y2 = 0.332 , total = 7 and h = 0.1
%-----------------------------------------------
%Implantacao da linguagem com interacao com o usuario
%=================================================
x = input('\n Informe o valor de etha (sugestao: zero) x : ');
y0 = input( 'Informe o valor inicial de f : ');
y1 = input( 'Informe o valor inicial de f1 : ');
y2 = input( 'informe o valor de f2 (sugestao: 0.332 : ');
total = input('Enter the value of total : ');
h = input( 'Informe o valor de analise (h) : ');
%==================================================
fprintf('\n Tamanho do passo no método =%5.3f é:',h);
fprintf('\n etha f f1 f2');
fprintf('\n%16.3f%16.3f%16.3f%16.3f',x,y0,y1,y2);
%=================================================
%Implantacao da solucao numerica por Runge Kutta
%================================================
 
for i = 1:(total/h)
 %Criacao do metodo
ak1y0 = func1(x,y0,y1,y2);
ak1y1 = func2(x,y0,y1,y2);
ak1y2 = func3(x,y0,y1,y2);
xx = x + h/2.;
yy0 = y0 + h*ak1y0/2.;
yy1 = y1 + h*ak1y1/2.;
yy2 = y2 + h*ak1y2/2.;
ak2y0 = func1(xx,yy0,yy1,yy2);
ak2y1 = func2(xx,yy0,yy1,yy2);
ak2y2 = func3(xx,yy0,yy1,yy2);
yy0 = y0 + h*ak2y0/2.;
yy1 = y1 + h*ak2y1/2.;
yy2 = y2 + h*ak2y2/2.;
ak3y0 = func1(xx,yy0,yy1,yy2);
ak3y1 = func2(xx,yy0,yy1,yy2);
ak3y2 = func3(xx,yy0,yy1,yy2);
xx = x + h;
yy0 = y0 + h*ak3y0;
yy1 = y1 + h*ak3y1;
yy2 = y2 + h*ak3y2;
ak4y0 = func1(xx,yy0,yy1,yy2);
ak4y1 = func2(xx,yy0,yy1,yy2);
ak4y2 = func3(xx,yy0,yy1,yy2);
y0 = y0 + (ak1y0 + 2.*ak2y0 + 2.*ak3y0 + ak4y0)*h/6.;
y1 = y1 + (ak1y1 + 2.*ak2y1 + 2.*ak3y1 + ak4y1)*h/6.;
y2 = y2 + (ak1y2 + 2.*ak2y2 + 2.*ak3y2 + ak4y2)*h/6.;
x = x + h;
%Vetores dinâmicos
v([i])=x;
v1([i])=y1;
fprintf('\n%16.3f%16.3f%16.3f%16.3f', x,y0,y1,y2);
end
%Imprime o grafico de f' em função de etha
plot(v,v1);grid, title('Camada-Limite laminar ao longo de uma placa plana. By: Juliano Lima'), xlabel('eta'); ylabel('df1')
Resultados
 Informe o valor de etha (sugestao: zero) x : 0
Informe o valor inicial de f : 0
Informe o valor inicial de f1 : 0
informe o valor de f2 (sugestao: 0.332) : 0.332
Enter the value of total : 10
Informe o valor de analise (h) : 0.5
 Tamanho do passo no método =0.500 é:
 etha f f1 f2
 0.000 0.000 0.000 0.332
 0.500 0.042 0.166 0.331
 1.000 0.166 0.330 0.323
 1.500 0.370 0.487 0.303
 2.000 0.650 0.630 0.267
 2.500 0.996 0.751 0.217
 3.000 1.397 0.846 0.161
 3.500 1.837 0.913 0.108
 4.000 2.305 0.955 0.064
 4.500 2.790 0.979 0.034
 5.000 3.283 0.991 0.016
 5.500 3.780 0.997 0.007
 6.000 4.279 0.999 0.002
 6.500 4.779 1.000 0.001
 7.000 5.279 1.000 0.000
 7.500 5.778 1.000 0.000
 8.000 6.278 1.000 0.000
 8.500 6.778 1.000 0.000
 9.000 7.278 1.000 0.000
 9.500 7.778 1.000 0.000
 10.000 8.278 1.000 0.000