Nono Simpósio de Mecânica Computacional
26 a 28 de maio de 2010
Universidade Federal de São João Del-Rei – MG
Associação Brasileira de Métodos Computacionais em Engenharia
Aplicação do Método dos Elementos de Contorno com Dupla
Reciprocidade em Problemas Difusivos-Advectivos Não Lineares
C. F. Loeffler1; F. P. Neves.2; P. C. Oliveira3
1
Programa de Pós-Graduação em Engenharia Mecânica
Universidade Federal do Espírito Santo, Vitória, ES – Brasil
CEP 29075-910
e-mail:[email protected]
2
Corpo de Bombeiros Militar do Espírito Santo, Vitória, ES – Brasil
e-mail: [email protected]
3
Universidade Federal do Espírito Santo, Alegre, ES – Brasil
e-mail: [email protected]
Resumo: Neste trabalho desenvolve-se um modelo numérico para simular
computacionalmente a distribuição de pressões, velocidades, temperaturas e fluxos de
calor estacionários em volumes de controle bidimensionais. O campo de temperaturas e
velocidades é governado pela equação da Difusão-Advecção, resolvida através da
formulação com Dupla Reciprocidade do Método dos Elementos de Contorno. Admite-se a
lei de Darcy para representar a Equação de Momentum, associando pressão e velocidade,
que no caso linear consiste num modelo matemático dado pela Equação de Laplace. A
análise não-linear é caracterizada pela dependência do campo de velocidade com a
temperatura, resultando num modelo expresso pela Equação de Poisson. Os resultados
para a velocidade são, então, implementados na Equação da Difusão-Advecção, gerando
temperaturas e fluxos de calor. Estes valores são re-introduzidos na Equação de Poisson
para gerar novas pressões e campos de velocidade. O processo é repetido até que haja
convergência dos resultados para valores estáveis.
Palavras chaves: Método dos Elementos de Contorno; Problemas Difusivos-Advectivos;
Formulação com Dupla Reciprocidade; Modelos Darcianos em Escoamentos.
Nono Simpósio de Mecânica Computacional
1
Universidade Federal de São João Del-Rei – MG – ABMEC
INTRODUÇÃO
São muito freqüentes na engenharia problemas físicos que envolvem transporte de
massa ou energia através de difusão associada à advecção, seja em escoamentos, misturas
ou deslocamento de particulados. Os casos mais conhecidos consistem na transferência de
calor junto à camada limite do fluido em escoamentos laminares, típico em aletas de
trocadores de calor e em aerofólios, e do transporte de fluidos incompressíveis com baixa
viscosidade em tubulações sujeitas a diferenciais de temperatura. Existem ainda situações
menos tradicionais, mas não menos importantes, como a dispersão superficial de poluentes
ou misturas em meio aquosos homogêneos, a secagem de produtos agrícolas, a absorção de
líquidos em região não saturada ou a extração de óleo em meio poroso, este último
problema atualmente de grande interesse na indústria de petróleo. Outro caso interessante
consiste dos problemas de aeração em ambientes fechados, onde o conhecimento das
regiões de estagnação e os principais pontos de insuflamento são importantes.
No caso geral, à complexidade dos fenômenos físicos associados a esses problemas
corresponde uma modelagem matemática igualmente elaborada, composta por equações
diferenciais parciais não-lineares, que requerem necessariamente o emprego de métodos
numéricos para sua solução aproximada. No entanto, o tratamento desses casos usualmente
é substituído pela adoção de modelos simplificados, que viabilizam sua solução e ao
mesmo tempo atendam às necessidades práticas. Nesse trabalho, adota-se um modelo nãolinear mais acessível para a equação de momentum linear, gerado a partir da Lei de Darcy,
que resulta numa relação escalar entre velocidades e o campo de pressões atuantes.
2
PROBLEMA FÍSICO E EQUAÇÕES BÁSICAS
O domínio físico compõe-se de um volume de controle bidimensional Ω(x,y), dentro
do qual atua um campo de velocidades e pressões e em cujas fronteiras Γ(x,y) são
prescritas temperaturas ou fluxos. A figura 1 ilustra tais características:
Figura 1: Volume de Controle atravessado por campo de velocidades (vx,vy)
Considerando-se escoamentos incompressíveis homogêneos, em regime permanente e
com propriedades isotrópicas, a Equação da Continuidade (White, 1986) se expressa por:
∂v x ∂v y
+
=0
∂x
∂y
(1)
Na equação anterior, vx e vy são as componentes do vetor velocidade V do
escoamento, enquanto x e y caracterizam as coordenadas globais do volume de controle. A
Equação da Energia (White, 1986), no caso é dada pela Equação da Difusão-Advecção:
ρC v [
∂T
∂T
∂ 2T ∂ 2T
vx +
v y ] = K[ 2 + 2 ]
∂x
∂y
∂x
∂y
(2)
Nono Simpósio de Mecânica Computacional
Universidade Federal de São João Del-Rei – MG – ABMEC
Onde Cv é o calor específico a volume constante e K é a condutividade térmica do
meio fluido. A Equação do Momentum (Özisik, 1968), admitindo-se presença de forças de
corpo bx e by, escreve-se como:
∂V
∂V
ρ[
vx +
v y ] = −∇p + µ∇ 2 V + ρb
∂x
∂y
(3)
Onde ρ é a massa específica, p é o campo de pressão no fluido e µ é a viscosidade. A
derivada total presente no lado esquerdo da equação (3) é proporcional à força requerida
para sobrepujar a inércia do elemento fluido e esse termo pode ser omitido das equações
quando se admite o fluxo viscoso, ou seja, escoamentos com Número de Reynolds (Ra)
bastante baixos (Bejan, 1993). Nessas condições, a mencionada equação para um fluido
incompressível se simplifica e é dada por:
∇p = ρb + µ∇ 2 V
(4)
Tal hipótese se aplica na descrição do escoamento de partículas com massa
desprezível através de um meio fluido. Essa hipótese também é pertinente no
equacionamento dos problemas de escoamento em meios porosos, nos quais pode se
aplicar a lei de Darcy.
3
MODELO DARCIANO
Admitindo-se que a porosidade seja relativamente baixa no meio, vigora a Lei de
Darcy (Bennett and Myers, 1982). Ressalta-se que tal lei é válida somente para fluxo lento
e viscoso. Nessas condições, pode-se escrever:
µΦ V = ρb − ∇p
χ
(5)
Onde χ é a permeabilidade do meio e Φ é a porosidade. Assim, nas condições de
escoamento laminar, da comparação entre o valor do campo de velocidade dado pela Lei
de Darcy e a expressão correlata dada pela Equação de Momentum para o escoamento
viscoso incompressível, tem-se:
φ ∇2 V = − V
χ
4
(6)
APROXIMAÇÃO DE BOUSISINESQ
Para poros rochosos saturados, a aproximação de Boussinesq (Zhao et al, 1996) pode
ser empregada para descrever a mudança na massa específica dos poros devido à mudança
na temperatura. Tal aproximação estabelece que nos casos em que atua a força
gravitacional, as alterações na massa específica não são pequenas suficientemente para
serem desprezadas. Expressando a massa específica como uma função da temperatura:
ρF = ρO [1 −β(T − TO )]
(7)
Onde ρF é a massa específica a uma temperatura TF qualquer, enquanto ρO é a massa
específica na temperatura TO de referência do meio e β é o coeficiente de expansão térmica
do meio poroso. Considerando-se agora a variação da massa específica, tem-se:
µΦ V = ρo [1 − β(T − To )]b − ∇p
χ
Adimensionalizando-se a equação anterior, seus componentes reescrevem-se como:
(8)
Nono Simpósio de Mecânica Computacional
v*X = R AX T * −
∂p *
,
∂x *
Universidade Federal de São João Del-Rei – MG – ABMEC
v*Y = R AY T * −
∂p *
∂y *
(9)
A equação de Momentum fica, então, dada por:
∇ 2 p* = −(R XA
5
∂T *
∂T *
+ R YA
)
∂x
∂y
(10)
SOLUÇÃO DO PROBLEMA ACOPLADO
Devido à aproximação de Boussinesq, a equação da Difusão-Advecção está acoplada à
Equação de Poisson, uma vez que ambas dependem da temperatura. Isto significa que o
problema como um todo não é linear. Aqui se utiliza um processo iterativo para resolver
esse acoplamento, de modo que para iniciar o processo, toma-se o lado direito da Eq. (10)
como um termo fonte, em que os fluxos iniciais são identificados de acordo com alguma
referência. Resolvendo essa equação, determina-se o campo de pressões e velocidades
necessário para compor a equação da Difusão-Advecção. Obtidos novos valores de
temperaturas e fluxos, o termo de fonte previamente arbitrado é reajustado e a Equação de
Poisson é resolvida novamente, fornecendo novos valores de pressões e velocidades. O
processo continua até que as magnitudes das grandezas envolvidas já não sofram alterações
significativas entre duas iterações sucessivas.
6
FORMULAÇÃO DO MEC PARA A EQUAÇÃO DE MOMENTUM
O tratamento matemático tradicional com o Método dos Elementos de Contorno
(MEC) para o lado esquerdo da Equação de Momentum é bem conhecido na literatura
(Brebbia, 1978). Consiste do estabelecimento do modelo numa forma integral, usando-se
uma função auxiliar, u*(ξ;x,y), a solução fundamental. Por ser um operador auto-adjunto,
o Laplaciano permite a aplicação dupla do esquema de integração por partes, ou seja:
∫[
Ω
∂2p * ∂2p *
∂2u * ∂2u *
+
]u * dΩ = ∫ [ 2 +
]p * dΩ − ∫ {p * q*}dΓ + ∫ {v * u*}dΓ
2
2
∂x
∂y
∂x
∂y 2
Ω
Γ
Γ
(11)
Onde foi aplicado devidamente o Teorema da Divergência. Nessa última equação
figuram:
q* =
∂u * ∂x ∂u * ∂y ∂u *
+
=
,
∂x ∂n ∂y ∂n ∂n
v* =
∂p * ∂x ∂p * ∂y ∂p *
+
=
∂x ∂n ∂y ∂n ∂n
(12)
Devido às propriedades da solução fundamental u*(ξ;x,y), sua substituição no lado
direito da equação (13) permite transformá-la numa expressão constituída de integrais de
contorno e uma função de ponto, tal que a Equação de Poisson na forma integral fica:
C(ξ)p * (ξ) + ∫ {p * q*}dΓ − ∫ {v * u*}dΓ = ∫ qu *dΩ
Γ
Γ
(13)
Ω
O coeficiente C(ξ) está ligado ao posicionamento do ponto fonte ξ com relação ao
domínio físico Ω(x,y). Para contornos suaves, se o ponto fonte se situa sobre o contorno,
C(ξ)=0,5. Já o lado direito da equação (17) pode ser interpretado como um termo fonte e
para esse existem algumas alternativas de abordagem. A mais elementar é o uso de células,
muito similar ao utilizado pelo Método dos Elementos Finitos, mas que requer a
discretização do domínio. Uma opção mais elegante e já bem consolidada com o MEC é o
uso da Formulação com Dupla Reciprocidade (FDR) (Partridge et al, 1992), que interpola
a função distribuída no domínio empregando usualmente funções radiais de base e
Nono Simpósio de Mecânica Computacional
Universidade Federal de São João Del-Rei – MG – ABMEC
escolhendo adequadamente tais funções de modo a se transformar todas as integrais de
domínio em integrais de contorno. A FDR propõe substituir a ação de domínio q(x,y) por
uma soma finita composta por novas funções auxiliares αjFj(x,y), onde αj são coeficientes
inicialmente desconhecidos e as Fj são funções de aproximação, cuja quantificação
depende de pontos Xj escolhidos tanto no contorno, quanto no interior do domínio. Assim,
o lado direito da equação (17) pode ser reescrito como se segue:
 j *

