124
SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE TRANSFERÊNCIA DE CALOR E MASSA Ana Cláudia Magalhães Pimentel Tese de Doutorado apresentada ao Programa de Pós- graduação em Engenharia Mecânica, COPPE, da Universidade Federal do Rio de Janeiro, como parte dos requisitos necessários à obtenção do título de Doutor em Engenharia Mecânica. Orientador(es): Helcio Rangel Barreto Orlande Marcelo José Colaço Rio de Janeiro Abril de 2014

SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

  • Upload
    others

  • View
    3

  • Download
    0

Embed Size (px)

Citation preview

Page 1: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE TRANSFERÊNCIA DE

CALOR E MASSA

Ana Cláudia Magalhães Pimentel

Tese de Doutorado apresentada ao Programa de Pós-

graduação em Engenharia Mecânica, COPPE, da

Universidade Federal do Rio de Janeiro, como parte

dos requisitos necessários à obtenção do título de

Doutor em Engenharia Mecânica.

Orientador(es): Helcio Rangel Barreto Orlande

Marcelo José Colaço

Rio de Janeiro

Abril de 2014

Page 2: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE TRANSFERÊNCIA DE

CALOR E MASSA

Ana Cláudia Magalhães Pimentel

TESE SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO LUIZ

COIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA (COPPE) DA

UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOS

REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE DOUTOR EM

CIÊNCIAS EM ENGENHARIA MECÂNICA.

Examinada por:

____________________________________

Prof. Helcio Rangel Barreto Orlande, PhD

____________________________________

Prof. Marcelo José Colaço, D.Sc.

____________________________________

Prof. Fernando Alves Rochinha, D.Sc.

____________________________________

Prof. Nilson Costa Roberty, D.Sc.

____________________________________

Prof. Carlos Frederico Trotta Matt, D.Sc.

RIO DE JANEIRO, RJ - BRASIL

ABRIL DE 2014

Page 3: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

iii

Pimentel, Ana Cláudia Magalhães

Simulações sob incertezas em problemas de transferência de

calor e massa/Ana Cláudia Magalhães Pimentel. – Rio de

Janeiro: UFRJ/COPPE, 2014.

IX, 115 p.: il.; 29,7cm.

Orientadores: Helcio Rangel Barreto Orlande

Marcelo José Colaço

Tese (doutorado) – UFRJ/COPPE/Programa de Engenharia

Mecânica, 2014.

Referências Bibliográficas: p. 107-113.

1. Polinômio caos. 2. Método de colocação estocástica. 3.

Função de base radial. I. Orlande, Helcio Rangel Barreto et al.

II. Universidade Federal do Rio de Janeiro, COPPE, Programa

de Engenharia Mecânica. III. Título.

Page 4: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

iv

Aos meus pais, marido e afilhado.

Page 5: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

v

Agradecimentos

A Deus por ter chegado aqui!

Aos meus queridos orientadores, Helcio Rangel Barreto Orlande e Marcelo José

Colaço por me orientarem em todos os passos desta tese, por disponibilizarem seu

tempo quando precisei. Agradeço pela atenção, paciência e amizade!

Aos meus pais e irmã por me incentivar a ampliar os meus estudos.

Ao meu marido, Thiago Pimentel, pelo companheirismo e incentivo.

Aos amigos que conquistei na COPPE/UFRJ, em especial aos amigos

Wellington Betencurte da Silva, Camila Lacerda, Karolina Lopes, Marcus Costa e

Gabriel Guerra (por sua ajuda).

Aos professores e funcionários do Programa de Pós-Graduação da Engenharia

Mecânica da COPPE, em especial aos funcionários do LTTC e LMT.

Ao CNPQ pela ajuda financeira, através da bolsa de doutorado.

As minhas amigas Aurelina Maria Medeiros e Maria das Graças Cerqueira por

me apoiar nesta jornada.

A todos que de alguma forma me ajudaram para a conclusão deste trabalho.

Page 6: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

vi

Resumo da Tese apresentada à COPPE/UFRJ como parte dos requisitos necessários

para a obtenção do grau de Doutor em Ciências (D.Sc.)

SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE TRANSFERÊNCIA DE

CALOR E MASSA

Ana Cláudia Magalhães Pimentel

Abril/2014

Orientador: Helcio Rangel Barreto Orlande

Marcelo José Colaço

Programa: Engenharia Mecânica

Esta tese tem como objetivo aplicar técnicas de soluções estocásticas a

problemas de transferência de calor e massa. Em particular foram estudadas aplicações

do método de polinômio caos com projeção de Galerkin e com colocação através de

funções de interpolação de Lagrange. Inicialmente estes métodos foram aplicados a dois

casos particulares, envolvendo uma equação diferencial ordinária e, em seguida, uma

equação diferencial parcial representativa de um fenômeno de condução de calor

unidimensional transiente. Depois, um novo problema que trata da liberação controlada

de fertilizantes em uma esfera, contendo incertezas em alguns parâmetros, cuja solução

sob incertezas não foi encontrada na literatura, foi apresentado. A equação do balanço

de massa foi resolvida e foram estudados casos com um e dois parâmetros com

incertezas. Ao final, um problema de Magnetohidrodinâmica foi abordado, também

contendo incertezas em alguns parâmetros. Este problema acopla a abordagem do

método de colocação estocástica já referenciado à casos anteriores junto ao método de

colocação baseado na função de base radial.

Page 7: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

vii

Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Doctor of Science (D.Sc.)

SIMULATIONS UNDER UNCERTAINTY IN PROBLEMS OF HEAT AND MASS

TRANSFER

Ana Cláudia Magalhães Pimentel

April/2014

Advisors: Helcio Rangel Barreto Orlande

Marcelo José Colaço

Department: Mechanical Engineering

This thesis aims to apply techniques of stochastic solutions to problems of heat

and mass transfer. In particular, applications of the polynomial chaos method with

Galerkin projection and collocation through functions of Lagrange interpolation were

studied. Initially these methods were applied to two particular cases involving an

ordinary differential equation, and then a partial differential equation representing a

phenomenon one-dimensional transient heat conduction. Then a new problem that

comes from controlled release of fertilizer in a sphere, containing uncertainty in some

parameters, whose solution under uncertainty was not found in the literature, was

presented. The mass balance equation was solved and cases with one and two

parameters with uncertainties were studied. At the end, a Magnetohydrodynamics

problem was approached, also containing uncertainty in some parameters. This problem

couples the stochastic collocation method already applied to the earlier cases with the

collocation method based on the radial basis function.

Page 8: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

viii

Sumário

1. Introdução 1

2. Revisão Bibliográfica 4

3. Fundamentação Teórica 19

3.1 Polinômios Caos 19

3.1.1 Polinômios Caos Generalizado................................................ 20

3.2 Aplicações de Polinômio Caos a equações diferenciais 22

3.3 Método de Colocação Estocástica 24

3.3.1 Colocação via Interpolação de Lagrange................................. 25

3.3.2 Cálculo dos pontos da quadratura de Gauss e seus

respectivos pesos......................................................................

26

3.4 Método de Colocação Estocástica para variáveis estocásticas

multidimensionais

29

3.4.1 Método de Colocação com malhas esparsas (Algoritmo de

Smolyak)..................................................................................

30

3.4.2 Seleção dos pontos de colocação............................................ 31

3.5 Função de Base Radial (RBF) 34

3.5.1 Parâmetro de forma (c) e distribuição de centros no domínio 35

3.5.2 Expansão com RBF.................................................................. 36

3.5.3 Aplicação do Método de colocação por RBF à equações

diferenciais............................................................................... 36

4 Problemas Abordados 38

4.1 Problema 1: Modelagem da liberação controlada de fertilizantes 39

4.1.1 Método de Colocação aplicado à equação de balanço de

massa (Dimensão estocástica 1N )....................................... 41

Page 9: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

ix

4.1.2 Método de Colocação aplicado à equação de balanço de

massa (Dimensão estocástica 2N )...................................... 43

4.2 Problema 2: Problema de Magnetohidrodinâmica 46

4.2.1 Método de Colocação estocástica aplicado ao problema de

Magnetohidrodinâmica (Dimensão estocástica 1N )............ 50

4.2.2 Método de Colocação estocástica aplicado ao problema de

Magnetohidrodinâmica (Dimensão estocástica 2N ).......... 53

4.2.3 Método quasi-Newton.............................................................. 56

5 Resultados 57

5.1 Modelagem da liberação controlada de fertilizantes 59

5.1.1 Variável estocástica Unidimensional....................................... 62

5.1.2 Variável estocástica Bidimensional......................................... 68

5.2 Problema de Magnetohidrodinâmica 76

5.2.1 Variável estocástica Unidimensional....................................... 77

5.2.1.1 Análise de convergência do método de Monte

Carlo.......................................................................... 78

5.2.1.2 Análise de convergência do método de Colocação

Estocástica.................................................................... 83

5.2.1.3 Perfil da velocidade e temperatura............................. 86

5.2.1.4 Linhas de corrente e isotermas.................................. 88

5.2.2 Variável estocástica Bidimensional......................................... 91

5.2.2.1 Análise de convergência do método de Monte

Carlo.......................................................................... 91

5.2.2.2 Análise de convergência do método de Colocação

Estocástica.................................................................... 96

5.2.2.3 Perfil da velocidade e temperatura............................. 100

5.2.2.4 Linhas de corrente e isotermas.................................. 101

6 Conclusões 105

Referências Bibliográficas 107

Apêndice 114

Page 10: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

1

Capítulo 1

Introdução

Em engenharia, a análise de modelos matemáticos determinísticos, sem

considerar incertezas nos dados de entrada, é muito comum. Entretanto, na prática, essas

situações ideais são raramente encontradas, deparando-se com a necessidade de

enfrentar incertezas nesses dados. Essas incertezas podem ser caracterizadas como

propriedades do material, condição inicial, variação nos dados experimentais e

dificuldade em descrever o sistema físico, podendo causar erros na geometria e

variações nas condições de contorno. As incertezas também podem ocorrer devido à

representação matemática do sistema físico como, por exemplo, erros em aproximações

matemáticas e nas discretizações.

Para calcular com precisão o desempenho do sistema é essencial incluir o efeito

dessas incertezas no modelo analisado e entender como elas se propagam e alteram a

solução final. A presença de incertezas pode ser modelada no sistema através da

reformulação das equações governantes como, por exemplo, equações diferenciais

parciais estocásticas. Para resolver essas equações, os métodos probabilísticos vêm

sendo desenvolvidos [1, 2].

Os métodos probabilísticos podem ser divididos em duas classes: estatísticos e

não-estatísticos. A primeira classe abrange simulações de Monte Carlo [3, 4],

amostragem estratificada, amostragem por “Latin Hypercube”, etc. Estes métodos

envolvem amostragem que, na maioria dos casos, são de fácil aplicação. No entanto, a

precisão depende do tamanho da amostra que pode tornar o método caro,

principalmente para problemas onde o caso determinístico já é complicado. Entre os

métodos não-estatísticos, o mais popular é o método da perturbação [5, 6], onde o

campo aleatório é expandido através de séries de Taylor em torno de sua média e

Page 11: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

2

truncado na ordem certa. É empregado normalmente em sistemas de equações de

segunda ordem, pois para além da segunda ordem, torna-se extremamente complicado.

Uma limitação desse método está na magnitude das incertezas, não podendo ser muito

grandes em relação ao seu valor médio (em torno de 10%). Outra abordagem é baseada

na manipulação dos operadores estocásticos nas equações governantes, que inclui a

expansão de Neumann [7, 8] e método de ponderação integral que também são restritos

a pequenas incertezas.

Outra metodologia não estatística desenvolvida recentemente é o método de

Galerkin baseado em polinômios caos generalizado [9-16], uma generalização dos

polinômios caos clássico, onde a solução estocástica é expressa como polinômios

ortogonais dos parâmetros aleatórios de entrada. Esses polinômios podem ser

escolhidos, de acordo com a distribuição adotada para as incertezas em parâmetros de

entrada, para se chegar a uma boa convergência. Dentro dessa metodologia,

fundamentada em polinômios caos, o método de colocação [17-23], também se destaca

por apresentar bom desempenho. Este método utiliza polinômios caos baseado na

interpolação de Lagrange entre os pontos de colocação, e esses pontos coincidem com

os pontos da quadratura de Gauss, que resulta em um conjunto de equações

desacopladas.

Em um dos problemas analisados, o problema de Magnetohidrodinâmica, será

utilizado também um método sem malha baseado na função de base radial para a

solução do problema determinístico. O parâmetro contendo incerteza será expandido

por polinômio caos por interpolação de Lagrange e a solução das equações será

expandida por funções de base radial. Aplicando-se a função de base radial nas

equações governantes e condições de contorno obtem-se um sistema algébrico não

linear que é resolvido através de um método quasi-Newton.

Neste trabalho, a classe não-estatística será aplicada através de métodos

baseados em polinômios ortogonais, mais precisamente, polinômios caos.

De início, na fase de qualificação foi analisada uma equação diferencial

ordinária estocástica para algumas distribuições de um parâmetro de entrada e um

problema de condução de calor unidimensional transiente em que algumas variáveis

foram consideradas estocásticas. Como continuidade, a metodologia será utilizada na

modelagem da liberação controlada de fertilizantes e também no problema de

magnetohidrodinâmica configurando-se como problemas inéditos. A presente tese de

doutorado está estruturada da seguinte maneira:

Page 12: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

3

No Capítulo 2 é realizada uma revisão bibliográfica, na qual métodos não

estatísticos foram empregados e trabalhos referentes a temas abordados nesta tese foram

utilizados. Também referente ao método sem malha baseado em função de base radial

serão apresentadas referências que darão o embasamento a sua aplicação.

O Capítulo 3 aborda os métodos baseados em polinômios ortogonais

(unidimensionais e multidimensionais), métodos gPC (polinômios caos generalizado),

método de colocação, método de colocação com malhas esparsas, métodos para resolver

o sistema gerado por esses polinômios e o método baseado em funções de base radial.

No Capítulo 4 descrevem-se os problemas propostos e suas formulações

matemáticas.

No Capítulo 5 são apresentados os resultados obtidos e a comparação com

aqueles da literatura e aqueles obtidos pelo método de Monte Carlo.

No Capítulo 6 apresentam-se conclusões e possíveis trabalhos de continuação

desta tese.

A contribuição inédita deste trabalho se configura na aplicação de variáveis

estocásticas ao problema de magnetohidrodinâmica, que caracteriza incertezas quanto

aos valores de certas propriedades termofísicas. Também a modelagem da liberação

controlada de fertilizantes é inovadora, visto que a abordagem estocástica deste

problema não foi encontrada na literatura.

Page 13: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

4

Capítulo 2

Revisão Bibliográfica

Uma revisão bibliográfica foi feita sobre os métodos e problemas abordados nesta

tese. De início, apresenta-se uma revisão de trabalhos que abordam os métodos baseados

em polinômios caos generalizados e suas aplicações, a fim de aprimorar o conhecimento de

tais métodos. Um dos métodos baseado em polinômios caos, a colocação estocástica por

interpolação de Lagrange que utiliza polinômios ortogonais, tem apresentado boas soluções

em muitos problemas. Além disso, em problemas com variáveis estocásticas

multidimensionais, o método de colocação estocástica com malhas esparsas ofereceu bons

resultados.

GHANEM e SPANOS [1] consideraram sistemas de engenharia com parâmetros

randômicos, os quais foram modelados como processos estocásticos de segunda ordem,

definidos por sua média e funções de covariância. Eles foram os primeiros a combinar os

polinômios caos de Wiener com o método de elementos finitos modelando incertezas em

aplicações de mecânica dos sólidos.

XIU e KARNIADAKIS [10] definiram um algoritmo formado pela expansão de

polinômios caos generalizado para a solução de equações diferenciais parciais elípticas

estocásticas em regime permanente sujeitas a entradas com incertezas (equação de Poisson

com condições de contorno, termo fonte e difusividade randômica). A projeção de Galerkin

foi aplicada a fim gerar um conjunto de equações determinísticas resolvidas por técnicas

Page 14: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

5

iterativas (método de Gauss-Seidel). Distribuições randômicas discretas e contínuas foram

consideradas e a convergência foi verificada no problema modelo e por simulações de

Monte Carlo. A decomposição Karhunen-Loève (uma maneira de representar um processo

aleatório, baseada na expansão espectral da função de correlação do referido processo e que

proporciona uma redução da dimensão do espaço aleatório) foi usada para reduzir a

dimensionalidade do espaço randômico. A escolha adequada de uma expansão de

polinômios caos de acordo com a entrada randômica gerou soluções com convergência

exponencial.

XIU e KARNIADAKIS [11] apresentaram um novo método para resolver equações

diferenciais estocásticas baseado na projeção de Galerkin e em extensões de polinômios

caos de Wiener. Os processos estocásticos foram representados por uma base da família

Askey de polinômios ortogonais que reduz a dimensionalidade do sistema e leva à

convergência exponencial do erro. Vários processos contínuos e discretos foram tratados.

Uma equação diferencial ordinária estocástica específica com diferentes tipos de entradas

aleatórias foi resolvida e as taxas de convergência para as expansões foram mostradas,

comparando os resultados numéricos com a solução exata correspondente.

XIU e KARNIADAKIS [12] apresentaram um algoritmo de polinômios caos

generalizado para a solução de condução de calor transiente sujeita a propriedades com

incertezas, ou seja, condutividade e capacidade térmicas randômicas. Tais propriedades e a

solução estocástica foram representadas espectralmente pelas funções polinomiais

ortogonais do esquema Askey, como uma generalização da idéia do polinômio caos

original de Wiener. O conjunto resultante de equações determinísticas foi posteriormente

discretizado pelo método do elemento espectral no espaço físico e integrado no tempo.

Exemplos numéricos foram apresentados, mostrando que, quando a expansão caos

apropriada é escolhida de acordo com a entrada aleatória, a solução do polinômios caos

generalizado converge exponencialmente para o problema modelo.

XIU e KARNIADAKIS [13] desenvolveram um algoritmo para modelar incertezas

nos dados de entrada e sua propagação em simulações de escoamento incompressível.

Exemplos numéricos foram apresentados para condições de contorno com incerteza. O

método também pode ser aplicado a modelos com incerteza no domínio da fronteira, como

por exemplo uma superfície rugosa, no coeficiente de transporte (viscosidade turbulenta e

Page 15: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

6

outros modelos de transporte) ou em forças de interação para os problemas acoplados. A

entrada estocástica foi representada espectralmente empregando o polinômio ortogonal do

esquema Askey como base para representar o espaço randômico. A projeção de Galerkin

foi aplicada na dimensão randômica para obter as equações resultantes. O sistema de

equações determinísticas foi então resolvido com métodos padrão a fim de se obter a

solução para cada modo randômico. A eficiência e a convergência foram estudadas por

comparação com soluções exatas, bem como soluções numéricas obtidas por simulações de

Monte Carlo. Foi mostrado que o método de polinômio caos generalizado consegue uma

boa aceleração em comparação com o método de Monte Carlo. A utilização de diferentes

tipos de polinômios ortogonais do esquema Askey também fornece uma maneira mais

eficiente de representar processos gerais não Gaussianos em comparação com as expansões

originais Wiener-Hermite.

WAN et al. [14] resolveram uma equação de advecção-difusão bidimensional com

velocidade de transporte randômica. Diferentes tipos de distribuições randômicas foram

consideradas, onde a taxa de convergência foi examinada tomando como base a solução

exata. Nas várias distribuições utilizadas, através da expansão por polinômios caos, foi

observado que os erros para a média e para a variância diminuíram à medida que o grau do

polinômio aumentou. Outro problema também foi tratado onde a velocidade de transporte

randômica foi considerada bidimensional.

WAN et al. [15] apresentaram o método ME-gPC (polinômios caos generalizado

multi-elemento). A principal idéia do ME-gPC é decompor o espaço de entradas

randômicas quando o erro relativo para a variância se torna maior que um valor limite

estipulado. Primeiro o espaço de entradas randômicas é decomposto em pequenos

elementos. Em cada elemento é gerada uma nova variável randômica e o gPC (Polinômios

caos generalizado) é aplicado novamente. O grau de perturbação em cada elemento é

reduzido proporcionalmente ao tamanho do elemento randômico, obtendo-se um polinômio

(gPC) de baixa ordem em cada elemento. Vários problemas foram resolvidos pelos

métodos gPC e ME-gPC e comparados com o método Monte Carlo, onde foi observado o

melhor desempenho para o ME-gPC. O ME-gPC com polinômio de grau 3 apresentou

melhor desempenho se comparado com o gPC de grau 10 quando aumentou-se a

decomposição do espaço em elementos disjuntos.

Page 16: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

7

GAUTSCHI [16] mostrou um conjunto de programas no Matlab (parte do seu livro

"Polinômios Ortogonais: Computação e Aproximação”, Universidade de Oxford, Oxford,

2004). O pacote contém rotinas para geração de polinômios ortogonais, bem como rotinas

com aplicações. O objetivo desse trabalho foi apresentar vários procedimentos para gerar os

coeficientes da relação de recorrência de três termos para polinômios ortogonais e mais

relações de recorrência gerais de polinômios ortogonais de Sobolev. Os métodos baseados

no momento e os métodos de discretização, bem como sua implementação em Matlab,

estão entre os principais temas discutidos.

LOEVEN et al. [17] mostraram que a eficiência da quantificação de incerteza para

múltiplos parâmetros é aumentada quando se segue uma abordagem de dois passos com

colocação caos. No primeiro passo, uma análise de sensibilidade é utilizada para identificar

o parâmetro mais importante do problema. No segundo, a solução estocástica da solução é

obtida para o parâmetro mais importante usando o método de colocação caos. A eficiência

do método de colocação caos foi comparada com o método de Monte Carlo, o método

polinômios caos e o método de colocação estocástica. A eficiência dos métodos foi

comparada pela convergência do erro com respeito ao esforço computacional. A abordagem

de 2 passos com colocação caos é demonstrada para um pistão linear com condições de

contorno instáveis.

GANAPATHYSUBRAMANIAN e ZABARAS [18] investigaram os efeitos de

incertezas nos dados de entrada em problemas de convecção natural. Os seguintes

problemas foram analisados: convecção natural com condição de contorno randômica,

convecção natural com topologia de contorno randômico (através da variação da ondulação

e rugosidade) e convecção em meios porosos heterogêneos. Tais problemas foram

resolvidos por meio de malhas adaptativas e malhas convencionais (para o método de

colocação), pelos polinômios caos e pelo método Monte de Carlo. Para altas dimensões foi

mostrado que malhas adaptativas têm uma redução de erro maior do que malhas

convencionais.

FOO et al. [19] investigaram problemas tridimensionais em mecânica dos sólidos,

com incertezas nas propriedades do material ou com cargas estocásticas. Através dos

métodos baseados em polinômios caos generalizado e colocação estocástica, tais problemas

foram reduzidos a um sistema de EDPs determinísticas de alta dimensão, formuladas e

Page 17: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

8

resolvidas numericamente usando o método de elementos finitos para a discretização

espacial. Foi demonstrado que o método gPC para sistemas de elasticidade linear fornece

resultados precisos e eficientes bem mais rápido se comparado ao tradicional método de

Monte-Carlo. O método de colocação estocástica foi aplicado a problemas lineares e não

lineares. Para o problema linear de elasticidade com carga de pressão estocástica, os

métodos de colocação estocástica e gPC mostraram-se comparáveis em termos de

eficiência.

WAN e KARNIADAKIS [20] estudaram um problema elíptico com coeficientes

randômicos não-Gaussianos (um modelo estocástico típico em física envolvendo meios

porosos). Em particular, re-analisaram os métodos com polinômios caos acoplados com a

expansão K-L (Karhunen-Loève) em que, dada uma função de correlação, a expansão K-L

é empregada para reduzir a dimensionalidade das entradas aleatórias. A idéia central é

baseada na observação de que o método do polinômio caos proporciona um modelo

aproximado de EDP estocástica, onde sua solução é utilizada no método de Monte Carlo.

