Embed Size (px)
Citation preview
FENÔMENOS DE TRANSPORTE E HIDRÁULICA: EXEMPLOS DE
MÉTODOS NUMÉRICOS E COMPUTACIONAIS PARA A
RESOLUÇÃO DE PROBLEMAS
André Luiz Andrade Simões1, Romualdo José Romão Brito*
2,
Harry Edmar Schulz3, Rodrigo de Melo Porto
4
1,3,4Escola de Engenharia de São Carlos, Universidade de São Paulo, Brasil.
2,Instituto Superior Politécnico de Tete, Moçambique
*E-mail: [email protected]
RESUMO
São apresentados diferentes exemplos de resolução de problemas normalmente apresentados
aos estudantes, profissionais e pesquisadores ligados às áreas dos fenômenos de transporte,
mecânica dos fluidos e hidráulica. Os exemplos fornecem códigos computacionais livres, que
permitem o uso imediato em exemplos ou na aplicação prática. Essa postura decorre do fato
de o engenheiro lidar, nessa área, com problemas cujas soluções não são triviais, sendo
necessário recorrer, com certa frequência, ao uso de códigos gratuitos e comerciais que podem
ser utilizados em calculadoras programáveis ou computadores de diferentes tipos. Assim,
objetiva-se colaborar nos aspectos vinculados ao ensino de engenharia associados ao uso de
recursos computacionais e métodos numéricos, que, conduzidos de forma adequada,
certamente trarão bons frutos à formação de profissionais e pesquisadores.
Palavras-chave: Hidráulica didática, métodos numéricos em fenômenos de transporte.
INTRODUÇÃO E OBJETIVOS
O mundo moderno demanda projetos mais seguros e sustentáveis em diferentes áreas das
atividades humanas, que, associados às novas técnicas de produção e construção, e aos novos
materiais e técnicas de produção, implica em maior detalhamento e rapidez na preparação
desses projetos. Sendo assim, a Engenharia também deve ampliar seus métodos, apoiando-se
em técnicas cada vez mais aprimoradas para auxiliar os profissionais. Neste contexto, o
engenheiro necessita efetuar cálculos cada vez mais complexos e diversificados, o que
resultou em desenvolvimentos acelerados em áreas como a modelagem matemática,
fluidodinâmica computacional, simulação de processos etc. Para recorrer aos métodos e
recursos desses ramos é necessário dominar aspectos básicos de métodos numéricos, além de
ter à disposição programadores, softwares livres ou licenças de softwares, o que pode limitar a
difusão de determinadas ferramentas devido aos custos elevados.
O objetivo geral deste trabalho pode ser enunciado como oferecer conhecimentos básicos de
forma estruturada para auxiliar no uso e disseminação de ferramentas computacionais em
fenômenos de transporte, mecânica dos fluidos aplicada e hidráulica. O objetivo específico
pode ser enunciado como apresentar diferentes aplicações de pacotes e códigos
computacionais elaborados para a resolução de problemas nas áreas mencionadas.
Os exemplos propostos e resolvidos neste trabalho tiveram como objetivo destacar alguns
aspectos essenciais vinculados à física-matemática dos mesmos, aos métodos numéricos e ao uso de códigos livres e comerciais. A planilha eletrônica Excel
® é empregada para resolver
problemas de determinação da vazão escoada entre dois reservatórios ligados por uma
tubulação, para o ajuste de curvas e para a obtenção de soluções para equações diferenciais
ordinárias e parciais. Foram escritos códigos equivalentes que funcionam no software
Matlab® ou, com poucas adaptações, no software livre Octave. A seguir são apresentadas
sínteses sobre os problemas de aplicação propostos.
CÁLCULO DA VAZÃO ESCOADA ENTRE DOIS RESERVATÓRIOS
Um problema básico de hidráulica, que consta provavelmente em todos os programas de
Engenharia Civil, consiste em determinar a vazão escoada através de uma tubulação entre
dois reservatórios com nível de água constante, localizados em cotas diferentes. Este
problema é modelado com a primeira lei da termodinâmica escrita convenientemente e com a
equação de Darcy-Weisbach para o cálculo da energia dissipada (embora hajam alternativas
com maior base empírica, como a equação de Hazen-Williams). Uma vez que o fator de
resistência de Darcy-Weisbach é calculado através do diagrama de Moody, traduzido por
equações não-lineares e com partes logarítmicas, a solução explícita para o problema não é
viável. Sendo assim, é necessário utilizar um método numérico para calcular as raízes da
função objetivo definida a partir das equações mencionadas (método de Newton-Raphson,
entre outros). O uso do recurso solver do Excel® é uma alternativa interessante para solução
deste tipo de problema, sendo explorado aqui.
A Figura 1, apresentada a seguir, ilustra de forma esquemática o problema em questão. O
reservatório 1 (da esquerda) alimenta o reservatório 2 por meio de uma tubulação com
diâmetro D, comprimento L e rugosidade absoluta equivalente .
Figura 1 – Dois reservatórios de níveis constantes interligados por uma tubulação
Considera-se que as dimensões dos reservatórios são suficientemente grandes para que os
níveis z1 e z2 das superfícies livres sejam considerados constantes. Para tais condições, a
equação da energia assume a seguinte forma:
N
i
ig
VK
g
V
D
Lfzz
1
22
2122
(1)
em que: z = cota correspondente à energia potencial gravitacional, f = fator de resistência de
Darcy-Weisbach, L = comprimento da tubulação, V = velocidade média do líquido na
tubulação, D = diâmetro interno do conduto, g = aceleração da gravidade, Ki = coeficientes de
perda de carga localizada, com i = 1, 2,...,N, sendo N = número de perdas localizadas.
O problema modelado com a equação 1 exige o conhecimento do fator de resistência de
Darcy-Weisbach. Com o uso de análise dimensional (com o teorema de Vaschy-
Buckingham), sabe-se que este adimensional é uma função do número de Reynolds (Re) e da
rugosidade relativa (/D). A aplicação do referido teorema não resulta em uma solução, mas
em uma função que requer alguma teoria adicional ou a realização de experimentos para que
ela possa ser utilizada. O trabalho de Nikuradse é um exemplo clássico que relaciona f com
Re e /D. Atualmente, o cálculo do fator de resistência é realizado com equações ajustadas ao diagrama de Moody que, assim como o diagrama de Nikuradse, abrange os regimes laminar e
turbulento. As equações para tubos e condutos com seções genéricas podem ser deduzidas a
partir das equações de Navier-Stokes para escoamentos laminares, sendo a solução de Hagen-
Poiseuille um exemplo interessante e bastante conhecido (f = 64/Re, válida para tubos de base
circular). Outros casos com seções transversais genéricas podem ser vistos em Simões et al.
(2010), inclusive para o regime variável. Considerando o escoamento turbulento, há equações
que podem ser deduzidas a partir de hipóteses e informações experimentais. Entre elas, pode-
se citar a equação de von Kármán, de 1930 (Schulz, 2003), e a equação de Blasius que, a
despeito da sua simplicidade, ajusta-se satisfatoriamente a resultados experimentais em tubos
lisos, como P.V.C. para 3000<Re<105 (Porto, 2006). A possibilidade de generalizar o ajuste
aos dados do diagrama de Moody é de interesse prático e levou ao desenvolvimento de
algumas equações, como as de Churchill (1977), Chue (1984), Pereira e Almeida (1986) e
Swamee (1993). Neste trabalho é utilizada a equação de Swamee (1993), escrita a seguir:
125.016
6
9.0
8
Re
2500
Re
74.5
7.3ln5.9
Re
64
Df
(2)
As informações relacionadas às perdas de carga localizadas podem ser obtidas em
tabelas ou gráficos desenvolvidos experimentalmente e, em poucos casos, por meio de
formulações teóricas (um exemplo é a equação de Borda-Carnott). Deste modo, escolhidos os
valores de Ki, para valores de z1, z2, L, D, , (viscosidade cinemática) conhecidos, tem-se
um sistema fechado de equações que possibilita calcular, por exemplo, a velocidade média e a
vazão escoada Q, sabendo-se que Q = VA e A = área da seção transversal do conduto. De
forma similar, se for dado o valor de Q é possível calcular o diâmetro e, em uma terceira
situação, pode-se determinar o valor de para Q e D conhecidos (além das demais grandezas). O cálculo de z1, z2 ou L não requer o uso de métodos numéricos iterativos. No presente texto
foi elaborada uma planilha de cálculo (nesse caso utilizou-se a planilha Excell®) que
possibilita o cálculo das grandezas V, D ou . O passo inicial para a sua construção consiste
na definição de uma função objetivo “F” com a seguinte forma:
022
),,(
1
22
12
N
i
ig
VK
g
V
D
LfzzDVF
(3)
Para duas variáveis fixas, a determinação da terceira requer o uso de um método
iterativo. Uma alternativa clássica e robusta é o método de Newton, definido a seguir tomando
como exemplo o cálculo da velocidade média:
)('
)(1
k
kkk
VF
VFVV
(4)
Trata-se de um método iterativo, isto é, são necessárias iterações (repetições) até que a
solução seja encontrada. Vk+1 é o valor obtido com base em um valor Vk (k = 1, 2, ..., n, sendo
n = número total de iterações). Para k = 1, é necessário estabelecer uma estimativa inicial,
formalmente denominada de “semente” e, no Brasil, vulgarmente denominada de “chute
inicial” (analogamente a uma partida de futebol, cujo resultado depende da tática dos times e
não do chute inicial, que pode ser de qualquer intensidade e direção). Nota-se que o uso do
método de Newton exige a derivação da função objetivo. Nem sempre esta é uma tarefa
trivial, o que conduziu à elaboração de métodos denominados quase-Newton, que empregam
derivadas aproximadas.
A planilha desenvolvida utiliza o recurso solver (Excel®) que, entre os métodos
numéricos, possui o método de Newton (ou o método quase-Newton). Considerando o
problema de determinação da vazão, os dados conhecidos e demais informações para
construção da planilha podem ser vistos na Figura 2a. Nessa figura, as células A3 a H3
possuem dados de entrada para o problema (neste caso foram desprezadas as perdas
localizadas). A célula A8 apresenta o resultado correspondente à velocidade. Nesta célula é
necessário inserir a estimativa inicial para a velocidade média, ou seja, Vk com k = 1. A célula
B8 calcula o número de Reynolds, Re = VD/. A célula C8 calcula o fator de resistência por meio da equação 2. A célula D8 contém a função objetivo (equação 3) e a célula E8 calcula a
vazão com a velocidade média e a área da tubulação. A Figura 2b indica quais células devem
ser utilizadas para resolver o problema. A célula de destino corresponde à função objetivo,
que deve ser igual a zero. A célula variável, para este problema de determinação da vazão, é a
que contém a velocidade média (A8). No item “opções” é possível escolher o método de
Newton e outros detalhes relacionados à solução numérica, como número máximo de
iterações, tempo máximo de cálculo, etc. Seguindo o mesmo procedimento, pode-se utilizar a
mesma planilha para calcular uma das outras variáveis. Para D, por exemplo, fixa-se V e
atribui-se uma estimativa inicial para D, ou seja, Dk com k = 1. A célula de destino é a mesma
e a célula variável passa a ser D3. Procedimento similar pode ser feito para determinar a
rugosidade .
(a)
(b)
Figura 2 – Planilha para o cálculo da vazão: (a) Dados de entrada e células com as equações necessárias para
resolução do problema; (b) Parâmetros do Solver.
Pode ser interessante calcular vazões para diferentes diâmetros concomitantemente.
Para realizar esse tipo de cálculo a célula de destino continua sendo a função objetivo,
enfatizando-se que se usa apenas uma célula porque não é permitido que este campo seja
preenchido com mais de uma célula. As células variáveis são as diferentes velocidades para k
= 1. Para que os cálculos sejam realizados concomitantemente, pode-se estabelecer uma
restrição no botão “Adicionar”. A restrição deve incluir todas as células que contém a função
objetivo, igualando-as a zero. A Figura 3 ilustra um exemplo deste tipo de procedimento.
Nesta planilha (Figura 3a) as células D8 a D27 contêm a função objetivo. A Figura 3b mostra
a forma utilizada para encontrar as raízes da função para diferentes diâmetros
simultaneamente no item “Submeter às restrições”.
(a)
(b)
Figura 3 – Determinação de vazões para diferentes diâmetros: (a) Dados de entrada e células com as equações
para diferentes diâmetros; (b) Parâmetros do Solver.
ALTURAS CONJUGADAS DE UM RESSALTO HIDRÁULICO
O ressalto hidráulico é uma onda estacionária que ocorre freqüentemente na transição entre
um escoamento supercrítico e um escoamento subcrítico, como ilustrado na Figura 4. A
determinação da relação entre as alturas conjugadas e o número de Froude para canais
retangulares é realizada com a aplicação das equações de conservação de quantidade de
movimento e de conservação de massa. A equação resultante é simples e estabelece uma
relação explícita entre as variáveis adimensionais. A determinação das alturas conjugadas em
um canal com seção transversal trapezoidal, entretanto, é um pouco mais complicada. O uso
dos mesmos princípios básicos resulta na equação 5 (Hager e Wanoschek, 1987; Porto, 2006,
p.341), que relaciona Y, Fr1 e M para um canal trapezoidal com taludes com a mesma
inclinação 1V:ZH, em que: Y = h2/h1 (ver Fig. 4), Fr1 é o número de Froude calculado com h1
para uma seção transversal trapeziforme e a partir da sua definição com a altura média (área
molhada dividida pela largura de topo). M é uma razão de aspecto da seção transversal, sendo
definida como M = Zh1/b.
1
1
1
21
11
31
2
12
21
32
MYY
M
M
MFrY
MY
(5)
Figura 4 – Ressalto hidráulico: Fr = número de Froude, x = coordenada longitudinal, h1 = conjugado
supercrítico, h2 = conjugado subcrítico, hc = profundidade crítica.
Para um dado M, pode-se calcular Fr1 para diferentes valores de Y>1 por meio da raiz
positiva da equação 5 (é uma equação de segundo grau simples em Fr1). Desta forma, torna-se
simples construir um gráfico que relacione as três variáveis adimensionais. A mesma
simplicidade não existe para determinar o valor de Y para valores dados de M e Fr1, uma vez
que se trata de uma equação de 5º grau em Y, não havendo solução geral conhecida para
equações algébricas desse grau. Sendo assim, este exemplo de aplicação tem como objetivo
ilustrar a possibilidade de programação do método de Newton em linguagem C com o intuito
de obter soluções de Y para a equação 5. Os resultados podem ser comparados com gráficos
construídos com o procedimento mencionado ou com resultados calculados através do Solver.
Para a elaboração do código com o método de Newton é necessário derivar a função objetivo
definida a partir da equação 5 em relação a Y:
22
22
12
22
132
)(
)21)(1(
21
)1()('
011
1
21
11
31
2
1)(
MYY
MYM
M
MFrMYYYF
MYY
M
M
MFrY
MYYF
Código em C
O código escrito para este exemplo de aplicação pode ser compilado gratuitamente no
programa Dev-C++. A sua estrutura é apresentada na Tabela 1. Ao longo do código são
apresentados alguns comentários com o formato /*comentário*/. O resultado são os valores
de Y que resolvem o problema posto.
Tabela 1 – Código em C com o método de Newton para resolução da equação 5
#include <stdio.h> /*biblioteca*/
#include <math.h> /*biblioteca*/
int k, decisao; /*int declara variáveis inteiras, como k e decisão*/
float Q, yt, b, Z, Yo; /*dados de entrada, i.e., fornecidos pelo usuário, com ponto flutuante*/
float M, YY, fY, Beta, y2, Fr1, Fr2;
double Y[100000000], e; /*Y é um vetor, no sentido computacional*/
main(){
printf("\n");/*\n separa as linhas*/
printf(" PROGRAMA PARA O CALCULO DO CONJUGADO SUBCRITICO DE UM\n");
printf(" RESSALTO ESTABELECIDO EM UM CANAL TRAPEZOIDAL SIMETRICO\n");
printf(" E DE FUNDO HORIZONTAL: SOLUCAO DA EQUACAO 11.15 DO LIVRO\n");
printf(" HIDRAULICA BASICA: PORTO, R. M. (2006, p.341) \n\n");
printf(" - ATENCAO:\n\n");
printf(" ADOTAR SISTEMA INTERNACIONAL DE UNIDADES (SI)\n\n");
printf(" METODO NUMERICO DE NEWTON-RAPHSON\n");
printf("\n\n");
printf("\n\n");
printf("\n PARA DADOS DE ENTRADA: Q, b, y1, Z - DIGITE 1\n");
printf("\n PARA DADOS DE ENTRADA: M e Fr1 - DIGITE 2\n");
scanf("%d", &decisao); /*scanf armazena o valor digitado*/
printf("\n\n");
printf(" DIGITE OS DADOS DO PROBLEMA\n\n");
printf("\n");
if(decisao==1){
printf(" VAZAO (Q)\n\n");
scanf("%f", &Q);
printf("\n LARGURA DE FUNDO (b)\n\n");
scanf("%f", &b);
printf("\n PROFUNDIDADE SUPERCRITICA DO RESSALTO (y1)\n\n");
scanf("%f", &yt);
printf("\n INCLINACAO DO TALUDE (Z)\n\n");
scanf("%f", &Z);
printf("\n CHUTE INICIAL (Sugestao: 20)\n");
scanf("%f", &Yo);
}
if(decisao==2){
printf(" M\n\n");
scanf("%f", &M);
printf("\n Fr1\n\n");
scanf("%f", &Fr1);
printf("\n CHUTE INICIAL (Sugestao: 20)\n");
scanf("%f", &Yo);
}
printf("\n-------------------------Iteracoes----------------------------\n\n");
if(decisao==1){
M = Z*yt/b;
Fr1 = Q*pow((b+2*Z*yt), 0.5)/(pow(9.8, 0.5)*pow((b*yt+Z*yt*yt), 1.5));
}
printf("\n Y[k] f(Y[k]) = 0\n\n");
e = 1;
Beta = (pow(Fr1, 2)*pow((1+M), 2)/(1+2*M));
k=0;
Y[0]=Yo;
while(e > 0.0000001){
fY = 0.5-0.5*Y[k]*Y[k] +M/3-(M/3)*pow(Y[k], 3) - Beta*((1+M)/(Y[k]+M*Y[k]*Y[k])) + Beta;
YY = -Y[k]-M*pow(Y[k], 2)+Beta*(1+M)*(1+2*M*Y[k])/pow((Y[k]+M*Y[k]*Y[k]), 2);
Y[k+1] = Y[k] - fY/YY;
e = 0.5-0.5*Y[k+1]*Y[k+1] +M/3-(M/3)*pow(Y[k+1], 3) - Beta*((1+M)/(Y[k+1]+M*Y[k+1]*Y[k+1])) + Beta;
if(e<0){
e = e*(-1); /*Apenas uma forma de considerar o módulo de e*/
}
printf(" %5f %5f", Y[k], e); /*Mostra na tela os valores de cada iteração*/
printf("\n");
k = k+1;
}
if(decisao==1){
y2 = Y[k]*yt;
Fr2 = Q*pow((b+2*Z*y2), 0.5)/(pow(9.8, 0.5)*pow((b*y2+Z*y2*y2), 1.5));
}
printf("\n-------------------------Resultados----------------------------\n\n");
printf(" Y = %10f \n\n", Y[k]);
if(decisao==1){
printf(" M = %10f\n\n", M);
printf(" Fr1 = %10f\n\n", Fr1);
printf(" Fr2 = %10f\n\n", Fr2);
printf(" y2 = %10f m\n\n", y2);
}
printf("\n");
printf("\n\n\n Para finalizar o programa pressione ENTER \n\n");
getchar();
getchar(); }
Download gratuito em http://stoa.usp.br/hidraulica/files/
AJUSTE DE MODELOS MATEMÁTICOS
O estudo e compreensão de fenômenos físicos normalmente têm como resultado a construção
de modelos matemáticos que os representem quantitativamente. Evidentemente os fenômenos
físicos também podem levar a modelos não-quantitativos, como é o caso do ciclo da água, por
exemplo. Os modelos quantitativos freqüentemente carregam em sua dedução uma ou mais
constantes que devem ser calculadas a partir de informações experimentais e, de acordo com a
forma matemática do modelo, pode ser necessário o uso de métodos numéricos para obtenção
de soluções (ex. Brito et al., 2010; Schulz et al., 2010). Neste exemplo o solver é utilizado
para calcular a constante de resistência oferecida pelo ar presente em um modelo matemático
que representa a queda livre de um objeto.
As forças atuantes em um corpo em queda livre são: (1) força gravitacional=mg, em que
m=massa do corpo e g=aceleração da gravidade; (2) resistência oferecida pelo ar, modelada
aqui como sendo proporcional ao quadrado da velocidade: Fresistência=kV2, em que k é
supostamente constante; (3) empuxo de Arquimedes (desprezado). Com estas forças, a
segunda lei de Newton resulta na seguinte equação diferencial ordinária separável:
2Vm
kg
dt
dV
(6)
Esta equação pode ser resolvida analiticamente resultando em:
1
1
1.
/
1
1./
0
0
0
0
./..2
./..20
./..2
./..2
0
ttmgk
ttmgk
ttmgk
ttmgk
e
e
kmg
V
e
ekmgV
V
(7)
Nesta equação V0 é a velocidade no instante inicial t0. Considerando que em t = t0=0, V0 = 0,
a equação anterior assume a seguinte forma:
mkgt
mkgt
e
e
k
mgV
/2
/2
1
1
(8)
Se for realizado um experimento para medir a velocidade de um corpo em queda livre
em função do tempo, pode-se comparar a teoria representada pela equação 8 com os dados
experimentais. Para um problema genérico, se a forma assumida pelos pontos experimentais
Vexp(t) for próxima de uma reta e a forma prevista pela teoria para V(t) for uma parábola, por
exemplo, conclui-se que a teoria não é capaz de descrever o fenômeno ou os pontos
experimentais estão errados. Neste exemplo são apresentados dados de um experimento
hipotético correspondente à queda livre com o objetivo de comparar dados experimentais com
modelo físico-matemático. Sendo Vexp(ti) um valor da velocidade medida e V(ti) um valor da
velocidade calculada no mesmo instante ti, pretende-se minimizar a soma dos erros (Vexp(ti)-
V(ti))2. Como alternativa, sugere-se o uso do recurso Solver da planilha Excel
® (Chapra e
Canale, 2008, p.472). A forma de organização para realizar estes cálculos é apresentada na
Figura 5a. Nesta figura mostram-se uma coluna para o tempo, uma para a velocidade
experimental e uma terceira para a velocidade calculada. A quarta coluna calcula o erro
err=(Vexp(ti)-V(ti))2 para cada instante ti (i = 1,2,...j, em que j = número de pontos medidos).
Na última célula abaixo da quarta coluna é realizada a soma dos erros. Com o uso do solver, a
célula D88 é definida como célula de destino. Escolhe-se a opção “min igual a zero” e a
célula variável é aquela que contém o valor de k (célula B3), que pode ser estimado por
tentativas antes de utilizar o solver. O aspecto geral após o ajuste é ilustrado na Figura 5b.
Figura 5 – Cálculo da constante k presente na equação 8: (a) Uma possível forma de construir a planilha; (b)
Resultados experimentais (hipotéticos) e solução analítica.
SOLUÇÕES DE EQUAÇÕES DIFERENCIAIS ORDINÁRIAS (EDO)
O modelo matemático citado no exemplo sobre queda livre é obtido a partir da solução
analítica de uma equação diferencial ordinária, a equação 6. Neste exemplo essa EDO é
resolvida com os métodos de Runge-Kutta até 5ª ordem (isto é, cinco métodos distintos).
Além de ilustrar o uso dos métodos numéricos de resolução de equações diferenciais, este
exemplo serve para demonstrar o significado da ordem de um método numérico através da
comparação da solução numérica com a solução analítica. É utilizado o software Matlab® para
resolver a EDO com os cinco esquemas numéricos. Os dados utilizados foram: g = 9,806
m/s2, k = 0,3 [N.s
2/m
2], m = 50 kg, tt(tempo total) = 60 s e V(0) = 0.
Optou-se por apresentar aqui, por economia de espaço, a parte do programa escrita para
o método de Runge-Kutta de 5ª ordem. Os demais, desde o método de Euler (primeira ordem)
até o método de Runge-Kutta de 4ª ordem são mais simples e podem ser programados
seguindo procedimentos similares aos apresentados aqui. As primeiras linhas do programa
definem os valores das grandezas físicas que são constantes ao longo da solução, como g, k,
m e tt. Além disto, define-se também o espaçamento para a malha temporal, escrito como dt e
não t no programa. Feito isto, o número total de pontos nesta malha é calculado por tt/dt+1. %Queda livre de um corpo (Unidades de acordo com o S.I.) %método de Runge-Kutta de quinta ordem
0
20
40
60
0 2.5 5 7.5 10
V [m
/s]
t [s]
Experimental
Teórico
clear
g=9.806; k=0.3;%input('Coeficiente de arrasto: '); m=50;%input('Massa do corpo: '); tt=60;%input('Tempo total de simulação: '); Vo=0;%input('Velocidade inicial: '); dt=0.05;%input('Espaçamento da malha temporal: '); Nt=tt/dt+1;%Número de nós na malha temporal.
O próximo passo consiste em definir vetores para armazenar os resultados. Ao definir
V=zeros(1,Nt), o programa cria uma matriz com uma linha e Nt colunas com elementos iguais
a zero que, ao longo da solução, são substituídos pelos valores calculados. A matriz Va
armazena os valores calculados com a solução analítica. Para substituir em V e Va o primeiro
elemento por V(0), utiliza-se V(1,1)=Vo;. Finalmente, n=1 corresponde ao instante t = 0.
V=zeros(1,Nt); Va=zeros(1,Nt); V(1,1)=Vo; Va(1,1)=Vo; Vt=sqrt(m*g/k); n=1;
As linhas apresentadas a seguir correspondem ao laço necessário para que seja calculada
a solução para dtttt. Uma vez calculados os valores de m1, m2,...,m6, é possível calcular V(1,n), que é a velocidade no instante dt, correspondente a n = 2, pois t0 = 0, que corresponde
a n = 1. A solução analítica também é computada dentro deste laço.
for t=dt:dt:tt n=n+1; m1=dt*(g-(k/m)*V(1,n-1)^2); m2=dt*(g-(k/m)*(V(1,n-1)+m1/4)^2); m3=dt*(g-(k/m)*(V(1,n-1)+m1/8+m2/8)^2); m4=dt*(g-(k/m)*(V(1,n-1)-m2/2+m3)^2); m5=dt*(g-(k/m)*(V(1,n-1)+(3/16)*m1+(9/16)*m4)^2); m6=dt*(g-(k/m)*(V(1,n-1)-(3/7)*m1+(2/7)*m2+(12/7)*m3... -(12/7)*m4+(8/7)*m5)^2); V(1,n)=V(1,n-1)+(1/90)*(7*m1+32*m3+12*m4+32*m5+7*m6); Va(1,n)=(Vo+Vt*(exp(2*t*(k*g/m)^0.5)-1)/(exp(2*t*(k*g/m)^0.5)+1))/... (1+(Vo/Vt)*(exp(2*t*(k*g/m)^0.5)-
1)/(exp(2*t*(k*g/m)^0.5)+1));%Solução analítica end
O código a seguir é utilizado para visualizar os resultados. O comando plot constrói, no
mesmo plano, a solução numérica e a analítica. Com o uso deste método não é possível
vislumbrar as diferenças entre as curvas a partir dos valores absolutos de velocidade lançados
em gráfico com as amplitudes totais da variação da velocidade. Para perceber tal diferença,
sugere-se o uso do recurso “Zoom In”. O segundo plot mostra o erro absoluto para cada
instante. Finalmente, a variável err armazena o valor máximo do erro absoluto. É este valor
que deve ser utilizado para o cálculo da ordem de precisão. A Figura 6a ilustra a solução
numérica e a solução analítica (com a amplitude total de variação) e a Figura 6b contém os
erros absolutos para cada instante. A Figura 6c demonstra a ordem de convergência dos cinco
métodos por meio dos expoentes das potências ajustadas às relações entre erros absolutos e
espaçamento da malha temporal.
%Gráficos e análise da ordem de precisão: plot(0:dt:tt,V(1,:),0:dt:tt,Va(1,:)) xlabel('t [s]') ylabel('V [m/s]') pause plot(0:dt:tt,abs(V-Va)) xlabel('t [s]') ylabel('|V-Va| [m/s]') pause err=max(abs(V(1,:)-Va(1,:)));
(a) (b)
(c)
Figura 6 – Comparação entre valores medidos e experimentais: (a) Soluções analíticas e numéricas pra a
velocidade; (b) Erro absoluto em função do tempo para um valor de t; (c) Erro absoluto em função de t para 5
métodos de discretização (Euler – 1ª ordem e Runge-Kutta, de 2ª a 5ª ordens)
SOLUÇÕES DE EQUAÇÕES DIFERENCIAIS PARCIAIS
Equações diferenciais parciais são propostas com grande freqüência como modelos
representativos de fenômenos físicos. Neste exemplo é calculada a distribuição bidimensional
de temperaturas em uma placa com as bordas aquecidas a 120 oC e com um núcleo a uma
temperatura igual a 160 oC. As soluções obtidas incluem o uso da planilha Excel
® e do
Matlab.
0 10 20 30 40 50 600
5
10
15
20
25
30
35
40
45
t [s]
V [
m/s
]
V
Va
0 10 20 30 40 50 600
0.2
0.4
0.6
0.8
1
1.2
1.4x 10
-11
t [s]
|V-V
a| [m
/s]
y = 1,9247x1,0224
R² = 1
y = 0,3543x2,0747
R² = 0,9999
y = 0,0272x3,0342
R² = 1
y = 0,0025x4,0728
R² = 1
y = 6E-05x5,0963
R² = 1
1,E-11
1,E-09
1,E-07
1,E-05
1,E-03
1,E-01
1,E+01
0,05 0,5
y=
max|V
-Va|
x=t
Euler RK2 RK3 RK4 RK5
A equação de conservação de energia pode ser simplificada de tal maneira que o modelo
resultante é uma equação de Laplace com a temperatura como variável dependente, do tipo:
T=0 ( = laplaciano). Para tanto, o regime deve ser permanente, não deve haver fontes ou
sumidouros e o meio deve estar parado, ou seja, Vi = 0. Considerando o caso bidimensional, a
equação T=0 assume a seguinte forma:
02
2
2
2
y
T
x
T
(9)
Utilizando uma aproximação por diferenças finitas centradas, a forma discretizada da
equação passa a ser aquela apresentada a seguir. A malha correspondente à discretização pode
ser vista na Figura 7.
022
2
1,,1,
2
,1,,1
y
TTT
x
TTT jijijijijiji
(10)
Figura 7 – Malha bidimensional
Para calcular a solução numérica T(x,y), a equação pode ser reescrita com a seguinte
forma:
1,1,2,1,12
1
22,1111
2 jijijijiji TTy
TTxyx
T
(11)
A equação 11 é uma forma conveniente para uso na planilha Excel®, por exemplo. No
presente exercício considera-se uma chapa metálica com a geometria indicada na Figura 8.
Para as temperaturas apontadas nos contornos, o objetivo é calcular as temperaturas nos nós
internos utilizando um recurso computacional.
Figura 8 – Geometria da placa metálica (T em oC)
Solução com Excel®
Para resolver o problema com esta planilha eletrônica é utilizada a equação 11 para os
nós internos. Não é necessário digitar a equação para cada nó, mas apenas para um nó e, em
seguida, os recursos copiar e colar podem ser utilizados. O leitor pode acessar a planilha
T(x,y)_Laplace em http://stoa.usp.br/hidraulica/files/ e ver um exemplo resolvido. Para
construir esta planilha, sugere-se que sejam seguidos os passos: (1) Adotar valores para os
lados e os espaçamentos para a malha; (2) Desenhar a geometria impondo as condições de
contorno; (3) Acessar “Opções do Excel”, “Fórmulas” e marcar a opção “Habilitar cálculo
iterativo”; (4) Digitar a equação 11 em um dos nós internos; (5) Copiar e colar a célula que
contém a equação 11 nos demais nós internos. É mais rápido arrastar a célula; (6) Finalmente,
a tecla F9 deve ser pressionada até que não sejam vistas mais alterações nas temperaturas.
A visualização dos resultados com os recursos desta planilha é limitada. Sugere-se o uso
de outro recurso, como o Matlab® ou o Octave. Para utilizar o Matlab
®, pode-se copiar e colar
os valores das temperaturas em outra planilha (apenas valores e não fórmulas). Em seguida,
esses valores devem ser copiados e colados no Notepad, assegurando-se de utilizar pontos
como separadores decimais. Ao abrir o Matlab®, acessar File, Import Data... e seguir as
instruções. O software deve criar uma matriz com os dados e com o nome do arquivo de
texto. O comando contourf(T) pode ser empregado para gerar a Figura 9a e
mesh(T,'FaceColor','c','EdgeColor','none'); camlight right; lighting phong produz o resultado
da Figura 9b. Além da limitação gráfica, se o número de células for muito grande, é possível
que a planilha não funcione adequadamente.
(a) (b)
Figura 9 – Visualização dos resultados (as cores em (a) correspondem a T [oC])
10 20 30 40 50
5
10
15
20
25
30
35
40
45
50
120
125
130
135
140
145
150
155
160
Solução com os softwares Octave ou Matlab®
A equação de Laplace é linear e a sua forma discreta pode ser escrita como Ax=b. O uso
destes softwares é interessante porque não exige que seja programado o algoritmo para a
solução do sistema Ax=b, sendo necessário apenas utilizar o comando x = A\b. O código
desenvolvido, denominado POISSON_MLAB_ex7, tem como função principal escrever as
matrizes A e b de acordo com a geometria do problema. Para que ele funcione é necessário
que o código ij2n esteja salvo na mesma pasta. Este código contém uma função capaz de
alocar adequadamente os elementos das matrizes em função da geometria. O nome parece não
ser adequado, pois a equação resolvida aqui é a de Laplace não de Poisson. Entretanto, para
que a equação de Poisson seja reduzida à equação de Laplace basta que a função f(x,y), de
f(x,y)=T, seja igual a zero. Esta função no código é denominada gradp, em referência ao
programa original, destinado ao cálculo de perfis de velocidades completamente
desenvolvidos (sendo gradp o gradiente de pressões no código original).
Para utilizar o programa é preciso inserir as seguintes informações: NI e NJ (número de
pontos nos eixos x e y), Lx e Ly (lados do domínio retangular), boundary-value = 120
(temperatura nos contornos externos), casemask = 2 (relacionado à geometria) e solid_value =
160 (temperatura do sólido interno). A geometria da chapa é desenhada impondo os valores
para a temperatura nos contornos do domínio. Isto é realizado com o uso de uma máscara,
escrita da seguinte maneira:
for i=2:NI-1 x = (i-1)*dx; for j=2:NJ-1 y = (j-1)*dy; if(y>=0.0 && y<=0.066 && x>=0.0 && x<=0.066) mask(i,j) = -1; end if((y >=0.234 && y<=Ly) && (x>=0.0 && x<=0.066)) mask(i,j) = -1; end if((y >=0 && y<=0.066) && (x>=0.234 && x<=Lx)) mask(i,j) = -1; end if((y >=0.234 && y<=Ly) && (x>=0.234 && x<=Lx)) mask(i,j) = -1; end if((y >=18*dy && y<=32*dy) && (x>=21*dx && x<=29*dx)) mask(i,j) = 0; end end end
Os dois laços do tipo “for” percorrem todos os nós internos. Os “if’s” são utilizados
para escolher posições específicas no domínio, ou seja, para escolher certos nós nos quais
serão fixados valores para a temperatura. Valores lógicos são utilizados para identificar
contornos internos, como -1 e 0. Se mask(i,j)=-1, a temperatura será igual a 120oC. Se for
zero, T=160oC. Estas informações não estão presentes neste trecho do código, mas sim no
trecho apresentado a seguir:
for i=1:NI
for j=1:NJ nij = ij2n(i,j,NI,NJ); % Aqui é utilizada a função ij2n if(mask(i,j) == 1) %Este é um nó interno no qual a velocidade é
calculada. nim1j = ij2n(i-1,j,NI,NJ); nip1j = ij2n(i+1,j,NI,NJ); nijm1 = ij2n(i,j-1,NI,NJ); nijp1 = ij2n(i,j+1,NI,NJ); Af(nij,nij) = val_c; Af(nij,nim1j) = val_l; Af(nij,nip1j) = val_r; Af(nij,nijm1) = val_b; Af(nij,nijp1) = val_t; elseif (mask(i,j) == 0) %T=160 Af(nij,nij) = 1.0; b(nij) = solid_value; elseif (mask(i,j) == -1) %T=120 nos cantos Af(nij,nij) = 1.0; b(nij) = boundary_value; else end end end
O aspecto da distribuição de temperaturas obtidas com este código é semelhante ao da
Figura 9 e, os erros absolutos entre os resultados obtidos no Excel® e com este código
apresentaram valor máximo igual a 0,0477 oC, com a distribuição mostrada na Figura 10.
Figura 10 – Erro absoluto calculado com os resultados do Excel e Matlab
®.
Com o intuito de ilustrar o uso do programa para uma geometria diferente da anterior,
foi simulada a distribuição de temperaturas em um triângulo com comprimento da base igual
a 1 m, lado esquerdo definido por y = 2x e lado direito por y = -2x+2. Inicialmente, define-se
os lados do domínio com Lx=Ly=1. A temperatura no contorno inferior (base do triângulo) é
definida em boundary_value (para este exemplo, igual a 20 oC). Os primeiros quatro laços
após esta definição devem ser alterados para a seguinte forma:
%Lado esquerdo for j=1:NJ mask(1,j) = 0;%-1; ndir=j+1; end %Lado direito
for j=1:NJ mask(NI,j) = 0;%-1; ndir=j+1; end %Lado inferior for i=1:NI mask(i,1) = -1; ndir=i+1; end %Lado superior for i=1:NI mask(i,NJ) = 0;%-1; ndir=i+1; end
As definições anteriores garantem que as temperaturas nos contornos esquerdo,
superior e direito sejam iguais à temperatura nos lados inclinados do triângulo. Escolhendo
casemask=3 e inserindo as informações dadas, obtém-se a distribuição de temperaturas
apresentada na Figura 12a, para uma malha com NI=NJ=91. É interessante notar que, mesmo
para este refinamento da malha, os lados inclinados do triângulo não são ajustados
perfeitamente pela malha. A partir do campo escalar de temperaturas é possível calcular o
campo vetorial do fluxo de calor cujas componentes são qx=-aT/x e qy=-aT/y. Para tanto, utiliza-se os seguintes comandos:
[DX,DY]=gradient(T,dx,dy); quiver(X,Y,-a*DX,-a*DY) xlim([0 1]) ylim([0 1])
Com este recurso, o software utiliza diferença finita centrada para os nós internos e
avançada para os nós dos contornos. A imagem (b) da Figura 12 mostra o campo de fluxo, a
imagem 12c contém um detalhe deste campo e a (d) apresenta o detalhe junto com o gráfico
de contornos.
(a) (b)
X= 0.67778
Y= 0.92222
Level= 200
X= 0.46667
Y= 0.88889
Level= 199.9132
0 0.2 0.4 0.6 0.8 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
20
40
60
80
100
120
140
160
180
200
0 0.2 0.4 0.6 0.8 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(c) (d)
Figura 12 – T(x,y) para uma geometria triangular e campo de fluxos. a) campo de temperaturas geral, b) campo
de fluxo geral, c) campo de fluxo no vértice esquerdo, d) campo de fluxo e de temperaturas no vértice esquerdo.
Dados: NI=NJ=91 e as cores indicam a temperatura.
CONCLUSÕES
Este trabalho teve como objetivo reunir problemas de engenharia relacionados aos fenômenos
de transporte e à hidráulica passíveis de resolução com o uso de recursos computacionais e
métodos numéricos de fácil alcance, e, portanto, com alto potencial didático. Explorou-se o
uso de softwares comerciais de custo relativamente baixo e o uso de softwares gratuitos. O
problema voltado ao cálculo de vazões entre dois reservatórios ilustrou a aplicação de uma
planilha eletrônica para obtenção de raízes de equações não lineares. Embora seja um caso
simples pertencente à engenharia hidráulica, a determinação da vazão escoada entre dois
reservatórios e das vazões para diferentes diâmetros mostra um caminho simples e prático
para resolução de sistemas não lineares. Um problema matematicamente semelhante
relacionado às características hidráulicas de ressaltos hidráulicos possibilitou o uso do método
de Newton em um código escrito em linguagem C. Em seguida, foi apresentado o uso de um
método de ajuste que possibilita minimizar os erros quadráticos entre valores medidos e
experimentais. O exemplo relacionado ao uso de discretizações de EDO’s pelos métodos de
Runge-Kutta ilustrou o uso do Matlab® e, o que é considerado mais importante nesse nível
didático, o significado da ordem de precisão de um método numérico de discretização. Ao
abordar o problema de transferência de calor por meio da solução de uma equação diferencial
parcial, exploraram-se aspectos básicos de programação numérica para equações diferenciais
parciais e formas de visualização normalmente empregadas em recursos computacionais
avançados.
REFERÊNCIAS
Chapra, S.C.; Canale, R.P. (2008) Métodos numéricos para engenharia. McGrawHill.
Chue, S.H. (1984) A pipe skin friction Law of universal applicability. Proc. Instn. Civ. Engrs.,
Part 2, 77, Mar., p.43-84.
Churchill, S.W. (1977) Friction factor equation spans all fluid regimes, Chem. Engng., 84, 7
Nov., p.91-92.
Brito, R.J.R., Simões, A.L.A., Schulz, H.E., Porto, R.M., Puche, A.A.S. (2010). Incorporação
de ar por aeradores de fundo em escoamentos de alta velocidade em vertedores: análise de
0.1 0.15 0.2 0.25 0.3
0.15
0.2
0.25
0.3
0.35
0 0.05 0.1 0.15 0.2 0.250
0.05
0.1
0.15
0.2
0.25
dados e modelos. In: Congresso Latinoamericano de Hidráulica, 24., Punta del Este,
Memorias... Punta del Este, Uruguay.
Hager, Willi H. and Wanoschek, Robert (1987). Hydraulic jump in triangular channel, Journal
of Hydraulic Research, 25:5, 549-564.
Pereira, A.J.S.; Almeida, A.B. (1986) Formulação explícita e universal da resistência em
tubos, XII Congresso Latino-Americano de Hidráulica, São Paulo.
Porto, R.M. (2006) Hidráulica básica. Projeto Reenge – EESC – USP.
Schulz, H.E. (2003) O essencial em fenômenos de transporte. Projeto Reenge – EESC – USP.
Schulz, H.E., Brito, R.J.R., Simões, A.L.A., Porto, R.M., Puche, A.A.S., Lobosco, R.J.
(2010). Modelos teóricos para análise de aeração em aeradores de fundo utilizando princípios
físicos. In: Congresso Latinoamericano de Hidráulica, 24., Punta del Este, Memorias... Punta
del Este, Uruguay.
Simões, A.L.A.; Schulz, H.E.; Porto, R.M. (2011) Equações e princípios básicos de mecânica
dos fluidos e hidráulica (Livro ainda não publicado).
Swamee, P.K. (1993) Design of a submarine pipeline. J.Transp. Eng. ASCE, 119(1), p.159-
170.
