Integraçao reduzida para problemas advectivos-difusivos escalares discretizados pela formulaçao Supg com operador de captura de descontinuidade

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 14,2, 145-166(1993) INTEGRACAO REDUZIDA PARA PROBLEMAS ADVECTIVOS-DIFUSIVOS ESCALARES DISCRETIZADOS PELA FORMULACÁO SUPG COM OPERADOR DE CAPTURA DE DESCONTINUIDADES C.M. DIAS e A.L.G.A. COUTINHO Progranza de El~gelzharia Civil - COPPE/UFRJ Cx. Postal 68506, Rio de Janeiro, RJ 21945-970, Brasil Tel.: + 55-21-280 9993 , Faz: + 55-21-280 9545 E-nzail: [email protected]. br Este traba1110 apresenta iiina técnica de controle das oscilacoes esgúrias que podein surgir quando as integrais proveiiieiites da discretiza~iio de problemas escalares advectivos-difusivos em regiine permanente pela foriiiulaciio SUPG com operador de captura de descontii~uidades, sao avaliadas utilizailclo-se uin ponto de quadratura de Gauss. Os termos de correcao siio obtidos & partir do priilcípio variacional de Hu-Wasliizu. Diversos exeinplos numéricos siio apresentados de forina a validar a técnica proposta. REDUCED INTEGRATION FOR STEADY ADVECTIVE-DIFFUSIVE PROBLEMS DISCRETIZED BY THE SUPG FORMULATION WITH DISCONTINUITY-CAPTURING SUMMARY This worlc presents a tecliiiique to control tlie spurious oscillations that inay appear when the integrals arisiilg froin the discretization of scalar steady-state advective-diffusive problems by the SUPG forinulatioii witli discontinuity-capturing are evaluated usiilg one-point Gaussiai~ quadrature. Tlie correction terms are clerived froin a Hu-Wasl-iizu variational forinulation. Severa1 numerical esamples a.re slzown, validating the proposed technique. Recibido: Septiembre 1996 OUniversitat Politecnica de Cataluriya (España) ISSN 0213-1315 146 C.M. DIAS E A.L.G.A. COUTINHO Na soluc5o de problemas advectivos-difusivos pelo Método dos Elementos Finitos surgein integrais que podem ser avaliadas utilizando-se uma regra de integras50 numérica, geralmente a quadratura de Gauss. O custo computacional da avaliaqáo destas integrais é proporcional ao número de pontos de integras20 utilizados. Assiin, a utilizac50 da integras50 com um ponto de quadratura vem reduziisignificativamente o custo destas avaliacóes. Atualmente muitos programas utilizam esta Integraq50 Reduzida com sucesso em problemas de mecanica dos sólidos, cascas e difusao, reduzindo o número de avaliacóes do operador semi-discreto e das equaqóes constitutivas. No entanto, tal procedin~ento pode provocar o aparecimento de oscilasóes espúrias conhecidas como modos hourglass, que podem comprometer os resultados. Desta forma, torna-se necessária a utilizasáo de uma técnica de controle destas oscilacóes. Existem diversas técnicas capazes de permitir a utilizas50 da Integras50 Reduzida com um ponto de quadratura em malhas de elementos quadriláteros e hexaedros lineares para estas classes de problemas. Dentre elas destaca-se a desenvolvida por Liu e Belytsclikol0 e Belytschko et al.l, onde os modos singulares espúrios que podem surgir com a Integra550 Reduzida s5o eliminados com a adic5o de uma parcela de estabilizaciio & matriz de elemento. Mallet et al. l1 estendeu tal metodologia de estabilizas50 & sistemas advectivos-difusivos. O presente trabalho voltado para a forinulasao SUPG (Brooks e Hughes2), acrescida de um operador de captura de descontinuidades (Codina4, Gale5o e Dutra do Carmo8), utiliza a Integrac50 Reduzida com um ponto de quadratura, estabilizada a partir do princípio variacional de Hu-Washizu, em problemas advectivos-difusivos escalares em regiine permanente, discretizados por elementos isoparamétricos quadriláteros bilineares. O artigo é organizado da seguinte forma. Primeiramente é feita uma rápida revisa0 das equacóes envolvidas e sua discretizas50 pelo Método dos Elementos Finitos para a forrnulaq50 SUPG acrescida do operador de captura de choque. É apresentada a forma matricial do problema. O item seguinte discute a metodologia de estabilizaciio empregada. Em seguida, apresentam-se alguns experimentos numéricos de validac50, e por último s5o feitos comentários e conclusóes sobre o trabalho. Sendo R um domínio em R2 com contorno dr = TI U F2, a equasao de advecc5odifusiio escalar em reginle permanente é dada por com condicóes de contorno +=$* emri (2) n.kV$=$, emT2 (3) onde f, +* e 4, sáo funsóes conhecidas, ,8 é a velocidade do fluido tal que V - B = 0, k é o tensor de difusividade, considerado muito pequena, n é o vetor unitário normal i fronteira e orientado para o exterior e q5 a incógnita do problema. Utilizando-se o método SUPG (Streamline Upwind Petrov-Galerltin) proposto por Brooks e Hughes2, acrescido de um operador de captura de descontinuidades (Galeáo e Dutra do Carmo8, Codina4), chega-se ii seguinte formulaciio de Elementos Finitos para a equacáo (1) nnos onde q5h = N;+; é a forma discreta de 4; q5;, i = 1,. . .iznos sáo as incógnitas iiodais; i N; siio as fuizc6es de interpolaciio usuais de elementos finitos e wh é a fuiiciio discreta de ponderacáo. Assume-se que a inalha possua nnos pontos nodais. As func6es de interpolacáo para elementos quadriláteros bilineares, restritas aos elementos sáo dadas por para o domínio de referencia, [E = - 1, E = 11 x [q = - 1, q = 11, coizforine mostrado na Figura 1. Na equacáo (4) o priineiro somatório corresponde ii parcela de estabilizaciio da formulacáo SUPG, eiiquaiito que o segundo somatório representa A parcela do terizio de captura de descontinuidades. A expressáo (5) correspoiide i parcela de Galerkin e o operador de resíduo discreto da equacáo de adveccáo-difusáo é dado pela equacáo (6). Plaiio físico (x, y) Plano de referencia (E, 11) Figura 1. Elemento quadrilátero bilinear - plano físico e plailo de referencia Na equacáo (4) nel é o número de elementos na malha, tais que izel R=URe e @=R;nR;, ifj e i,j=l, ...) 7zel (8) Considera-se que a contribuicáo descontínua 2i inontante das linhas de corrente é calculada pela expressáo assintótica (Brool<s e Hughes2) sendo Pe o número de Péclet do elemento, ou seja O comprimento característico do elemento he, pode ser calculado segundo Codina et aL3, isto é onde 11 .II indica a norma euclidiana e é a projecáo de P no plano de referencia. Para a parcela de captura de descontinuidades, as grandezas correspondentes sáo definidas como lzir P~II 7-2 = - min (?; O, 7) 2 onde o subíndice 11 indica a projecáo de P na direcáo paralela a Vq5. A expressáo (12) é uma aproximas50 assintótica obtida a partir da relacáo dada ein Codiiia4. Substituindose as aproximac6es de Elemento Finitos ein (4) A (6), chega-se a forma inatricial do problema expressa por K e F sáo obtidos de forma usual, a partir das contribuicGes dos elementos nel K = ~(KZ + K:, t K:, t KL,) e=l nel Considerando-se as fuiic6es de iiiterpolacáo bilineares dadas pela expressáo (7) e o eleinento quadrilátero da Figura 1, as parcelas de K presentes em (16) podem ser explicit adas como r Matriz de difusiio de Galerltin: [Kz] = /Re N6,jkGN,,idR, a,b = 1,. . .,4; i, j = 1,2 (l8) e Matriz de adveccáo de Galerltiii: [K:,] = / NbP VN,dR, a, b = 1,. . . ,4 Re (19) r Matriz de advecciío de Petrov-Galerkin: [K:,,] = / -Nb,jpi @@;N,,;, a,b = 1,. . .,4; i, j = 1,2 Re IlPll (20) r Operador de captura de descoiltiiluidades: onde A é a área do eleinento dada por 1 A = -(.2'311/42 + 224~31) 2 para Tais integrais podem ser avaliadas através da quadratura de Gauss com um ponto, conduzindo a [~;](~)=~b;l;.;~~b', ~j=1,2 ern [=?=O (24) [K:,](') = ~~~bj&p;bT, i = 1,2 ern [ = q = 0 (27) onde, segundo Codiila4 tein-se 150 C.M. DIAS E A.L.G.A. COUTINHO Observa-se que para o operador de captura do tipo CAU (Galeáo e Dutra do Carino8), P = 1. Nas expressóes aciina, bl e b2 sáo as parcelas de VN = [bl b2], dadas por e o vetor t represeiita un1 inoviineiito de corpo rígido As inatrizes iiitegradas iio poiito (E = 0, 11 = O ) so deficieiites cle posto, o que pode ocasioiiar oscilacóes espúrias (Hayesg). Desta forma, se faz iiecessária a utilizacáo de urna inetodologia capaz de aiiular os efeitos llegativos da deficieiicia de posto das inatrizes. Neste traba1110 os termos de correcáo das oscilacóes espúrias seráo deseiivolvidos coin base iios estudos de Liu e Belytsclilto10 e Belytscliko et al. para a equacáo de difusáo, e os trabalhos de Mallet et al. l1 e Dias6 para problemas predominaiiteineiite coiivectivos. Tais termos sáo obtidos a partir de uin priiicípio variacional do tipo Hu-Ilasliizu, como será visto a seguir. A formulacáo fraca e111 cada eleineiito, pelo fuiicioilal de Hu-Wasliizu é onde S iiidica a variacáo e k é o tensor de difusividade soinado i difusáo artificial oriuiida da forinulacáo SUPG coin operador de captura de descoiltiiluidades. A iiiterpolacáo de 4 pocle ser escrita coi110 45 = 450 + 4512 + 452~ + 4llh oiide a fuiicáo 1% é dada por Pode-se provar que sáo válidas as seguiiites propriedades de ortogoiializacáo sendo o vetor h, ideiitificado como o moclo lzozirglciss, dado por Utilizarido tais proprieclades, juiitaiileiite coin os valores iiodais de 4, toinados a partir de (35) cliega-se a espressáo oiide o superíiidice e indica que os valores siio tomados ein relacáo ao eleineiito e e o vetor y é calculado seguiido Belytscliko et al.', de forina que o posto do gradieiite discreto VN, seja 3. Defiiiiiido-se as aproxiinacoes onde onde g e q si50 as parcelas cle estabilizacáo do gradiente e do fluxo respectivaineiite. 152 C.M. DIAS E A.L.G.A. COUTINHO Substituindo-se (41) ein (34), e tainbéin as aproximacóes (44) e (45), vein S,. 64eT(cT + bjxj + 7h)(pi(gj1) + gh,i) - f)dfl6qFT(bi + -yhji)(p{') + óh,i)dfl+ S,. (6ql1lT + 6qT1x,;)(y{i) + yh,i - bTg5e - 7Tq5eh,i)dQ+ (1)T (1) + Le(6gi + 6~~h,i)(q, + rjh,; + kijgjl) t kijgh,j)dQ = 0, i = 1,2 A partir da primeira integral da equacáo (49) clzega-se A expressao do fl.uxo advectivo em cada elemento, e a partir da segunda, do fluxo difusivo observa-se que as duas primeiras parcelas correspondein ao termo de adveccáo de Galerkiii e as demais aos terinos de difusáo de Galerlcin e difusáo artificial oriunda dos termos do SUPG e do operador de captura de descontinuidades. Estas quando desmembradas, toriiain possível a identificacao das inatrizes de correcáo da deficiencia de posto provocada pela integracáo reduzida. Desta forma, tem-se Matriz de correcáo da difusáo: Matriz de correcáo da adveccáo: Matriz de correcao da adveccáo de Petrov-Galerkin: Matriz de correcáo para o operador de captura de descontinuidades: Assim, as parcelas da izzatriz K en1 (16) ser50 dadas pela soma dos termos provenientes da integragáo de Gauss con1 uin ponto, equacóes (24) a (27), e as matrizes de estabilizacao, equacóes (51) a (54)) isto é A foriiia esplícita destas iiiatrizes eiicoiitra-se 1x0 Apeiidice. e É iinportaiite notar que para o cálculo do resíduo iio interior do eleineiito, IL(C$~)~, toma-se os termos avaliados iio poiito de coorcleiiadas ([ = O, q = O). Logo, este terino tambéiii deverá ser estabilizado, utilizaiido-se uina correcáo seiiielliaiite a empregada iios casos aiiteriores, oiide e Nas equacóes (51) a (54), as integrais coiistituem os cliainados parametros de estabilizagiio Él = E zih,;clR, i, j = 1,2 (pa,ra aclvecciio de Galerliiii) S,. (61) & = E S, l~,;h,~dR, i = 1,2 (para os deinais casos) (62) e O coeficieiite livre E que aparece iias equacóes (61) e (62) iios dá a possibilidade de ajustar a precisáo dos resultados. Assiin, quaiido E = 1 verifica-se que [Ke] será exata para inalhas de eleineiitos retaiigulares e paralelograinos. Coin o objetivo de salieiitar o papel da escollia do coeficieiite E iia qualidade das solucóes obtidas, Liii e Belyts~lilio~~ apreseiitarain uina tabela coinparaiido as forinulag6es de Eleineiitos Fiiiitos para diversos valores de E coin as correspoiideiites moléculas de Difereiigas Fiiiitas para a equacáo de difusáo, coiisideraiido uina inallia de 4 eleineiitos quadrados de lado An: e tomaiido-se os resultados eiii toriio do izó central. Aquí, apreselita-se iia Tabela 1 resultaclos seinelliaiite aos citados aiiteriorilieiite, desta vez tomaiido-se as moléculas correspoiicleiites a iiioiitageiii global en1 toriio clo lió central da mesilla inallia de 4 eleineiitos citada aciina para as parcelas da advecg5o obtidas coin a forinulacáo aclui apreseiitada, coin E = O e E = 1, para aiigulos de escoaiiieiito de 45" e O". Verifica-se iia Tabela 1 que as moléculas siio identicas as apreseiitadas por Swamiiiatliaii e Voller13 para o iiiétoclo SUPG coi11 integra550 2 x 2. Observa-se iia Tabela 11 que as solucóes obtidas coi11 E = O, isto 6, sein estabilizacáo apreseiltain oscilacGes espiírias. Já as solucóes obtidas com a inetodologia de estabilizacáo, consideraiido-se E = 1, siio similares aos resultados obtidos com a iiitegracáo com 2 x 2 poiitos de quadratura. Ein seguida, o mesino exemplo foi solucionado para os casos aiiteriores, porém com a iilallia iiáo-uiliforine da Figura 3. Os resiiltados correspoiideiites eiicoiitram-se iia Tabela 111. Estes resultaclos sáo bastaiite aceitáveis cluaiido coizzparados coln os obtidos coi11 2 x 2 pontos de quadratura. 160 C.M. DIAS E A.L.G.A. COUTINHO Figura 3. Probleilla de adveccao pura com riialha 11.50-uiiiforine Problema de adveccao-difusa0 com operador de captura de descoi~tinuidades Consideraiido-se o mesillo doiníiiio do exeinplo anterior, com le = lo-' e velocidade coiistaiite de inódulo uiiitário, para uin iingulo O = 26,5", o probleina foi solucioiiado com uma malha regular de 20 x 20 elementos através das forinulacóes estabilizadas com operador de captura de descoiitinuidades. Os resultados obtidos eiicoiitram-se iia Figura 4, oilde pocle-se visualizar a solucáo da Sormulacáo SUPG, coi11 o operador de captura de Codiiia4 e com o operador CAU (Galeáo e Dutra clo Carmo8). Em todas as soluc6es considerou-se E = 1. As solucóes aciina difereiiz eiitre si basicaineilte juilto A camada limite. Nota-se que a iiitroducáo do operador de captura diminui sigiiificativaineiite as oscilacóes próximas As descoiltiiiuidacles. Observa-se tainbém que o operador CAU 6 visivelineiite mais dihsivo do que o operador dado coin a equaciio (28). Ein seguida, o l~roblema foi solucionado para outros valores do parainetro de estabilizacáo E. Os resultados para E = O e E = 2 sáo plotados iia Figura 5, ambos para o operador CAU. Como pode ser visto os resultados con1 E = O apreseiitain oscilacóes espiírias, eiiquarito que os resultados con1 E = 2 sáo bastante satisfatórios. SUPG coin E = 1 SUPG + operador de Codiiia SUPG + operador CAU COlll E = 1 coin E = 1 Figura 4. Testes con1 operador de captura E = O Figura 5. Testes com E = O e E 164 C.M. DIAS E A.L.G.A. COUTINHO 1. T. Belytschlto, J.S.J. Oiig, W.I<. Liu e J.M. Keililecly, "I-Iourglass Coiltrol iii Liilear aild Non-Linear Probleius", Conzp. Meth. Appl. Mech. Eilg~zg., Vol. 43, pp. 251-276, (1984). 2. A.N. Broolís e T.J.R. Hiiglies, "Streamliile Up~vind Petrov-Galerlíiil Foriiiulatioil for Convection Doi~~iilatecl Flows with Particular Einphasis oii tlie Iilcoinpressible NavierStolíes Equatioils", Co~izp. Meth. Appl. Mech. Eizg~zg. , Vol. 32, pp. 199-259, (1982). 3. R. Codiila, E. Oñate e M. Cervera, "Tlie Intriilsic Time for the Streamliiie Upwiiid/PetrovGalerltiil Forinulatioii Usiilg Quadratic Eleineiits", Conzp. kleth. Appl. Mech. E~zgizg., Vol. 94, pp. 239-262, (1992). 4. R. Codiiia, "A Discoiltiiiuity-Capturiiig Crosswiiicl-Dissipatioil for the Finite Eleineilt Solutioil of the Coi~vectioii-Diffusioii Equatioii", Conzp. Metlz. Appl. Mech. Engng., Vol. 110, pp. 325-342, (1993). 5. A.L.G.A. Coutiiiho e R.B. Aiiioriin, "Non-Syinmetric CG-lilte Scheines for the Eleineilt-byElement Solutioil of Fiiiite Eleineilt Equatioils", Aiiais do V Encontro Nacioilal de Ciencias Térmicas, ABCM, SP, Bra.si1, pp. 309-311, (1994). 6. C.M. Dias, "Uin estudo sobre a utilizas50 da iiitegrac50 reduzida ein problemas aclvectivosdifusivos", Tese de MSc., Programa de Eilgeilharia. Civil, COPPE/UFRJ, (1995). 7. FLOTRAN-Beilcliiiia.rlí Manual, pp. 4.1-4.7, (1991). S. A.C. Galeao e E.G. Dutra do Ca.riiio, "A Consistent Approxiiiiate Upwiiid Petrov-Galerkin Metliod for Coilvection-Domiiiated Problems", Conzp. ~ldeth. Appl. Mech. Engng., Vol. 68, pp 83-95, (1988). 9. L.J. Hayes, "Practica1 Stability Test for Fiiiite Eleiiieilts witli Reclucecl Iiitegration", 111-1. J. Num. Meth. Engng., Vol. 17, pp. 1689-695, (1981). 10. W.K. Liu e T. Belyt~clilio, "Efficieilt Linear ancl Noiiliiiear Heat Coilduction with a Quadrilateral Element", I~zt. J. Nunz. Meth. E~zgng., Vol. 20, pp. 931-948, (1984). 11. M. Mallet, C. Poirier e F. Slialiib, "A New Fiiiite Eleiiieilt Forinulatioii for Coinputational Fluid Dyilainics: Developinent of aii Hourglass Coiitrol Operator for Multiclimeiisional Advective-Diffusive Systems", Comp. Meth. Appl. Mech. Azg~zg., Vol. 94, pp. 429-442, (1992). 12. Y. Saad e M.H. Schultz, "GMRES: A Geileralized Miiiimal Residual Algoritlim for Solviiig Non-Syminetric Linear Systeilis", ,S'IAkf J. Sci. ,Ctat. Conzp., Vol. 7, pp. 856-869, (1986). 13. C.R. Swamiilatliaii e V.R. Voller, "Streainline Upwind Scliei~ie for Control-Voluine Fiilite Elemeilts, Part 1 - Foril-~ulations", Nunz. Heat Transfer, Part B, Vol. 22, pp. 95-107, (1992). 14. C.R. Swainiiiatl-iail e V.R. Voller, "Streamline Upwiilcl Scheine for Control-Volume Finite Eleineilts, Part 11: Impleineiitation aild Coi~~parison with the SUPG Fiiiite Eleineilt Sclieil~e", NPLI)~. Beat Trnnsfer, Part B, IJol. 22, pp. 109-124, (1992). As inatrizes de eleiiieiito para a equacao de advecsao-difusa0 escalar consideranclose uin elemeilto quadrilátero bilinear, integradas con1 uin ponto e estabilizadas coin a metodologia do item 3, eiicoiitrain-se A seguir. De forma geral, a primeira matriz representa o resultaclo da iiltegracáo coin uin poiito, eilquanto que a segunda iiiatriz é o termo de estabilizacáo. e AdveccZo (Galeikin): ~4 = -2A - ct1y13 - ct2x31 Difusiio (Galerltin): (A.13) (A. 14) (A.15) 166 C.M. DIAS E A.L.G.A. COUTINHO clcl clc2 clc3 c1c4 1 C2C1 ~2~2 ~2~3 ~2~4 C3C1 C3Cz C3C3 C3C4 cqc1 C4C2 C4C3 C4C4 onde 1 12 = p,?/;1 + ptJx;3 + 2PzPyi~31.2'13 Operador de captura de descontinuidades: onde (A.22) (A.23) (A. 24) (A.25) (A.26) (A.27) (A.28)