No modelo preditor-corretor gPC, foi usado o método de Monte Carlo para aperfeiçoar as

estatísticas dadas pelo modelo preditor gPC, onde a solução do polinômio caos serviu como

uma variável de controle para redução da variância, acelerando eficientemente a

convergência do método de Monte Carlo.

BABUSKA et al. [21] aplicaram o método de colocação estocástica na solução de

equações diferenciais parciais elípticas com coeficientes aleatórios. O método consiste

numa aproximação de Galerkin junto à colocação com pontos de Gauss o que leva a

solução de problemas determinísticos desacoplados.

XIU e HESTHAVEN [22] aplicaram o método de colocação e o método de Monte

Carlo a uma equação elíptica com entrada aleatória. Para o método de colocação é

conhecido que a escolha de seus pontos interfere no custo computacional e para isso várias

opções de conjuntos de pontos de colocação foram investigadas. Para altas dimensões

estocásticas, verificou-se que a malha esparsa baseada no algoritmo Smolyak ofereceu alta

precisão e rápida convergência, com uma fraca dependência do número de dimensões

aleatórias. Experimentos numéricos demonstraram que a estimativa do erro pode ser

separada em duas partes: o erro de discretização espacial e o erro na interpolação espacial

Page 18: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

9

aleatória. Constatou-se também que para dimensões aleatórias próximas a 50 os métodos de

colocação estocásticos são mais eficientes do que o método de Monte Carlo.

PEDIRODA et al. [24] utilizaram uma técnica de modelagem estocástica baseada

no método de colocação caos para medir a média e o desvio padrão para incertezas nos

parâmetros. Para uma dada precisão, o método de colocação caos requer menos amostras a

serem avaliadas em relação ao MC (Monte Carlo). A boa avaliação da média e do desvio

padrão pelo método de colocação caos torna-o bastante atrativo para ser usado com

métodos RDO (Robust Design Optimization). A RDO de um projeto de motores

automotivos foi realizada empregando o método de colocação caos. A estratégia de solução

foi implementada na integração de processo comercial e otimização de projeto com o

software modeFRONTIER.

XIU [25] apresentou alguns tipos de métodos estocásticos, mostrando comparações

entre eles. A metodologia do polinômio caos generalizado foi analisada. Ao utilizar os

polinômios caos, estratégias para estimar os coeficientes da expansão são necessárias, e

duas alternativas foram discutidas: o método de Galerkin e o método de colocação. Para

ilustrar os métodos, exemplos foram resolvidos e comparados com simulações Monte

Carlo.

ELDRED e BURKARDT [26] compararam o desempenho de uma abordagem não-

intrusiva do método baseado em polinômios caos e do método por colocação estocástica

com um problema clássico de soluções conhecidas e analisaram as diferenças encontradas

nas soluções a partir da utilização desses métodos. A principal diferença encontrada é que

os polinômios caos definem uma formulação de expansão e uma estimação do coeficiente

correspondente e a colocação estocástica requer somente uma definição dos pontos de

colocação a partir dos quais a expansão polinomial é derivada com base na interpolação de

Lagrange, seguida pela estimação dos coeficientes. Nos experimentos computacionais

analisados foram obtidos resultados bem próximos para ambos os métodos. Quando uma

diferença foi observada entre eles, a colocação estocástica apresentou uma precisão mais

elevada. Essa diferença pode ser em grande parte atribuída à expansão ou questões de

integração do polinômio caos, o que motivaram o estudo de abordagens para adaptação das

expansões por polinômios caos. Os resultados foram comparados com o método de Monte

Carlo. Para o caso da quadratura pelo produto tensorial, a utilização dos polinômios caos

Page 19: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

10

tem execuções idênticas à colocação estocástica, com suas diferenças completamente

eliminadas. Ambos os métodos superaram consistentemente o polinômios caos tradicional.

No entanto, as abordagens da quadratura pelo produto tensorial apenas superam abordagens

de malhas esparsas para problemas unidimensionais. Para problemas com mais de duas

dimensões, as abordagens de malhas esparsas mostraram-se superiores às abordagens da

quadratura pelo produto tensorial.

PREMPRANEERACH et al. [27] apresentaram diversos algoritmos para resolver

equações diferenciais ordinárias estocásticas, com aplicações em sistemas de energia

eletromecânica. O conceito de polinômios caos foi aplicado através da projeção de Galerkin

e do método de colocação se mostrando eficaz para problemas de baixas dimensões. Para

problemas de grandes dimensões ou na presença de descontinuidades, tanto a projeção de

Galerkin quanto o método de colocação perdem a convergência exponencial. Para isso, a

decomposição do espaço randômico proporcionou bons resultados, especialmente quando

foi combinada com a decomposição do domínio adaptativa.

PARUSSINI et al. [28] utilizaram o TeCC (Tensorial-expanded Chaos Collocation

Method) juntamente com o método de domínios fictícios para resolver problemas de

mecânica dos fluidos com incertezas geométricas. Este método de domínios fictícios é

baseado no LSqSEM (Least Squares Spectral Element Method), onde problemas

formulados em um domínio complexo podem ser resolvidos em um domínio fictício de

forma simples contendo o problema original. Desta forma, o domínio computacional, ou

seja, o domínio fictício, não é influenciado por pequenas variações nos contornos do

domínio original (sujeito a incertezas), que agora estão concentradas no domínio

computacional.

ONORATO et al. [29] compararam o método polinômios caos intrusivo (método

onde há modificação do código computacional) e o método colocação probabilística não-

intrusivo através de dois problemas contendo incertezas no número de Mach e no ângulo de

ataque. O método intrusivo apresentou uma característica importante: para cada incerteza

na formulação do modelo matemático uma nova dimensão foi introduzida na solução,

sendo essa incerteza considerada dependente da dimensão. A expansão convergente ao

longo destas novas dimensões foi buscada em termos de funções de base ortogonais, cujos

coeficientes foram usados para quantificar as incertezas da solução. Já no método não-

Page 20: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

11

intrusivo, a adição de uma nova variável com incerteza influenciou a solução através do

cálculo dos pontos de colocação, que foram obtidos pela quadratura de Gauss através do

produto tensorial.

WAN e KARNIADAKIS [30] utilizaram os métodos ME-gPC e gPC em equações

diferenciais parciais com coeficientes estocásticos (equação de advecção unidimensional

com velocidade de transporte randômica uniforme e um escoamento com “ruído” passando

por um cilindro circular estacionário, com duas condições de contorno diferentes na

entrada). Os resultados mostraram que o erro para o ME-gPC aumentou na mesma

velocidade que para o gPC. No entanto, o ME-gPC foi muito mais preciso. Também foi

observado que os erros para o ME-gPC são sempre delimitados pelo erro do gPC. O ME-

gPC pode melhorar o desempenho do gPC para problemas relacionados à freqüência

randômica, mas não tende a um valor assintótico. Para isso, uma opção é aumentar

adaptativamente o número de elementos do ME-gPC para manter uma precisão razoável em

um intervalo de tempo de integração desejável. No entanto, para entradas randômicas de

grande dimensão a eficácia do ME-gPC será enfraquecida desde que o número de

elementos aumente rapidamente.

MA e ZABARAS [31] analisaram o efeito de incertezas em dados de entrada na

solução de equações diferenciais ordinárias e parciais. Neste trabalho, o método de

colocação estocástica com malhas esparsas se mostrou uma boa alternativa para o método

dos elementos finitos estocásticos espectrais, por usar a interpolação polinomial de

Lagrange que acarreta em uma característica não intrusiva. Mas, apesar deste fato atrativo,

problemas que apresentem descontinuidades no espaço estocástico podem gerar soluções

não convergentes. Para isso, foi utilizada uma estratégia de colocação com malhas esparsas

adaptativas que detecta automaticamente a região de descontinuidade no espaço estocástico

e adaptativamente refina pontos de colocação nessa região.

GANDER e KARP [32] estudaram um problema padrão de radiação térmica, onde a

equação de transferência de calor por radiação foi resolvida utilizando a regra da quadratura

de Gauss. Para isso, foram apresentados dois algoritmos clássicos para calcular a regra da

quadratura de Gauss desejada, o algoritmo Stieltjes e o método usando momentos, que são

neste caso instáveis. Em seu lugar, para superar as instabilidades numéricas desses

Page 21: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

12

métodos, foi apresentado um método numericamente estável para a regra da quadratura de

alta ordem arbitrária.

BUNGARTZ e GRIEBEL [33] fizeram um levantamento dos fundamentos,

propriedades e aplicações de malhas esparsas, com foco na solução de equações

diferenciais parciais. Mostraram que o método de colocação com malhas esparsas é

bastante adequado para os problemas de alta dimensionalidade, podendo ser estendido a

malhas esparsas adaptativas para problemas mais complexos que, por exemplo, tenham

descontinuidades.

GERSTNER e GRIEBEL [34] apresentaram algoritmos novos e revisões sobre a

integração numérica de funções multivariadas usando vários tipos de malhas esparsas

aplicadas ao método de colocação introduzido pela primeira vez por Smolyak. Nesta

abordagem, as fórmulas de quadratura multidimensional foram construídas utilizando

combinações adequadas do produto tensorial de fórmulas unidimensionais. Foram exibidas

várias construções para as fórmulas de quadratura multivariada com malhas esparsas com

base em Newton- Cotes, Clenshaw -Curtis, Gauss e fórmulas estendidas de Gauss.

Resultados relativos aos custos de computação e ao erro indicaram uma implementação

numericamente estável mostrando a confiabilidade das fórmulas de malhas esparsas.

NOBILE et al. [35] analisaram equações diferenciais parciais com coeficientes

randômicos e termo fonte através do método de colocação estocástica com malhas esparsas

do tipo Smolyak. O método utilizou soluções aproximadas obtidas por elementos finitos,

correspondendo a um conjunto de equações determinísticas desacopladas. Em se tratando

de muitas variáveis aleatórias, o algoritmo de Smolyak utiliza abscissas de Clenshaw-Curtis

ou Gaussianas. Os resultados encontrados foram comparados com o método de colocação

com produto tensorial completo e com o método de Monte Carlo, mostrando ser bastante

eficaz.

XIU [36] combinou o método de colocação estocástica de mais alta ordem com a

expansão de polinômios caos generalizado, resultando em uma abordagem do tipo pseudo-

espectral. Os parâmetros de incerteza foram modelados como variáveis aleatórias e as

equações governantes foram tratadas como estocásticas. As soluções foram expressas como

séries convergentes de expansões polinomiais ortogonais em termos dos parâmetros

randômicos de entrada. Uma estimativa de erro foi apresentada, juntamente com exemplos

Page 22: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

13

numéricos para problemas com formas relativamente complicadas de equações

governantes. Várias opções de pontos de integração foram abordadas e a construção de

malhas esparsas mostrou ser uma escolha razoável para os problemas com grande número

de variáveis aleatórias.

BUNGARTZ e DIRNSTORFER [37] estudaram as malhas esparsas adaptativas

para variáveis aleatórias com mais de três dimensões. Os métodos convencionais

normalmente sofrem com o aumento da dimensão do problema ou apresentam precisão

insatisfatória. Para isso, uma estratégia de adaptação foi utilizada, a qual foi capaz de

detectar diferenças na importância das dimensões automaticamente. O algoritmo

apresentado foi aplicado a alguns problemas testes e comparados com outros métodos

existentes, revelando ser um método promissor em cenários que podem usar um esquema

de refinamento adaptativo.

ELDRED [38] investigou o desempenho relativo ao polinômio caos generalizado

não intrusivo e o método de colocação estocástica baseado na interpolação de Lagrange

aplicados a vários problemas com soluções conhecidas. A principal distinção entre esses

métodos é que o polinômio caos generalizado deve estimar os coeficientes para uma base

conhecida de polinômios ortogonais (por amostragem, regressão linear, produto tensorial,

ou malhas esparsas Smolyak) enquanto que a colocação estocástica deve formar uma

interpolação para os pontos conhecidos (usando quadratura ou malhas esparsas). O

desempenho desses métodos mostra-se muito semelhante e ambos demonstram boa

eficiência em relação aos métodos de amostragem de Monte Carlo, com boa precisão em

relação a outros métodos encontrados na literatura.

HOSDER et al. [39] abordaram a precisão e a eficiência computacional do ponto de

colocação no método não intrusivo de polinômio caos (NIPC) aplicado a vários problemas

estocásticos com múltiplas incertezas nas variáveis de entrada. Dois problemas estocásticos

com múltiplas variáveis aleatórias uniformes foram estudados para determinar o efeito de

diferentes métodos de amostragem (Randômico, “Latin Hipercube”, e Hammersley) para a

seleção dos pontos de colocação. O primeiro problema estocástico teve duas variáveis

aleatórias uniformes e o segundo problema incluiu quatro variáveis aleatórias uniformes. O

efeito do número de pontos de colocação na precisão das expansões polinomiais também

foi investigado. Os resultados dos problemas estocásticos mostraram que todos os três

Page 23: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

14

métodos de amostragem apresentam um desempenho semelhante em termos de exatidão e

eficiência computacional das expansões polinomiais.

A modelagem de liberação controlada de fertilizantes foi analisada através de

trabalhos encontrados na literatura e nenhum resultado com estudos estocásticos de

abordagem não intrusiva do problema foi encontrado.

AL-ZAHRANI [40] desenvolveu modelos matemáticos que podem descrever a taxa

de liberação de fertilizantes através de uma membrana polimérica esférica em diferentes

condições operacionais. A espessura da membrana polimérica desempenha um papel

importante para retardar a transferência de fertilizantes, controlando a transferência destes

no solo. Soluções detalhadas e aproximadas foram desenvolvidas e comparadas com a

solução numérica para a taxa de liberação de fertilizantes, sendo capaz de determinar o

perfil de concentração de fertilizantes ao longo da membrana polimérica a qualquer

momento. Com estas soluções é possível quantificar o papel desempenhado pelos diversos

parâmetros que caracterizam o fertilizante e a membrana polimérica, tendo grande

importância na seleção do polímero necessário para preparar tal membrana.

SHAVIV et al. [41] construíram um modelo estatístico para a liberação de uma

pluralidade de grânulos CRFs (fertilizantes de liberação controlada) com base em um

mecanismo de difusão. Uma análise de sensibilidade para os principais fatores que afetam a

liberação de uma população de grânulos (influência da permeabilidade do soluto e da água)

e características estatísticas como valor médio, desvio padrão e tipo de distribuição foram

feitas. As variações do raio e da espessura do revestimento afetam a liberação de

fertilizantes quando apresentam desvio padrão elevado ou quando a permeabilidade à água

é reduzida sem afetar a permeabilidade de solutos. O modelo forneceu uma ferramenta

eficaz para projetar e melhorar a eficácia agronômica e ambiental de CRFs revestidos de

polímeros.

TONG et al. [42] testaram o efeito de hidrogéis superabsorventes como portadores

de uréia com liberação controlada para reduzir a perda de lixiviação e aumentar o efeito de

fertilização por uréia. Também foi determinada a relação entre a taxa de liberação de uréia

e propriedades aparentes de hidrogel, além de desenvolver um modelo matemático com

Page 24: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

15

parâmetros facilmente obtidos, que poderiam ajudar a preparar formulações de liberação

controlada de uréia, cujo lançamento fosse sincronizado com as exigências de absorção das

plantas. Os métodos experimentais e de modelagem utilizado nesse estudo puderam

fornecer um método conveniente para a concepção e preparação de hidrogéis

superabsorventes com taxas de liberação adequadas para atender às necessidades práticas

das plantas, e assim produzir benefícios econômicos e ambientais.

BASU et al. [43] estudaram o tempo de saturação e de liberação de nutrientes a

partir de um grânulo esférico de fertilizante revestido com diferentes áreas de contato com

o solo, usando técnicas de modelagem matemática e simulações determinísticas. Vários

fatores que interferem nessa liberação foram estudados: o efeito do raio do grânulo, o

coeficiente de difusão, a área de contato no tempo de saturação, o tempo de liberação dos

nutrientes e o efeito da perda por evaporação no tempo de saturação. O estudo sobre esses

parâmetros que interferem na liberação de nutrientes foi bastante útil para determinar o raio

e a largura do revestimento do grânulo em sua fabricação. As condições do solo para os

grânulos de fertilizantes devem ser tal que a quantidade máxima de nutrientes liberados a

partir dos grânulos atinga as raízes das plantas vizinhas.

BASU et al. [44] estudaram a liberação de três tipos de nutrientes (Cloreto de

Potássio, Fosfato de Diamônio e Uréia) utilizados normalmente em fertilizantes comerciais,

contidos em um grânulo de fertilizante esférico. A liberação de cada nutriente é diferente

devido às suas características (um eletrólito forte, um eletrólito fraco e um não eletrólito).

Para a liberação controlada dos nutrientes foram estudados dois tipos de superfícies de

contato com o solo: um ponto e o hemisfério inferior da esfera. Foram ainda verificadas a

influência do raio do grânulo, a taxa de liberação dos nutrientes, a constante de associação,

o pH e a temperatura no tempo de liberação de nutrientes. Tais estudos são úteis na

fabricação de grânulos esféricos revestidos de acordo com as condições variáveis do solo e

também para informar os agricultores sobre o uso ideal dos grânulos, podendo definir a

forma de aplicação de nutrientes de acordo com a necessidade.

Page 25: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

16

O método sem malha baseado nas funções de base radial será aplicado ao problema

de magnetohidrodinâmica, para a solução do problema determinístico, junto com o método

de colocação estocástica. Para trabalhar com tal método foi feita uma revisão da literatura

desta técnica. Também, sobre o problema de magnetohidrodinâmica, foi elaborada uma

revisão.

COLAÇO et al. [45] empregaram o método sem malha baseado nas funções de base

radial para resolver um problema magnetohidrodinâmico bidimensional em uma cavidade

quadrada com escoamento incompressível, laminar, estacionário e com campo magnético

constante. O sistema de equações resultante foi resolvido utilizando um método quasi-

Newton e seus resultados foram comparados ao método dos Volumes Finitos. Os resultados

foram obtidos para vários tipos de distribuição de pontos no domínio. As funções de base

radial apresentaram resultados bem próximos aos da referência com tempos de computação

inferiores.

AL-NAJEM et al. [46] analisaram a influência do campo magnético no processo de

transferência de calor em uma cavidade quadrada inclinada e diferentes números de

Grashof. O método dos volumes de controle foi aplicado para resolver o problema e os

resultados obtidos para a vorticidade, temperatura e função corrente puderam ser

comparados com resultados encontrados na literatura. A análise dos resultados mostrou que

a característica do escoamento depende fortemente do ângulo de inclinação como também

do campo magnético aplicado na cavidade.

RUDRAIAH et al. [47] apresentaram um estudo numérico para a convecção

natural de um fluido eletricamente condutor em uma cavidade retangular bidimensional na

presença de um campo magnético vertical. As duas paredes verticais da cavidade eram

sujeitas a diferenças de temperatura e as paredes horizontais eram adiabáticas. Um esquema

de diferenças finitas que consiste no método modificado ADI (Alternating Direction

Implicit) e no método SLOR (Successive Line Over Relaxation) foi usado para resolver a

formulação função corrente-vorticidade do problema. Os resultados numéricos mostraram

que o campo magnético reduz a taxa de transferência de calor por convecção e o número de

Nusselt médio diminui com o aumento de número de Hartmann. Para estes resultados,

Page 26: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

17

diferentes números de Grashof e Hartmann foram empregados mantendo o número de

Prandtl fixo.

MAGALHÃES et al. [48] descreveram a solução das equações de Navier-Stokes e

da energia bidimensionais pelo método baseado nas funções de base radial (RBF). A

aproximação por funções de base radial foi utilizada nas equações governantes, resultando

num sistema não linear, o qual foi resolvido por um método quasi-Newton. A aplicação de

RBF elimina a tarefa de discretizar as equações e possíveis dificuldades na utilização de

uma malha co-localizada. O problema estudado refere-se a uma cavidade com convecção

natural, onde foram encontrados bons resultados com apenas poucos pontos de colocação.

COLAÇO et al. [49] aplicaram o método de Kansa, também conhecido como

método de colocação de funções de base radial para resolver problemas de difusão, como

transferência de calor em regime permanente e também problemas de convecção-difusão,

como o de regime permanente incompressível das equações de Navier-Stokes. No primeiro

caso, as multiquádricas renderam melhores resultados e as soluções apresentadas foram

obtidas após uma boa escolha do parâmetro de forma c. No caso das equações de Navier-

Stokes, os resultados mais significativos foram encontrados com a aplicação das funções de

Wendland.

CHINCHAPATNAM et al. [50] apresentaram um método sem malha baseado nas

funções de base radial para resolver as equações de Navier-Stokes. Utilizaram a formulação

função corrente-vorticidade chegando a uma equação bi-harmônica não linear. O sistema

não-linear foi resolvido pelo método de Levenberg Marquardt. Problemas com condições

de contorno foram solucionados com a estratégia de “centros fantasmas”. O método foi

aplicado para a cavidade quadrada e retangular com tampa deslizante para vários números

de Reynolds. Foram obtidos bons resultados com apenas poucos pontos, os quais foram

comparados com resultados benchmark.

Da revisão da literatura, identificou-se que o método de polinômio caos tem sido

usado em vários problemas, em especial problemas de condução de calor unidimensional

em regime permanente ou transiente, bem como em alguns problemas de convecção

natural. A proposta desta tese é a aplicação deste método a problemas mais complexos,

como aquele envolvendo a convecção natural na presença de forças de corpo externas,

Page 27: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

18

como por exemplo um problema de magnetohidrodinâmica. Neste tipo de problema, podem

existir incertezas quanto aos valores de certas propriedades termofísicas, o que justifica sua

aplicação. Inicialmente, na qualificação, foi apresentada a aplicação do método de

polinômios caos a dois problemas já encontrados na literatura, que teve como objetivo o

aprendizado do método e a validação dos resultados frente a outros já existentes. Agora, a

metodologia será utilizada na modelagem da liberação controlada de fertilizantes e no

problema de magnetohidrodinâmica. A modelagem da liberação controlada de fertilizantes

envolvendo a abordagem estocástica é inovadora visto que os trabalhos encontrados na

literatura abordam somente o problema determinístico.

Page 28: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

19

Capítulo 3

Fundamentação Teórica

Neste capítulo será apresentada a teoria geral dos polinômios caos, tanto para a

formulação de Galerkin, quanto para o uso com pontos de colocação estocástica. Junto a

isto, o método de colocação utilizando a interpolação com funções de base radial (RBF)

será estudado a fim de aplicá-lo a problemas onde exista a necessidade de utilizar

métodos sem malha. Suas aplicações serão apresentadas no capítulo seguinte.

3.1 Polinômios Caos

Dentre os métodos não estatísticos, o polinômio caos generalizado baseia-se na

teoria do caos de Wiener [9] e tive seu início com o PC (polinômios caos clássico) por

Ghanem e Spanos em muitos trabalhos, especialmente em [1] onde foi combinada a

expansão de Hermite com o método de elementos finitos para modelar incertezas

encontradas em vários problemas de mecânica dos sólidos. Karniadakis e outros

aplicaram uma expansão de polinômios caos generalizados para modelar incertezas em

difusão e problemas de condução de calor [10-12], [26- 28].

Historicamente, uma expansão com base em polinômios ortogonais de Hermite

em termos de variáveis aleatórias Gaussianas foi empregada. Depois, um esquema geral,

chamado Polinômios Caos Generalizado [10-16] (gPC) ou Esquema Askey de

polinômios foi introduzido. Este esquema geral inclui uma família de polinômios

ortogonais que satisfaz o tipo de equação diferencial analisada. Isto significa que o tipo

de polinômio ortogonal e a função de ponderação são escolhidos de acordo com a

Page 29: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

20

entrada estocástica, o que determina o tipo de polinômio ortogonal que será usado como

base [11]. Isso permite ao método alcançar uma convergência ótima por escolher uma

base adequada. Neste esquema de Askey, as variáveis aleatórias não se restringem a

variáveis aleatórias Gaussianas, sendo os polinômios de Hermite um subconjunto deste