*
j
j *
j
j *
j
(14)
∫Ω qu dΩ = Ω∫ α Ψ , ii u dΩ =α  ∫Γ η u dΓ - ∫Γ Ψ q dΓ-C ( ξ ) Ψ ( ξ )
Ressalta-se que ηj e ψj são funções conhecidas dependentes apenas das funções F j ,
que podem ser escolhidas arbitrariamente. A escolha das funções interpolantes Fj é
conteúdo de muitas pesquisas, mas a classe de funções mais flexível é a das funções radiais
de base, devido à sua flexibilidade, invariância, entre outras propriedades.
7
FORMULAÇÃO DO MEC PARA A EQUAÇÃO DA ENERGIA
Existem formulações outras do MEC para o tratamento de problemas difusivosadvectivos (Ramachandran, 1994). No entanto, para o tratamento da Equação da Energia
também foi escolhida a FDR, pois além de não impor restrições ao campo de velocidades,
aproveita-se a simplicidade matemática da formulação difusiva para se resolver o termo
advectivo, conforme procedimento similar ao exposto anteriormente. Deste modo, admitese que o termo advectivo seja considerado tal como uma ação de domínio g(x,y), ou seja:
∫
Ω
8


gu ∗dΩ = γ j  ∫ η j X j ; X u* ( ξ; X ) − Ψ j X j ; X q* ( ξ; X ) dΓ − C ( ξ ) Ψ j ( ξ ) 
Γ

