íifc ANAIS COBEM 81 VI CONGRESSO BRASILEIRO DE ENGENHARIA MECÂNICA RIO DE JANEIRO, IS -18 de dezembro de 1981 PROCEEOINGS - PDF
íifc ANAIS COBEM 81 VI CONGRESSO BRASILEIRO DE ENGENHARIA MECÂNICA RIO DE JANEIRO, IS -18 de dezembro de 1981 PROCEEOINGS
Download "íifc ANAIS COBEM 81 VI CONGRESSO BRASILEIRO DE ENGENHARIA MECÂNICA RIO DE JANEIRO, IS -18 de dezembro de 1981 PROCEEOINGS"
Gabriella Conceição Fraga
1 ANAIS COBEM 81 VI CONGRESSO BRASILEIRO DE ENGENHARIA MECÂNICA RIO DE JANEIRO, IS -18 de dezembro de 1981 PROCEEOINGS íifc TRABALHO PAPES N. D- P.P PUC/RJ APLICAÇÃO DE "ALTERNATING DIRECTION METHODS" A SOLUÇÃO DA EQUAÇÃO DE CONDUÇÃO DE CALOR CCM FONTE E EM REGIME TRANSITÓRIO Antonio Carlos de Oliveira Barroso Antonio Carlos Marques Alvim CÍMISSAO NACIONAL DE ENERGIA NUCLEAR Anibal N. Gebrin Dept» de Engenharia Nuclear - COPPE - UFRJ Rubens S. dos Santos Seção de Engenharia Nuclear - IMJ: SUMARIO Diversos tipos e variantes dos "alternating direction methods", (ADM) foram aplicados a solução da equação de condução de calor, com fon te, regime transitório, em regiões axissimétricas. Os resultados obtidos, mostram as vantagens de algumas variantes em relação ao método tradicional de Crank-Nicolson. 0 "alternating direction explicit", (ADK), foi en tre os AIM básicos, o que mais se destacou. Entre todos, o ADI: acoplado a una transformação exponencial parece ser a variante com maior potencial, principamente, quando usada com intervalo de tempo variável. SUMMARY Various types and variants of alternating direction methods, (ABM), were applied to the solution of the time-dependent heat conduction equation, with source, in regions with axial simmetry. Among the basic AlW's, the alternating direction explicit was lhe one which performed better. An exponential transformation coupled to the ADf: seems to be the variai* with gtcater potencial, especially if used with a variable time step scheme.
2 f Introdução li» trabalhos anteriores CQ'CO e DO» m étodo "altern.it ing d rection implicit" (ADT), proposto por Peaceman e Rachford \jq e modificado por Douglas [XJi foi aplicado a solução temporal de equação da con dução de calor em duas e três dimensões. Nesses trabalhos foi constatada uma clara vantagem desso método sobre o processo tradicional ilc Crnnk - Nicolson para a maioria das aplicações estudadas. Também em cinética de neutrons, os "alternating direction methods" (ADM) foram aplicados com sucesso ã solução da equação da difusão considerando vários grupos de energia QTJ, JJÊT] e J20* Nesta área contudo, para os problemas de real interesse o ADI mostrou-se apenas equivalente {jq a uma variante otimizada do método de Crank - Nicolson, implementada no programa TWIGL P ]. Tal comportamento atribuído à "stiffness" das equações de cinética, pode ser em parte compensado com a transformação exponencial [TQ, DO* Desta forma, o ADI aparecia como uma boa opção, contudo foi um outro tipo de AÜM, proposto por Larkin \j[j c denominado "alternating direction explicit" (ADE), que quando aplicado ãs equações de cinética QTJ, rjftj, mostrou ser a melhor opção. Nesse trabalho, as alternativas apresentadas na literatura consultada, bem como algumas variantes investigadas pelos autores, for.uii testa das para uma classe reduzida de problemas. Foram tratados problemas bidi_ mensionais com simetria axial e geometria semelhante a de uma vareta com bustível de um reator PWR. Os testes realizados foram todos muito simples, pois o objetivo do trabalho nesta fase é apenas selecionar as variantes mais promissoras que seriam então implementadas num código para análise da distribuição de temperaturas numa vareta de combustível duran te um acidente de perda de refrigerante primário. 2. A Equação da Difusão de Calor Espacialmente Discretizada Em regiões heterogêneas, mas que apresentam isotropia c simetria axial, a condução de calor ó descrita pela seguinte equação: onde p, c, k e q"' são respectivamente massa específica, calor específico, condutividade térmica e densidade volumétrica da fonte de calor. Usando-so o processo de "box integration" QO], obtém-se uma discretização espacial de diferenças finitas do tipo "corner mesh", com a (D
3 r seguinte forma: 2B7 a d,p onde os coeficientes a's e 0's'dependem dos parâmetros k, p, c c do espa çumente da malha eni volta do ioiito (i,j). A equação (2), juntamente com as formas discretizadas das condições de contorno podem ser escritas na seguinte forma compacta: (3) onde: T é um vetor coluna cujos componentes são as temperaturas nodais. A é a matriz que traduz o acoplamcnto entre os pontos nodais, de acordo com a discretização efetuada. S é o vetor com os termos de fonte, seus componentes são do tipo q'."./f Integração Temporal da Equação Scmi-discrctizada 0 método de Crank - Nicolson generalizado, aproxima a solução da e quação (3) do seguinte modo: + os(n + l) + 1_ Q) s (n) (1) ou seja. T ín)»0s tntl) *ll-oj s ín. onde j[ 5 a nntriz identidade, 0 é um fator peso, (0-0- I), entre as par to» implícitas o explícitas i* os sobio-csfritos usados indicam o inst;uite de tempo no qual a grandeza está sendo avaliada. Com a escolha de 8 = 0,5, a equação (5) recairia no método original de Crank - Nicolson. A idéia básica do método acima ao avançar um intervalo de tempo (át), é fazê-ío de forma parcialmente úiiplícita e parcialmente explícita, concomitantemente em todos os pontos nodais. Desta forma, as equações alftêbricas resultantes têm que ser resolvidas simultaneamente. Dentro da motivação de se obter sistemas uigébricos mais simples, surgem os métodos ADM. Nesses métodos, para problemas bidimensionais, por
4 288 exemplo, o intervalo de tcjnjx> (Atyéfracionado em duas partes, na primeira metade o avanço ê feito de forma implícita para um conjunto de pontos e explícita para os -outros pontos nodais. Na segunda parte 6 feita uma inversão quanto aos»ntos que foram calculados explicita c implicitamente. Considere as matrizes, (6) Lançando mão desta notação, os métodos ADM, para duas dimensões, podem ser descritos de uma forma geral. Na primeira parte do intervalo At, a equação (3) seria assim integrada, T(n+l/2)_ T (n) =t _ At A T (n+l/2) +A ^ )T (n i + * S (7) na segunda parte, a integração seria feita da seguinte forma: + f A,) 1 [(I- ^ A 3 )T (8) e combinando-se as equações (7) e (8) obtem-se: 3.1. Os Métodos ADM Básicos Considere as matrizes L. e L^, tais que L, contenha os elementos de A, que correspondem a discretizaçao do operador diferencial da equação (1) com relação ã direção re Lj aqueles correspondentes ã direção z. Note que l^ + Ir * A, O ADI de Peaceman e Rachford t Neste caso, as seguintes partições são usadas: 42 "
5 O ADI de Do»glas ÇQ Neste caso toma-se: 2 89 (6) 0 ADE Proposto por Larkin \jfj Considere as matrizes E, D e C, tais qiie D contenha os termos da diagonal principal de A, e E e C, respectivaiicnte, os termos abaixo e a- cima desta diagonal. Neste caso A é fracionada da seguinte forma: (7) 0 "Alternating Direction Checkerboard" (ADC), flfj Considere as matrizes: Io I* ' Io k X Jo 4 - lx V (8) escolhendo-se, matriz diagonal < 1, 0,1,0... A In m matriz diagonal <0,l,0,l,..A, f o método ADC será descrito pela seguinte partição: 0 Método de Linhas Alternadas,'line Hopscotch," (ADL), Considere as matrizes: Neste caso, I, e I 2 serão Matrizes bloco-diagonais cujos blocos são aitemadamente matrizes identidade e matrizes nula. Pcsta forma, este método c descrito pela seguinte partição: âl ' à 3 ' li à & - ÒA " Í2 è
6 290 Vale a pena notar, cjiie coin esta partição, ao se avançar At/2 são resolvidos 2 sistemas de equações separados, um totalmente explícito e outro equivalente a um ADI. Tais sistemas têm como incógnitas as tempera turas cm linhas alternadas da malha de discretizaçao Consistência, Estabilidade e Erro de Truncamento As propriedades dos métodos acima quanto ã estabilidade a consistência, não serão aqui discutidas, pois já foram claramente aborüaci-ts em trabalhos anteriores. A referência ptj por exemplo, apresenta uma prova da estabilidade e consistência do ADI e ADE, valida também paru o caso de regiões heterogêneas. Quanto ao erro de truncamento, apresentaremos uma derivarão da for ma de seu termo principal, posto que esta sorvirii de base ã proposição de algumas variantes dos ADM básicos. Em uma forma mais compacta, a equação (1) pode ser escrita como. - MT + q' (10) onde M 5 um operador diferencial linear que pode ser definido por comparação entre as equações (1) e (10). A equação (3), obtida de uma discretização espacial, é uma aproximação da equação (10), onde ao se aproximar o operador M pela matriz A, 2 2 introduziu-se um erro de ordem 0 (Ar ) + 0 (Az ). Por outro lado, a solução formal da equação (3), ao se avançar um intervalo de tempo At, é dada por: Note-se contudo, que em relação 5 solução exata das equações (1) ou (10), os valores dados pela equação (11) estão afetados de um erro de triuicaiiiento de 3- ordem, ou seja 0 (At Ar 2 ) + 0 (At Az ). Por meio de algumas expansões em série, a equação (11) pode ser pos_ ta na seguinte forma: (12) H-At 2 A+2(At 3 )] sh AtJ-O(At 3 fj s ín+l >
7 291 Procedendo-se de forma análoga e considerando-se apenas partições simétricas (âj * 3 e A3 «A^), a equação ( 9 ) pode ser reescrita. a tf 1 (13) Comparãndo-se as expressões (12) e (13) e lembrando o erro jâ implícito na equação (11), obtem-se uma expressão para o erro de truncamen to, válida p;ira os AfiM simétricos. Definindo-se o erro de truncamento (fi.t.) como a diferença entre os valores aproximados e exatos, pode-se escrever:! : -. íí E.T. =- At* Í (n+1) -S (nj l+o(at 3 )+0(ôtAr 2 )+0(AtAz 2 ) (14) 3.3. Algumas Variantes dos ADM Básicos Os ADM tós-corrigidos (AD1>C) A idéia básica é se adicionar um "termo corretor" de modo a cancelar o termo de ordem At na equação (14). Desta forma, após se computar T(n+1) u s a n c] o a equação (9) t adiciona-se a este o termo ^ A^L S^n+1 ' - " S J, para melhorar a ordem do erro de truncamento do método. Note que com exceção do ADI de Douglas, esse tipo de variante pode ser aplicado a qualquer um dos métodos básicos. A Transformação Exponencial Utilizando-se a seguinte mudança de variável, í IL" I T (t) = exp(ftt)8 (t), (15) ondo g ô uma nuitriz diagonal, a equação (3) tran?forma-se em: -è- 6 «expí-ílt) exp(-í)t) 0 (16) A motivação desta transformação ê atenuar possíveis variações de ordem exponencial na variável original, obtendo-so desta forma, uma variável transformada de comportamento temporal mais suave. I:spera-se portanto, que a equação (16) possa ser integrada com intervalos de tempo
8 maiores do que a equação (3). 292 Na pratica, para se usar esta transformação é necessário os ciumentos da matriz Q a cadu intervalo de tempo. Vários algoritmos foram tentados com essa finalidade fít], DO e os melhores resultados furam obu dos com a expressão abaixo: M - 1 (n) " Et In r i.j TV < t < t., n+1 (17) Métodos de Partição Variável A motivação aqui é também melhorar o erro de truncamcntu, purún ao invés de se adicionar iun termo corretor, tenta-se anular o termo principal através de uma partição adequada. Desta forma, os termos das diagonais de A, e A, seriam escolhidos, a cada intervalo de tempo, de modo que ) -s Cn) ]= o 0 Método "Locally One-üimensional" (LOD) Este método surgiu da classe de métodos de intervalos frac Lonados, propostos na literatura russa. Usando o mesmo formalism) empregado na descrição do ADI, tem-se: A 4 Como pode ser notado, este método não ê um ADM no sentido estrito, mas foi testado por haver haver boas indicações na literatura, quanto no seu potencial. 4. Resultados c Conclusões Foram utilizados seis problemas-teste, todos com a mesma geometria, mas com diferentes combinações de condições de contorno e fonte térmica. A respeito das variantes descritas na seção 3.3, observou-se que a adição do termo de correção, testada com o ADB e com o ADC, não acarretou mudanças sensíveis nos resultados. Das possíveis variantes com partição variável, foi testada apenas uma, com partições de configuração semelhan tesa do ADE. Como os resultados obtidos não foram consistentes, desistiu-sc dessas variantes.
9 293 Quanto ao LOD,. foi abandoaido por apresentar dcscnpenho muito inferior ao dos ADM básicos. Exceto para as variantes citadas acima, são apresentados resultados obtidos com os vários métodos testados, para um cilindro de 12,85 ft de comprimento e 0,21 in de diâmetro, con as seguintes propriedades: p lbm/ft 3, k» 1,54 Btu / hr.ft. F e c = 0.0S9 Btu/lh» F. As condições de contorno correspondem as de convecção para a água ã temperatura de 55.0 F e coeficiente de transferência de calor h = 5364 Btu/hr.ft F. As tabelas 1,2 e 3 a seguir resumem os resultados obtidos para tran sitórios cem fonte em rampa, fonte simulando calor de decaimento e fonte em degrau. O erro máximo E, referido nestas tabelas, foi definido cano: i G * max ij x 100 (I), onde Ti. 6 o valor da temperatura no ponto (i,j),dada por uma convergida. solução Tabela 1. Fonte em Rampa - Incremento de Tempo (At)para um Erro Máximo de 11 Método CN ADE AI)F ADC ADL AI)T a'" q * < l^1 ' Incremento de Tempo, At (seg) < 32 < 8 < A t < o c ' = 2.12xlO Btu/ft? hr rt,t>o Y =i.27xlo Btu/ft? hr 2 gli Pi ill
10 r 29U Tabela 2. Calor de Decaimento - Incremento de Tempo para um Erro Máximo de lt Método CN ADE ADF ADC ADL ADI yq»t Incremento de Tempo, At (seg) < 1, 0 < 0, 5 < o. 5 < 0, 5 < o. 5 < 0, 5 t < o q" 1 s 2, 65x10 7 Btu/ft? hr - 0,26 ^ JC 095 o, Tabela 3. Fonte em Degrau - Incremento de Tempo para um Erro Máximo de lt i'! Mé q'" todo CN ADE ADF ADC ADL ADI í "iff t t Incremento de Tempo, < 8 At (seg) < o >o ^'-2,65x10'Btu/ítr hr Observando os resultados obtidos pode-se verificar que, para a mes^ ma precisão, o CN admite, em geral, intervalos de tempo de 2 a 8 vezes maiores que os necessários aos demais métodos. O pior desempenho desse me todo foi para a perturbação cm degrau. Entretanto, o tempo de processamento é cerca de IS vezes maior, que, por exemplo, o ADE. A transformação exponencial parece ajudar, quando usada para modificar o ADE. 0 método resultante foi designado ADF. Este método mostrou ser apenas murginalmente superior ao ADE quanto à precisão, tendo porém, a desvantagem de um tempo de computação 2 vezes maior que o ADE. Deve-se observar, contudo, que algumas melhorias estão sendo introduzidas e que poderão mudar a conclusão acima. Por exemplo, as "freqüências" ft^ poderiam ser avuliadus pura o conqiortumento médio de um grupo de pontos c não para a totalidade dos mesmos. Outra Modificação seria a introdução de um
11 a mes vezes.seme essanodi- I i I ; ;trou ' u-cin, j fe-se! que ' >odc- i c não 1 e urn 1! r 295 controle de intervalo de tempo, mantendo o erro abaixo de determinado l'± mite. 0 ADF mostrou ser o método mais sensível a este controle, uma vez quo no início dos transientos, ao so usar un At maior do que o indicado para o ADE simples, "injeta-se" un erro que será realimentado pelo próprio calculo das freqüências. Isto foi verificado ao se testar o caso em degrau, onde se utilizou intervalos de tempo variáveis para o ADF, tendo este desempenhado bem melhor que o ADE. 0 ADC admitiu intervalos de tempo quase do mesmo tamanho que o ADE, sendo tão eficiente computacionalmente quanto este. 0 ADL é mais rápido que o ADI, sendo ambos equivalentes an precisão ao ADE. Entretanto, ficam aquém do ADE ou do ADC em rapidez de computação. Com base nos casos analisados, concluímos pela seguinte hierarquia de métodos, considerando precisão e tempo de processamento: ADE, ADF ou ADC e ADI. m ra DO L 6 1 Ü1 REFERENCIAS Snith, J.H., "Engineering Applications of ADI Methods to Piecewise Linear Multidimensional Heat Transfer", J.CoEip.P'iys. vol.17 (197:0, pp Snith, J.H., "Survey of Three Dimensional Finite Differences of Heat Equation", Sandia Laboratories Report SC-M-7083 (1970). Douglas, J.Jr., "Alternating Direction Methods for Three Space Variables", Numerische Mathematik. Vol. 4 (1962), pp Peaceman, D.W., Racltford, , "The Numerical Solution of Parabolic and Elliptic Differential Equations", J. Soc.Ind.Appl.Math., vol. 3 (1955), Read, W.H., Hansen, K.F., "Alternating Direction Methods for the Reactor Kinetics Equations", Nucl.Sci.Eng.. vol.41 (1970), pp Fergunson, D.R., Hansen K.F., "Solution of tfte Space - Dependent Reactor Kinetics Equations in Three Dimensions", Nucl. Sci. Eng., vol. 51 (1973), pp S. Hageman, L.A., Yasinsky, J.B., "Comparison of Alternating-Direction Time - Differencing Methods with other Implicit Methods for the Solution of the Neutron Group Diffusion liquations", Nucl. Set. Iftg. vol. 38 (1969), pp Yansinsky, J.B., Natelson, M., Hageaan. L.A., "TWIGL - A Program to Solve the Tw - Dimensional, Two-Group, Space-Time Neutron Diffusion liquations with Temperature Feedback", WAI'ü-lM-743, Bcttis Atomic Power Laboratory (1968). ;: ;} hi i w ill f 111! m
12 296 larkin, U.K., "Somo Stable lixplicit Difference Approximations In the Diffusion Equation", Math. Comp., vol. 18 (1964), pp ?. 00] Varga, R.S., Matrix Iterative Analysis, ch. 6, Prentice - Rill, Englewood Cliffs, N.J. (1962). [Iff Gourlay, A.R., "ibpscotch: A Fast Second - Order Partial Differential Equation Solver", J. Tnst. Maths. Applies., vol. 6 (J970J, pp Q2] Gourlay, A.R., McGuire, G.R., "General Hopscotch Algoritlm for the Numerical Solution of Partial Differential Equations", J. Inst. Maths. Applies., vol. 7 (1971), pp