esquema.

O esquema de Askey, também chamado de Wiener-Askey, é apresentado mais a

frente (Tabela 1). Nele encontram-se a função densidade de probabilidade e o polinômio

ortogonal referente a cada tipo de variável aleatória. A seleção ótima da base polinomial

deriva da sua ortogonalidade com respeito às funções de ponderação que correspondem

às funções densidade de probabilidade (PDFs).

3.1.1 Polinômios Caos Generalizado

Os polinômios caos generalizado aproximam uma processo estocástico através

de polinômios ortogonais de variáveis aleatórias, sendo um meio de representar

processos estocásticos de segunda ordem ( , ; )X t x , é um evento randômico e t é o

tempo e x o espaço. De acordo com esta metodologia, pode-se escrever ( , ; )X t x

como [11]:

1 1

1

1

1 2 1 2

1 2

1 2

1 2 3 1 2 3

1 2 3

0 0 1

1

2

1 1

3

1 1 1

( , ; ) ( , ) ( , ) ( ( ))

( , ) ( ( ), ( ))

( , ) ( ( ), ( ), ( )) ...,

i i

i

i

i i i i

i i

i i

i i i i i i

i i i

X t a t a t

a t

a t

x x x

x

x

(3.1)

onde 1

( ,..., )nn i i denota o polinômio caos generalizado de grau n nas variáveis

1( ,..., )

ni i . Tais polinômios são ortogonais em termos da variável aleatória

multidimensional 1

( ,..., )ni i ξ , com 0 1 . Por conveniência, reescreve-se a equação

(3.1) da seguinte forma

0

ˆ( , ; ) ( , ) ( )i i

i

X t a t

x x ξ

(3.2)

Page 30: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

21

onde existe uma correspondência unívoca entre as funções 1

( ,..., )nn i i e ( )i ξ , bem

como entre seus respectivos coeficientes 1 2 3 ...i i ia e ˆ

ia . Os polinômios ortogonais

satisfazem a relação de ortogonalidade

2

i j i ij (3.3)

onde ij é a função delta de Kronecker e é o produto interno entre eles. Por

definição, o produto interno entre polinômios ortogonais é dado por

( ), ( ) ( ) ( ) ( )f g f g w d ξ ξ ξ ξ ξ ξ

(Caso contínuo)

(3.4)

( ), ( ) ( ) ( ) ( )f g f g w

ξ ξ ξ ξ ξ

(Caso discreto)

(3.5)

onde ( )w ξ denota a função densidade de probabilidade e o intervalo de integração ou

do somatório é o domínio onde a variável aleatória é definida.

O polinômio ortogonal é escolhido de acordo com a função densidade de

probabilidade da variável aleatória em análise. Junto a isso, também é escolhida a

função de ponderação, que é empregada na relação de ortogonalidade.

A equivalência entre polinômios e variáveis aleatórias é obtida na Tabela 1.

Em termos práticos, a equação (3.2) divide o campo estocástico ( , ; )X t x em

uma parte determinística, o coeficiente ˆ ( , )ia tx , e a parte estocástica, o polinômio caos

( )i ξ . As expansões (3.1) ou (3.2) são truncadas em um número finito de termos, por

exemplo P termos, onde a variável aleatória é um vetor aleatório n -dimensional.

0

ˆ( , ; ) ( , ) ( )P

i i

i

X t a t

x x ξ (3.6)

Page 31: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

22

Tabela 1: Correspondência de polinômios ortogonais e funções densidade de

probabilidade de diferentes variáveis aleatórias (N≥0 inteiro finito) [11]

Variáveis

Aleatórias

Função densidade

( )w

Polinômios

Ortogonais ( )

Domínio

Contí

nua

Gaussiana

2 2

2

xe

Hermite ( , )

Gamma ( 1)

xe x

Laguerre 0,

Beta 1

(1 ) (1 )

2 ( 1, 1)

x x

B

Jacobi 1,1

Uniforme 1

2 Legendre 1,1

Dis

cret

a

Poisson !

xe

x

Charlier 0,1,2,...

Binomial (1 )x N xN

p px

Krawtchouk 0,1,..., N

Binomial

negativa

(1 )

!

x

xc c

x

Meixner 0,1,2,...

Hipergeométrica x N x

N

Hahn 0,1,..., N

O número total de termos de uma expansão de polinômios caos é 1P (quando o

somatório se iniciar em zero), obtido pela dimensionalidade da variável aleatória ( n )

e pela ordem da expansão polinomial desejada ( p ), dada por [4]

( )!1

! !

n pP

n p

(3.7)

3.2 Aplicações de Polinômios caos a equações diferenciais

O método de polinômio caos, descrito brevemente na seção anterior, pode ser

aplicado à solução de equações diferenciais, tal como foi objeto deste trabalho. Apenas

a título ilustrativo, considere a equação diferencial estocástica com condição de

contorno:

Page 32: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

23

( , , ; ) ( , ; )L t u f t x x

(3.8.a)

( , , ; ) ( , ; )B t u g t x x

(3.8.b)

onde é o parâmetro randômico, ( , ; )u t x é a solução, ( , ; )f t x e ( , ; )g t x são

termos fontes, L é o operador diferencial no espaço/tempo definido no domínio, que

pode ser não linear e B é o operador diferencial definido no contorno. Considerando a

solução ( , ; )u t x como um processo aleatório, pode-se expandi-la por polinômios caos

0

( , , ) ( , ) ( ( ))P

i i

i

u t u t

x x ξ (3.9)

onde o coeficiente ( , )iu tx é determinístico e a base polinomial ( ( ))i ξ absorve o

parâmetro randômico, isto é, a aleatoriedade é transferida para ela. A expansão para

( , ; )u t x pode ser truncada em um número finito de termos como dado em (3.6).

Substituindo a expansão polinômios caos (3.9) na equação diferencial (3.8),

obtém-se:

0

, , ; ( , ; )P

i i

i

L t u f t

x x

(3.10)

A projeção de Galerkin da equação (3.10) para cada base polinomial i é

realizada pelo produto interno da equação com cada base. Usando a ortogonalidade da

base polinomial, pode-se obter um conjunto de 1P equações acopladas para cada

coeficiente determinístico ( , )iu x t onde {0,1,..., }i P .

0

, , ; , ( , ; ),P

i i l l

i

L t u f t

x x

0,1,...,l P .

(3.11)

Assim, as equações governantes para os coeficientes da expansão de ( , , )u t x

resultantes da equação (3.11) são determinísticas.

Page 33: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

24

As discretizações no espaço e no tempo podem ser feitas por qualquer técnica

determinística usual.

A expansão por polinômios caos fornece uma boa estrutura para problemas

estocásticos. Uma abordagem peculiar para resolver os coeficientes determinísticos da

expansão por esses polinômios é a realização de uma projeção de Galerkin, já obtida

pela expressão (3.11).

A aplicação do método de Galerkin resulta em um conjunto acoplado

(característica intrusiva) de equações determinísticas, em que todos os coeficientes da

expansão por polinômios caos são calculados simultaneamente e técnicas numéricas

padrões podem ser aplicadas. Esta abordagem de Galerkin tem sido aplicada desde os

primeiros trabalhos com polinômios caos clássico demonstrando resultados eficazes.

Entretanto, quando as equações estocásticas do problema admitem formas complexas

este método torna-se um desafio. Além disso, por ser um método intrusivo, a

modificação do código computacional torna-se necessária quando se modifica a

expansão de polinômios caos em um problema.

3.3 Método de Colocação Estocástica

O método de colocação estocástica [17-23] consiste em obter uma função de

interpolação multidimensional que representa a solução do sistema em função das

variáveis randômicas. Essas funções de interpolação são mutuamente ortogonais

tornando o sistema resultante desacoplado.

Nesta abordagem, a solução determinística é avaliada em vários pontos do

espaço estocástico e, em seguida, constrói-se uma função de interpolação que melhor

aproxima a solução.

Na análise de problemas, primeiramente define-se uma formulação de expansão

por polinômio de Lagrange para cada variável estocástica, depois substituem-se essas

aproximações na equação governante e resolve-se o problema para cada ponto de

colocação, estimando-se os coeficientes determinísticos de suas variáveis estocásticas

correspondentes. Esses pontos de colocação, para o caso de variável estocástica

unidimensional, correspondem aos pontos da quadratura de Gauss [51].

A escolha da função de ponderação para a regra da quadratura de Gauss é muito

importante. Para garantir a convergência, a função de ponderação é igual à função

densidade de probabilidade do parâmetro de incerteza.

Page 34: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

25

Este método se aproxima do método de Galerkin à medida que também utiliza

um polinômio ortogonal como base, com o diferencial inerente ao método de colocação

que aplica polinômios de Lagrange, gerando um sistema desacoplado. Este fato

determina a característica não intrusiva do método, necessitando apenas de realizações

repetitivas do código computacional para obter soluções.

A garantia de convergência está ligada à escolha de um bom conjunto de pontos

de colocação para obter uma convergência rápida e uma alta precisão do método de

colocação.

O método de colocação, ao contrário do método de Galerkin, requer um número

maior de equações a serem resolvidas, pois cada equação será resolvida para cada ponto

de colocação. Por isso, sua aplicabilidade é restrita ao menor número de variáveis

aleatórias, pois o número de pontos de colocação cresce exponencialmente à medida

que a dimensão estocástica aumenta.

Neste caso, uma técnica utilizada para o problema do crescimento exponencial

dos pontos necessários à solução é servir-se de malhas esparsas, a qual reduz o número

de pontos de colocação. Esses seriam pontos “ótimos” da solução. Esta técnica será

vista mais a frente.

Em resumo, o método de colocação é baseado na construção de polinômios de

interpolação em um espaço aleatório multidimensional. Esse método exige apenas a

resolução de um conjunto de equações determinísticas em cada ponto pré-determinado

no espaço aleatório. Usando esses resultados obtidos em cada ponto, a interpolação é

construída, fornecendo uma aproximação para a solução estocástica.

Ao empregar os polinômios de Lagrange é necessário que os pontos de

colocação sejam escolhidos de tal maneira que as integrais a serem aproximadas

possuam o menor erro possível. Este fato para apenas uma variável estocástica ( 1N )

pode ser obtido pelo uso dos pontos da quadratura de Gauss.

3.3.1 Colocação via Interpolação de Lagrange

Nesta abordagem, a solução do problema e cada variável dependente do

parâmetro de incerteza são expandidas como:

1

( , , ) ( , ) ( )PN

i i

i

u t x u x t h

(3.12)

Page 35: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

26

onde o coeficiente ( , )iu x t é a solução ( , ; )u t x para cada ponto de colocação i , com

1,..., Pi N e ( )ih é o polinômio caos por interpolação de Lagrange correspondente ao

ponto de colocação i . PN é o número de pontos de colocação. O polinômio por

interpolação de Lagrange é uma função em termos da variável aleatória i , a qual é

escolhida tal que o parâmetro de incerteza do problema seja uma transformação linear

dessa variável. Os polinômios caos por interpolação de Lagrange são polinômios ( )h

de ordem P que passam por 1PN P pontos de colocação, com ( )i j ijh ( ij é o

delta de Kronecker), onde:

0

( )P

ki

k i kk i

h

(3.13)

Os pontos de colocação correspondem aos pontos da quadratura de Gauss,

usados para integrar a função ( , ; )u t x no domínio da variável aleatória. A solução é

então integrada a fim de se obter sua média e variância.

Quando múltiplos parâmetros de incerteza estão envolvidos, os pontos de

colocação são obtidos pelo produto tensorial dos pontos unidimensionais. O número de

pontos de colocação é então dado por ( 1)N

PN P , onde P é a ordem do polinômio e

N é a dimensão do problema estocástico (número de parâmetros de incertezas).

Para se obter os pontos da quadratura de Gauss e os pesos adequados à solução

se estabelece o procedimento descrito na seção a seguir.

3.3.2 Cálculo dos pontos da quadratura de Gauss e seus respectivos pesos

Um algoritmo para gerar a regra da quadratura Gaussiana definida pela função

peso é dada pela relação de recorrência de três termos. Essa relação é definida para uma

sequência de polinômios que são ortogonais com respeito à função peso em um

intervalo de integração (algoritmo de Golub-Welsch [51]). Os coeficientes de

recorrência são calculados utilizando um método clássico para medidas contínuas da

função peso, o procedimento de discretização de Stieltjes [32], que consiste em uma

sequência de cálculos de coeficientes i e i iniciados com 1i indo até Pi N , onde

Page 36: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

27

PN é o número de pontos de colocação. Esse método usa o fato de que os coeficientes

de recorrência podem ser expressos em termos de polinômios ortogonais. Essa

discretização pela função peso leva a um algoritmo estável como pode ser visto em [26].

A Relação de Recorrência de Três Termos é dada por:

1 1( ) ( ) ( ) ( )i i i i i

1,2,..., Pi N

(3.14)

0( ) 0 1( ) 1

onde i e i são coeficientes de recorrência determinados por (fórmula de Darboux

[9]):

( ), ( )

( ), ( )

i i

i

i i

1,2,..., Pi N

(3.15)

1 1

( ), ( )

( ), ( )

i i

i

i i

2,3,..., Pi N

(3.16)

onde .,. denota o produto interno e o primeiro coeficiente 1 é obtido por

1 1( ), ( ) . O produto interno entre duas funções p e q é definido como:

1

( )

( ), ( ) ( ) ( ) ( ) ( ) ( )N

N j j j

jD w

p q p q d p q

(3.17)

onde 1

( ) ( )N

N j j

j

com 0j , D é o domínio do parâmetro de incerteza

, é a função delta Dirac e ( )w é a função de ponderação dada pela função

densidade de probabilidade.

A partir dos coeficientes de recorrência i e i , 1,2,..., Pi N , os pontos de

colocação i e seus correspondentes pesos iw são calculados usando o algoritmo de

Golub-Welsch. A matriz Jacobiana J que é obtida pela relação de recorrência de três

termos conforme mostrado no Apêndice A pode ser construída como:

Page 37: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

28

1 2

2 2 3

3 3

1

1 1

0 0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0 0

P

P P P

P P

N

N N N

N N

J

(3.18)

Os autovalores de J são os pontos de colocação i , 1,2,..., Pi N , os quais são

raízes dos polinômios de ordem PN de um conjunto de polinômios ortogonais. A

distribuição de é usada para mapear os pontos de colocação dos parâmetros de

entrada. Os pesos são calculados por 2

1 1,i iw v , 1,2,..., Pi N , onde 1,iv é o primeiro

componente do autovetor normalizado correspondente ao autovalor i .

A média e a variância da solução são calculadas por:

1

( , )PN

u i i

i

u x t w

(3.19)

2

2 2

1 1

( ( , )) ( , )P PN N

u i i i i

i i

u x t w u x t w

(3.20)

onde ui(x,t) é a solução obtida no ponto de colocação i com o respectivo peso wi.

Em resumo, distribuições específicas de entradas para todos os parâmetros

incertos são determinadas por sua função densidade de probabilidade, através da qual

podem ser determinados seus pontos de colocação e seus pesos correspondentes.

Com os pontos de colocação obtidos, a equação diferencial governante do

problema analisado pode ser resolvida. Primeiramente, cada variável aleatória presente

na equação governante é substituída pela expansão dada por (3.12). Depois, essa nova

equação é calculada para cada ponto de colocação, o que torna a equação governante

um sistema de equações desacopladas.

Page 38: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

29

3.4 Método de Colocação Estocástica para variáveis estocásticas

multidimensionais

No caso de mais de uma variável aleatória, uma metodologia empregada é a

construção de pontos de colocação e funções de interpolação pelo produto tensorial dos

pontos de colocação e funções de interpolação unidimensionais. O produto tensorial dos

nós unidimensionais generaliza as propriedades da interpolação unidimensional. Então,

na equação de expansão unidimensional (3.12), o coeficiente ( , )iu x t será substituído

pelo coeficiente ( , )i

ju x t e o polinômio de Lagrange ( )ih será representado pelo termo

( )i

jh , onde o índice i denota o nível de interpolação (corresponde a quantidade de

pontos de colocação) e o índice j denota a quantidade de pontos no determinado nível

i . Logo, a solução ( , ; )iu t x para cada ponto de colocação i

j , com 1,..., Pj N ( PN é o

número de pontos de colocação) em cada nível i pode ser expressa por:

1

( , , ) ( , ) ( )PN

i i i

j j

j

u t x u x t h

(3.21)

a fim de poder empregá-la como base do produto tensorial. A expansão (3.21) denota

somente o caso de uma variável estocástica ( 1N ) como também a expansão (3.12).

Para o caso multidimensional, a expansão (3.21) pode ser estendida a partir do produto

tensorial das funções de interpolação como:

1

1

1 1

1 1

1 1 1

( ( , , ) ... ( , , ))

... ( , ),..., ( , ) . ( ) ... ( )

N

PP N

N N

N N

N

ii

NN

i ii i

j j j j

j j

u t x u t x

u x t u x t h h

(3.22)

A desvantagem desta técnica é que o número de pontos necessários para uma

boa solução cresce rapidamente com altas dimensões do espaço estocástico. Por

exemplo, se são utilizados “ q ” pontos em cada dimensão, então o número total de

pontos em um espaço N-dimensional é NQ q . Para uma aproximação com apenas três

pontos ( 3q ) em cada dimensão, 3NQ para 1N (por exemplo, para 10N ,

10 43 ~ 6 10 ), tornando a técnica do produto tensorial pouco empregada.

Page 39: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

30

Malhas esparsas têm sido estudadas no contexto de interpolação e integração

multivariada. As malhas esparsas, baseadas no algoritmo de Smolyak, são um

subconjunto do produto tensorial da malha. O subconjunto é escolhido estrategicamente

de tal forma que a aproximação das propriedades para 1N são preservados para

1N , tanto quanto possível.

3.4.1 Método de Colocação com malhas esparsas (Algoritmo de Smolyak)

O algoritmo de Smolyak [33-35] fornece uma maneira de construir funções de

interpolação com base em uma quantidade mínima de pontos no espaço

multidimensional estendendo as funções de interpolação unidimensionais para o caso

multidimensional através de uma forma particular do produto tensorial. Esse algoritmo

é uma combinação linear da fórmula do produto (3.22) com a propriedade característica

de usar somente o produto com um número relativamente pequeno de pontos.

Pode-se então definir o interpolante esparso ,q NA com produtos de funções

unidimensionais da seguinte forma [33]:

1i

,

1 i

1( , , ) ( 1) . ...

iN

q ii

q N

q N q

NA t x u u

q

(3.23)

onde N é o número de dimensões estocásticas, q N é a ordem de interpolação e ki

representa o nível de interpolação na k -ésima dimensão com 1,..., Nk N .

O algoritmo de Smolyak constrói a função de interpolação adicionando uma

combinação de funções unidimensionais de ordem ki com a restrição de que i esteja

entre 1q N e q ( 1i ... Ni i ).

O método de colocação com malhas esparsas apresenta uma boa escolha de

pontos, pois o erro obtido na interpolação é da mesma ordem que o alcançado pelo

produto tensorial completo em um mesmo nível de interpolação com a vantagem de

utilizar um número mínimo de pontos [33].

Podemos também definir o interpolante esparso (3.23) de outra forma

considerando i como a diferença entre duas interpolações aproximadas como:

Page 40: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

31

1i i iu u , 0 0u (3.24)

Utilizando a interpolação incremental (3.24) no interpolante esparso (3.23),

temos:

1

1

,

1,

( , , ) ( ... )( , , )

( , , ) ( ... )( , , )

N

N

ii

q N

i q

ii

q N

i q

A t x x t

A t x x t

(3.25)

Para calcular a interpolação ,q NA é necessário avaliar a função em cada ponto da

malha esparsa ,q NH onde iY é o conjunto de pontos em cada nível de interpolação i :

1

,

1 i

( ... )Nii

q N

q N q

H Y Y

(3.26)

A construção do algoritmo permite utilizar todos os resultados gerados

anteriormente para melhorar a interpolação. Ao escolher pontos apropriados por

interpolação da função unidimensional, pode-se garantir que o conjunto de pontos

iY está contido no conjunto de pontos 1iY ( 1i iY Y ). Logo, para calcular a função em

um nível superior ( 1i ), somente deve-se avaliá-la nos pontos da malha exclusivos

daquele nível. Para expressar a malha esparsa através de pontos da malha exclusivos de

cada nível de interpolação (iY ):

1

,

i

( ... )Nii

q N

q

H Y Y

(3.27)

3.4.2 Seleção dos pontos de colocação

Esta seção abordará o principal elemento do método de colocação: a

localização de seus pontos. Como já foi falado, quando utilizamos a interpolação por

polinômios de Lagrange é necessário que os pontos de colocação sejam escolhidos de

tal maneira que as integrais a serem aproximadas possuam o menor erro possível. Para o

caso de variável estocástica unidimensional, isso pode ser obtido pelo uso dos pontos da

Page 41: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

32

quadratura de Gauss. Para o caso de variáveis estocásticas multidimensionais pode-se

aplicar abscissas de Clenshaw-Curtis [33-34], [52] entre outras, pois apresentam as

melhores propriedades para o polinômio de Lagrange.

Tomando o suporte da variável aleatória unidimensional i no intervalo [0,1]

para 1,...,i N , pode-se expressar, sem perda da generalidade, que para o caso de

variável aleatória multidimensional, o conjunto de pontos iY terá suporte no intervalo

[0,1]N (um hipercubo), no qual as variáveis aleatórias definidas em qualquer intervalo

[a, b] podem sempre ser mapeadas para o intervalo [0,1].

As abscissas de Clenshaw-Curtis são obtidas a partir dos extremos dos

polinômios de Chebyshev e para qualquer 1im , são dadas por:

i

jY

( 1)cos 0.5, 1,..., 1

2( 1)i i

i

jpara j m se m

m

0.5, 1 1ipara j se m

(3.28)

em que o número de abscissas im em cada nível é encontrado pela expressão:

1

1 1 2 1 1i

im e m para i

Esta escolha de distribuição de pontos de forma pré-determinada verifica a

propriedade ( 1i iY Y ) em que aproveita a vantagem da recorrência dos pontos de

colocação para níveis de interpolação superiores.

Abaixo são apresentados exemplos da distribuição de pontos obtidos pela

abscissa Clenshaw-Curtis para 2N , malha esparsa bidimensional, nos oito primeiros

níveis no intervalo [ 1,1] . O código para a geração dessas malhas esparsas e os

respectivos conjuntos de pesos referentes aos pontos de tais malhas foi obtido da página

pessoal do professor John Burkardt da Universidade Estadual da Flórida nos Estados

Unidos [67].

Page 42: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

33

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

-1-0.500.51-1

-0.5

0

0.5

1

(a) Nível 0. (b) Nível 1.

-1 -0.5 0 0.5 1-1

-0.5

0

0.5

1

-1 -0.5 0 0.5 1-1

-0.5

0

0.5

1

(c) Nível 2. (d) Nível 3.

-1 -0.5 0 0.5 1-1

-0.5

0

0.5

1

-1 -0.5 0 0.5 1-1

-0.5

0

0.5

1

(e) Nível 4. (f) Nível 5.

-1 -0.5 0 0.5 1-1

-0.5

0

0.5

1

-1 -0.5 0 0.5 1-1

-0.5

0

0.5

1

(g) Nível 6. (h) Nível 7.

Figura 3.1: Malhas esparsas 2D para abscissas de Clenshaw-Curtis. (a) Nível 0, (b)

Nível 1, (c) Nível 2, (d) Nível 3, (e) Nível 4, (f) Nível 5, (g) Nível 6 e (h) Nível 7.

Page 43: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

34

3.5 Funções de Base Radial (RBF)

Os métodos sem malha emergiram como técnicas numéricas eficazes para

resolver problemas da ciência e da engenharia. O número de métodos que têm sido

propostos recentemente é uma evidência do interesse crescente da engenharia nestes

tipos de técnicas numéricas. Apesar de métodos clássicos para solução numérica de

equações diferenciais, como o método dos Elementos Finitos, das Diferenças Finitas, ou