(
)
(
)
(15)
DISCRETIZAÇÃO E FORMAÇÃO DAS EQUAÇÕES MATRICIAIS
O procedimento de discretização com o MEC é bastante conhecido e muito simples.
Divide-se o contorno Γ(x) em N elementos discretos, nos quais a variável básica e derivada
normal e a forma geométrica dos elementos, todos são aproximados por funções típicas de
interpolação. Usualmente, toma-se o ponto fonte coincidente com os pontos nodais,
gerando um conjunto de equações resolvível. No caso da FDR, também as funções de
interpolação são aproximadas ao longo dos elementos de contorno. No caso da Equação de
Poisson, o lado esquerdo gera as tradicionais matrizes H e G referentes ao problema de
Laplace (Brebbia, 1978); no lado direito, graças à reciprocidade das integrais, tem-se:
Hp − Gv = {H Ψ − G η}α = {H Ψ − Gη}F −1q = Q
(16)
Onde foi escolhido um número de funções Fj igual ao número de nós de discretização
para substituição do vetor de coeficientes α. Na Equação da Energia, o procedimento é
similar:
HT − Gq =
ρC v
{HΨ − Gη}γ
K
(17)
No entanto, a explicitação dos coeficientes α agora é mais elaborada:
HT − Gq =
ρCv
ρC v
ρC v
{HΨ − Gη}γ =
{HΨ − Gη}F−1g =
[ HΨ − Gη] F−1Vi T,i
K
K
K
(18)
As derivadas espaciais da temperatura podem ser eliminadas simplesmente por
interpolação, de modo que o equacionamento matricial definitivo pode ser escrito como:
Nono Simpósio de Mecânica Computacional
HT - Gq =
Universidade Federal de São João Del-Rei – MG – ABMEC
ρCv
[ HΨ - GN] F-1Vi F,i F-1T = RT → (HT - RT) = Gq
K
(19)
Ressalta-se que o campo de velocidades deve ser completamente conhecido nas
fronteiras do volume de controle; considerando que a solução da Equação de Poisson
fornece apenas os valores das velocidades normais, as demais componentes precisam ser
calculadas para compor o campo de velocidades na Equação da Difusão-Advecção. Tais
velocidades podem ser obtidas usando-se funções de base radial, ou seja:
p,i = vi = F,i β = F,i F−1p
(20)
O mesmo procedimento é requerido para a determinação das derivadas tangenciais da
temperatura, que compõem o termo fonte da Equação de Poisson. Também os valores dos
fluxos em pontos internos necessitam ser calculados, especialmente de maneira a melhorar
a precisão do procedimento interpolativo da FDR na solução da Equação de Poisson.
9
EXEMPLO
O exemplo consiste de um silo vertical, no qual atuam campos de pressão, velocidade
e temperatura, conforme mostra a Figura 2 (a-c). Na composição da malha foram
empregados 40 elementos de contorno constantes, de mesmo tamanho. A numeração dos
elementos é crescente no sentido trigonométrico ao longo do contorno, a partir da origem.
Para melhor entendimento dos resultados, as faces do volume de controle são devidamente
numeradas, conforme destacado na Figura 2(d).
Figura 2: (a) características geométricas; (b) pressões prescritas no contorno;
(c) temperaturas prescritas no contorno; (d) especificação das diferentes faces do silo.
Conforme mencionado previamente, o problema é não-linear. Duas partes distintas
compõem a solução completa, na qual o campo de velocidades é obtido inicialmente e
então aplicado na equação da Difusão-advecção para determinação do campo de
temperaturas e fluxos. Novos valores dos fluxos integram o termo fonte da Equação de
Poisson. Para o problema proposto os valores de temperatura são adimensionais, sendo o
valor de referência igual a 100oC. O número de Rayleigh (RA) considerado é igual a 5.
As curvas presentes na Figura 3 mostram os valores da componente normal das
velocidades atuantes na face 3, calculadas para as três iterações realizadas. Da primeira
para a segunda iteração, é significativa a diferença dos valores obtidos, devido à inserção
do efeito gravitacional. No entanto, não foi necessário um maior número de interações com
o valor do RA adotado, pois nenhuma diferença significativa ocorreu da segunda para a
terceira iteração. Ainda sobre a Figura 3, pode ser percebido que a presença do efeito
gravitacional reduziu a intensidade das componentes normais da velocidade ao longo da
face, conforme esperado.
Nono Simpósio de Mecânica Computacional
Universidade Federal de São João Del-Rei – MG – ABMEC
Figura 3: Velocidades Normais nos pontos nodais localizados na face 3.
O efeito gravitacional também reduz a intensidade da pressão no fundo do silo (face 1),
conforme mostrado na Figura 4. A queda na magnitude dos valores da pressão ao longo da
direção vertical se deve ao distanciamento da entrada com maior valor de pressão,
localizada na face 4.
Figura 4: Valores de Pressão nos pontos nodais localizados na face 1.
Figura 5: Valores do Fluxo Difusivo nos pontos nodais localizados na face 3.
Diferentemente do cálculo das pressões e velocidades normais, a solução do problema
de temperaturas e fluxos praticamente não foi afetada pela ação gravitacional, tal como
mostra a figura 5, onde os resultados para as iterações realizadas praticamente coincidem.
Nono Simpósio de Mecânica Computacional
Universidade Federal de São João Del-Rei – MG – ABMEC
10 CONCLUSÕES
Muitos experimentos numéricos preliminares mostraram que o MEC se aplica com
precisão necessária em problemas governados separadamente pela Equação de Laplace,
Poisson e da Difusão-Advecção (Ramachandran, 1994), este último modelo considerando
casos em que o número de Peclet não é elevado, ou seja, quando o fenômeno difusivo não
é desprezível. Assim, esperava-se que o método poderia ser aplicado com sucesso em
casos onde houvesse a presença de não-linearidades como as produzidas pelo acoplamento
entre o modelo difusivo-advectivo e o campo de pressão, gerado pela inclusão da ação
gravitacional. Este é um problema que normalmente se manifesta no caso de escoamentos
em meios porosos.
De fato, a aplicação do FDR foi satisfatória, pois se conseguiram resultados
aproximados com facilidade e rapidez para as diversas variáveis de estado que compõem o
problema físico proposto, tais como os valores da velocidade normal e tangencial,
pressões, temperaturas e fluxos de calor, normal e tangencial, seja nas paredes do silo
como no seu interior. Por razões de espaço, apenas uma pequena parte deles pode ser aqui
apresentada. Pela mesma razão não pode ser mostrada a comparação dos resultados do
MEC com o Método dos Volumes Finitos, que tiveram excelente concordância, exceto
junto aos pontos nodais dos elementos dos cantos.
Destaca-se que para o modelo numérico proposto, o custo da solução computacional
foi muito reduzido. Ressalta-se ainda a facilidade de entrada de dados e interatividade entre
os resultados obtidos pela solução das equações de Poisson e Difusão-Advecção.
11 BIBLIOGRAFIA
Bejan, A., 1993, “Heat Transfer”, John Wiley & Sons, New York.
Bennet, C.O., Myers, J.E., 1982, “Momentum, Heat and Mass Transfer”, McGraw-Hill,
New York.
Brebbia, C.A., 1978, “The Boundary Element Method for Engineers”, Pentech Press,
London.
Özisik, M.N., 1968, “Boundary Value Problems of Heat Conduction”, Dover Publ., New
York.
Partridge P.W., Brebbia C.A., Wrobel L.C., 1992, “The Dual Reciprocity Boundary
Element Method”, Computational Mechanics Publications and Elsevier, London.
Ramachandran P.A.,1994, “Boundary Element Methods in Transport Phenomena”,
Computational Mechanics Publication and Elsevier Applied Science, London.
White F.M., 1986,“Fluid Mechanics”, McGraw-Hill Int., Singapore.
Zhao, C., Hobbs, B.E., Muhlhaus, H.B., 1997. Finite Element Modelling of Temperature
Gradient Driven Rock Alteration and Mineralization in Porous Rock Masses, Int. J.
Numer. Anal. Meth. Geomech., 21, 863-881.
Zhao, C., Hobbs, B.E., Hornby, P., Ord. A., Peng, S., 2006. Numerical Modelling of Fluids
Mixing, Heat Transfer and Non-equilibrium Redox Chemical Reactions in Fluidsaturated Porous Rocks, Int. Journal for Num. Methods in Engineering, 66, 10611076.
12
DIREITOS AUTORAIS
Os autores são os únicos responsáveis pelo conteúdo do material impresso incluído no
seu trabalho.