A matemática como ferramenta para compreender os fenômenos físicos: o caso da Camada Limite Laminar

     Certa vez resolvi uma lista de exercícios, na qual o professor pediu para deduzir e resolver numericamente a denominada Solução de Blasius da camada limite laminar sobre uma placa plana. O resultado foi o gráfico aparece à esquerda da figura, cuja abscissa representa a velocidade do escoamento e a ordenada representa a distância a partir da superfície da placa em contato com o escoamento. Até aí nada de mais, porém algum tempo atrás eu encontrei uma fotografia da camada limite laminar (infelizmente eu não tenho certeza do nome do livro onde achei a foto). Ela aparece no centro da figura. A figura da direita mostra a sobreposição da solução do modelo matemático com o resultado experimental.


    O código-fonte escrito em FORTRAN90 e usado para resolver a Equação de Blasius está listado a seguir. Para a sua implementação não foi empregado o processo de marcha denominado Shooting Method, sugerido em Bejan (2013, p.83), mas sim foi feita a discretização do modelo matemático com as condições de contorno em ambas as extremidades do domínio (veriável transformada eta) e uso do solver TDMA. O programa requer que esteja instalada a versão 3.5 do Wgnuplot.

program Equacao_de_Blasius

 

use portlib

implicit none

 

integer:: i,j

integer:: n

real*8:: etta

real*8:: d_etta

real*8:: a

real*8,allocatable,dimension(:):: f

real*8,allocatable,dimension(:):: fl

real*8,allocatable,dimension(:):: ap,aw,aww,awww,bp

integer :: jmax

real*8 :: tol

 

open(1,file='saida.txt')

1 format(i4, 5(1x, 1pe10.3))

tol = 1d-8      ! Default: 1.0d-7 ou 1.0d-8

jmax = 200

d_etta = 5.0d-4 ! Default: 1.0d-3 ou 5.0d-4

etta = 1.0d1

n = 1 + etta / d_etta

 

write(1,*)'etta:', etta

write(1,*)'d_etta:', d_etta

write(1,*)'n:', n

 

allocate(f(0:n-1),ap(0:n-1),aw(0:n-1),aww(0:n-1),awww(0:n-1),bp(0:n-1),fl(0:n-2))

f = 0.0d0

 

i = 0

ap(i) = 1.0d0

aw(i) = 0.0d0

aww(i) = 0.0d0

awww(i) = 0.0d0

bp(i) = 0.0d0

write(1,1) i, ap(i), aw(i), aww(i), awww(i), f(i)

 

i = 1

ap(i) = 1.0d0

aw(i) = 0.0d0

aww(i) = 0.0d0

awww(i) = 0.0d0

bp(i) = 0.0d0

write(1,1) i, ap(i), aw(i), aww(i), awww(i), f(i)

 

! Chute inicial para f(2):

f(2) = d_etta**2.0d0

 

do j = 1, jmax

 

    do i = 3, n-1

 

          ap(i)  =  2.0d0/(d_etta**3.0d0) + f(i)/(d_etta**2.0d0)

          aw(i)  = -6.0d0/(d_etta**3.0d0) - 2.0d0*f(i)/(d_etta**2.0d0)

          aww(i) =  6.0d0/(d_etta**3.0d0) + f(i)/(d_etta**2.0d0)

          awww(i)= -2.0d0/(d_etta**3.0d0)

          bp(i)  =  0.0d0

 

          f(i) = ( bp(i) - aw(i)*f(i-1) - aww(i)*f(i-2) - awww(i)*f(i-3)) / ap(i)

          if(j==jmax)write(1,1) i, ap(i), aw(i), aww(i), awww(i), f(i)

 

    end do

    !end do

    if(j==jmax) write(1,*)'=================='

    ! correção de f(2):

    a = (f(n-1) - f(n-2))/d_etta

!    a = (f(n-1) - 2*f(n-2) + f(n-3) )/d_etta**2

   

!    write(*,*) j, a !;pause

!    write(*,*) a;pause

!    f(2) = d_etta**2.0d0 * a**(-1.5d0)

    if(a>1.0d0) then

        f(2) = f(2) - tol*abs(1-a)

    elseif(a<1.d0) then

        f(2) = f(2) + tol*abs(1-a)

    elseif(a==1.d0) then

        exit

    endif

    if(j==jmax) write(1,*)'a:',a

    if(j==jmax) write(1,*)'f(2):',f(2)

 

    do i = 3, n-1

 

          ap(i)  =  2.0d0/(d_etta**3.0d0) + f(i)/(d_etta**2.0d0)

          aw(i)  = -6.0d0/(d_etta**3.0d0) - 2.0d0*f(i)/(d_etta**2.0d0)

          aww(i) =  6.0d0/(d_etta**3.0d0) + f(i)/(d_etta**2.0d0)

          awww(i)= -2.0d0/(d_etta**3.0d0)

          bp(i)  =  0.0d0

 

          f(i) = ( bp(i) - aw(i)*f(i-1) - aww(i)*f(i-2) - awww(i)*f(i-3)) / ap(i)

          if(j==jmax) write(1,1) i, ap(i), aw(i), aww(i), awww(i), f(i)

 

    end do