dos Volumes Finitos se mostrarem eficazes, algumas etapas podem gerar

inconvenientes, como a geração das malhas e a utilização de esquemas de acoplamento

pressão-velocidade em problemas incompressíveis. O problema do acoplamento

pressão-velocidade pode, em problemas bidimensionais, ser eliminado usando-se a

formulação função corrente-vorticidade. Contudo, tal transformação de variáveis produz

equações cuja solução por métodos numéricos tradicionais requer discretizações de alta

ordem a fim de minimizar erros de truncamento [53]. Pela característica de não requerer

nenhuma estrutura de malhas, métodos sem malha se destacam sobre os métodos

tradicionais.

Inseridas nesta abordagem, a técnica que faz uso das funções de base radial

capazes de construir um esquema de interpolação de boa qualidade para a resolução de

equações diferenciais parciais [54] é a mais utilizada. O conceito de função de base

radial (RBF) foi empregado inicialmente em aproximações multidimensionais por

Hardy [55] e depois em trabalhos baseados em RBFs do tipo multiquádrica, aplicadas à

derivadas parciais e a equações diferenciais parciais por Kansa [56, 57].

As funções de base radial podem ser definidas da seguinte forma:

Seja : R+ R uma função contínua que depende apenas da distância entre o

centro jx e um ponto qualquer x . Se cada centro jx pertence ao domínio, então

( )j jx x x , onde é a norma Euclidiana e j é chamada de função de base

radial.

As RBFs se dividem em funções locais e globais. As funções locais fornecem

uma resposta significativa apenas na vizinhança do centro e são definidas em uma parte

do domínio (suporte compacto), já as funções globais são definidas em todo o domínio

(suporte global) [58].

Page 44: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

35

Neste trabalho será utilizada a RBF de suporte global multiquádrica definida

como:

2

2

1

( )N

j j j

j

x x x c

(3.29)

onde c é chamado de parâmetro de forma.

3.5.1 Parâmetro de forma (c) e distribuição de centros no domínio

Será discutido aqui um aspecto importante para a implementação das RBFs.

Chamada de parâmetro de forma, a constante c determina a região de influência da

RBF de suporte global (GSRBF). A precisão da solução numérica depende do valor

escolhido para o parâmetro da forma. Neste trabalho, o parâmetro é estimado através da

expressão:

c d (3.30)

onde é um inteiro positivo e d é a menor distância entre dois pontos quaisquer do

domínio. Logo o parâmetro c é definido inicialmente como a menor distância entre dois

centros no domínio. A partir daí, aumenta-se o valor do parâmetro de acordo com

(3.30).

O valor da constante se inicia por 1 e vai aumentando discretamente. O valor

de d é sempre o mesmo para cada malha utilizando o mesmo critério regular de

distribuição de centros [50]. Neste trabalho, os sistemas não lineares, resultantes da

aplicação das funções de base radial nas equações diferenciais parciais, foram

resolvidos pelo método Broyden [70], que é um método quasi-Newton. Assim, o melhor

valor do fator de forma “ c ” foi aquele que gerou o menor resíduo na solução dos

sistemas não lineares. Em geral, quanto maior o parâmetro c , melhor a solução. Porém,

isso implica no fato do sistema ficar muito mal condicionado [50]. Assim, em geral o

melhor c é aquele que seja grande o suficiente para gerar uma boa aproximação para a

solução, porém não tão grande a ponto de tornar a solução do sistema de equações não

lineares inviável. Portanto, quanto mais próximo do mal condicionamento estiver a

matriz, melhor a solução, até o ponto em que qualquer novo aumento no valor de c

resultar em um sistema tão mal condicionado que nenhuma resposta significativa será

encontrada. Foi utilizado a variação da constante de 1 até 5, obtendo-se cinco valores

Page 45: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

36

para o parâmetro de forma. Destes cinco valores o que resultar no menor erro residual

do sistema resolvido será o melhor valor de c .

Quanto à distribuição de pontos, no presente trabalho foram usadas distribuições

regulares de centros, em que os centros foram uniformemente espaçados no domínio de

interesse.

3.5.2 Expansão com RBF

Seja 1 2, ,..., n

Nx x x um conjunto de centros. A função de base radial é a

função j , definida anteriormente. Cada função j é radialmente simétrica sobre o

centro jx . Dado ( )i if f x , 1,...,i N , a aproximação de interpolação com RBF é:

2

2

1 1

( ) ( )N N

i j j j i j j

j j

s x s x s x x c

(3.31)

onde N é o número de centros. Os coeficientes da expansão js , são escolhidos de modo

que ( ) ( )i is x f x . Esses coeficientes js , são obtidos pela resolução do sistema de N

equações lineares:

Hs f

(3.32)

onde 1 ,...,i i N ix x x x H é a matriz conhecida das RBFs,

1,...,T

Ns ss é o vetor de parâmetros desconhecidos, 1,...,T

Nf ff é o vetor de

termos independentes e 1,...,i N .

3.5.3 Aplicação do método de colocação por RBF à equações diferenciais

Utilizando a expansão em termos de funções de base radial, é possível

solucionar um sistema de equações diferenciais parciais. Assim como o método de

colocação estocástica já visto, basta aplicar os operadores diferenciais correspondentes à

função de base radial e depois usar colocação em um conjunto de centros no domínio e

no contorno. Em resumo, a colocação é a aplicação dos operadores diferenciais do

Page 46: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

37

domínio ( L ) e do contorno ( B ) a um conjunto de N M centros no domínio ( ) e

M centros no contorno ( ), respectivamente. Assim, pode-se obter o seguinte sistema

de equações lineares:

1

( )N

i j i j

j

Lu x s L x x

1,...,i N M

(3.33)

1

( )N

i j i j

j

Bu x s B x x

1,...,i N M N

(3.34)

Page 47: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

38

Capítulo 4

Problemas abordados

Neste capítulo descrevem-se os problemas analisados, as formulações

matemáticas, a implementação do método baseado em polinômio caos e colocação

estocástica, bem como uma extensão para variáveis aleatórias bidimensionais utilizando

malhas esparsas. Também será apresentada a implementação do método baseado em

função de base radial para o problema determinístico de magnetohidrodinâmica.

De início, na etapa de qualificação, foi estudada uma equação diferencial

ordinária estocástica encontrada na referência [11] e depois um problema

unidimensional transiente de transferência de calor encontrado em [12]. Os problemas

ora em análise foram escolhidos para validar a implementação de alguns métodos

abordados no capítulo anterior. Como continuidade, a modelagem da liberação

controlada de fertilizantes abordada em [40] foi estudada, alocando incertezas no

coeficiente de distribuição de fertilizantes e depois no coeficiente de difusão. No

problema de Magnetohidrodinâmica [45] foram consideradas incertezas no número de

Grashof e no número de Hartmann.

Page 48: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

39

4.1. Problema 1: Modelagem da liberação controlada de fertilizantes

A utilização de fertilizantes com liberação controlada (CRFs) é uma tentativa

para melhorar a sua eficácia de utilização, pois estima-se que a porcentagem de

fertilizantes recuperados pelas plantas quando aplicados nas formas convencionais

limita-se a 30-50%. Esse uso também satisfaz as necessidades do plantio, reduzindo

significativamente a poluição do meio ambiente. A liberação controlada de fertilizantes

é proposta desde 1962, ano em que o primeiro estudo sobre a aplicação dessa tecnologia

foi apresentada. O controle da liberação de fertilizantes mantém as suas concentrações

em níveis eficazes no solo e os libera quando a planta mais necessita dele. Essa

tecnologia resulta no controle máximo do uso de fertilizantes, diminuindo a quantidade

de fertilizantes empregados e reduzindo as perdas pelos efeitos da volatilização e

lixiviação [40-44].

A maioria dos CRFs são grânulos com nutrientes solúveis em água, os quais

são encapsulados (fertilizantes revestidos por uma membrana polimérica) e essa

liberação controlada se dá por difusão através do revestimento, ou seja, a liberação de

fertilizantes ocorre por toda a superfície da esfera. O estudo de CRFs é útil na

fabricação de grânulos esféricos revestidos de acordo com as condições variáveis de

solo e também para informar os agricultores sobre o uso ideal dos grânulos, podendo

definir a forma de aplicação de nutrientes de acordo com a necessidade.

Figura 4.1: Grânulo de fertilizante com recorte da membrana polimérica [72]

O processo de liberação de fertilizantes é mostrado na Figura 4.2. Primeiro

mostra o grânulo coberto pela membrana polimérica (1), em seguida ocorre a

penetração da água pela membrana (2). Em (3) acontece a dissolução dos nutrientes

formando uma solução concentrada e em (4) ocorre a liberação controlada dos

nutrientes. [72]

Page 49: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

40

(1) (2) (3) (4)

Figura 4.2: Mecanismo de liberação de fertilizantes [72].

O modelo matemático para a liberação controlada de um fertilizante é expresso

pela equação de balanço de massa em uma partícula esférica:

2 **

2

2C C D CD

t r r r

(4.1)

onde C é a concentração do fertilizante no grânulo e *D é o coeficiente de difusão. A

condição de contorno e as condições iniciais são dadas por:

2 *4s

r R

V C Cr R R D

k t r

em 0t

(4.2)

0 0C

rr

em 0t

(4.3)

0r R C C (4.4)

0r R C (4.5)

onde sV é o volume da solução dentro da esfera que se mistura com a concentração, k é

o coeficiente de distribuição de fertilizante e 0C é a concentração inicial. A condição de

contorno representa a taxa de perda de fertilizantes pela área da esfera.

As equações podem ser escritas na forma adimensional utilizando-se as

seguintes variáveis adimensionais:

0

C C

C C

(4.6)

Page 50: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

41

*

2

D t

R

(4.7)

r

R

(4.8)

3

0

4

3 s

C R k

C C V

(4.9)

Substituindo as variáveis adimensionais (4.6)-(4.9) na equação governante (4.1),

na condição de contorno (4.2) e nas condições iniciais (4.3)-(4.5), obtém-se:

2

2

2

(4.10)

1 3

em 0 (4.11)

0 0

em 0

(4.12)

1 1 (4.13)

1 (4.14)

4.1.1. Método de Colocação aplicado à equação de balanço de massa (Dimensão

estocástica 1N )

O parâmetro de incerteza inicialmente considerado é a variável k (coeficiente

de distribuição de fertilizantes). Considerando 34

13

e

s s

VR

V V

(volume da esfera eV

igual ao volume da solução sV ) na equação (4.9), tornamos o parâmetro de incerteza,

logo:

1 0( ) ( )k (4.15)

em que 1 é a variável estocástica com coeficiente de dispersão ( ) e 0 é a variável

de referência.

Page 51: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

42

Aplicando a expansão de colocação por interpolação de Lagrange na variável

( ) e na concentração adimensional , podemos escrever:

1

( ) ( )PN

i i

i

h

(4.16)

1

( , ; ) ( , ) ( )PN

i i

i

h

(4.17)

onde ( )ih é o polinômio de Lagrange.

Aplicando essas expansões na equação diferencial e nas condições de contorno e

iniciais, obtêm-se:

2

21 1 1

( , ) ( , ) ( , )2( ) ( ) ( )

P P PN N N

i i ii i i

i i i

h h h

(4.18)

1 1 1

( ,1) ( ,1)( ) 3 ( ) ( )

P P PN N N

i ii i i i

i i i

h h h

(4.19)

1

( ,0)( ) 0

PN

ii

i

h

(4.20)

1

( , 1) ( ) 1PN

i i

i

h

(4.21)

1 1

(0,1) ( ) ( )P PN N

i i i i

i i

h h

(4.22)

Usando a propriedade do polinômio de Lagrange ( ( )i j ijh ) e aplicando essa

propriedade nas equações (4.18)-(4.22), obtém-se:

Page 52: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

43

2

2

( , ) ( , ) ( , )2i i i

(4.23)

( ,1) ( ,1)3i i

i

(4.24)

( ,0)

0i

(4.25)

( , 1) 1i

(4.26)

(0,1)i i

(4.27)

para 1,..., Pi N , onde pN é o número de pontos de colocação. Com a aplicação das

expansões (4.16)-(4.17) nas equações (4.10)-(4.14), a solução do problema será dada

para cada ponto ( i ). A média e a variância da solução serão calculadas usando as

soluções obtidas em cada ponto de colocação.

4.1.2. Método de Colocação aplicado à equação de balanço de massa (Dimensão

estocástica 2N )

Neste caso, além do coeficiente de distribuição, trata-se o coeficiente de difusão

*D como sendo uma variável aleatória, ou seja,

2 0*( ) DD D (4.28)

em que o parâmetro *D que representa o coeficiente de difusão contém incertezas, 0D

é o coeficiente de difusão de referência e 2 é a variável estocástica com coeficiente de

dispersão ( D ). Neste caso, a expressão de adimensionalização que define a variável

independente adimensional é redefinida da seguinte maneira:

Page 53: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

44

0

2

D t

R

(4.29)

Substituindo as variáveis adimensionais (4.6), (4.8)-(4.9) e (4.29) na equação

governante (4.1), na condição de contorno (4.2) e nas condições iniciais (4.3)-(4.5),

obtem-se:

2

2

2D

(4.30)

1 3 D

em 0 (4.31)

0 0

em 0

(4.32)

1 1 (4.33)

1 (4.34)

Aplicando a expansão de colocação por interpolação de Lagrange nas variáveis

( ) , ( )D e na concentração adimensional , podemos escrever:

1

( ) ( )PN

i i

i

h

(4.35)

1

( ) ( )PN

i i

i

D D h

(4.36)

1

( , ; ) ( , ) ( )PN

i i

i

h

(4.37)

Aplicando essas expansões na equação diferencial e nas condições de contorno

e iniciais, obtém-se:

Page 54: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

45

2

21 1 1 1

( , ) ( , ) ( , )2( ) ( ) ( ) ( )

P P P PN N N N

i i ii i i i i

i i i i

h D h h h

(4.38)

1 1 1 1

( ,1) ( ,1)( ) 3 ( ) ( ) ( )

P P P PN N N N

i ii i i i i i

i i i i

h h D h h

(4.39)

1

( ,0)( ) 0

PN

ii

i

h

(4.40)

1

( , 1) ( ) 1PN

i i

i

h

(4.41)

1 1

(0,1) ( ) ( )P PN N

i i i i

i i

h h

(4.42)

Também como no caso de variável estocástica unidimensional, admite-se a

propriedade do polinômio de Lagrange ( ( )i j ijh ) e aplica-se essa propriedade nas

equações (4.38)-(4.42), obtendo-se:

2

2

( , ) ( , ) ( , )2i i iiD

(4.43)

( ,1) ( ,1)3i i

i iD

(4.44)

( ,0)

0i

(4.45)

( , 1) 1i

(4.46)

(0,1)i i

(4.47)

Page 55: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

46

para 1,..., Pi N , onde pN é o número de pontos de colocação. Com a aplicação das

expansões (4.35)-(4.37) nas equações, a solução do problema será dada para cada ponto

(i ). A média e a variância da solução serão calculadas usando as soluções obtidas em

cada ponto de colocação.

4.2. Problema 2: Problema de Magnetohidrodinâmica

O problema de Magnetohidrodinâmica bidimensional, incompressível, em

regime permanente, com escoamento laminar e campo magnético constante aplicado em

uma cavidade quadrada será estudado aqui. Nessa cavidade, as paredes superior e

inferior são isoladas e as paredes da esquerda e da direita estão submetidas a diferentes

temperaturas. Com essa diferença de temperatura passa a existir uma força de empuxo,

que é modelada pela hipótese de Boussinesq. O fluido é permeado por um campo

magnético constante que leva a uma força de empuxo adicional [45].

Figura 4.3: Problema de magnetohidrodinâmica

Page 56: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

47

O escoamento convectivo é regido pela combinação da força de empuxo e da

força eletromagnética. Neste trabalho, o número de Reynolds magnético é assumido

pequeno de maneira que o campo magnético induzido causado pelo movimento do

fluido eletricamente condutor pode ser desprezado em comparação com o campo

magnético 0B [46]. Além disso, podemos desprezar o aquecimento por efeito Joule e

dissipação viscosa por ser supostamente muito pequenos.

O problema é modelado pelas equações de conservação de massa, conservação

da quantidade de movimento em x e em y e conservação de energia, dadas

respectivamente por [45]:

0u v

x y

(4.48)

2 2

2 2

1u u p u uu v

x y x x y

(4.49)

22 2

0

2 2

1( ) E

T c

Bv v p v vu v g T T

x y y x y

(4.50)

2 2

2 2T

T T T Tu v

x y x y

(4.51)

Nestas equações, u e v são as componentes de velocidade nas direções x e y ,

respectivamente, p é a pressão, T é a temperatura, é a viscosidade cinemática do

fluido, é a massa específica do fluido, g é a aceleração da gravidade, T é o

coeficiente de expansão térmica, T é a difusividade térmica do fluido, E é a

condutividade elétrica, 0B é o campo magnético constante e cT é a temperatura da

parede da direita que é menor do que a temperatura da parede da esquerda.

As condições de contorno da cavidade, como já foram citadas, correspondem às

temperaturas prescritas em ambas as paredes verticais, sendo a temperatura da parede da

esquerda maior do que a temperatura da parede da direita. As paredes horizontais são

Page 57: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

48

termicamente isoladas. Todas as condições de contorno para a velocidade incluem a

condição de não deslizamento, bem como o seu valor prescrito igual à zero em todos os

contornos. Além disso, um campo magnético constante é aplicado na direção x (da

esquerda para a direita).

0 vu 0

y

T

0y (4.52)

0 vu HTT 0x (4.53)

0 vu 0

y

T

Dy (4.54)

0 vu CTT Dx (4.55)

As equações governantes são reescritas em termos de função corrente e

combinadas em uma equação bi-harmônica, a fim de eliminar o gradiente de pressão e a

necessidade de um esquema de acoplamento pressão-velocidade, além de reduzir o

número de equações a serem resolvidas. Dessa forma, pode-se definir a função corrente

como:

uy

(4.56)

vx

(4.57)

As equações (4.48)-(4.51) juntamente com a aplicação da formulação (4.56)-

(4.57) podem ainda ser escritas na forma adimensional utilizando-se as seguintes

variáveis adimensionais:

D

xx

(4.58)

D

yy

(4.59)

Page 58: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

49

(4.60)

CH

C

TT

TTT

(4.61)

onde D é a largura da cavidade, HT é a temperatura da parece da esquerda, CT é a

temperatura da parede direita e ,, yx e T são variáveis dimensionais.

As variáveis Ha é o número de Hartmann, Gr é o número de Grashof e Pr é o

número de Prandtl, definidos por:

0Ha=B ED

3

T

2

g ( - )Gr= h cT T D

T

Pr

(4.62)-(4.64)

O número de Hartmann ( Ha ) é a relação entre a força eletromagnética e as

forças viscosas e é usado para controlar a intensidade do campo magnético. O número

de Grashof ( Gr ) é a relação entre o empuxo e as forças viscosas e é usado para

controlar a intensidade de convecção natural. O número de Prandtl ( Pr ) é uma

propriedade do fluido, representando a relação entre a difusão de quantidade de

movimento e a difusão térmica.

Substituindo as variáveis adimensionais (4.58)-(4.64) e a definição de função

corrente (4.56)-(4.57) nas equações governantes (4.48)-(4.51) e nas condições de

contorno (4.52)-(4.55), obtemos:

3 3 3 3

2 3 2 3

4 4 4 22

4 2 2 4 22 Ha Gr

y x y x x y x y

T

x x y y x x

(4.65) (4.4.2.e)

2 2

2 2

1

Pr

T T T T

y x x y x y

(4.66) (4.4.2.f)

Page 59: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

50

0 0

x

1T em 0x (4.67)

0 0

x

0T em 1x (4.68)

0 0

y

0

T

y

em 0y (4.69)

0 0

y

0

T

y

em 1y (4.70)

A combinação das equações de conservação de massa e conservação das

quantidades de movimento origina uma equação bi-harmônica a qual necessita de duas

condições de contorno para cada parede.

Os números de Hartmann e Grashof podem ser assumidos como variáveis

estocásticas devido a incertezas nas propriedades físicas como, por exemplo,

condutividade elétrica ( E ) e coeficiente de expansão térmica ( T ), dentre outras que

as compõem.

4.2.1 Método de Colocação estocástica aplicado ao problema de

Magnetohidrodinâmica (Dimensão estocástica 1N )

O parâmetro de incerteza inicialmente considerado estocástico é o número de

Grashof ( Gr ):

0 1Gr( ) Gr (1 )Gr (4.71)

em que 1 é a variável estocástica com coeficiente de dispersão dos dados ( Gr ) e 0Gr é

o Grashof de referência.

Aplicando a expansão de colocação estocástica por interpolação de Lagrange à

variável estocástica Grashof ( Gr ), temos:

1

Gr( ) Gr ( )PN

k k

k

h

(4.72)

Page 60: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

51

em que ( )kh é o polinômio de Lagrange.

Aplicando essa expansão (4.72) na equação (4.65) obtêm-se:

3 3 3 3

2 3 2 3

4 4 4 22

4 2 2 4 21

2 Ha Gr ( )PN

k k

k

y x y x x y x y

Th

x x y y x x

(4.73) (4.4.2.e)

Admitindo a propriedade do polinômio de Lagrange ( ( )k j kjh ) e aplicando

essa propriedade na equação (4.73), obtem-se:

3 3 3 3

2 3 2 3

4 4 4 22

4 2 2 4 22 Ha Gr

k k k k k k

k k k k kk

y x y x x y x y

T

x x y y x x

1,..., pk N

(4.74) (4.4.2.e)

onde pN é o número de pontos de colocação.

A equação (4.74) será resolvida para cada ponto de colocação e junto à equação

(4.66) resulta em um sistema não-linear. Neste sistema será aplicada a expansão em

RBF à solução para T e :

,

1

( , ) ( )M

k j k j

j

T x y T r

(4.75)

,

1

( , ) ( )N

k i k i

i

x y r

(4.76)

em que ,Tj k e ,i k são coeficientes desconhecidos, é o mesmo para as duas

expansões, M e N são os números de centros utilizados em cada expansão por RBF.

Substituindo as expansões (4.75)-(4.76) nas equações governantes (4.74) e

(4.65) e nas condições de contorno (4.67)-(4.70), obtêm-se o sistema:

Page 61: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

52

3 3

, , ,2 31 1 1

3 3

, , ,2 31 1 1

4 4

, ,4 2 21 1

42

, ,41

( ) ( ) ( )

( ) ( ) ( )

( ) ( )2

( )Ha

N N Ni i i

i k i k i k

i i i

N N Ni i i

i k i k i k

i i i

N Ni i

i k i k

i i

Ni

i k i

i

r r r

y x y x

r r r

x y x y

r r

x x y

r

y

2

,21 1

( )( )Gr

N Mji

k k j k

i j

rrT

x x

(4.77)

, , , ,

1 1 1 1

2 2

, ,2 21 1

( ) ( )( ) ( )

( ) ( )1

Pr

N M N Mj ji i

i k j k i k j k

i j i j

M Mj j

j k j k

j j

r rr rT T

y x x y

r rT T

x y

(4.78)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

x

,

1

( ) 1M

j k j

j

T r

em 0x (4.79)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

x

,

1

( ) 0M

j k j

j

T r

em 1x (4.80)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

y

,

1

( )0

Mj

j k

j

rT

y

em 0y (4.81)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

y

,

1

( )0

Mj

j k

j

rT

y

em 1y (4.82)

que será resolvido k vezes com 1,..., pk N . Estes sistemas de equações algébricas

serão resolvidos pelo método de Broyden [70]. Logo, para cada ponto de colocação

estocástica ( k ) será solucionado um sistema não-linear, para se achar a solução do

problema determinístico associado.

No problema determinístico, as equações (4.77)-(4.78) devem ser escritas para

cada centro dentro do domínio, e as equações (4.79)-(4.82) para cada centro da

Page 62: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

53

fronteira. Assim, para L centros no domínio tem-se L equações (4.77) e L equações

