Construindo o Diagrama de Moody

    Um tempo atrás implementei uma rotina para resolver a equação de Colebrook-White, que descreve o fator de atrito em escoamentos turbulentos. Sendo uma equação transcendental, não é possível isolar o fator de atrito diretamente. Por isso usei o método da bisseção para resolvê-la numericamente.

    Depois tive uma ideia: construir o Diagrama de Moody que aparece nos livros de Mecânica dos Fluidos. Tomei os valores típicos de rugosidade relativa e/D e resolvi para cada valor, e pelo menos 1000 valores do número de Reynolds. Também inclui a solução para escoamento laminar, obtendo o seguinte diagrama:


A rotina de solução da Equação de Colebrook-White é mostrada a seguir. Ela é escrita em FORTRAN90.

subroutine Fator_de_Atrito_Moody(ro,um,ni,e_D,Re_D,fator_de_atrito)

 implicit none

integer :: tipo_esc ! Tipo de escoamento (1 p/ laminar; 2 p/ transição e 3 p/ turbulento)

real(8),intent(in) :: ro      ! Massa específica do fluido [kg/m3]

real(8),intent(in) :: um      ! Velocidade média do escoamento [m/s]

real(8),intent(in) :: ni      ! Viscosidade cinemática do fluido [m2/s]

real(8),intent(in) :: Re_D  ! Número de Reynolds

real(8),intent(out) :: fator_de_atrito ! Fator de atrito de Moody ou Fator de atrito de Darcy [s/ dimensão]

real(8),intent(in) :: e_D   ! e/D = Rugosidade relativa (ver valores típicos de 'e' na tabela abaixo)

                          !==============================

                         !Processo de   |

                         !de fabricação | e[micrometros]

                         !                    |

                         !Extrudado     | 1.5d0

                         !Aço comercial | 46.0d0

                         !Ferro fundido | 260.0d0

                         !Concreto      | 3.0d2 à 3.0d3

                         !==============================

! Auxiliary variables

integer :: it

integer :: itmax

real(8) :: x

real(8) :: raiz_f0

real(8) :: raiz_f1

real(8) :: raiz_f, f  !Raiz quadrada do fator de atrito de Moody, f**5.d-1

 

! Equação transcendental de Colebrook para escoamento turbulento:

f(x) = 1.0d+0/x + 2.0d+0*dlog10(e_D/3.7d+0 + 2.51d+0/(Re_D*x))

 

! Inicialização das variáveis

itmax = 10000

if(Re_D<=2.0d+3)                 tipo_esc = 1 ! Escoamento laminar

if(Re_D >2.0d+3.and.Re_D<4.0d+3) tipo_esc = 2 ! Escoamento de transição

if(Re_D>=4.0d+3)                 tipo_esc = 3 ! Escoamento turbulento

select case (tipo_esc)

      case (1) ! Escoamento Laminar

            fator_de_atrito = 6.4d+1/Re_D

      case (2) ! Escoamento de Transição (calcula tanto o fator laminar quanto o turbulento, mas mantém o turbulento)

        fator_de_atrito = 6.4d+1/Re_D

        raiz_f0 = sqrt(8.0d-3)

        raiz_f1 = sqrt(1.0d-1)

        do it = 1, itmax

                  raiz_f = 0.5d0 * ( raiz_f0 + raiz_f1 )

                  if ( f(raiz_f) * f(raiz_f0) > 0.0d+0 ) then

                        raiz_f0 = raiz_f

                  else

                        raiz_f1 = raiz_f

                  end if

                  if ( abs(raiz_f1-raiz_f0) < 1.d-15 ) exit

            end do

            if (it > itmax) write(*,*)"O processo de calculo do fator de atrito nao convergiu com", itmax, "iteracoes!"

            fator_de_atrito = raiz_f*raiz_f

            !write(*,*)"Fator de atrito Turbulento:",fator_de_atrito

    case (3) ! Escoamento Turbulento

        raiz_f0 = sqrt(8.0d-3)

        raiz_f1 = sqrt(1.0d-1)

        do it = 1, itmax

                  raiz_f = 0.5d0 * ( raiz_f0 + raiz_f1 )

                  if ( f(raiz_f) * f(raiz_f0) > 0.0d+0 ) then

                        raiz_f0 = raiz_f

                  else

                        raiz_f1 = raiz_f

                  end if

                  if ( abs(raiz_f1-raiz_f0) < 1.d-15 ) exit

            end do

            if (it > itmax) write(*,*)"O processo de calculo do fator de atrito nao convergiu com", itmax, "iteracoes!"

            fator_de_atrito = raiz_f*raiz_f

end select

end subroutine Fator_de_Atrito_Moody

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