end do

 

close(1)


write(*,*) f(2)/d_etta**2

 

open(2,file='saida2.txt')

2 format(2(1x, 1pe25.10))

! Cálculo da derivada de primeria ordem usando DDS:

do i = 0, n-2

      fl(i) = (f(i+1)-f(i))/d_etta

      write(2,2) fl(i),i*d_etta

end do

close(2)

 

open(3,file='grafico_fl_x_eta.txt')

    write(3,*) "set grid"

    write(3,*) "set xlabel 'Derivada de f [u/U_inf]'"

    write(3,*) "set ylabel 'etta = (y/x)Re**(1/2)'"

    write(3,*) "set xrange[ ", -0.01, ": ", 1.01, " ]"

    write(3,*) "set yrange[ ", -.01, ": ", etta, " ]"

    write(3,*) "set title 'Perfil de velocidades - solução similaridade'"

      write(3,*) "plot 'saida2.txt' using 1:2 with lines"

close(3)

i=system('wgnuplot grafico_fl_x_eta.txt')

!i=system('saida.txt')

end program Equacao_de_Blasius

Comentários

Postar um comentário

Postagens mais visitadas deste blog

Engenharia Reversa do "Boundary Layer Module" - Parte 1

Imagem
Introdução O BLM ( Boundary Layer Module ) é um módulo de cálculo da camada limite que se desenvolve sobre as paredes internas da câmara de empuxo de motores-foguete, entretanto pode ser aplicado em diversas aplicações envolvendo aerodinâmica e transferência de calor. Em sua concepção original o seu objetivo é predizer a espessura da camada limite em cada posição discreta ao longo do eixo longitudinal do motor-foguete. Dentro de um processo iterativo, essa espessura é adicionada à geometria do motor, permitindo que a região não afetada pela camada limite ( i.e. núcleo potencial do escoamento) mantenha a geometria prevista para ótimo desempenho ( e.g. o núcleo potencial terá a geometria calculada com o Método das Características, otimizando o empuxo do motor). O BLM é um dos módulos que compõe um programa maior denominado TDK ( Two-Dimensional Kinetic ), desenvolvido sob contrato da NASA na década de 80. Além do BLM há um módulo para cálculo das propriedades químicas dos produtos de c...
Leia mais

Engenharia Reversa do "Boundary Layer Module" - Parte 2

Imagem
Nesta segunda postagem da série sobre a engenharia reversa do programa Boundary Layer Module (BLM) é detalhada a transformação de coordenadas da equação da QML (balanço da quantidade de movimento linear) do domínio físico para um domínio transformado. O objetivo dessa transformação é não resolver o problema no domínio físico, em geral de geometria complexa, mas sim em um domínio transformado, mais simples e de dimensões unitárias. A estratégia consiste em encontrar a equação no sistema transformado que é análoga àquela do domínio físico. Parte da estratégia é selecionar uma função corrente que permite resolver a equação da QML simultaneamente com a equação da conservação da massa. A primeira parte da postagem pode ser lida neste link , por isso aqui a equação da QML será apresentada sem maiores detalhes e concentrar-se-á no processo de transformação de coordenadas. A nomenclatura das variáveis empregadas aqui é aquela tradicionalmente utilizada na literatura e pode ser encontrada no l...
Leia mais

Ordenação Lexicográfica: Aumento de Complexidade versus Redução de Recurso Computacional

Imagem
 Durante a implementação de um código para modelar a transferência de calor por radiação em um domínio bidimensional, me deparei com o seguinte problema: mesmo levando em conta as simplificações, a grandeza fundamental do problema, a intensidade da radiação, precisaria ser descrita por quatro dimensões, ou seja, é uma variável quadridimensional. As dimensões são: duas coordenadas espaciais: r e z (o problema é em coordenadas cilíndricas) e duas coordenadas que definem a direção de propagação da radiação: o ângulo polar theta e o ângulo azimutal phi . Não é incomum que o modelo numérico de problemas de CFD e CHT produzam matrizes esparsas, ou seja, matrizes onde apenas a diagonal principal e as diagonais adjacentes são formadas por elementos não nulos. Neste caso, o uso de matrizes para computar a relação entre um ente discreto (ponto, volume ou elemento) e os vizinhos faz com que seja necessário alocar muita memória que não será usada. Além de limitar o tamanho do problema que o...
Leia mais