(4.78). Para P centros na fronteira tem-se 2P equações de fronteira para a equação bi-

harmônica, por ser tratar de uma equação de 4a ordem. Já para a equação da energia

tem-se P equações de fronteira. Logo, a equação bi-harmônica (4.77) tem L P

variáveis e 2L P equações e a equação da energia (4.78) tem L P variáveis e L P

equações.

Nota-se que para a equação bi-harmônica, o número de equações é superior ao

número de variáveis. Para igualar esses números, faltam P centros nas fronteiras.

Assim, para resolver esse problema, foi empregada a estratégia de “centros fantasmas”

[50] em que P centros extras são colocados fora do domínio e perto das condições de

contorno, igualando o número de equações e de variáveis, obtendo 2N L P e

M L P .

4.2.2 Método de Colocação estocástica aplicado ao problema de

Magnetohidrodinâmica (Dimensão estocástica 2N )

Neste caso, os números de Grashof ( Gr ) e Hartmann ( Ha ) serão considerados

variáveis aleatórias assumidas como:

0 1Gr( ) Gr (1 )Gr (4.83)

0 2Ha( ) Ha (1 )Ha (4.84)

onde Gr e Ha são os coeficientes de dispersão das variáveis aleatórias 1 e 2

respectivamente e 0Ha e 0Gr são números de Hartmann e Grashof de referência.

Aplicando a expansão de colocação estocástica por interpolação de Lagrange nas

variáveis estocásticas Grashof ( Gr ) e Hartmann ( Ha ), podemos escrever:

1

Gr( ) Gr ( )PN

k k

k

h

(4.85)

(4.86)

Page 63: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

54

1

Ha( ) Ha ( )PN

k k

k

h

Aplicando as expansões (4.85)-(4.86) na equação (4.62) obtêm-se:

3 3 3 3

2 3 2 3

24 4 4 2

4 2 2 4 21 1

2 ( ) Gr ( )P PN N

k k k k

k k

y x y x x y x y

THa h h

x x y y x x

(4.87) (4.4.2.e)

Admitindo a propriedade do polinômio de Lagrange ( ( )k j kjh ) e aplicando

essa propriedade na equação (4.87), obtem-se:

3 3 3 3

2 3 2 3

4 4 4 22

4 2 2 4 22 Gr

k k k k k k

k k k k kk k

y x y x x y x y

THa

x x y y x x

1,..., pk N

(4.88) (4.4.2.e)

onde pN é o número de pontos de colocação.

A equação (4.88) será resolvida para cada ponto de colocação e junto à equação

(4.66) resulta em um sistema não-linear. Neste sistema, e também nas condições de

contorno (4.67)-(4.70), serão aplicadas as expansões em RBF (4.75)-(4.76) à solução

para T e resultando no sistema de equações abaixo:

3 3

, , ,2 31 1 1

3 3

, , ,2 31 1 1

4 4

, ,4 2 21 1

42

, 41

( ) ( ) ( )

( ) ( ) ( )

( ) ( )2

( )Ha

N N Ni i i

i k i k i k

i i i

N N Ni i i

i k i k i k

i i i

N Ni i

i k i k

i i

Ni

i k k i

i

r r r

y x y x

r r r

x y x y

r r

x x y

r

y

2

, ,21 1

( )( )Gr

N Mji

k k j k

i j

rrT

x x

(4.89)

Page 64: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

55

, , , ,

1 1 1 1

2 2

, ,2 21 1

( ) ( )( ) ( )

( ) ( )1

Pr

N M N Mj ji i

i k j k i k j k

i j i j

M Mj j

j k j k

j j

r rr rT T

y x x y

r rT T

x y

(4.90)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

x

,

1

( ) 1M

j k j

j

T r

em 0x (4.91)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

x

,

1

( ) 0M

j k j

j

T r

em 1x (4.92)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

y

,

1

( )0

Mj

j k

j

rT

y

em 0y (4.93)

,

1

( ) 0N

i k i

i

r

,

1

( )0

Ni

i k

i

r

y

,

1

( )0

Mj

j k

j

rT

y

em 1y (4.94)

que será resolvido k vezes com 1,..., pk N . Estes sistemas de equações algébricas

serão resolvidos pelo método de Broyden [70]. Logo, para cada ponto de colocação

estocástica ( k ) será solucionado um sistema não-linear, encontrando a solução do

problema determinístico associado. Analogamente, assim como na seção 4.2.1, será

utilizada a estratégia de “centros fantasmas” para que o sistema não-linear seja possível

e determinado.

Page 65: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

56

4.2.3 Método quasi-Newton

O problema de magnetohidrodinâmica resulta em um sistema não linear de

equações algébricas que pode ser muito mal condicionado, dependendo do valor do

fator de forma c . Para a solução deste sistema existem vários métodos numéricos. Um

deles, o método Broyden [70], se mostrou mais robusto. Tal método é classificado como

um método quasi-Newton. No método de Newton tradicional, o cálculo da inversa da

matriz Hessiana pode ser custoso podendo demandar elevado tempo computacional,

dependendo do problema. Assim, o método quasi-Newton fornece uma aproximação da

inversa da matriz Hessiana, usando apenas informações do gradiente. A matriz

aproximação é atualizada em cada iteração usando a informação das derivadas

(gradiente) da iteração anterior. Dessa forma, o método torna-se rápido, mas não muito

custoso.

Page 66: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

57

Capítulo 5

Resultados

Nesta seção são apresentados resultados numéricos para os dois problemas

encontrados no capítulo 4. Anteriormente, no trabalho de qualificação, uma equação

diferencial ordinária estocástica aproximada pela expansão de polinômio caos foi

resolvida para quatro tipos de distribuições do parâmetro k (Gaussiana, Beta, Gamma e

Poisson). A EDO foi reproduzida a fim de ampliar o conhecimento dos métodos

baseados em polinômio caos. O problema de condução de calor unidimensional

transiente foi resolvido para a distribuição uniforme da variável aleatória. As variáveis

do problema foram aproximadas pela expansão de Legendre do conjunto de polinômio

caos Wiener-Askey. O problema unidimensional foi resolvido para polinômio caos nos

graus 1, 2, 3 e 4 e também através do método de colocação estocástica. Através destes

dois primeiros problemas pode-se constatar que o método de Galerkin gerou resultados

pouco satisfatórios quando comparados ao método de colocação estocástica, devido sua

característica inerente de obter resultados através de equações acopladas. Com os

problemas anteriores pude ampliar a confiabilidade na metodologia de colocação frente

à de Galerkin. Sendo assim, nos problemas enunciados no capítulo 4 será somente

aplicado o método de colocação estocástica que será verificado pelo método de Monte

Carlo.

A modelagem da liberação controlada de fertilizantes em uma esfera, que se

expressa por uma equação de balanço de massa, foi determinada para a distribuição

uniforme do parâmetro estocástico. Neste caso foi imposto que o coeficiente de

distribuição de fertilizante ( k ) e o coeficiente de difusão ( *D ) contêm incertezas.

Page 67: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

58

Primeiramente, a solução determinística da taxa de liberação para vários valores de

(variável relacionada ao coeficiente de distribuição de fertilizantes) foi obtida e

comparada ao artigo de referência [40]. Também uma análise de convergência de

pontos do domínio foi feita, mostrando uma convergência da solução média com 1000

pontos. Num segundo momento, adotando parâmetros estocásticos, as variáveis do

problema foram aproximadas pelo método de colocação estocástica por interpolação de

Lagrange utilizando 3, 5, 10 e 20 pontos de colocação (no caso unidimensional) e 5, 13

e 29 pontos de colocação ou níveis de interpolação 1, 2 e 3, respectivamente, no caso

bidimensional, utilizando as abscissas de Clenshaw-Curtis. Suas soluções, para a taxa

de liberação e concentração, foram comparadas com o método de Monte Carlo. As

soluções foram obtidas para três valores do coeficiente de dispersão (0.1, 0.5 e 0.9) a

fim de mostrar a influência da amplitude da incerteza na solução. Além disso, a média e

o desvio padrão da solução podem ser avaliados em conjunto, identificando a dispersão

da solução.

O segundo problema, de Magnetohidrodinâmica, exprime um caráter estocástico

por apresentar o número de Grashof e o número de Hartman, como variáveis aleatórias.

As incertezas associadas a estas variáveis, por sua vez, foram adotadas com

distribuições uniformes. Para o problema estocástico, com variáveis estocásticas

unidimensional e bidimensional, nos resultados obtidos pelo método de colocação

estocástica utilizou-se a regra de Clenshaw-Curtis para alguns níveis de interpolação.

Para comparação, foi utilizado o método de Monte Carlo em que a análise de

convergência da solução foi verificada pela variância do número de Nusselt nas paredes

verticais da cavidade. As soluções foram obtidas para dois valores do coeficiente de

dispersão (0.2 e 0.7), a fim de mostrar a influência da amplitude da incerteza na solução.

Nos resultados obtidos em cada problema com a distribuição uniforme, não é

possível definir um intervalo de confiança da mesma forma como definido na

distribuição Gaussiana. Assim, será definido um intervalo de credibilidade como sendo

uma variação de 2.576em torno da média, onde é dado pela equação (3.20) e a

média pela equação (3.19).

Page 68: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

59

5.1 Modelagem da liberação controlada de fertilizantes

Nesta seção serão mostrados resultados do problema de liberação controlada de

fertilizantes, apresentado nas equações (4.23)-(4.27). Para a verificação de alguns dos

resultados foi utilizado o artigo de Al-Zahrani [40]. Neste trabalho, o autor desenvolveu

modelos matemáticos para descrever a taxa de liberação de fertilizantes na membrana

polimérica de uma esfera de fertilizantes. A taxa de liberação de fertilizantes é obtida

pela equação:

( ,1)1T

(5.1)

onde a taxa de liberação de fertilizantes (T ) depende da concentração ( ) na casca da

esfera ao passar do tempo e do coeficiente de distribuição ( k ), que está relacionado à

variável ( ).

A Figura 5.1 apresenta a solução da concentração de fertilizantes na casca da

esfera pelo tempo para cada conjunto de pontos de discretização por diferenças finitas

(50, 100, 200, 500, 1000, 10000 pontos) distribuídos no domínio sem incerteza no

parâmetro , exibindo uma análise de convergência. Através dela, adotamos um

conjunto de 1000 pontos no domínio, o que já se pode garantir resultados confiáveis

quanto à discretização espacial.

Figura 5.1: Análise da convergência para um conjunto de pontos ( p ) no domínio.

t

Co

nce

ntr

açã

o

0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

Page 69: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

60

A Figura 5.2 apresenta a taxa de liberação obtida por Al-Zahrani [40] e no

presente trabalho, obtida por diferenças finitas. As taxas foram comparadas para

diversos valores de . Neste resultado não foram consideradas incertezas, apenas

verificou-se o programa computacional em estudo com um resultado da literatura, visto

que o artigo para verificação apenas disponibiliza o resultado para a variável taxa de

liberação.

Figura 5.2: Comparação da taxa de liberação de fertilizantes com resultado da literatura.

A Figura 5.3(a) mostra a concentração de fertilizantes pelo tempo para três

valores de (que é diretamente proporcional ao coeficiente de distribuição k )

diferentes. Pode-se constatar que aumentando-se o valor , ocorre a diminuição da

concentração em qualquer instante de tempo. Este comportamento exibe a concentração

de fertilizantes dentro da esfera que, com o tempo, libera o fertilizante tornando sua

concentração interna igual à zero. Analogamente, a Figura 5.3(b) mostra a taxa de

liberação pela raiz quadrada do tempo para três valores de diferentes, em que

aumentando-se o valor de ocorre o aumento da taxa de liberação.

t

Ta

xa

0 0.1 0.2 0.3 0.4 0.5 0.60

0.2

0.4

0.6

0.8

1

alfa 0.1,Al-Zahrani

alfa 0.5,Al-Zahrani

alfa 1.0,Al-Zahrani

alfa 0.1,presente trabalho

alfa 0.5,presente trabalho

alfa 1.0,presente trabalho

Page 70: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

61

(a) (b)

Figura 5.3: (a) Concentração x tempo para três valores de . (b) Taxa de liberação x

raiz quadrada do tempo para três valores de .

A modelagem da liberação controlada de fertilizantes abordada em [40] será

estudada na referida tese, alocando incertezas no coeficiente de distribuição ( k ) e

depois também no coeficiente de difusão ( *D ). Para solucionar o problema, foram

adotados os seguintes valores de referência: 0 1 e 0 1D . A incerteza com variável

estocástica bidimensional necessita do produto tensorial dos pontos de colocação

unidimensionais, o que pode acarretar no crescimento exagerado do número de pontos

para se obter uma boa solução. Neste caso, uma malha esparsa será empregada a fim de

se garantir soluções ótimas com apenas poucos pontos.

As soluções do presente trabalho por colocação, para uma variável estocástica,

foram obtidas para 3, 5, 10 e 20 pontos de colocação gerados pela quadratura de Gauss.

O sistema de equações diferenciais parciais determinístico resultante da aplicação do

método de colocação estocástica foi resolvido pelo método de diferenças finitas. Para o

caso de duas variáveis aleatórias (dimensão estocástica N=2), foram utilizados 5, 13 e

29 pontos de colocação gerados pela abscissa Clenshaw-Curtis [34]. As soluções foram

obtidas para três valores do coeficiente de dispersão (0.1, 0.5 e 0.9) da variável

aleatória, a fim de avaliar a solução frente à dispersão dos dados de entrada. Para a

verificação do método foi utilizado o método de Monte Carlo.

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

t

Ta

xa

0.1 0.2 0.3 0.4 0.5 0.60

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Page 71: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

62

5.1.1 Variável estocástica Unidimensional – Coeficiente de distribuição ( k )

O parâmetro contendo incerteza é dado pela equação (4.15). A distribuição de

é uniforme no intervalo [ 1,1] e para o coeficiente de dispersão de dados ( ) foram

utilizados os seguintes valores: 0.1 , 0.5 e 0.9 . A concentração ( C ) e a

taxa de liberação (T ) são dadas pelas expansões:

1

( , ; ) ( , )PN

i i

i

C w

(5.2)

1

( ,1)

( , ; ) 1

PN

i i

i

w

T

(5.3)

Para verificar os resultados obtidos pelo método de colocação foi utilizado o

método de Monte Carlo, em que foi aplicado um conjunto de 100, 500, 1000 e 5000

amostras, a fim de se observar a convergência do número adequado de amostras para

variável estocástica com 10% de dispersão (Figura 5.4(a)-(b)). Com pequena

incerteza nesta variável de entrada, o gráfico da variância da solução para a

concentração e taxa de liberação, obtidos para vários conjuntos de amostras, mostrou-se

com pequenas alterações. Pode-se considerar com isso a convergência da variância com

o método de Monte Carlo com um conjunto de 5000 amostras, a qual será considerada

na comparação abaixo.

(a) (b)

Figura 5.4: Convergência da variância com o método de Monte Carlo para 5000

amostras. (a) Concentração (b) Taxa de liberação.

t

0.2 0.4 0.6 0.8 1

10-24

10-19

10-14

10-9

10-4

MC 100

MC 500

MC 1000

MC 5000

2

t

0.2 0.4 0.6 0.8 1

10-13

10-11

10-9

10-7

10-5

10-3

MC 100

MC 500

MC 1000

MC 5000

2

Page 72: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

63

Para o método de colocação utilizou-se 3, 5 e 10 pontos para 0.1 . Um

desvio MC entre a solução média obtida por colocação e a solução média obtida pelo

método de Monte Carlo foi utilizada, para verificar a convergência da solução por

colocação.

( ) ( )( )

( )

coloc MCMC

MC

sol t sol tt

sol t

(5.4)

em que colocsol é a solução média obtida pelo método de colocação e MCsol é a solução

média obtida pelo método de Monte Carlo com um número de amostras convergentes.

A Figura 5.5(a)-(b) apresenta o desvio entre a solução média obtida pelo método

de Monte Carlo e a solução média obtida pelo método de colocação para diferentes

valores de NP (número de pontos de colocação). Pode-se verificar que não houve

mudança neste desvio, portanto pode-se dizer que o método de colocação converge com

3 pontos.

(a) (b)

Figura 5.5: (a) Desvio entre soluções para a concentração em alguns valores de tempo.

(b) Desvio entre soluções para a taxa de liberação em alguns valores de tempo.

A Figura 5.6(a)-(b) exibe a comparação entre as soluções obtidas pelo método de

Monte Carlo com 5000 amostras e pelo método de colocação com 3 pontos, ambas

convergidas com seus respectivos intervalos de credibilidade (IC). Nele pode-se

t

0 0.2 0.4 0.6 0.8 110

-4

10-3

10-2

coloc3

coloc5

coloc10

t

0 0.1 0.2 0.3 0.4 0.5 0.610

-6

10-5

10-4

10-3

coloc3

coloc5

coloc10

Page 73: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

64

verificar a influência do parâmetro de incerteza.

(a) (b)

Figura 5.6: Comparação do método de colocação com o método de Monte Carlo: (a)

Concentração. (b) Taxa de liberação.

Observou-se que para a dispersão da variável estocástica de 10% , com um

conjunto de 3 pontos de colocação, foi alcançada a convergência da solução.

Depois, a variável estocástica foi admitida com 50% ( 0.5 ) de dispersão.

Novamente para verificar os resultados obtidos pelo método de colocação foi utilizado o

método de Monte Carlo, em que foi aplicado um conjunto de 500, 1000, 5000 e 10000

amostras a fim de se observar a convergência do número adequado de amostras para

este caso (Figura 5.7(a)-(b)). Aumentando a variabilidade na incerteza, o gráfico para a

variância da solução para a concentração e taxa de liberação obtida para vários

conjuntos de amostras é possível notar uma pequena discrepância. Através da análise

gráfica pode-se garantir que, para um conjunto de 10000 amostras, uma convergência

foi alcançada.

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

coloc 3

coloc 3-IC

MC 5000

MC 5000-IC

t

Ta

xa

0.2 0.4 0.6 0.8 1

0.2

0.4

0.6

0.8

1

coloc 3

coloc 3-IC

MC 5000

MC 5000-IC

Page 74: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

65

(a) (b)

Figura 5.7: Convergência da variância com o método de Monte Carlo para 10000

amostras. (a) Concentração (b) Taxa de liberação.

Para o método de colocação utilizou-se 3, 5, 10 e 20 pontos de colocação para

0.5 . A Figura 5.8(a)-(b) exibe o desvio MC entre as soluções obtidas por Monte

Carlo e por colocação. Com este desvio verifica-se a convergência da solução por

colocação com 3 pontos.

(a) (b)

Figura 5.8: (a) Desvio entre soluções para a concentração em alguns valores de tempo.

(b) Desvio entre soluções para a taxa de liberação em alguns valores de tempo.

A Figura 5.9(a)-(b) apresenta a comparação dos métodos em análise com seus

respectivos intervalos de credibilidade. Pode-se notar a convergência da solução através

de seus intervalos de credibilidade que se apresentam idênticos. Também é possível

t

0.2 0.4 0.6 0.8 1

10-23

10-18

10-13

10-8

10-3

MC 500

MC 1000

MC 5000

MC 100002

t

0.2 0.4 0.6 0.8 1

10-11

10-9

10-7

10-5

10-3

10-1

MC 500

MC 1000

MC 5000

MC 10000

2

t

0 0.2 0.4 0.6 0.8 110

-4

10-3

10-2

coloc 3

coloc 5

coloc 10

coloc 20

t

0 0.1 0.2 0.3 0.4 0.5 0.610

-5

10-4

10-3

10-2

coloc 3

coloc 5

coloc 10

coloc 20

Page 75: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

66

notar que uma maior dispersão na variável estocástica leva ao aumento do intervalo de

credibilidade. Conclui-se, que para a dispersão de dados da variável estocástica de

50%, com um conjunto de 3 pontos de colocação, foi alcançada a convergência da

solução.

(a) (b)

Figura 5.9: Comparação do método de colocação com o método de Monte Carlo: (a)

Concentração. (b) Taxa de liberação.

Por fim, a variável estocástica foi admitida com 90% ( 0.9 ) de dispersão.

A verificação dos resultados obtidos pelo método de colocação foi feita pelo método de

Monte Carlo em que foi aplicado um conjunto de 500, 1000, 5000 e 10000 amostras a

fim de se observar a convergência do número adequado de amostras para este caso

(Figura 5.10(a)-(b)). Com uma variabilidade na incerteza de entrada de quase 100% , é

possível observar uma leve discrepância no gráfico para resultados obtidos para vários

conjuntos de amostras, em se tratando da variância das soluções obtidas por Monte

Carlo. Com isso, para uma convergência aceitável pode-se adotar um conjunto de 10000

amostras.

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

coloc 3

coloc 3-IC

coloc 5

coloc 5-IC

coloc 10

coloc 10-IC

MC 10000

MC 10000-IC

t

Ta

xa

0.2 0.4 0.6 0.8 1

0.2

0.4

0.6

0.8

1

coloc 3

coloc 3-IC

coloc 5

coloc 5-IC

coloc 10

coloc 10-IC

MC 10000

MC 10000-IC

Page 76: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

67

(a) (b)

Figura 5.10: Convergência da variância com o método de Monte Carlo para 10000

amostras. (a) Concentração (b) Taxa de liberação.

Para o método de colocação utilizou-se 3, 5, 10 e 20 pontos de colocação para

0.9 . Através da Figura 5.11(a)-(b) pode-se perceber a convergência do método de

colocação com 5 pontos, já que a partir deste conjunto de pontos o desvio entre soluções

não se modifica.

(a) (b)

Figura 5.11: (a) Desvio entre soluções para a concentração em alguns valores de tempo.

(b) Desvio entre soluções para a taxa de liberação em alguns valores de tempo.

Na Figura 5.12(a)-(b) que apresenta uma comparação com o método de Monte

Carlo, observa-se que, para a dispersão de dados da variável estocástica de 90% ,

com um conjunto de 5 pontos de colocação, foi alcançada a convergência da solução.

t

0.2 0.4 0.6 0.8 1

0

0.0005

0.001

0.0015

0.002 MC 500

MC 1000

MC 5000

MC 10000

2

t

0 0.2 0.4 0.6 0.8 1

0

0.001

0.002

0.003

0.004

0.005

0.006

MC 500

MC 1000

MC 5000

MC 10000

2

t

0.2 0.4 0.6 0.810

-4

10-3

10-2

coloc 3

coloc 5

coloc 10

coloc 20

t

0 0.1 0.2 0.3 0.4 0.5 0.610

-6

10-5

10-4

10-3

coloc 3

coloc 5

coloc 10

coloc 20

Page 77: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

68

(a) (b)

Figura 5.12: Comparação do método de colocação com o método de Monte Carlo: (a)

Concentração. (b) Taxa de liberação.

Pode-se notar nos resultados apresentados acima que, para uma variabilidade na

variável estocástica de quase 100%, apenas poucos pontos de colocação foram

necessários para uma boa solução.

5.1.2 Variável estocástica Bidimensional – Coeficiente de distribuição ( k ) e

Coeficiente de difusão ( *D )

Os parâmetros contendo incertezas são dados pela equação (4.15) para o

coeficiente de distribuição e pela equação (4.28) para o coeficiente de difusão. A

distribuição da variável estocástica é uniforme no intervalo [ 1,1] e para o coeficiente

de dispersão de dados ( ) foram utilizados os seguintes valores: 0.1 , 0.5 e

0.9 . A concentração ( C ) e taxa de liberação (T ) são dadas pelas expansões (5.2)-

(5.3).

Inicialmente, apresenta-se o gráfico com a média da solução estocástica obtida

pelo método de colocação para 5, 13 e 29 pontos de colocação, para a concentração de

fertilizantes e para a taxa de liberação, obtida por (5.2)-(5.3) juntamente com os

intervalos de credibilidade (IC) gerado pela variância (3.20), associados a essa média

(Figura 5.13(a)-(b)). Ambos, média e variância, são funções dependentes do tempo. As

variáveis estocásticas e *D foram admitidas com 10% ( 0.1 ) de dispersão.

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

sem_incerteza

coloc 5

coloc 5-IC

MC 10000

MC 10000-IC

t

Ta

xa

0.2 0.4 0.6 0.8 1

0.2

0.4

0.6

0.8

1

sem_incerteza

coloc 5

coloc 5-IC

MC 10000

MC 10000-IC

Page 78: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

69

Pode-se verificar a diminuição do intervalo de credibilidade com o aumento de pontos

de colocação.

(a) (b)

Figura 5.13: (a) Concentração x tempo para 0.1 . (b) Taxa de liberação x raiz

quadrada do tempo para 0.1 .

Para validar os resultados obtidos pelo método de colocação foi utilizado o

método de Monte Carlo em que foi aplicado um conjunto de 1000, 3000 e 5000

amostras a fim de se observar a convergência para este caso (Figura 5.14(a)-(b)). Os

gráficos apresentam resultados para a variância da solução para a concentração e taxa de

liberação e pode-se verificar a convergência da variância com o método de Monte Carlo

com um conjunto de 5000 amostras.

(a) (b)

Figura 5.14: Convergência da variância com o método de Monte Carlo para 5000

amostras. (a) Concentração (b) Taxa de liberação.

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

coloc 5

coloc 5-IC

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

t

Ta

xa

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

coloc 5

coloc 5-IC

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

t

0.2 0.4 0.6 0.8 1

0

0.0005

0.001

0.0015

0.002

MC 1000

MC 3000

MC 5000

2

t

0.2 0.4 0.6 0.8 1

0

5E-05

0.0001

0.00015

0.0002

MC 1000

MC 3000

MC 5000

2

Page 79: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

70

Para o método de colocação utilizaram-se 5, 13 e 29 pontos de colocação para

0.1 . Na Figura 5.15(a)-(b) verifica-se a convergência da solução média pelo

método de colocação com 13 pontos, quando nenhuma mudança significante no desvio

MC foi constatada em relação à solução com 29 pontos.

(a) (b)

Figura 5.15: (a) Desvio entre soluções para a concentração em alguns valores de tempo.

(b) Desvio entre soluções para a taxa de liberação em alguns valores de tempo.

A Figura 5.16(a)-(b) exibe a comparação entre os métodos de colocação e de

Monte Carlo em que é possível verificar a convergência através do intervalo de

credibilidade. Observou-se que para a dispersão de dados das variáveis estocásticas e

*D de 10%, com um conjunto de 13 pontos de colocação, foi alcançada a convergência

da solução.

t

0 0.2 0.4 0.6 0.8 110

-4

10-3

10-2

10-1

100

101

coloc 5

coloc 13

coloc 29

t

0 0.1 0.2 0.3 0.4 0.5 0.610

-6

10-5

10-4

10-3

10-2

10-1

coloc 5

coloc 13

coloc 29

Page 80: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

71

(a) (b)

Figura 5.16: Comparação do método de colocação com o método de Monte Carlo: (a)

Concentração. (b) Taxa de liberação.

A Figura 5.17(a)-(b) apresenta o gráfico com a média da solução estocástica

alcançada pelo método de colocação para 5, 13 e 29 pontos de colocação, para a

concentração de fertilizantes e para a taxa de liberação respectivamente, obtida por

(5.2)-(5.3) juntamente com o intervalo de credibilidade (IC) gerado pela variância

(3.20), associado a essa média. As variáveis estocásticas e *D foram admitidas com

50% ( 0.5 ) de dispersão. Verifica-se que o aumento da dispersão dos dados

estocásticos influencia a solução, o que é claramente notado no intervalo de

credibilidade em comparação à Figura 5.13 para 0.1 .

(a) (b)

Figura 5.17: (a) Concentração x tempo para 0.5 . (b) Taxa de liberação x raiz

quadrada do tempo para 0.5 .

t

Co

nce

ntr

açã

o

0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

MC 10000

MC 10000-IC

t

Ta

xa

0.2 0.4 0.6 0.8 1

0.2

0.4

0.6

0.8

1

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

MC 10000

MC 10000-IC

t

C

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

coloc 5

coloc 5-IC

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

t

Ta

xa

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

coloc 5

coloc 5-IC

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

Page 81: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

72

Para verificar os resultados obtidos pelo método de colocação foi utilizado o

método de Monte Carlo em que foi aplicado um conjunto de 1000, 3000, 5000 e 10000

amostras (Figura 5.18(a)-(b)). Aumentando o grau de incerteza, não é verificado

discrepâncias acentuadas no gráfico para a variância da solução para a concentração e

taxa de liberação obtida para vários conjuntos de amostras. O que se pode concluir que

para estes conjuntos de amostras houve convergência da solução média obtida por

Monte Carlo. Para comparação com o método de colocação, considerando incertezas

com coeficientes de dispersão de dados de 50%, será adotada a solução com o método

de Monte Carlo obtida com 10000 amostras.

(a) (b)

Figura 5.18: Convergência da variância com o método de Monte Carlo para 10000

amostras. (a) Concentração (b) Taxa de liberação.

Para o método de colocação utilizou-se 5, 13 e 29 pontos de colocação para

0.5 , resultados já vistos na Figura 5.17(a)-(b).

A Figura 5.19(a)-(b) mostra o desvio MC entre as soluções obtidas com alguns

conjuntos de pontos de colocação, comparadas com o método de Monte Carlo com

10000 amostras. Pode-se constatar uma diminuição no desvio para um conjunto de 29

pontos de colocação para a concentração. Por outro lado, não se observa um

comportamento característico do desvio MC para a taxa de liberação.

t

0.2 0.4 0.6 0.8 1

0

0.01

0.02

0.03

0.04

MC 1000

MC 3000

MC 5000

MC 10000

2

t

0.2 0.4 0.6 0.8 10

0.001

0.002

0.003

0.004

0.005

MC 1000

MC 3000

MC 5000

MC 10000

2

Page 82: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

73

(a) (b)

Figura 5.19: (a) Desvio entre soluções para a concentração em alguns valores de tempo.

(b) Desvio entre soluções para a taxa de liberação em alguns valores de tempo.

A Figura 5.20(a)-(b) apresenta a comparação entre os métodos utilizados em que

se pode verificar uma tendente convergência entre os métodos. Essa convergência para

o método de colocação foi alcançada com 29 pontos.

(a) (b)

Figura 5.20: Comparação do método de colocação com o método de Monte Carlo: (a)

Concentração. (b) Taxa de liberação.

Por último, as variáveis estocásticas e *D foram admitidas com 90%(

0.9 ) de dispersão. A Figura 5.21(a)-(b) apresenta o gráfico com a média da solução

estocástica alcançada pelo método de colocação para 5, 13 e 29 pontos de colocação,

para a concentração de fertilizantes e para a taxa de liberação, obtida por (5.2)-(5.3)

t

0 0.2 0.4 0.6 0.8 110

-4

10-3

10-2

10-1

100

coloc 5

coloc 13

coloc 29

t

0 0.1 0.2 0.3 0.4 0.5 0.610

-5

10-4

10-3

10-2

coloc 5

coloc 13

coloc 29

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

MC 10000

MC 10000-IC

t

Ta

xa

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

MC 10000

MC 10000-IC

Page 83: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

74

juntamente com o intervalo de credibilidade (IC) gerado pela variância (3.20), associado

a essa média. Podemos observar a influência da dispersão das variáveis estocásticas nos

intervalos de credibilidade, em que é visível a discrepância entre estes intervalos para

vários conjuntos de pontos de colocação.

(a) (b)

Figura 5.21: (a) Concentração x tempo para 0.9 . (b) Taxa de liberação x raiz

quadrada do tempo para 0.9 .

A verificação dos resultados obtidos pelo método de colocação foi feita pelo

método de Monte Carlo em que foi aplicado um conjunto de 1000, 3000, 5000 e 10000

amostras (Figura 5.22(a)-(b)). Estes resultados mostram a convergência da variância da

solução para a concentração e taxa de liberação. Com isso pode-se adotar um conjunto

de 10000 amostras para a comparação que se segue.

(a) (b)

Figura 5.22: Convergência da variância com o método de Monte Carlo para 10000

amostras. (a) Concentração (b) Taxa de liberação.

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

coloc 5

coloc 5-IC

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

t

Ta

xa

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

coloc 5

coloc 5-IC

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

t

0.2 0.4 0.6 0.8 1

0

0.02

0.04

0.06

0.08

0.1

MC 1000

MC 3000

MC 5000

MC 10000

2

t

0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

MC 1000

MC 3000

MC 5000

MC 10000

2

Page 84: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

75

Para o método de colocação foi aplicado 5, 13 e 29 pontos para 0.9 . A

Figura 5.23(a)-(b) exibe estes desvios. Pode-se verificar o seu declínio para 29 pontos

de colocação. Este declínio garante a convergência da solução obtida pelo método de

colocação com 29 pontos como pode ser vista na Figura 5.24(a)-(b).

(a) (b)

Figura 5.23: (a) Desvio entre soluções para a concentração em alguns valores de tempo.

(b) Desvio entre soluções para a taxa de liberação em alguns valores de tempo.

(a) (b)

Figura 5.24: Comparação do método de colocação com o método de Monte Carlo: (a)

Concentração. (b) Taxa de liberação.

Os resultados obtidos quando duas variáveis de entrada são aleatórias indicam

que as incertezas na solução final podem ser significativas.

t

0 0.2 0.4 0.6 0.8 110

-6

10-5

10-4

10-3

10-2

10-1

100

coloc 5

coloc 13

coloc 29

t

0 0.1 0.2 0.3 0.4 0.5 0.610

-5

10-4

10-3

10-2

10-1

coloc 5

coloc 13

coloc 29

t

Co

nce

ntr

açã

o

0.2 0.4 0.6 0.8 10

0.5

1

1.5

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

MC10000

MC 10000-IC

t

Ta

xa

0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

coloc 13

coloc 13-IC

coloc 29

coloc 29-IC

MC 10000

MC 10000-IC

Page 85: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

76

5.2 Problema de Magnetohidrodinâmica

Nesta seção serão mostrados resultados do problema estocástico de

magnetohidrodinâmica, apresentado nas equações (4.62)-(4.67). O problema

determinístico foi estudado por Colaço et al [45], onde os autores empregaram o

método sem malha baseado nas funções de base radial para resolver um problema de

magnetohidrodinâmica (MHD) bidimensional em uma cavidade quadrada com

escoamento incompressível, laminar, estacionário e com campo magnético constante. O

sistema de equações resultante foi resolvido utilizando um método quasi-Newton e seus

resultados foram comparados ao método dos Volumes Finitos.

No presente trabalho, o método sem malha baseado nas funções de base radial,

juntamente com o método de colocação estocástica, foi aplicado ao mesmo problema. A

aplicação do método de colocação estocástica se dá pelo emprego de incertezas em

alguns parâmetros de entrada.

O problema de magnetohidrodinâmica abordado em [45] será estudado na

referida tese, alocando incertezas no número de Grashof ( Gr ) e depois também no

número de Hartman ( Ha ). Para solucionar o problema, foram adotados os seguintes

valores de referência: 4

0Gr 10 e 0Ha 25 [45]. No caso de variável estocástica

unidimensional, para obter os pontos de colocação, foi testada a regra da quadratura de

Gauss e a regra de Clenshaw-Curtis [34]. Devido à necessidade de gerar muitos pontos

para alcançar boa precisão na solução do problema de magnetohidrodinâmica, com a

utilização da primeira regra não foi possível gerar pontos suficientes, pois a sua

composição necessitaria construir uma matriz Jacobiana (3.18) de coeficientes de

recorrência muito grande, o que dificultaria a sua aplicação. Neste caso, o uso da regra

de Clenshaw-Curtis gerou resultados mais rápidos e sem esta dificuldade de geração de

pontos. Tanto para variável estocástica unidimensional como para variável estocástica

bidimensional, a malha esparsa, gerada pela regra de Clenshaw-Curtis, será empregada

a fim de garantir boas soluções com menos pontos do que uma malha completa.

As soluções do presente trabalho por colocação estocástica, para uma e duas

variáveis estocásticas, foram obtidas para quatro níveis diferentes das abscissas de

Clenshaw-Curtis, em que pelo menos em um dos níveis utilizados foi verificada a

convergência da solução. Os níveis de interpolação utilizados em cada dimensão

estocástica são apresentados na Tabela 5.1. O sistema de equações diferenciais parciais

Page 86: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

77

obtido pelo método utilizado foi resolvido pelo método baseado nas funções de base

radial. As soluções foram obtidas para dois valores do coeficiente de dispersão (0.2 e

0.7) da variável aleatória a fim de avaliar a solução frente à dispersão dos dados de

entrada. Para verificação do método de colocação foi utilizado o método de Monte

Carlo.

Tabela 5.1: Número de pontos de colocação de acordo com o nível de interpolação e

dimensão estocástica.

Dimensão

estocástica

( )N

Nível de interpolação ( i ) Número de pontos de colocação ( PN )

1

7 129

8 257

9 513

10 1025

2

6 321

7 705

8 1537

9 3329

5.2.1 Variável estocástica Unidimensional – Número de Grashof ( Gr )

O parâmetro contendo incerteza é dado pela equação (4.71). A distribuição de

é uniforme no intervalo [ 1,1] e para o coeficiente de dispersão de dados ( ) foram

utilizados os seguintes valores: 0.2 e 0.7 .

Na solução do problema determinístico, para cada ponto de colocação

estocástica, foi aplicada a técnica das funções de base radial com 15x15 centros

uniformemente espaçados. O número de Prandlt é uma variável adimensional do

problema e seu valor foi tomado como 0.71. O número de Hartman foi escolhido como

zero, ou seja, a solução do problema será obtida sem considerar campo magnético.

Como este trabalho faz referência principalmente ao método de colocação estocástica, a

quantidade de centros empregados na expansão por RBF será mantida fixa a fim de

Page 87: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

78

investigar primordialmente os efeitos do método de colocação estocástica através de

cada conjunto de pontos de colocação (níveis de interpolação).

A fim de ponderar a dispersão da variável estocástica na solução do problema

de magnetohidrodinâmica serão apresentados os resultados para os dois graus de

dispersão ( 0.2 e 0.7 ) em conjunto para cada variável analisada.

5.2.1.1 Análise de convergência do método de Monte Carlo

Inicialmente, foi realizada uma análise de convergência do método de Monte

Carlo. Esta análise foi feita através da comparação da média e da variância do número

de Nusselt nas paredes verticais da cavidade, por exprimir a taxa de transferência de

calor.

A Figura 5.25(a) apresenta a média da solução do número de Nusselt em 0x ,

ao longo da parede vertical, obtida por simulações de Monte Carlo com vários

conjuntos de amostras para um coeficiente de dispersão da variável estocástica de 20%.

Da mesma forma, a Figura 5.25(b) apresenta resultados análogos para um coeficiente de

dispersão da variável estocástica de 70%. As Figuras 5.25(a)-(b) não aparentam

discrepâncias da média da solução do número de Nusselt para o número de amostras

utilizadas, mas o aumento da dispersão da variável estocástica gera uma leve alteração

na média do número de Nusselt, como é possível verificar para 0.2 , Nu0=3.23 e

para 0.7 , Nu0=3.17 em 0y .

Page 88: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

79

(a) (b)

Figura 5.25: Média do número de Nusselt em 0x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

A Figura 5.26(a) apresenta a média da solução do número de Nusselt em 1x ,

ao longo da parede vertical, obtida por simulações de Monte Carlo com vários

conjuntos de amostras para um coeficiente de dispersão da variável estocástica de 20%

e a Figura 5.26(b) para um coeficiente de dispersão da variável estocástica de 70%.

Novamente, em conformidade com o número de Nusselt em 0x , a média da solução

para o número de Nusselt não é modificada para o conjunto de amostras usadas, mas

uma leve discrepância (em 1y ) pode ser verificada quando se aumenta o grau de

incerteza do parâmetro estocástico.

(a) (b)

Figura 5.26: Média do número de Nusselt em 1x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

+ + + + + ++

++

++

++

++

++

++

++

++

++

++

+ + +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

+ + + ++

++

++

+

+

++

++

+

+

++

++

+ +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

+ + + ++

++

++

++

++

++

++

++

++

++

++

+ + + + +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

+ + ++

++

+

+

+

++

++

+

++

++

++

+ + +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

Page 89: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

80

A Figura 5.27(a) apresenta a variância do número de Nusselt em 0x , ao longo

do eixo vertical, obtida por simulações de Monte Carlo com vários conjuntos de

amostras para um coeficiente de dispersão da variável estocástica de 20% e a Figura

5.27(b) para um coeficiente de dispersão da variável estocástica de 70%. Para a

variância do número de Nusselt é possível verificar a alteração de acordo com o número

de amostras. Também com o aumento do grau de incerteza no número de Grashof é

visível a amplitude da variância. Na Figura 5.27 verifica-se o intervalo onde ela está

definida, com o coeficiente de dispersão de 20%, 20 0.02 e com o coeficiente de

dispersão de 70%, 20 0.25 .

(a) (b)

Figura 5.27: Variância do número de Nusselt em 0x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

A Figura 5.28(a) apresenta a variância do número de Nusselt em 1x , ao longo

do eixo vertical, obtida por simulações de Monte Carlo com vários conjuntos de

amostras para um coeficiente de dispersão da variável estocástica de 20% e a Figura

5.28(b) para um coeficiente de dispersão da variável estocástica de 70%. Do mesmo

modo, as características percebidas para variância do número de Nusselt em 0x são

encontradas em 1x .

++

++

++

++

++

++

++

++

++

++

+ + + +

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

++

++

++

++

++

++

++

++

+

++

++

++ + +

y

0 0.2 0.4 0.6 0.8 10

0.05

0.1

0.15

0.2

0.25

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

Page 90: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

81

(a) (b)

Figura 5.28: Variância do número de Nusselt em 1x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

Nas Figuras 5.27(b) e 5.28(b) percebe-se uma discrepância da variância obtida

com 600 amostras em relação aos outros conjuntos de amostras, indicando a

necessidade de maiores amostras para se realizar a convergência desta medida.

A fim de verificar com mais clareza o número de amostras que define uma

convergência aceitável, foi realizada uma análise dada pela diferença ( ) entre a

variância do número de Nusselt nas paredes verticais da cavidade obtida pelo método de

Monte Carlo para o conjunto de amostras próximas, por exemplo, a diferença entre a

variância obtida com 1200 amostras e 600 amostras ( 1200 600 ). A Figura 5.29 mostra

esta diferença para o coeficiente de dispersão de dados de 20% e a Figura 5.30 mostra

esta diferença para o coeficiente de dispersão de dados de 70%.

+ + + ++

++

++

++

+

+

+

++

++

++

++

+

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

+ + + ++

++

+

+

++

++

++

++

++

+

++

++

+

y

0 0.2 0.4 0.6 0.8 10

0.05

0.1

0.15

0.2

0.25

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

Page 91: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

82

(a) (b)

Figura 5.29: Diferença entre a variância do número de Nusselt em 0x (a) e em 1x

(b) ao longo do domínio, obtida por Monte Carlo para coeficiente de dispersão de 20%.

(a) (b)

Figura 5.30: Diferença entre a variância do número de Nusselt em 0x (a) e em 1x

(b) ao longo do domínio, obtida por Monte Carlo para coeficiente de dispersão de 70%.

Pode-se concluir que, apresentando pequenas variações nas amostras utilizadas,

pode-se considerar um conjunto de 7200 amostras do método de Monte Carlo, tanto

para o coeficiente de dispersão de 20% quanto para o coeficiente de dispersão de 70%,

como convergida.

y

0 0.2 0.4 0.6 0.8 1

10-6

10-5

10-4

10-3

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

y

0 0.2 0.4 0.6 0.8 1

10-7

10-6

10-5

10-4

10-3

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

y

0 0.2 0.4 0.6 0.8 1

10-5

10-4

10-3

10-2

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

y

0 0.2 0.4 0.6 0.8 1

10-7

10-6

10-5

10-4

10-3

10-2

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

Page 92: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

83

5.2.1.2 Análise de convergência do método de Colocação Estocástica

Com os resultados obtidos pelo método de Monte Carlo, pode-se verificar como

se comporta o problema de MHD quando se aplica o método de colocação para

diferentes níveis de interpolação. Assim como no método de Monte Carlo, foi analisada

a média e variância do número de Nusselt nas paredes verticais da cavidade, com o

método de colocação estocástica. Inicialmente será apresentada a média da solução do

número de Nusselt para vários níveis de interpolação, juntamente com a solução obtida

pelo método de Monte Carlo com o número de amostras convergentes.

A Figura 5.31(a) mostra a média da solução do número de Nusselt em 0x , ao

longo da parede vertical, obtida pelo método de colocação em vários níveis de

interpolação e obtida pelo método de Monte Carlo para um coeficiente de dispersão da

variável estocástica de 20% e a Figura 5.31(b) para um coeficiente de dispersão da

variável estocástica de 70%. Para o coeficiente de dispersão de 20%, a média da solução

do número de Nusselt com nível de interpolação 7 apresenta diferenças observáveis em

comparação com outros níveis, indicando a necessidade de um nível de interpolação

maior. Em níveis de interpolação superiores, o resultado para a média do número de

Nusselt se mostrou convergente nos dois tipos de graus de incerteza.

(a) (b)

Figura 5.31: Média do número de Nusselt em 0x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

Da mesma forma, a Figura 5.32(a)-(b) apresenta a média da solução para o

número de Nusselt em 1x com 20% e 70% de dispersão da variável estocástica.

+ + + + + ++

++

++

++

++

++

++

++

++

++

++

+ + +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

Nível 7

Nível 8

Nível 9

Nível 10

MC 7200+

+ + + + + ++

++

++

++

++

++

++

++

++

++

++

+ + +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

Nível 7

Nível 8

Nível 9

Nível 10

MC 7200+

Page 93: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

84

Novamente, com o nível de interpolação 7, não foi possível obter um resultado

satisfatório com 0.2 . No entanto, para 0.7 não foi percebido discrepâncias da

solução entre os níveis.

(a) (b)

Figura 5.32: Média do número de Nusselt em 1x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

A Figura 5.33 apresenta a variância do número de Nusselt em 0x , ao longo

do eixo vertical, obtida pelo método de colocação para vários níveis de interpolação em

comparação com o método de Monte Carlo, para coeficiente de dispersão da variável

estocástica de 20% (a) e para coeficiente de dispersão da variável estocástica de 70%

(b). Como já detectada pela média do número de Nusselt, o nível de interpolação 7

gerou resultados não convergidos, o que também foi detectado pela variância, de forma

tão significativa, que não foi possível inseri-la na figura 5.33(a). Pela variância do

número de Nusselt é possível identificar sua magnitude crescente quando se aumenta o

grau de incerteza, ou seja, para grau de incerteza de 20%, 20 0.02 e para grau de

incerteza de 70%, 20 0.25 .

+ + + ++

++

++

++

++

++

++

++

++

++

++

+ + + + +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

Nível 7

Nível 8

Nível 9

Nível 10

MC 7200+

+ + + ++

++

++

++

++

++

++

++

++

++

++

+ + + + +

y

0 0.2 0.4 0.6 0.8 1

3.5

3

2.5

2

1.5

1

0.5

Nível 7

Nível 8

Nível 9

Nível 10

MC 7200+

Page 94: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

85

(a) (b)

Figura 5.33: Variância do número de Nusselt em 0x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

A Figura 5.34 apresenta a variância do número de Nusselt em 1x , ao longo do

eixo vertical, obtida pelo método de colocação para vários níveis de interpolação em

comparação com o método de Monte Carlo, para coeficiente de dispersão da variável

estocástica de 20% (a) e para coeficiente de dispersão da variável estocástica de 70%

(b).

(a) (b)

Figura 5.34: Variância do número de Nusselt em 1x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

A tabela 5.2 apresenta o número de Nusselt médio em 0x para os dois graus

de incertezas aplicadas obtidos pelo método de colocação e pelo método de Monte

++++

++

++

++

++

++

++

++

++

++

++

++

+ + + + + +

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

Nível 8

Nível 9

Nível 10

MC 7200+

2+ + + +

++

++

++

++

++

++

++

++

++

++

++ + + + +

y

0 0.2 0.4 0.6 0.8 10

0.05

0.1

0.15

0.2

0.25

Nível 7

Nível 8

Nível 9

Nível 10

MC 7200+

2

+ + + ++

++

++

++

+

+

+

++

++

++

++

+

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

Nível 8

Nível 9

Nível 10

MC 7200+

2

+ + + + + ++

++

++

++

++

++

++

++

++

++

++

+ + +

y

0 0.2 0.4 0.6 0.8 10

0.05

0.1

0.15

0.2

0.25

Nível 7

Nível 8

Nível 9

Nível 10

MC 7200+

2

Page 95: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

86

Carlo. Percebe-se que a aplicação de 20% de dispersão na variável estocástica não

implica em mudanças aparentes no número de Nusselt médio, pois ele preserva o valor

do problema determinístico. Verifica-se que para 0.2 com o nível 8 de interpolação

o número de Nusselt médio foi capturado.

Tabela 5.2: Número de Nusselt médio em 0x

Nusselt Médio em 0x 410Gr Nível 7 Nível 8 Nível 9 Nível 10 Monte Carlo

0Ha

20% 1.94 2.01 2.01 2.01 2.01

70% 1.94 1.94 1.94 1.94 1.97

Determinístico 2.01

Assim como na Figura 5.33(a), a Figura 5.34(a) indica resultados convergentes

com nível de interpolação 8. Através da Figura 5.33(b) e da Figura 5.34(b) pode-se

verificar resultados bem próximos em todos os níveis de interpolação e pode-se admitir

a convergência a partir do nível 7.

5.2.1.3 Perfil da velocidade e temperatura

A Figura 5.35(a) mostra o perfil médio da velocidade para os pontos localizados

ao longo da direção vertical para 0.5x , com o coeficiente de dispersão da variável

estocástica de 20% e a Figura 5.35(b) mostra o mesmo perfil da velocidade para

coeficiente de dispersão da variável estocástica de 70%. Estas figuras apresentam

também o intervalo de credibilidade da solução de acordo com o nível de interpolação

utilizado no método de colocação. Percebe-se que com o aumento do nível de

interpolação ocorre a diminuição do intervalo de credibilidade. Também, à medida que

aumenta a dispersão da variável estocástica, este intervalo de credibilidade, no caso da

velocidade, foi levemente ampliado para um nível de interpolação convergente (Figura

5.35(c)-(d)).

Page 96: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

87

(a) (b)

(c) (d)

Figura 5.35: Perfil da velocidade longo do eixo vertical para 0.5x . (a) 0.2 , (b)

0.7 , (c) Ampliação de 0.2 e (d) Ampliação de 0.7 .

A Figura 5.36(a) apresenta a média da temperatura para os pontos localizados ao

longo da direção horizontal quando 0.5y , para coeficiente de dispersão da variável

estocástica de 20% e a Figura 5.36(b) mostra a média da temperatura para coeficiente de

dispersão da variável estocástica de 70%. Também são apresentados os intervalos de

credibilidade associados com as soluções obtidas com diferentes números de pontos de

colocação. Nas Figuras 5.36(a)-(b) a dispersão da variável estocástica é evidente,

mostando que mesmo em um nível de interpolação convergente, ela se apresenta com

maior amplitude.

++

+ + + + + + + + + + + + + + + + + + + + + ++

+

y

U

0 0.2 0.4 0.6 0.8 1-1.5

-1

-0.5

0

0.5

1

1.5

2Nível 7Nível 7 - ICNível 8Nível 8 - ICNível 9Nível 9 - ICNível 10Nível 10 - ICMC 7200+

++

+ + + + + + + + + + + + + + + + + + + + + + ++

y

U

0 0.2 0.4 0.6 0.8 1-1.5

-1

-0.5

0

0.5

1

1.5

2

Nível 7Nível 7 - ICNível 8Nível 8 - ICNível 9Nível 9 - ICNível 10Nível 10 - ICMC 7200+

+ + + + + + + + + + + + + + + + + + + + + + + + +

y

U

0.2 0.4

-0.5

0

0.5

Nível 7

Nível 7 - IC

Nível 8

Nível 8 - IC

Nível 9

Nível 9 - IC

Nível 10

Nível 10 - IC

MC 7200+

+ + + + + + + + + + + + + + + + + + + + + + + + +

y

U

0.2 0.4

-1

-0.5

0

0.5

Nível 7Nível 7 - ICNível 8Nível 8 - ICNível 9Nível 9 - ICNível 10Nível 10 - ICMC 7200+

Page 97: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

88

(a) (b)

Figura 5.36: Perfil da temperatura longo do eixo horizontal para 0.5y . (a) 0.2 e

(b) 0.7

Para 0.2 pode-se verificar a convergência para o nível de interpolação 9,

considerando soluções convergentes para o par velocidade e temperatura. Os resultados

para 0.7 são análogos. Para o perfil da velocidade é verificado a diminuição do

intervalo de credibilidade com o aumento do nível de interpolação ao contrário do perfil

de temperatura o qual não revela esta característica de convergência.

5.2.1.4 Linhas de corrente e isotermas

A Figura 5.37 apresenta as linhas de corrente médias obtidas pelo método de

colocação para alguns níveis de interpolação em comparação com o método de Monte

Carlo. As linhas de corrente obtidas para cada nível de interpolação aproximaram-se do

resultado obtido pelo método de Monte Carlo tanto para 20% como para 70% de

dispersão da variável estocástica. A Figura 5.38 mostra as isotermas médias obtidas

pelo método de colocação para alguns níveis de interpolação em comparação com o

método de Monte Carlo. As isotermas obtidas para cada nível de interpolação mostram

boa comparação com o método de Monte Carlo tanto para 20% como para 70% de

dispersão da variável estocástica. Para efeito de comparação, as soluções

determinísticas, obtidas com as simulações computacionais sem incertezas no

parâmetro, também são incluídas nas Figuras 5.37 e 5.38.

+

+

+

+

++

++

+

++

++

+ + + + + + + + ++

++

+

+

++

++

+

+

+

x

T

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1Nível 7

Nível 7 - IC

Nível 8

Nível 8 - IC

Nível 9

Nível 9 - IC

Nível 10

Nível 10 - IC

MC 7200+

+

+

+

++

++

+

++

++

++ + + + + + + + +

++

++

+

++

++

+

+

+

x

T

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

Nível 7

Nível 7 - IC

Nível 8

Nível 8 - IC

Nível 9

Nível 9 - IC

Nível 10

Nível 10 - IC

MC 7200+

Page 98: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

89

Linhas de Corrente

0Ha 20% 70% M

on

te C

arlo

7200

Nív

el 7

Nív

el 8

Nív

el 9

Nív

el 1

0

Det

erm

inís

tico

Figura 5.37: Linhas de corrente para 20% e 70% de dispersão na variável estocástica

obtidas para diferentes níveis de interpolação do método de colocação e pelo método de

Monte Carlo.

Page 99: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

90

Isotermas

0Ha 20% 70% M

on

te C

arlo

7200

Nív

el 7

Nív

el 8

Nív

el 9

Nív

el 1

0

Det

erm

inís

tico

Figura 5.38: Isolinhas para 20% e 70% de dispersão na variável estocástica obtidas para

diferentes níveis de interpolação do método de colocação e pelo método de Monte

Carlo.

Page 100: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

91

Tanto para as linhas de corrente quanto para as isotermas não foram verificadas

diferenças entre os níveis de interpolação, mesmo para diferentes graus de dispersão da

variável estocástica. A incidência de incerteza no número de Grashof indica a

modificação do vórtice central expandindo-o (Figura 5.37) no caso das linhas de

corrente. Já para as isotermas foi identificada uma boa concordância dos resultados

médios e da simulação determinística.

Levando em consideração toda a análise feita para as variáveis do problema de

MHD para 1 parâmetro contendo incerteza, pode-se concluir que foi atingida a

convergência da solução para o nível de interpolação 9 tanto para 20% como para 70%

de incerteza.

O crescimento do grau de incerteza da variável estocástica não acarretou, de

maneira geral, no aumento do nível de interpolação. No caso desta análise, o método de

colocação estocástica mostrou acurácia equivalente ao método de Monte Carlo com

bem menos realizações do problema direto.

5.2.2 Variável estocástica Bidimensional – Número de Grashof ( Gr ) e

Número de Hartman ( Ha )

Agora os parâmetros contendo incertezas são os números de Grashof e de

Hartman, dados pelas equações (4.83) e (4.84). A distribuição de 1 e 2 são uniformes

no intervalo [ 1,1] e para os coeficientes de dispersão de dados ( Gr e Ha ) foram

utilizados os seguintes valores: 0.2Gr com 0.2Ha e 0.7Gr com 0.7Ha .

Considera-se novamente na solução do problema determinístico, para cada

ponto de colocação estocástica, 15 15 (centros) uniformemente espaçados no domínio

espacial e o número de Prandlt igual 0.71.

5.2.2.1 Análise de convergência do método de Monte Carlo

Também para o caso de variável estocástica bidimensional, foi realizada uma

análise de convergência do Método de Monte Carlo, utilizando-se os números de

Nusselt nas paredes verticais da cavidade.

A Figura 5.39(a) apresenta a média da solução do número de Nusselt em 0x ,

obtida por simulações de Monte Carlo com vários números de amostras, para um

Page 101: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

92

coeficiente de dispersão das variáveis estocásticas de 20%. Da mesma forma, a Figura

5.39(b) apresenta a média da solução do número de Nusselt em 0x , obtida por

simulações de Monte Carlo com vários conjuntos de amostras para um coeficiente de

dispersão das variáveis estocásticas de 70%.

A Figura 5.39(a) não aparenta discrepância da média da solução do número de

Nusselt para o número de amostras analisadas, mas o aumento da dispersão das

variáveis estocásticas (Figura 5.39(b)) mostra a necessidade de mais amostras, como

pode-se perceber em 0 0.3y e também gera uma leve alteração na média da

solução do número de Nusselt, como é possível verificar em 0y .

(a) (b)

Figura 5.39: Média do número de Nusselt em 0x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

A Figura 5.40(a) apresenta a média da solução do número de Nusselt em 1x ,

ao longo do eixo vertical, obtida por simulações de Monte Carlo com vários conjuntos

de amostras para um coeficiente de dispersão das variáveis estocásticas de 20% e a

Figura 5.40(b) para um coeficiente de dispersão das variáveis estocásticas de 70%.

Novamente, em conformidade com o número de Nusselt em 0x , percebe-se que a

medida que se aumenta o grau da incerteza nos parâmetros estocásticos ocorre à

necessidade de mais amostras. Também, o intervalo em que a média é definida se

amplia com o aumento da incerteza aplicada ao problema.

+ ++

++

+

++

++

++

+

+

+

++

++

++

++ +

y

0 0.2 0.4 0.6 0.8 1

1.8

1.6

1.4

1.2

1

0.8

0.6

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

+ + + + ++

+++++++++++++++++++++++++++++++

++ + + + + +

y

0 0.2 0.4 0.6 0.8 1

2

1.8

1.6

1.4

1.2

1

0.8

0.6

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

Page 102: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

93

(a) (b)

Figura 5.40: Média do número de Nusselt em 1x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

A Figura 5.41(a) apresenta a variância do número de Nusselt em 0x , obtida

por simulações de Monte Carlo com vários conjuntos de amostras para um coeficiente

de dispersão das variáveis estocásticas de 20% e a Figura 5.41(b) para um coeficiente de

dispersão das variáveis estocásticas de 70%. Para variância do número de Nusselt é

possível verificar a alteração de acordo com o número de amostras. Também com o

aumento do grau de incerteza no número de Grashof e Hartman é visível a amplitude da

variância. Na Figura 5.41 verifica-se o intervalo onde ela está definida, com o

coeficiente de dispersão de 20%, 20 0.03 e com o coeficiente de dispersão de

70%, 20 0.42 .

+ + ++

++

++

++

+

+

+

+

+

+

+

++

++

+ +

y

0 0.2 0.4 0.6 0.8 1

1.8

1.6

1.4

1.2

1

0.8

0.6

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

+ + + + + ++

++++++++++++++++++++++++++++++++

+ + + + +

y

0 0.2 0.4 0.6 0.8 1

2

1.8

1.6

1.4

1.2

1

0.8

0.6

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

Page 103: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

94

(a) (b)

Figura 5.41: Variância do número de Nusselt em 0x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

A Figura 5.42(a) apresenta a variância do número de Nusselt em 1x , obtida

por simulações de Monte Carlo com vários conjuntos de amostras para um coeficiente

de dispersão das variáveis estocásticas de 20% e a Figura 5.42(b) para um coeficiente de

dispersão das variáveis estocásticas de 70%. Do mesmo modo, as características

percebidas para variância do número de Nusselt em 0x são encontradas em 1x .

(a) (b)

Figura 5.42: Variância do número de Nusselt em 1x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

Nas Figuras 5.42(b) e 5.43(b) percebe-se uma significativa discrepância da

variância obtida com 600 amostras em relação aos outros conjuntos de amostras,

++

+

+

+

+

+++++++

+

+

++

++

+ + + + + + +

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

0.025

0.03

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2+ + + +

+++++++++++++++++++++++++

+ + + + + + + + + + + + +

y

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

+ + + + + + ++

++

++

+

+

+

+++++

+

+

++

++

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

0.025

0.03MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

+ + + + + + + + + + + + + ++++++++++++++++++++++++++ + + +

y

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

Page 104: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

95

indicando a necessidade de maior número de amostras para se realizar a convergência

desta variável.

(a) (b)

Figura 5.43: Variância do número de Nusselt em 1x ao longo do domínio, obtida por

Monte Carlo. (a) 0.2 e (b) 0.7

A fim de avaliar a variação da variância entre amostras do método de Monte

Carlo, novamente será realizada uma análise dada pela diferença ( ) entre a variância

do número de Nusselt nas paredes verticais da cavidade. A Figura 5.44 mostra esta

diferença para o coeficiente de dispersão de dados de 20% e a Figura 5.45 mostra esta

diferença para o coeficiente de dispersão de dados de 70%.

(a) (b)

Figura 5.44: Diferença entre a variância do número de Nusselt em 0x (a) e em 1x

(b) ao longo do domínio, obtida por Monte Carlo para coeficiente de dispersão de 20%.

+ + + + + + ++

++

++

+

+

+

+++++

+

+

++

++

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

0.025

0.03MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

+ + + + + + + + + + + + + ++++++++++++++++++++++++++ + + +

y

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

MC 600

MC 1200

MC 3600

MC 7200

MC 8400+

2

y

0 0.2 0.4 0.6 0.8 110

-8

10-7

10-6

10-5

10-4

10-3

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

y

0 0.2 0.4 0.6 0.8 1

10-8

10-7

10-6

10-5

10-4

10-3

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

Page 105: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

96

(a) (b)

Figura 5.45: Diferença entre a variância do número de Nusselt em 0x (a) e em 1x

(b) ao longo do domínio, obtida por Monte Carlo para coeficiente de dispersão de 70%.

A Figura 5.44 apresenta uma diferença entre a variância do número de Nusselt

obtida por 8400 amostras e a variância do número de Nusselt obtida por 7200 amostras

menores que 410, para grau de incerteza de 20%. Já a Figura 5.45 apresenta uma

diferença entre a variância do número de Nusselt obtida por 8400 amostras e a variância

do número de Nusselt obtida por 7200 amostras em torno de 410 para grau de incerteza

de 70%. Assim, também para o caso de variável estocástica bidimensional, pode-se

considerar a solução com um conjunto de 7200 amostras do método de Monte Carlo,

tanto para o coeficiente de dispersão de 20% quanto para o coeficiente de dispersão de

70%, como convergida.

5.2.2.2 Análise de convergência do método de Colocação Estocástica

O método de colocação estocástica será verificado utilizando-se a solução

convergida do método de Monte Carlo. Foi analisada a média e também a variância do

número de Nusselt nas paredes verticais da cavidade, obtidas com o método de

colocação estocástica para alguns níveis de interpolação.

A Figura 5.46(a) mostra a média da solução do número de Nusselt em 0x ,

obtida pelo método de colocação em vários níveis de interpolação e obtida pelo método

de Monte Carlo para um coeficiente de dispersão das variáveis estocásticas de 20% e, a

Figura 5.46(b) mostra resultados análogos para um coeficiente de dispersão das

y

0 0.2 0.4 0.6 0.8 1

10-6

10-5

10-4

10-3

10-2

10-1

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

y

0 0.2 0.4 0.6 0.8 1

10-6

10-5

10-4

10-3

10-2

10-1

1200 - 600

3600 - 1200

7200 - 3600

8400 - 7200

Page 106: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

97

variáveis estocásticas de 70%. Para os dois graus de incerteza nos parâmetros não foram

identificadas discrepâncias das soluções entre os diferentes níveis de interpolação,

apenas foi verificado que a amplitude da média do número de Nusselt.

(a) (b)

Figura 5.46: Média do número de Nusselt em 0x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

Da mesma forma, a Figura 5.47(a)-(b) apresenta a média da solução para o

número de Nusselt em 1x com 20% e 70% de dispersão das variáveis estocásticas.

Outra vez, quando aumenta-se o coeficiente de dispersão, foi identificado o aumento do

intervalo da média.

(a) (b)

Figura 5.47: Média do número de Nusselt em 1x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

+ + + + +++

+++++++++++++++++++++++++++++++

+ + + + + + +

y

0 0.2 0.4 0.6 0.8 1

1.8

1.6

1.4

1.2

1

0.8

0.6

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

+ + + + ++

+++++++++++++++++++++++++++++++

++ + + + + +

y

0 0.2 0.4 0.6 0.8 1

2

1.8

1.6

1.4

1.2

1

0.8

0.6

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

+ + + + + ++

+++++++++++++++++++++++++++++++++

+ + + + +

y

0 0.2 0.4 0.6 0.8 1

1.8

1.6

1.4

1.2

1

0.8

0.6

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

+ + ++

++

++

++

++

++

++

++

++

++

++

++

+ + + +

y

0 0.2 0.4 0.6 0.8 1

2

1.8

1.6

1.4

1.2

1

0.8

0.6

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

Page 107: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

98

A Figura 5.48 apresenta a variância do número de Nusselt em 0x , obtida pelo

método de colocação para vários níveis de interpolação, em comparação com o método

de Monte Carlo, para coeficientes de dispersão das variáveis estocásticas de 20% (a) e

para coeficiente de dispersão das variáveis estocásticas de 70% (b). A Figura 5.48(b) já

detecta a necessidade de maiores níveis de interpolação, com o aumento do grau de

incerteza. Da mesma forma que para a média, amplia-se o intervalo da variância: com

20% de dispersão das variáveis estocásticas 20 0.03 e com 70% de dispersão das

variáveis estocásticas 20 0.4 , sinalizando alterações quando inclui ainda mais

incertezas.

(a) (b)

Figura 5.48: Variância do número de Nusselt em 0x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

A Figura 5.49 apresenta a variância do número de Nusselt em 1x , obtida pelo

método de colocação para vários níveis de interpolação em comparação com o método

de Monte Carlo, para coeficiente de dispersão das variáveis estocásticas de 20% (a) e

para coeficiente de dispersão das variáveis estocásticas de 70% (b). As Figuras 5.49(a)-

(b) também apresentam comportamentos próximos a variância do número de Nusselt

em 0x , no que diz respeito ao grau de incerteza.

+ ++

+++++++++++++++++

++

++

+ + + + + + + + + + +

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

0.025

0.03

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

2

+ + ++

+++++++++++++++++

++

++

+ + + + + + + + + +

y

0 0.2 0.4 0.6 0.8 1

0.1

0.2

0.3

0.4

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

2

Page 108: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

99

(a) (b)

Figura 5.49: Variância do número de Nusselt em 1x ao longo do domínio, obtida por

colocação. (a) 0.2 e (b) 0.7

A tabela 5.3 apresenta o número de Nusselt médio em 0x para os dois graus

de incertezas, obtidos pelo método de colocação e pelo método de Monte Carlo. O

resultado para o caso determinístico também é incluído nesta tabela. Percebe-se,

novamente, que a aplicação de 20% de dispersão nas variáveis estocásticas não implica

mudanças aparentes no número de Nusselt médio, pois ele mantém o valor do problema

sem incerteza.

Tabela 5.3: Número de Nusselt médio em 0x

Nusselt Médio em 0x 410Gr Nível 6 Nível 7 Nível 8 Nível 9 Monte Carlo

25Ha

20% 1.17 1.17 1.17 1.17 1.17

70% 1.28 1.28 1.28 1.28 1.25

Determinístico 1.17

Pelas Figuras analisadas nesta seção, pode-se verificar a convergência para o

coeficiente de dispersão de dados de 20% com o nível de interpolação 8 e para o

coeficiente de dispersão de 70% com o nível de interpolação 9. Mas, a fim de examinar

outras grandezas do problema, como por exemplo, o perfil de velocidade, será

apresentado nas seguintes seções como se apresentam todos os níveis de interpolação

utilizados neste caso de variável estocástica bidimensional.

+ + + + + + + + + + ++

++

+++++++++++++++++

++

+ +

y

0 0.2 0.4 0.6 0.8 10

0.005

0.01

0.015

0.02

0.025

0.03

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

2

+ + + + + + + + + ++

++

+++++++++++++++++

++

++ +

y

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

Nível 6

Nível 7

Nível 8

Nível 9

MC 7200+

2

Page 109: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

100

5.2.2.3 Perfil da velocidade e temperatura

A Figura 5.50(a) mostra o perfil da velocidade para os pontos localizados ao

longo do eixo vertical quando 0.5x para coeficiente de dispersão da variável

estocástica de 20% e a Figura 5.50(b) mostra o perfil da velocidade para os pontos

localizados ao longo do eixo vertical quando 0.5x para coeficiente de dispersão da

variável estocástica de 70%. Estas figuras também apresentam o intervalo de

credibilidade da solução de acordo com o nível de interpolação utilizado no método de

colocação. Percebe-se, novamente, que com o aumento do nível de interpolação ocorre

a diminuição do intervalo de credibilidade. Também, à medida que aumenta a dispersão

da variável estocástica, este intervalo de credibilidade, no caso da velocidade, foi

levemente ampliado para um nível de interpolação convergente.

(a) (b)

Figura 5.50: Perfil da velocidade longo do eixo vertical para 0.5x . (a) 0.2 e (b)

0.7

A Figura 5.51(a) apresenta a temperatura para os pontos localizados ao longo do

eixo horizontal quando 0.5y para coeficiente de dispersão da variável estocástica de

20% e a Figura 5.51(b) mostra a temperatura para os pontos localizados ao longo do

eixo horizontal quando 0.5y para coeficiente de dispersão da variável estocástica de

70%. Nas Figuras 5.51(a)-(b) a dispersão da variável estocástica é evidente, mostando

que mesmo em um nível de interpolação convergente, ela se apresenta com maior

amplitude para maiores níveis de incerteza nos números de Grashof e Hartman.

++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+

y

U

20 40 60 80 100-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8Nível 6

Nível 6 - IC

Nível 7

Nível 7 - IC

Nível 8

Nível 8 - IC

Nível 9

Nível 8 - IC

MC 7200+

++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+

y

U

20 40 60 80 100-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8Nível 6

Nível 6 - IC

Nível 7

Nível 7 - IC

Nível 8

Nível 8 - IC

Nível 9

Nível 9 - IC

MC 7200+

Page 110: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

101

(a) (b)

Figura 5.51: Perfil da temperatura longo do eixo horizontal para 0.5y . (a) 0.2 e

(b) 0.7

Para 0.2 pode-se verificar a convergência para o nível de interpolação 8,

considerando soluções convergentes para o par velocidade e temperatura, assim como

para 0.7 . Para o perfil da velocidade é verificado que em todos os níveis de

interpolação utilizados foi chegado à convergência, visto que o intervalo de

credibilidade não foi modificado.

5.2.2.4 Linhas de corrente e isotermas

A Figura 5.52 apresenta as médias das linhas de corrente obtidas pelo método de

colocação para alguns níveis de interpolação em comparação com o método de Monte

Carlo. As linhas de corrente obtidas para cada nível de interpolação aproximou-se do

resultado obtido pelo método de Monte Carlo tanto para 20% como para 70% de

dispersão das variáveis estocásticas.

A Figura 5.53 mostra as médias das isotermas obtidas pelo método de colocação

para alguns níveis de interpolação em comparação com o método de Monte Carlo e com

a solução determinística.

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

+

x

T

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1Nível 6

Nível 6 - IC

Nível 7

Nível 7 - IC

Nível 8

Nível 8 - IC

Nível 9

Nível 9 - IC

MC 7200+

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

++

x

T

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

Nível 6

Nível 6 - IC

Nível 7

Nível 7 - IC

Nível 8

Nível 8 - IC

Nível 9

Nível 9 - IC

MC 7200+

Page 111: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

102

Linhas de Corrente

25Ha 20% 70% M

on

te C

arlo

7200

Nív

el 6

Nív

el 7

Nív

el 8

Nív

el 9

Det

erm

inís

tico

Figura 5.52: Linhas de corrente para 20% e 70% de dispersão nas variáveis estocásticas

obtidas para diferentes níveis de interpolação do método de colocação e pelo método de

Monte Carlo.

Page 112: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

103

Isotermas

25Ha 20% 70% M

on

te C

arlo

7200

Nív

el 6

Nív

el 7

Nív

el 8

Nív

el 9

Det

erm

inís

tico

Figura 5.53: Isolinhas para 20% e 70% de dispersão na variável estocástica obtidas para

diferentes níveis de interpolação do método de colocação e pelo método de Monte

Carlo.

Page 113: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

104

As isotermas obtidas para cada nível de interpolação mostraram-se satisfatórias

em comparação com o método de Monte Carlo tanto para 20% como para 70% de

dispersão das variáveis estocásticas. Para ambas, linhas de corrente e isotermas, não foi

verificado diferenças entre cada nível de interpolação mesmo para diferentes graus de

dispersão da variável estocástica. A presença de campo magnético detecta nas linhas de

corrente uma modificação no vórtice central e diminui a distorção das isolinhas de

temperatura. A incidência de incerteza no Grashof e Hartman indica a modificação do

vórtice central expandindo-o, no caso das linhas de corrente. Novamente, para as

isotermas não foram verificadas diferenças significativas em relação ao caso

determinístico.

Analisando o problema de MHD para 2 parâmetros contendo incertezas, pode-se

considerar que aumentando o grau de dispersão das variáveis estocásticas necessitou de

um nível maior de interpolação (pontos de colocação). Em comparação com a variável

estocástica unidimensional, o acúmulo de incertezas (bidimensional) gera resultados

com maior variância.

Para o método de Monte Carlo, mesmo aplicando um maior grau de incerteza

em um parâmetro de entrada, ou em dois parâmetros de entrada não é verificada a

necessidade de mais amostras, no caso de convergência. Para o método de colocação

estocástica, em geral, isso também foi constatado, no caso do problema de MHD. O

método de colocação mostrou-se robusto, pois se compara de forma positiva com o

método de Monte Carlo com bem menos “amostras” (níveis de interpolação).

Page 114: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

105

Capítulo 6

Conclusões e trabalhos futuros

O estudo atual aplicou um método não-intrusivo, baseado em polinômio caos,

para a modelagem de liberação controlada de fertilizantes em uma esfera e um problema

de magnetohidrodinâmica . Inicialmente, na fase de qualificação, dois problemas mais

simples foram examinados com o objetivo de se aprimorar o método em estudo.

Na modelagem de liberação controlada de fertilizantes em uma esfera, para

incerteza alocada no coeficiente de distribuição de fertilizantes, foi utilizada a

quadratura de Gauss apresentando bons resultados. Quando também se aplica incerteza

no coeficiente de difusão, o uso de malhas esparsas (no caso as abscissas de Clenshaw-

Curtis) é indicado, pois reduz significativamente o número de pontos necessários para

uma solução precisa. A robustez do método pode ser constatada quando se aplica ao

mesmo problema o método de Monte Carlo e notam-se soluções convergentes com

menos pontos de colocação, no caso do método de Monte Carlo, número de amostras.

No método de Monte Carlo, em geral, foram utilizadas 10000 amostras para uma

solução convergente enquanto que no método de colocação estocástica foi utilizado no

máximo 29 pontos de colocação, também chamado de nível de interpolação 3.

No problema de magnetohidrodinâmica, a incerteza analisada ocasionou

perturbações em todas as medidas avaliadas (número de Nusselt, temperatura,

velocidade, linhas de Corrente). Estas incertezas propagadas nas soluções do problema

permitem prever a realidade de sistemas experimentais, em que é importante considerar

perturbações que indicam indefinições em propriedades do material. Neste problema,

inicialmente, foi considerada incerteza no número de Grashof. Para obter resultados

para a variação neste parâmetro, o método de colocação estocástica foi aplicado com as

Page 115: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

106

abscissas de Clenshaw-Curtis, por necessitar de muitos pontos de colocação devido à

dificuldade do problema e a regra da quadratura inibir este tipo de aplicação. A malha

esparsa, neste caso, gera menos pontos escolhidos estrategicamente. Além deste

parâmetro, também foi considerada incerteza no número de Hartman. Os resultados

apresentados para estes casos mostraram, para níveis baixos de interpolação

divergências em relação aos resultados obtidos pelo método de Monte Carlo,

necessitando-se de maiores níveis de interpolação para se obter bons resultados. Mesmo

assim, o método se mostrou robusto, pois necessitou de menos pontos do que o método

de Monte Carlo.

De modo geral, os resultados apresentados mostraram uma boa concordância em

relação aos obtidos através do método de Monte Carlo. Também, o método de

colocação estocástica com malhas esparsas ofereceu bons resultados e ainda reduziu o

custo computacional expressivamente se comparado com o produto tensorial completo

dos pontos de colocação.

Nos problemas analisados nesta tese não foi contemplada a adaptatividade com

relação a malhas esparsas, ou seja, todas as dimensões estocásticas foram tratadas de

forma semelhante, isto é, para cada dimensão o número de pontos de colocação foi o

mesmo. Para pesquisas futuras, como continuidade deste trabalho, podem-se tratar

problemas de alta dimensão, onde é interessante a aplicação de uma adaptatividade de

malhas esparsas. Neste tipo de abordagem é possível reduzir ainda mais o trabalho

computacional através de processos que identificam automaticamente as dimensões que

são mais ou menos importantes.

Page 116: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

107

Referência Bibliográfica

[1] GHANEM, R.G., SPANOS, P. D. Stochastic Finite Elements: a Spectral

Approach, Springer-Verlag, 1991.

[2] BHARUCHA-REID, A. T. “Probabilistic methods in applied

mathematics”, Vol. 1, Academic Press, New York, 1968.

[3] WINKLER, R. “An Introduction to Bayesian Inference and Decision”,

Probabilistic Publishing, Gainsville, Florida, USA, 2003.

[4] LEE, P. M. “Bayesian Statistics”, Oxford University Press, London,

2004.

[5] NAYFEH, A. H. “Perturbation Methods”, John Wiley and Sons, New

York, 1973.

[6] JORDAN, D. W., SMITH, P. “Nonlinear Ordinary Differential

Equations”, Oxford University Press, Oxford, 1977.

[7] FREDHOLM, I. “Sur une class d’equations fonctionnelles”, Acta

Mathematica, Vol. 27, pp.365-390, 1903.

[8] MIKHLIN, S. G. “Integral Equations”, Pergamon Press, Oxford, 1957.

[9] WIENER, N. “The homogeneous chaos”, Amer. J. Math., 60 (1938), pp.

897-936.

[10] XIU, D., KARNIADAKIS, G.E. “Modeling uncertainty in steady state

diffusion problems via generalized polynomial chaos”, Comput. Methods Appl.

Mech. Engrg., vol. 191, pp. 4927-4948, 2002.

[11] XIU, D., KARNIADAKIS, G.E. “The Wiener-Askey Polynomial Chaos

for stochastic differential equations”, SIAM Journal of Scientific Computing, vol

24, no. 2, pp. 619-644, 2002.

Page 117: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

108

[12] XIU, D., KARNIADAKIS, G.E. “A new stochastic approach to transient

heat conduction modeling with uncertainty”, Int. J. Heat & Mass Transfer, vol.

46, pp. 4681-4693, 2003.

[13] XIU, D., KARNIADAKIS, G.E. “Modeling uncertainty in flow

simulations via Generalized Polynomial Chaos”, J. Comp. Phys., vol. 187, pp.

137-167, 2003.

[14] WAN, X., XIU, D., KARNIADAKIS, G.E. “Stochastic solutions for the

two-dimensional advection-diffusion equation”, SIAM J. Sci. Comput., vol.

26(2), pp. 578-590, 2004.

[15] WAN, X., KARNIADAKIS, G. E. “An adaptive multi-element

generalized polynomial chaos method for stochastic differential equations”, J.

Comp. Phys., vol. 209(2), pp. 617-642, 2005.

[16] GAUTSCHI, W. “Orthogonal Polynomials (in Matlab)”, Journal of

Computational and Applied Mathematics, vol. 178, pp. 215-234, 2005.

[17] LOEVEN, A., WITTEVEEN, J., BIJL, H. “Efficient uncertainty

quantification using a two-step approach with chaos collocation”, European

Conference on Computational Fluid Dynamics, 2006.

[18] GANAPATHYSUBRAMANIAN, B., ZABARAS, N. “Sparse grid

collocation schemes for stochastic natural convection problems”, J. Comp. Phys.

225, pp. 652-685, 2007.

[19] FOO, J., YOSIBASH, Z., KARNIADAKIS, G. E. “Stochastic simulation

of riser-sections with uncertain measured pressure loads and/or uncertain

material properties”, Comput. Methods Appl. Mech. Engr., vol. 196, pp. 4250-

4271, 2007.

[20] WAN, X., KARNIADAKIS, G. E. “Solving elliptic problems with non-

Gaussian spatially-dependent random coefficients: algorithms, error analysis and

applications”, Comput. Methods Appl. Mech. Engr., invited, 2008.

[21] BABUSKA, M., NOBILE, F., TEMPONE, R. “A stochastic collocation

method for elliptic partial differential equations with random input data”. SIAM

J. Numer. Anal, 43(3):1005–1034, 2007.

[22] XIU, D., HESTHAVEN, J. S. “High-order collocation methods for

differential equations with random inputs”. SIAM J. Sci. Comput., 27(3):1118–

1139, 2005.

Page 118: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

109

[23] MATHELIN, L., HUSSAINI, M. Y. “A Stochastic Collocation

Algorithm for Uncertainty Analysis”, NASA, 2003.

[24] PEDIRODA, V., PARUSSINI, L., POLONI, C., PARASHAR, S.,

FATEH, N., POIAN, M. “Efficient Stochastic Optimization using Chaos

Collocation Method with modeFRONTIER”, SAE International, 2008.

[25] XIU, D. “Fast Numerical Methods for Stochastic Computations: A

Review”, J. Comp. Phys., vol. 5,Nº 2-4, pp. 242-272, 2009.

[26] ELDRED, M. S., BURKARDT, J. “Comparison of Non-Intrusive

Polynomial Chaos and Stochastic Collocation Methods for Uncertaint

Quantification”, American Institute of Aeronautics and Astronautics, 2009.

[27] PREMPRANEERACH, P., HOVER, F. S., TRIANTAFYLLOU, M. S.,

KARNIADAKIS, G. E. “Uncertainty quantification in simulations of power

systems: Multi-element polynomial chaos methods”, Reliability Engineering

and System Safety 95, pp. 632–646, 2010.

[28] PARUSSINI, L., PEDIRODA, V., POLONI, C. “Effects of geometric

tolerance in fluid dynamics”, ECCOMAS CFD 2010.

[29] ONORATO, G., LOEVEN, G. J. A., GHORBANIASL, G., BIJL, H.,

LACOR, C. “Comparison of intrusive and non-intrusive polynomial chaos

methods for CFD applications in aeronautics”, ECCOMAS CFD 2010.

[30] WAN, X., KARNIADAKIS, G. E. “Long-term behavior of polynomial

chaos in stochastic flow simulations”, Comput. Methods Appl. Mech. Engrg.,

vol. 195, pp. 5528-5596, 2006.

[31] MA, X., ZABARAS, N. “An adaptive hierarchical sparse grid

collocation algorithm for the solution of stochastic differential equations”, J.

Comp.Phys., vol. 228, pp. 3084–3113, 2009.

[32] GANDER, M. J., KARP, A. H. “Stable Computation of High Order

Gauss Quadrature Rules using Discretization for Measures in Radiation

Transfer”, Journal of Quantitative Spectroscopy Radiative Transfer, vol. 68, pp.

213-223, 2001.

[33] BUNGARTZ, H., GRIEBEL, M. “Sparse grids”, Acta Numerica, v. 13,

pp. 147–270, 2004.

[34] GERSTNER, T., GRIEBEL, M. “Numerical integration using sparse

grids”, Numerical Algorithms, v. 18, pp. 209–232, 1998. ISSN: 1017-1398.

Page 119: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

110

[35] NOBILE, F., TEMPONE, R., WEBSTER, C. G. “A Sparse Grid

Stochastic Collocation Method for Partial Differential Equations with Random

Input Data”, SIAM J. Numer. Anal., v. 46, pp. 2309–2345, May 2008.

[36] XIU, D. “Efficient Collocational Approach for Parametric Uncertainty

Analysis”, Commun. Comput. Phys., vol. 2, pp. 293-309,2007

[37] BUNGARTZ, H. J., DIRNSTORFER, S. “Multivariate quadrature on

adaptive sparse grids”, Computing, pp. 89-114, 2003.

[38] ELDRED, M. S. “Recent advances in non-intrusive polynomial chaos

and stochastic collocation methods for uncertainty analysis and design”,

Structures, Structural Dynamics, and Materials Conference, 2009.

[39] HOSDER. S., WALTERS, R. W., BALCH, M. “Efficient sampling for

non-intrusive polynomial chaos applications with multiple uncertain input

variables”, Structures, Structural Dynamics, and Materials Conference, 2007.

[40] AL-ZAHRANI, S. M. “Controlled-release of fertilizers: modelling and

simulation”, Int. J. of Eng. Science, vol. 37, pp.1299-1307, 1999.

[41] SHAVIV, A., RABAN, S., ZAIDEL, E. “Modeling Controlled Nutrient

Release from a Population of Polymer Coated Fertilizers: Statistically Based

Model for Diffusion Release”, Environ. Sci. Technol., vol.37, pp. 2257-2261,

2003.

[42] TONG, Z., YUHAI, L., SHIHUO, Y., ZHONGYI, H. “Superabsorbent

hydrogels as carriers for the controlled-release of urea: Experiments and a

mathematical model describing the release rate”, Biosystems Engineering,

pp.44-50, 2009.

[43] BASU, S. K., KUMAR, N. “Mathematical model and computer

simulation for release of nutrients from coated fertilizer granules”, Mathematics

and Computers in Simulation, vol. 79, pp. 634–646, 2008.

[44] BASU, S. K., KUMAR, N., SRIVASTAVA, J. P. “Modeling NPK

release from spherically coated fertilizer granules”, Simulation Modelling

Practice and Theory, vol. 18, pp.820-835, 2010.

[45] COLAÇO, M. J., DULIKRAVICH, G. S., ORLANDE, H. R. B.

“Magnetohydrodynamic Simulations Using Radial Basis Functions",

International Journal of Heat and Mass Transfer, v. 52, n. 25-26, pp. 5932-

5939,2009.

Page 120: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

111

[46] AL-NAJEM, N., KHANAFER, K., EL-REFAEE, M. “Numerical Study

of Laminar Natural Convection in Tilted Enclosure with Transverse Magnetic

Field”, International Journal of Numerical Methods for Heat and Fluid Flow, v.

8, n. 6, pp. 651-672, 1998.

[47] RUDRAIAH, N., BARRON, R. M., VENKATACHALAPPA, M.,

SUBBARAYA, C. K. “Effect of a magnetic field on free convection in a

rectangular enclosure”, Int. J. of Eng. Science, vol. 33, pp.1075-1084, 1995.

[48] MAGALHÃES, A. C., Métodos sem malha aplicados a problemas

difusivos-convectivos, Dissertação de M.Sc., IME, Rio de Janeiro, RJ, 2008.

[49] COLAÇO, M. J., DULIKRAVICH G. S. “A multilevel hybrid

optimization of magnetohydrodynamic problems in double-diffusive fluid flow”,

J. Phys.Chem. Solids, vol. 67, pp. 1965–1972, 2006.

[50] CHINCHAPATNAM, P. P., DJIDJELI, K., NAIR, P. B. “Meshless RBF

Collocation for Steady Incompressible Viscous Flows", 36th AIAA Fluid

Dynamics Conference and Exhibit, 2006.

[51] GOLUB, G. H., WELSCH, J. H. “Calculation of Gauss Quadrature

Rules”, Mathematics of Computation, vol. 23(106), pp. 221-230, 1969.

[52] TREFETHEN, L. N. “Is Gauss Quadrature better than Clenshaw-

Curtis?”, SIAM Rev., v. 50, pp. 67–87, February 2008.

[53] ANDERSON, D. A., TANNEHILL, J. C., PLETCHER, R. H.

“Computational fluid mechanics and heat transfer”, Hemisphere, 1984.

[54] COLAÇO, M. J., ORLANDE, H. R. B., ROBERTY, N. C., ALVES, C.,

LEITÃO, V. “On the solution of difusion-convection problems by means of

RBF Approximations", Proceedings of the 11th Brazilian Congress of Thermal

Sciences and Engineering - ENCIT, 2006.

[55] HARDY, R. “Multiquadric equations of topography and other irregular

surfaces”, J. Geophys. Res., v. 176, pp. 1905-1930, 1971.

[56] KANSA, E. J. “Multiquadrics, a scattered data approximation scheme

with applications to computational fluid-dynamics I - Surface approximations

and partial derivative estimates”, Computers Mathematics with Applications, v.

19, n. 8-9, pp. 127-145, 1990.

[57] KANSA, E. J. “Multiquadrics, a scattered data approximation scheme

with applications to computational fluid-dynamics II - Solutions to parabolic,

Page 121: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

112

hyperbolic and elliptic partial differential equations”, Computers Mathematics

with Applications, v. 19, n. 8-9, pp. 147-161, 1990.

[58] BUHMANN, M. D. “Radial basis functions: Theory and

Implementations”, Cambridge University, 2004.

[59] SOBOL, I., LEVITAN, Y. L. “The Production of Points Uniformly

Distributed in a Multidimensional Cube", USSR Academy of Sciences, v. 40,

1976.

[60] DULIKRAVICH, G. S., LYNN, S. R. “Unified electro-magneto-fluid

dynamics (EMFD): a survey of mathematical models”, Int. J. Non-Linear Mech.

vol. 32, n. 5, pp. 923–932, 1997.

[61] GARANDET, J. P., ALBOUSSIERE, T., MOREAU, R. “Buoyancy

driven convection in a rectangular enclosure with a transverse magnetic field”,

Int. J. Heat Mass Transfer, v. 35,n. 4, pp. 741–748,1992.

[62] ASKEY, R., WILSON, J. “Some basic hypergeometric polynomials that

generalize Jacobi polynomials”, Memoirs Amer. Math. Soc., AMS, Providence

RI, 319 (1985).

[63] ALBUQUERQUE, E. C. M. C. “Preparação, caracterização de

Lipossomas e Microesferas para aplicação na terapia por dessensibilização de

reação alérgica.”, UNICAMP, 2005.

[64] BERNADÁ, G. M. G., Quantificação de incertezas em problemas de

interação fluido estrutura via método de colocação estocástica, Tese de D.Sc.,

COPPE/UFRJ, Rio de Janeiro, RJ, 2011.

[65] KLIMKE, A. Uncertainty Modeling using Fuzzy Arithmetic and Sparse

Grids, PhD Thesis, Universitt Stuttgart, Shaker Verlag, Aachen, 2006.

[66] LACERDA, C. R., Solução de problemas convectivos-difusivos na

presença de campos magnéticos usando métodos sem malha, Dissertação de

M.Sc., COPPE/UFRJ, Rio de Janeiro, RJ, 2013.

[67] BURKARDT, J. “FORTRAN90 software”. http://people.sc.fsu.edu/

jburkardt/, 2000.

[68] NOBILE, F., TEMPONE, R., WEBSTER, C. G. “An Anisotropic Sparse

Grid Stochastic Collocation Method for Partial Differential Equations with

Random input Data”, SIAM J. Numer. Anal., v. 46, pp. 2411–2442, June 2008.

Page 122: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

113

[69] BARTHELMANN, V., NOVAK, E., RITTER, K. “High dimensional

polynomial interpolation on sparse grids”, Adv. Compu. Math. 12 (2000) 273–

288.

[70] PRESS, W., FLANNERY, B., TEUKOLSKY, S. “Numerical Recipes

FORTRAN", Cambridge Universaty Press, 1992.

[71] FOO, J., KARNIADAKIS, G. E. “Multi-element probabilistic collocation

method in high dimensions”, J. Comp. Phys., vol. 229, pp. 1536-1557, 2010.

[72] COMPO, do Brasil Ltda. Disponível em:http://www.compo-expert.com/br.

Acesso em: 8 nov. 2011, 13:30.

[73] MAGALHÃES, A. C., COLAÇO, M. J., ORLANDE, H. R. B. “Radial

Basis Functions Applied to the Solution of the Natural Convection Problem in

Square Cavities", Proceedings of the 12th Brazilian Congress of Thermal

Sciences and Engineering - ENCIT, 2008.

[74] WERTZ, J., KANSA, E. J., LING, L. “The role of the multiquadric shape

parameters in solving elliptic partial differential equations", Computers and

Mathematics with Applications, v. 51, n. 8, pp. 1335-1348, 2006.

[75] FISHMAN, G. S. Monte Carlo: Concepts, Algorithms and Applications.

Third ed. Department of Operations Research, Stanford University, Stanford,

CA 94305, Springer-Verlag, 1999.

Page 123: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

114

Apêndice A

Obtenção da matriz Jacobiana ( J )

Através da Relação de Recorrência de Três Termos (3.4.2.a) podemos descrever

o seguinte sistema de equações [51]:

1i 1 1 2( ) ( ) ( ) 0

2i 2 1 2 2 3( ) ( ) ( ) ( ) 0

3i 3 2 3 3 4( ) ( ) ( ) ( ) 0

1Pi N 1 2 1 1( ) ( ) ( ) ( ) 0P P P P PN N N N N

Pi N 1( ) ( ) ( ) 0P P P PN N N N

Esse sistema pode ser escrito na forma matricial:

1 1

2 2 2

3 3 3

1 1 1

( ) 1 0 0 0 0

( ) 1 0 0 0

0 ( ) 0 00

0 0 1 0

0 0 0 ( ) 1

0 0 0 0 ( )

P P P

P P P

N N N

N N N

Page 124: SIMULAÇÕES SOB INCERTEZAS EM PROBLEMAS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/1321/pemufrj2014dsc… · equação diferencial parcial representativa de um fenômeno

115

Ou na forma:

1 1 1

2 2 2 2

3 3 3 3

1 1 1 1

A

1 0 0 0 0

1 0 0 0

0 0 0

0 0 1 0

0 0 0 1

0 0 0 0

P P P P

P P P P

N N N N

N N N N

Logo, as raízes desses polinômios ortogonais são também os autovalores dessa

matriz. Aplicando a transformação diagonal de semelhança, a qual garante que a matriz

inicial A têm os mesmos autovalores que a matriz encontrada J (matriz jacobiana

simétrica) [23], encontramos:

1 2

2 2 3

3 3

1

1 1

0 0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0 0

P

P P P

P P

N

N N N

N N

J

Os pontos de colocação e os pesos da regra da quadratura podem ser calculados

através desse autosistema da matriz tridiagonal simétrica.