View
225
Download
0
Category
Preview:
Citation preview
Mauro Artemio Carrión Pachás
Análise Limite com Otimizador de Grande Escala e Análise de Confiabilidade
TESE DE DOUTORADO
Tese apresentada ao Programa de Pós-graduação em Engenharia Civil da PUC-Rio como requisito parcial para obtenção do título de Doutor em Engenharia Civil. Ênfase: Geotécnica
Orientador: Eurípedes do Amaral Vargas Júnior
Co-orientadores: Luíz Eloy Vaz
José Herskovits Norman
Rio de Janeiro, março de 2009
Mauro Artemio Carrión Pachás
Análise Limite com Otimizador de Grande Escala e Análise de Confiabilidade
Tese apresentada ao Programa de Pós-Graduação em Engenharia Civil da PUC-Rio como requisito parcial para obtenção do título de Doutor em Engenharia Civil. Aprovada pela Comissão Examinadora abaixo assinada.
Prof. Eurípedes do Amaral Vargas Júnior Presidente/Orientador
Departamento de Engenharia Civil – PUC-Rio
Prof. Luiz Eloy Vaz Co-Orientador-UFRJ
Prof. José Herskovits Norman Co-Orientador-UFRJ
Prof. Luiz Fernando Marta Departamento de Engenharia Civil – PUC-Rio
Prof. Ivan Menezes PUC- Rio
Prof. Aldo Duran Farfán UENF
Prof. Silvia Almeida UFG
Prof. José Eugênio Leal Coordenador Setorial do Centro Técnico Científico – PUC–Rio
Rio de Janeiro, março de 2009
Todos os direitos reservados. É proibida a reprodução total ou parcial do trabalho sem autorização da universidade, do autor e do orientador.
Mauro Artemio Carrión Pachás Graduou-se em Engenharia Civil na Universidade Nacional de Engenharia (UNI-PERU) em 1996. Trabalhou como pesquisador no Centro Peruano Japonês de Investigações Sísmicas e Mitigação de Desastres CISMID em Lima-Perú no período de 1997 a 2002. Estudou mestrado na PUC-Rio, em Engenharia Civil, na área de Geotecnia, no período 2002.2-2004.1. Ingressou no curso de doutorado na PUC-Rio no período 2004.2, atuando na linha de Pesquisa Numérica. Análise Limite com Otimizador de Grande Escala e Análise de Confiabilidade.
Ficha Catalográfica
Carrión Pachás, Mauro Artemio
Análise Limite com Otimizador de Grande Escala e Análise de Confiabilidade / Mauro Artemio Carrión Pachás; orientador: Eurípedes do Amaral Vargas Júnior; co-orientadores: Luiz Eloy Vaz, José Herskovits Norman. - 2009.
188 f. : il. ; 30 cm
Tese (Doutorado em Engenharia Civil) – Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro, 2009.
Incluí referências bibliográficas.
1. Engenharia civil – Teses. 2. Análise limite. 3. Métdo de elementos finito. 4. Otimização. 5. Escoamento. 6. GEOLIMA. 7. Confiabilidade. 8. Função de falha. 9. FORM. I. Vargas Júnior, Eurípedes do Amaral. II. Vaz, Luiz Eloy. III. Norman, José Herskovits. IV. Pontifícia Universidade Católica do Rio de Janeiro. Departamento de Engenharia Civil. V. Título.
CDD: 624
Para meus pais, Juan Carrión e Marina Pachás,
pelo grande amor, confiança e exemplo..
Agradecimentos
Desejo expressar minha gratidão ao professor Eurípedes do Amaral Vargas Júnior,
por ter-me convencido para fazer o curso de doutorado.
Aos professores Eurípedes do Amaral Vargas Júnior, Luíz Eloy Vaz e José
Herskovits Norman pela orientação durante a realização deste trabalho.
Ao Professor Luiz Fernando Martha pela confiança e apoio.
Aos professores dos departamentos de Civil, Informática e Elétrica, pelos
conhecimentos transmitidos em cada uma das disciplinas que cursei.
A meus pais, irmão, irmãs, meu primo Felix e a toda minha família, que sempre
me apoiaram e incentivaram para a realização deste curso de doutorado.
À Mishel, pelo grande amor e compreensão, muito obrigado.
À Rita, secretária da pós-graduação, por sua atenção e disponibilidade.
À PUC-Rio, CAPES e CNPq pelos auxílios financeiros concedidos e ao
TECGRAF/PUC-Rio pela oportunidade de poder trabalhar e estudar, sem a qual
não teria sido possível realizar este estudo.
A todos os colegas do trabalho, do estudo e das peladas, muito obrigado pela
convivência.
A Deus, porque sem a ajuda d’Ele nada acontece.
Resumo
Pachás, Mauro Artemio Carrión; Vargas Júnior, Eurípedes do Amaral; Vaz, Luíz Eloy e Herskovits, José Norman. Análise Limite com Otimizador de Grande Escala e Análise de Confiabilidade. Rio de Janeiro, 2008. 188p. Tese de Doutorado - Departamento de Engenharia Civil, Pontifícia Universidade Católica do Rio de Janeiro.
O presente trabalho tem por objetivo desenvolver um otimizador eficiente
de grande escala, que permita a aplicabilidade prática da Análise Limite Numérica
pelo MEF, para resolver problemas reais da Engenharia Geotécnica. Para isto, foi
desenvolvido um otimizador para o programa GEOLIMA (GEOtechnical LIMit
Analysis) (Carrión, 2004) baseado no algoritmo de Pontos Interiores,
computacionalmente mais eficiente que os otimizadores comerciais existentes.
Pelo fato das propriedades do solo serem de natureza aleatória, a possibilidade de
aplicar Análise de Confiabilidade com a Análise Limite pelo método FORM em
problemas geotécnicos é pesquisada também. Sendo a grande vantagem do
método FORM a possibilidade de se aplicar para funções de falha quaisquer e
variáveis com distribuição quaisquer. Inicialmente, são apresentados os
fundamentos da teoria de Análise Limite e sua formulação numérica pelo MEF
(Método dos Elementos Finitos). A seguir, é investigada a possibilidade de se usar
otimizadores comerciais para resolver o problema matemático resultante da
aplicação de Análise Limite com o MEF e são descritos os fundamentos teóricos
do otimizador implementado baseado no algoritmo de Pontos Interiores. Um
resumo dos fundamentos teóricos da Análise de Confiabilidade é apresentado. É
descrito o processo de cálculo pelo método FORM e dois exemplos de aplicação
são realizados. Finalmente, análises de diferentes problemas resolvidos com o
otimizador implementado são apresentados indicando o grande potencial da
Análise Limite Numérica, na solução de problemas reais da Engenharia
Geotécnica.
Palavras-chave Análise limite; método de elementos finitos; otimização; função de
escoamento; GEOLIMA; confiabilidade; função de falha; FORM.
Abstract
Pachás, Mauro Artemio Carrión; Vargas Júnior, Eurípedes do Amaral; Vaz, Luíz Eloy e Herskovits, José Norman (Advisors). Limit Analysis with Large Scale Optimizer and Reliability Analysis. Rio de Janeiro, 2008. 188p. DrSc. These - Department of Civil Engineering, Pontifícia Universidade Católica do Rio de Janeiro.
This work has, as its main objective, the development of an efficient and
large scale optimizer, that allows the practical application of Numerical Limit
Analysis (NLA) with Finite Element Method (FEM) to solve real problems in
Geotechnical Engineering. For that purpose, an optimizer was developed for
GEOLIMA (GEOtechnical LIMit Analysis) program (Carrión, 2004), based on
Interior Points algorithm, computationally more efficient than the existing
commercial optimizers. Due to the fact that soils have random properties, the
possibility to apply Reliability Analysis with Limit Analysis using the FORM
method was also investigated. Initially, Limit Analysis theory was presented
together with its numerical formulation using the FEM. In sequence, the use of
commercial optimizers was investigated in order to solve the resulting
mathematical problem. Subsequently, the theorical foundations of the developed
optimizer, based on the Interior Points algorithm were described. A summary of
Reliability Analysis was also presented together with a description of
computational procedures using FORM and two examples were developed.
Finally, analyses of different problems solved with developed optimizer were
presented. The obtained results demonstrated the great potential of Numerical
Limit Analysis (NLA), in the solution of real problems in Geotechnical
Engineering.
Keywords Limit analysis; finite element method; optimization; yield function;
GEOLIMA; reliability; failure function; FORM.
Sumário
1 INTRODUÇÃO 21
2 ANÁLISE LIMITE NUMÉRICA 23
2.1. Teoremas da Análise Limite 24
2.1.1. Campos de Tensões Estaticamente Admissíveis 24
2.1.2. Campos de Velocidades Cinematicamente Admissíveis 24
2.1.3. Teorema de Limite Inferior 25
2.1.4. Teorema de Limite Superior 25
2.2. Considerações na Análise Limite 25
2.2.1. Consideração de Plasticidade Perfeita 25
2.2.2. Considerações sobre Escoamento 26
2.2.3. Considerações sobre a Lei de Fluxo 28
2.3. Principio dos Trabalhos Virtuais 30
2.4. Critérios de Escoamento 30
2.4.1. Critério de Mohr-Coulomb 31
2.4.2. Critério de Drucker-Prager 33
2.5. Formulação Numérica da Análise Limite pelo MEF 35
2.5.1. Condição de Equilíbrio 36
2.5.2. Condições de Contorno 38
2.5.3. Condição de Escoamento 38
2.5.4. Problema de Análise Limite 39
2.5.5. Elementos Finitos Implementados 39
3 SOLUÇÃO NUMÉRICA DA ANÁLISE LIMITE 41
3.1. Otimizadores Matemáticos Testados 42
3.1.1. Otimizador Lingo 43
3.1.2. Otimizador Minos 47
3.1.3. Otimizador Lancelot 50
3.1.4. Comparação de Desempenho de Otimizadores Testados 52
3.2. Otimizador Implementado 54
3.2.1. Condições de Otimalidade 54
3.2.2. Algoritmo de Pontos Interiores 55
3.2.3. Inicialização 55
3.2.4. Direção de Busca 56
3.2.4.1. Técnica de deflexão 57
3.2.4.2. Técnica de relaxação-contração 58
3.2.4.3. Técnica Vetorial - Proposto 59
3.2.5. Comprimento de Passo 60
3.2.6. Atualização das Variáveis 61
3.2.7. Teste de Desempenho dos Algoritmos 62
3.2.8. Manipulação de sistemas lineares a serem resolvidos 63
3.2.8.1. Manipulação matricial global 64
3.2.8.2. Manipulação matricial por elementos 65
3.2.8.3. Solução direta sem manipulação - proposta 67
3.2.8.4. Teste de Desempenho das Manipulações 68
3.2.9. Resolvedores Implementados 75
3.2.9.1. Método dos gradientes conjugados (CG) 76
3.2.9.2. Teste de desempenho de resolvedores implementados 78
3.2.10. Resolvedor SAMG Testado 80
3.3. Melhora do desempenho 80
3.3.1. Tratamento de Matriz Esparsa 80
3.3.1.1. Formato CSR(Compressed Sparse Row) 80
3.3.1.2. Teste de desempenho de CG com tratamento da matriz esparsa81
3.3.2. Precondicionamento 82
3.3.2.1. Escala Diagonal (DS) 83
3.3.2.2. Escala Simétrica (SS) 83
3.3.2.3. Fatoração Incompleta de Cholesky (ICF) 84
3.3.2.4. Pré-condicionadores mistos - proposto 85
3.3.2.5. Teste de Desempenho de Pré-condicionadores 85
3.4. Teste de desempenho do Otimizador Implementado 87
3.4.1. Teste com problema em 2D 87
3.4.2. Teste com problema em 3D 90
4 ANÁLISE DE CONFIABILIDADE COM ANÁLISE LIMITE 93
4.1. Conceitos Fundamentais da Análise de Confiabilidade 95
4.1.1. Incertezas 95
4.1.2. Função de Falha 95
4.1.3. Função Densidade de Probabilidade Conjunta 97
4.1.4. Probabilidade de Falha 99
4.1.5. Confiabilidade 100
4.1.6. Índice de Confiabilidade 101
4.1.7. Espaço Reduzido 103
4.1.8. Distribuição Normal Equivalente 105
4.2. Métodos de Cálculo 107
4.2.1. Método FORM (First-Order Reliability Method) 112
4.2.1.1. Transformação de Variáveis 112
4.2.1.2. Pesquisa de Ponto de Projeto 115
4.2.1.3. Processo de Cálculo 117
4.3. Exemplos de Aplicação 118
4.3.1. Talude 2D 118
4.3.2. Talude Confinado 3D 122
5 APLICAÇÕES 125
5.1. Análise 2D - Talude Infinito Homogêneo 126
5.2. Análise 2D - Talude Infinito Heterogêneo 134
5.3. Análise 2D - Talude com Percolação 142
5.4. Análise 2D - Barragem de Terra 146
5.5. Análise 3D - Talude Confinado 151
5.6. Análise 3D - Depósito de Rejeito 155
6 CONCLUSÕES E SUGESTÕES 161
6.1. Conclusões 161
6.2. Sugestões para Futuras Pesquisas 165
7 REFERÊNCIAS BIBLIOGRÁFICAS 166
A CONCEITOS ESTATÍSTICOS 171
A.1. Variáveis Determinísticas e Aleatórias 171
A.2. Espaço Amostral, Evento e Valor Observado 171
A.3. Medidas de Tendência Central 171
A.3.1. Momento estatístico de ordem m 172
A.3.2. Média aritmética 172
A.3.3. Características da média aritmética 173
A.4. Medidas de Dispersão da Variável Aleatória 173
A.4.1. Desvio 173
A.4.2. Momento estatístico central de ordem m 174
A.4.3. Desvio absoluto médio 174
A.4.4. Variância 174
A.4.5. Desvio padrão 175
A.4.6. Coeficiente de variação 176
A.5. Medidas de Correlação de Variáveis Aleatórias 177
A.5.1. Covariância 177
A.5.2. Coeficiente de correlação 178
A.6. Caracterização de Variáveis Aleatórias 180
A.6.1. Função Densidade de Probabilidade (PDF) 180
A.6.2. Função Distribuição de Probabilidade 182
A.6.3. Coeficiente de Inclinação de uma Distribuição 185
A.6.4. Coeficiente de Curtose 186
A.7. Esperança Matemática de uma Variável Aleatória 186
Lista de figuras
Figura 2.1 – Relação tensão deformação para solo real e ideal (Chen, 1975). 26
Figura 2.2 – Superfície de escoamento no espaço de tensões principais. 27
Figura 2.3 – Superfície de escoamento e vetor de deformação plástica. 29
Figura 2.4 – Superfície de escoamento, critério de Mohr-Coulomb. 32
Figura 2.5 – Critério de escoamento de Mohr-Coulomb 2D. 33
Figura 2.6 – Critério de escoamento de Drucker & Prager. 34
Figura 2.7 – Elemento finito: (a) quadrilateral (2D), (b) hexaédrico (3D). 40
Figura 3.1 – Problema para teste de otimizadores. 42
Figura 3.2 – Malhas: (a) 28, (b) 64, (c) 126, (d) 225, (e) 360, (f) 500 e (g) 750
elementos. 43
Figura 3.3a – Variação da memória usada pelo otimizador LINGO. 45
Figura 3.3b – Variação de número de iterações do otimizador LINGO. 45
Figura 3.3c – Desempenho do otimizador LINGO. 46
Figura 3.3d – Variação do fator de colapso obtido com LINGO. 46
Figura 3.4a – Variação da memória usada pelo otimizador MINOS. 48
Figura 3.4b – variação de número de iterações do otimizador MINOS. 48
Figura 3.4c – Desempenho do otimizador MINOS. 49
Figura 3.4d – Variação do fator de colapso obtido com MINOS. 49
Figura 3.5a – Comparação de uso da memória entre Lingo e Minos 52
Figura 3.5b – Comparação de número de Iterações entre Lingo e Minos 52
Figura 3.5c – Comparação de desempenho entre Lingo e Minos 53
Figura 3.5d – Comparação de variação de fator de colapso, obtidos com Lingo e
Minos. 53
Figura 3.6 – Técnica de deflexão. 57
Figura 3.7 – Técnica de Relaxação-Contração. 59
Figura 3.8– Técnica Vetorial - Proposta. 59
Figura 3.9 – Comprimento de passo s . 61
Figura 3.10 – Atualização da variável x . 61
Figura 3.11 Geometria do problema para teste de algoritmos (malha de 25
elementos). 62
Figura 3.12 – Memória requerida pela manipulação matricial global. 69
Figura 3.13 – Tempo empregado na manipulação e solução do sistema. 70
Figura 3.14 – Memória requerida para a solução do sistema. 71
Figura 3.15 - Tempo empregado na manipulação e solução do sistema. 71
Figura 3.16 - Memória requerida para a solução do sistema. 72
Figura 3.17 – Tempo empregado na solução do sistema. 73
Figura 3.18 – Memória usada pelas técnicas. 74
Figura 3.19 – Tempo empregado pelas técnicas. 75
Figura 3.20 – Algoritmo de Gradiente Conjugado pré-condicionado (Sandoval,
2006). 77
Figura 3.21 – Malhas: (a) 8, (b) 25, (c) 64, (d) 150, (e) 400 e (f) 676 elementos. 78
Figura 3.22 – Comparação do desempenho dos métodos implementados. 79
Figura 3.23a – Armazenamento da matriz esparsa (SMAILBEGOVIC et all,
2006). 81
Figura 3.23b – Produto de uma matriz esparsa A por um vetor d . 81
Figura 3.24 – Desempenho do método CG com tratamento da matriz esparsa CSR.
82
Figura 3.25 – Desempenho de CG com os pré-condicionadores implementados. 86
Figura 3.26 – Comparação de uso da memória pelos otimizadores. 88
Figura 3.27 – Comparação de número de iterações. 88
Figura 3.28 – Comparação do desempenho dos otimizadores. 89
Figura 3.29 – Variação de fator de colapso. 89
Figura 3.30 – Geometria da estrutura a ser analisada. 91
Figura 3.31 – Malha de elementos finitos (676 elementos e 945 nós). 91
Figura 3.32 – Mecanismo de colapso da estrutura obtido pelo programa
GEOLIMA 2.0. 92
Figura 3.33 – Mecanismo de ruptura obtido a partir de ensaios em modelo físico
em escala reduzida (Sterpi,1996) 92
Figura 4.1 – Função de falha. 96
Figura 4.2 – Função densidade de probabilidade conjunta (Melchers, 2002). 99
Figura 4.3 – Probabilidade de falha. 100
Figura 4.4 – Função densidade de probabilidade. 102
Figura 4.5 – Espaço original e espaço reduzido para uma variável. 104
Figura 4.6 – Espaço original e reduzido para duas variáveis. 105
Figura 4.7 – Funções densidade de probabilidade PDF. 105
Figura 4.8 – Funções distribuição de probabilidade CDF. 106
Figura 4.9 – Variáveis em coordenadas reduzidas: função de falha não linear. 110
Figura 4.10 – Espaço original e reduzido para duas variáveis. 116
Figura 4.11 – Fluxograma da Análise de Confiabilidade pelo método FORM. 117
Figura 4.12 – Malha de elementos finitos 2D (25 elementos com 36 nós). 119
Figura 4.13 – Zonas de plastificação no MPP(Most Probable Point). 120
Figura 4.14 – Superfície de falha no MPP. 120
Figura 4.15 – Malha de elementos finitos 3D (27 elementos com 64 nós). 123
Figura 4.16 – Zonas de plastificação no MPP. 123
Figura 4.17 – Superfície de falha no MPP. 124
Figura 5.1 – Talude Infinito homogêneo. 127
Figura 5.2 – Malha de elementos finitos (500 elementos e 561 nós). 127
Figura 5.3(a) – Zonas de plastificação (talude com material coesivo). 128
Figura 5.3(b) – Vetor de velocidades (material coesivo). 128
Figura 5.3(c) – Superfície de falha (talude com material coesivo). 129
Figura 5.3(d) – Mecanismo de colapso (talude com material coesivo). 129
Figura 5.4(a) – Zonas de plastificação (talude com material com atrito). 130
Figura 5.4(b) – Vetor de veolicidades (talude com material com atrito). 130
Figura 5.4(c) – Superfície de falha (talude com material com atrito). 131
Figura 5.4(d) – Mecanismo de colapso (talude com material com atrito). 131
Figura 5.5(a) – Zonas de plastificação (talude com material maciço rochoso). 132
Figura 5.5(b) – Vetor de velocidades (talude com material maciço rochoso). 132
Figura 5.5(b) – Superfície de falha (talude com material maciço rochoso). 133
Figura 5.5(c) – Mecanismo de colapso (talude com material maciço rochoso). 133
Figura 5.6 – Talude infinito heterogêneo. 136
Figura 5.7(a) – Malha de elementos finitos 1 (100 elementos e 126 nós) . 136
Figura 5.7(b) – Zonas de plastificação (malha: 100 elementos e 126 nós). 137
Figura 5.7(c) – Superfície de falha (malha: 100 elementos e 126 nós). 137
Figura 5.8(a) – Malha de elementos finitos 2 (500 elementos e 561 nós). 138
Figura 5.8(b) – Zonas de plstificação (malha: 500 elementos e 561 nós). 138
Figura 5.8(c) – Superfície de falha (malha: 500 elementos e 561 nós). 139
Figura 5.9(a) – Malha de elementos finitos 3 (4500 elementos e 4681 nós). 139
Figura 5.9(b) – Zonas de plastificação (malha: 4500 elementos e 4681 nós). 140
Figura 5.9(c) – Vetor de velocidades (malha: 4500 elementos e 4681 nós). 140
Figura 5.9(d) – Superfície de falha (malha: 4500 elementos e 4681 nós). 141
Figura 5.10 – Geometria de talude. 143
Figura 5.11 – Malha de elementos finitos (200 elementos e 231 nós). 143
Figura 5.12(a) – Zonas de plastificação, sem considerar percolação. 144
Figura 5.12(b) – Superfície de falha sem considerar percolação. 144
Figura 5.12(c) – Vetor de velocidades sem considerar percolação. 144
Figura 5.13(a) – Zonas de plastificação, considerando percolação. 145
Figura 5.13(b) – Superfície de falha, considerando percolação. 145
Figura 5.13(c) – Vetor de velocidades, considerando percolação. 145
Figura 5.14 – Plano em planta do projeto (CISMID, 2002). 147
Figura 5.15 – Secção 3-3 Barragem de Terra (CISMID, 2002). 148
Figura 5.16 – Malha de elementos finitos (624 elementos e 689 nós). 148
Figura 5.17(a) – Zonas de plastificação. 149
Figura 5.17(b) – Superfície de falha. 149
Figura 5.17(c) – Vetores de velocidade. 149
Figura 5.18(a) – Zonas de plastificação. 150
Figura 5.18(b) – Superfície de falha. 150
Figura 5.18(c) – Vetor de velocidade. 150
Figura 5.19 – Malha de elementos finitos (8000 elementos e 9261 nós). 152
Figura 5.20(a) – Zonas que plastificam. 152
Figura 5.20(b) – Vetor de veolocidades. 153
Figura 5.20(c) – Distribuição de campo de velocidades. 153
Figura 5.20(d) – Distribuição de campo de velocidades na seção longitudinal. 154
Figura 5.20 (e) – Mecanismo de colapso do talude. 154
Figura 5.21 – Topografia e limite do material de rejeito (CISMID, 1998). 156
Figura 5.22 – Modelo 3D da bacia sem material de rejeito. 157
Figura 5.23 – Material de rejeito. 158
Figura 5.24 – Malha de elementos finitos (2000 elementos e 2541 nós). 158
Figura 5.25(a) – Zonas de plastificação. 159
Figura 5.25(b) – Zonas de plastificação na seção longitudinal. 159
Figura 5.25(c) –Vetor de velocidades. 159
Figura 5.25(d) – Distribuição de velocidades. 160
Figura 5.25(e) – Distribuição de velocidades na seção longitudinal. 160
Figura A.1 – Interpretação gráfica de coeficiente de correlação. 179
Figura A.2 – Função densidade de probabilidade normal. 182
Figura A.3 – Função distribuição de probabilidade normal. 183
Figura A.4 – (a) Função de variável aleatória; (b) Função de ponderação; (c)
Esperança matemática da função de variável aleatória. 188
Lista de tabelas
Tabela 3.1 – Testes realizados com o programa Lingo. 44
Tabela 3.2 – Testes realizadas com o otimizador Minos. 47
Tabela 3.3 – Comparação da eficiência entre Lancelot e Minos (Bongartz et all,
1997). 51
Tabela 3.4 – Desempenho dos algoritmos implementados 63
Tabela 3.5 – Teste da manipulação matricial global. 69
Tabela 3.6 – Teste da manipulação matricial por elementos. 70
Tabela 3.7 – Teste da solução direta sem manipulação. 72
Tabela 3.8 – Comparação de resultados com as três técnicas de solução. 74
Tabela 3.9 – Comparação do desempenho dos métodos implementados. 79
Tabela 3.10 – Desempenho dos métodos com o tratamento da matriz esparsa. 82
Tabela 3.11 – Teste de pré-condicionadores implementados. 86
Tabela 3.12 – Resultados do teste em 2D do otimizador implementado. 87
Tabela 4.1 – Índice de confiabilidade e probabilidade de falha. 103
Tabela 4.2 – Coeficiente de correlação equivalente. 113
Tabela 5.1 – Propriedades de materiais e resultados da Análise Limite. 126
Tabela 5.2 – Problema de otimização e resultados da Análise Limite. 135
Tabela 5.3 – Propriedades de materiais e resultados da Análise Limite. 143
Tabela 5.4 – Propriedades dos materiais (CISMID, 2002). 148
Tabela 5.4 – Resultados da Análise Limite. 148
Tabela 5.5 – Propriedades dos materiais (CISMID, 1998). 157
Tabela A.1 – Grão de dependência das variáveis aleatórias. 179
Tabela A.2 – Funções densidade probabilidade (PDF) e de distribuição de
probabilidade (CDF) ( Lopes, 2007). 184
Lista de Símbolos
Na Análise Limite (AL) )(F ijσ Função de escoamento em termos de tensões
ijσ Campo de tensões
I1 Primeiro invariante de tensor de tensões
J2 Segundo invariante desviador
J3 Terceiro invariante desviador
ijε& Velocidade de deformação total eijε& Velocidade de deformação elástica
pijε& Velocidade de deformação plástica
λ& Fator de proporcionalidade (escalar).
σ1 Tensão principal maior
σ2 Tensão principal intermédia
σ3 Tensão principal menor
θ Ângulo de Lode
φ Ângulo de atrito do material
C Coesão do material
α,k Parâmetros do material (critério de Drucker-Prager)
oc Resistência à compressão uniaxial
No MEF
Ni Função de interpolação
r,s,t Coordenadas locais
ri,si,ti Coordenadas locais nodais
σ Vetor campo de tensões
σ Vetor campo de tensões nodais
u& Vetor campo de velocidades
u& Vetor campo de velocidades nodais
ε& Vetor campo de velocidades de deformação
0F Força de volume inicial
0T Força de superfície inicial
a Fator de colapso
B Matriz de equilíbrio
b Vetor de carregamentos nodais
Otimização x Vetor das variáveis
µλ, Multiplicadores de Lagrange
n Número de variáveis
m Número de restrições de desigualdade
p Número de restrições de igualdade
)(αf Função objetivo
)( xg Restrições de desigualdade
G Matriz de escoamento (diagonal)
S Matriz gradiente
H Matriz Hessiana
θ Ângulo de deflexão
Confiabilidade Xi Variável aleatória i
ij x j-ésimo valor observado da variável aleatória Xi
mi
s X Momento estatístico de ordem m da amostra
is X Média da amostra da variável aleatória Xi
mi
s X Momento estatístico central de ordem m da amostra
is X Desvio absoluto médio da amostra
)x(Var is Variância da variável aleatória Xi
ixσ Desvio padrão da variável aleatória Xi
σ Matriz desvio padrão
ixδ Coeficiente de variação da variável aleatória Xi
Cov(X1,X2) Covariância das variáveis aleatórias X1 e X2
S Matriz covariância
kj x,xρ Coeficiente de correlação das variáveis aleatórias Xj e Xk
ρ Matriz coeficiente de correlação p(x) Função densidade de probabilidade
)( axP Função distribuição de probabilidade
]Pr[ axX ≤
Probabilidade de que a variável X seja menor ou igual a xa
X Vetor de variáveis aleatórias
)(XF Função de falha
F Média da função de falha
Fσ Desvio padrão da função de falha
fP Probabilidade de falha
C Confiabilidade
β Índice de confiabilidade
Y Variável aleatória reduzida
φ Função distribuição de probabilidade normal padrão
Φ Função cumulativa normal padrão
J Jacobiano
L Matriz triangular inferior
∇ Operador diferencial *y Ponto de pesquisa
1 INTRODUÇÃO
A avaliação de segurança ou estabilidade de Estruturas Geotécnicas é
resolvida, na prática, por métodos simples como Equilíbrio Limite, que é um
método de tentativas e erros, baseado na comparação de forças atuantes e
resistentes sobre uma superfície de ruptura predeterminada pelo usuário ou por
algum programa de computador. Tais métodos se resumem a encontrar uma
superfície com menor fator de segurança, definida normalmente como o quociente
entre a força atuante e a força resistente.
Por outro lado existem na teoria de plasticidade, métodos como a Análise
Limite, baseados nos teoremas do limite inferior e superior. Esses métodos, apesar
de oferecerem boas perspectivas para a análise de estabilidade de problemas
Geotécnicos, encontram-se atualmente em nível de pesquisa no âmbito
acadêmico, não sendo aplicados na prática, por falta de algoritmos ou métodos
eficientes para resolver o problema matemático resultante da aplicação da teoria
de Análise Limite (AL) e do Método de Elementos Finitos (MEF). Na formulação
mista fraca, empregada no presente trabalho, o problema consiste em encontrar
um campo de tensões correspondente à carga que leva ao colapso uma estrutura
(carga de colapso). Este problema pode ser resolvido por alguma técnica de
otimização, onde a função objetivo é o fator de carga e as restrições do problema
são expressas na forma de um sistema de equações que correspondem às equações
de equilíbrio com número de variáveis muito maior que o número de equações e
um outro sistema de inequações desacopladas não lineares que satisfazem um
critério de escoamento do material.
O objetivo principal do presente trabalho é aplicar a Análise Limite na
análise de problemas reais de Geotecnia de maneira computacionalmente
eficiente, e pelo fato de testes realizados com otimizadores comerciais existentes
ter-se mostrado ineficientes, o grande desafio é desenvolver um otimizador
eficiente para resolver o problema descrito no parágrafo anterior.
22
Este trabalho é uma continuação da pesquisa realizada na dissertação de
mestrado do autor, onde foi desenvolvido o programa computacional GEOLIMA
(GEOtechnical LIMit Analysis), versão 1.0, e foi feita também a validação do
algoritmo da Análise Limite implementado. Naquela versão, o programa fazia uso
do otimizador Minos 5.5 para resolver o problema de otimização, onde somente
pequenos problemas puderam ser analisados. Neste trabalho, apresenta-se
GEOLIMA versão 2.0 com algumas melhoras, mas com um otimizador próprio,
mais eficiente, baseado no algoritmo de Pontos Interiores. A validação do
otimizador implementado é feita comparando-se os resultados obtidos com os
otimizadores comerciais LINGO(Lindo system inc., 1997) e MINOS (Murtagh e
Saunders, 1998).
Muitas pesquisas e publicações já foram realizadas na área da Análise
Limite Numérica, a grande maioria usando otimizadores comerciais como LINGO
(Lindo system inc., 1997), MINOS (Murtagh e Saunders, 1998) e LANCELOT
(Conn, Gould e Toin, 1992). Existem também vários trabalhos sobre Análise
Limite com uso do algoritmo de Pontos Interiores como Borges (1991), Zouain
(1993), Lyamin e Sloan (1993) e Farfan (2000); mas não se encontraram
aplicações com malhas suficientemente refinadas ou aplicações reais. Neste
trabalho propõe-se uma implementação alternativa do algoritmo de Pontos
Interiores, computacionalmente mais eficiente.
Sendo as propriedades dos materiais geotécnicos de natureza aleatória,
como complemento ao objetivo principal, neste trabalho aplica-se a Análise de
Confiabilidade com a Análise Limite pelo método FORM (First Order Reliability
Method).
Não sendo o objetivo deste trabalho expor as teorias da Análise Limite
(AL), Método de Elementos Finitos (MEF), Programação Matemática (PM) e da
Confiabilidade, no Capítulo 2, é apresentado somente um resumo de alguns
conceitos sobre a Análise Limite e sua formulação pelo MEF. No Capítulo 3, é
apresentada a solução numérica do problema da Análise Limite. No Capítulo 4,
são apresentados os fundamentos da Análise de Confiabilidade. No Capítulo 5,
são apresentados exemplos de aplicação da Análise Limite, realizados com o
programa GEOLIMA 2.0 e finalmente, no Capítulo 6, são apresentadas as
conclusões e sugestões do presente trabalho.
2 ANÁLISE LIMITE NUMÉRICA
O objetivo principal da Análise Limite é determinar a carga que leva uma
estrutura ao colapso (carga de colapso). As formulações existentes na Análise
Limite para o cálculo da carga de colapso são baseadas nos teoremas de Análise
Limite (inferior ou superior). A formulação pelo limite inferior (formulação
estática) deve satisfazer as condições de admissibilidade estática dos campos de
tensões; a formulação pelo limite superior (formulação cinemática) deve satisfazer
as condições de admissibilidade cinemática dos campos de velocidades; e,
finalmente, a formulação mista tem a forma da formulação pelo limite inferior,
com variáveis em tensões, mas o campo de tensões não é estaticamente
admissível, a não ser de forma aproximada.
A solução de problemas geotécnicos, então, pode ser feita por meio da
aplicação de qualquer uma das formulações estática, cinemática ou mista. Estas
formulações podem também ser classificadas em formulação forte e formulação
fraca. A formulação é forte quando satisfaz explicitamente as condições dos
teoremas da Análise Limite, e é fraca quando as condições são satisfeitas por meio
do princípio dos trabalhos virtuais.
No presente trabalho é usada a formulação mista fraca baseada no princípio
dos trabalhos virtuais para representar as equações de equilíbrio do sistema, o qual
estabelece que, para qualquer pequeno movimento cinematicamente admissível, o
trabalho das forças externas será igual ao trabalho das forças internas.
Não é o propósito deste trabalho descrever em detalhe toda a teoria da
Análise Limite, portanto, neste Capítulo é apresentado somente um resumo de
alguns conceitos importantes que permitem a compreensão das formulações
estabelecidas no Capítulo seguinte.
24
2.1. Teoremas da Análise Limite
Os teoremas fundamentais da Análise Limite foram apresentados por
Gvosdev (1960); Drucker et. al. (1952); Drucker, Greenberg and Prager (1952).
Chen (1975) apresentou estes teoremas para materiais com comportamento
plástico perfeito (Figura 2.1). Neste tipo de comportamento, as deformações
plásticas ocorrem sem que haja mudança nas tensões. Apesar de representar
apenas uma aproximação do comportamento real dos materiais, a hipótese pode
ser utilizada com segurança porque as deformações elásticas são desprezíveis em
relação às plásticas.
Os teoremas da Análise Limite são formulados em termos de
admissibilidade dos campos, conforme descrito a seguir.
2.1.1. Campos de Tensões Estaticamente Admissíveis
Define-se que um campo de tensões em um corpo tridimensional é
estaticamente admissível quando satisfaz as seguintes condições:
• Equilíbrio no volume.
• Condições de contorno.
• Critério de escoamento.
2.1.2. Campos de Velocidades Cinematicamente Admissíveis
Define-se que um campo de velocidades é cinematicamente admissível ou
geometricamente compatível quando:
• Satisfaz as condições de contorno em termos de velocidades.
• Satisfaz as condições de compatibilidade cinemática em termos de
deformações.
• Trabalho externo é igual à dissipação de energia interna.
25
2.1.3. Teorema de Limite Inferior
Se um campo de tensões correspondente a uma carga atuante (ou externa) é
estaticamente admissível, então a carga atuante é menor ou igual à carga de
colapso. A carga ou fator de colapso assim obtido será um limite inferior do fator
de colapso real.
2.1.4. Teorema de Limite Superior
Se um campo de velocidades é cinematicamente admissível, o colapso
ocorre quando o trabalho externo (forças externas) for igual ou maior ao trabalho
interno (forças internas). A carga ou fator de colapso assim obtido será um limite
superior ao fator de colapso real.
2.2. Considerações na Análise Limite
Algumas considerações importantes são apresentadas a seguir, para um
melhor entendimento dos conceitos utilizados na formulação pelo MEF.
2.2.1. Consideração de Plasticidade Perfeita
O solo exibe comportamento elástico para pequenas deformações.
Evidências de experimentos em laboratório mostram que deformações
permanentes ou irrecuperáveis ocorrem quando elas ultrapassam uma deformação
de tolerância. Este tipo de comportamento pode ser descrito por meio da Teoria da
Plasticidade.
A Análise Limite faz a consideração de que o material, no colapso,
comporta-se como perfeitamente plástico (Figura 2.1). Neste tipo de
comportamento, as características de endurecimento e amolecimento são
ignoradas.
26
σ
ε
Perfeitamente Plástico (Ideal)
Real
σ
ε
Perfeitamente Plástico (Ideal)
Real
Figura 2.1 – Relação tensão deformação para solo real e ideal (Chen, 1975).
2.2.2. Considerações sobre Escoamento
Para se caracterizar o comportamento de um material submetido a um
estado de tensões complexo, a Teoria de Plasticidade define uma superfície de
escoamento (Figura 2.2), caracterizada por uma função de escoamento. A
superfície de escoamento pode ser interpretada como: para um estado de tensões
dentro da superfície, só acontecem deformações elásticas; se o estado de tensões
alcança a superfície, ocorrem deformações plásticas e finalmente, os estados de
tensões acima da superfície de escoamento são excluídos.
Assim, os estados de tensões para o qual 0)( >ijF σ são excluídos,
enquanto que, 0)( <ijF σ implica em comportamento elástico e 0)( =ijF σ indica
que ocorre fluxo plástico.
27
Figura 2.2 – Superfície de escoamento no espaço de tensões principais.
Para um material perfeitamente plástico, a função de escoamento F
depende somente do conjunto de componentes de tensões σij e não de
componentes de deformação εij. Portanto, a função de escoamento é fixa no espaço
de tensões e o fluxo plástico ocorre quando a função de escoamento é igual a zero,
ou seja:
0)(F ij =σ (2.1)
Se o material é isotrópico, o escoamento plástico depende apenas da
magnitude das três tensões principais e não de suas direções. Neste caso, de uma
maneira geral, a superfície de escoamento pode ser expressa matematicamente em
função dos invariantes de tensões, ou seja:
0321 =)J,J,I(F (2.2)
onde:
I1 Primeiro invariante do tensor de tensões
J2 Segundo invariante do tensor desviador
J3 Terceiro invariante desviador
28
2.2.3. Considerações sobre a Lei de Fluxo
Diz-se que um fluxo plástico acontece quando o estado de tensões no espaço
das tensões (Figura 2.3) alcança a superfície de escoamento. A magnitude do
fluxo plástico total ( pijε ) nestas condições é ilimitado. Portanto é evidente que não
se pode falar nada sobre a deformação plástica total e por isso a descrição do
processo de fluxo plástico é feita pela cinemática de fluxo plástico, ou seja, por
meio das velocidades de deformação plástica ( pijε& ) em vez da deformação plástica
total ( pijε ). A velocidade de deformação total ( ijε& ) é composta por uma parte
elástica e outra plástica, que pode ser expressa da seguinte forma (Chen, 1975):
pij
eijij εεε &&& += (2.3)
onde:
ijε& Velocidade de deformação total
eijε& Velocidade de deformação elástica
pijε& Velocidade de deformação plástica
As velocidades de deformação elástica e as tensões estão relacionadas por
meio da lei de Hooke. Mas, como as velocidades de deformação elástica são
pequenas em relação às plásticas, elas são consideradas desprezíveis. Não entanto,
as velocidades de deformação plástica dependem do estado de tensões (Chen e
Liu, 1990).
Para a discussão sobre as velocidades de deformação plástica, é necessário
definir suas direções. Os eixos das coordenadas no espaço de tensões já
mencionados para a superfície de escoamento podem também ser usados para
representar simultaneamente velocidades de deformação plástica. Sendo assim,
cada eixo de tensões representará também um eixo de velocidades de deformação
plástica (Figura 2.3).
29
Figura 2.3 – Superfície de escoamento e vetor de deformação plástica.
Para materiais estáveis, nos quais o incremento de tensões e deformações
realiza um trabalho positivo, o vetor de velocidades de deformação plástica tem
direção e sentido das normais à superfície de escoamento (Figura 2.3). Portanto, a
velocidade de deformação plástica pode ser expressa por (Lancellotta, 1995):
ij
ijpij
)(Fσσ
λε∂
∂= && (2.4)
onde:
σij Campo de tensões
F(σij) Função de escoamento
λ& Fator de proporcionalidade (escalar) p
ijε& Velocidade de deformação plástica
A Equação 2.4 representa a lei de fluxo associada ou lei da normalidade, por
ser associada com a normal à superfície de escoamento do material.
30
2.3. Principio dos Trabalhos Virtuais
O princípio dos trabalhos virtuais pode ser usado para tratar problemas de
colapso de estruturas em materiais geotécnicos. Este princípio é uma expressão de
trabalho equilibrado e pode ser aplicado para materiais deformáveis como o solo.
Para a análise limite, o princípio dos trabalhos virtuais deve ser escrito em termos
do campo de velocidades virtuais e campo de tensões reais:
“se uma estrutura está em equilíbrio, o trabalho das forças externas sobre
um campo de velocidades virtuais dado, compatível com as condições de
fronteira, deve ser igual ao trabalho interno realizado pelas tensões sobre as
velocidades de deformações virtuais, compatíveis com o campo de velocidades
virtuais dado” (Lancellotta, 1995).
2.4. Critérios de Escoamento
Sendo de interesse conhecer as tensões de escoamento para definir um
campo de tensões estaticamente admissível (condição do teorema de limite
inferior da Análise Limite), neste trabalho são mostrados dois critérios de
escoamento que, mediante equações matemáticas, definem o lugar geométrico
onde se dá o escoamento do material. De uma forma geral estas funções de
escoamento podem ser definidas em termos dos invariantes de tensões, ou seja:
3211 σσσ ++=I (2.5)
[ ])()()(61
1332212 σσσσσσ −+−+−=J (2.6)
[ ] 321212331
2232
213213 9
4)()()(91)(
272 σσσσσσσσσσσσσσσ ++++++−++=J (2.7)
onde:
I1 Primeiro invariante de tensor de tensões
J2 Segundo invariante do tensor de tensões desviador
J3 Terceiro invariante desviador
σ1 Tensão principal maior
31
σ2 Tensão principal intermediaria
σ3 Tensão principal menor
2.4.1. Critério de Mohr-Coulomb
De acordo com o critério de ruptura de Mohr-Coulomb, a função de
escoamento F que define a superfície de escoamento é expressa em sua forma
mais geral pela Equação 2.8 (Chen e Liu, 1990).
)cos()sin()3
cos(3
131)
3sin(),,( 21221 φφπθπθθ CJIJJIF +⎥⎦
⎤⎢⎣⎡ ++−+= (2.8)
e
⎥⎦
⎤⎢⎣
⎡= −
2/32
31
233sin
31
JJ
θ (2.9)
onde:
θ Ângulo de Lode
φ Ângulo de atrito
C Coesão
No espaço de tensões principais, o critério de escoamento de Mohr-
Coulomb é representado por uma pirâmide hexagonal irregular (Fig. 2.4).
32
Figura 2.4 – Superfície de escoamento, critério de Mohr-Coulomb.
Este critério apresenta algumas deficiências, como por exemplo:
• Em 2D, a envoltória no diagrama de Mohr-Coulomb é composta de linhas
retas, o que indica que o parâmetro φ (ângulo de atrito) não muda com a
pressão de confinamento. Esta aproximação é razoável somente para um
intervalo limitado de pressões de confinamento.
• Em 3D, a superfície de escoamento tem descontinuidades, as quais dificultam
seu uso em análise numérica. É, portanto, inadequada matematicamente para
seu uso na análise tridimensional de problemas geotécnicos (Chen e Liu,
1990).
Pela deficiência descrita no parágrafo anterior, no presente trabalho, este
critério é somente usado para análise de problemas em 2D.
Para o caso bidimensional (2D), a Equação 2.8 fica simplificada para:
)cos()sin(21),( 1221 φφ CIJJIF −−= (2.12)
O caso bidimensional (2D) do criterio de Mohr-Coulomb é representado no
espaço de tensões principais pela Figura 2.5.
33
ψ
3σ
1σ
0c )sin(1)cos(2
0 φφ
−= Cc
)sin(1)sin(1)tan( φ
φψ −+=
ψ
3σ
1σ
0c )sin(1)cos(2
0 φφ
−= Cc
)sin(1)sin(1)tan( φ
φψ −+=
Figura 2.5 – Critério de escoamento de Mohr-Coulomb 2D.
2.4.2. Critério de Drucker-Prager
Drucker e Prager (1952) propuseram uma extensão do critério de Von
Mises. Este critério é muito usado em inúmeras aplicações práticas. Em termos
dos invariantes das tensões, a função de escoamento F que define a superfície de
escoamento (Figura 2.6) é expressa por:
kIJkJIF −−= 1221 ),,,( αα (2.11)
onde:
1I Primeiro invariante do tensor de tensões
2J Segundo invariante do tensor de tensões desviador
k,α Parâmetros do material
34
Figura 2.6 – Critério de escoamento de Drucker & Prager.
As principais características deste critério são:
• O critério de ruptura é relativamente simples.
• Possui somente dois parâmetros, k e α . Estes parâmetros podem ser
facilmente aproximados a partir dos parâmetros de resistência (C e φ ) no caso
de materiais geotécnicos.
• A superfície é contínua e, por isso, matematicamente apropiada para uso em
aplicações tridimensionais (3D).
Os parâmetros k e α podem ser obtidos a partir de uma aproximação do
critério de Mohr-Coulomb para o qual são ajustados dois círculos diferentes, um
deles passando pelos vértices de compressão e outro passando pelos vértices de
tração.
A primeira aproximação é obtida fazendo o ângulo de Lode de °= 0θ no
critério de escoamento geral de Mohr-Coublomb (Equação 2.8), obtendo-se os
parâmetros k e α do círculo superior ou de compressão como mostradas pelas
seguintes expressões (Chen e Liu, 1990):
35
))sin((
)sin(φ
φα−
=33
2 (2.12)
))sin((
)cos(Ckφ
φ−
=33
6 (2.13)
A segunda aproximação é obtida para um valor de ângulo de Lode de
°= 60θ no critério de escoamento geral de Mohr-Coulomb, obtendo-se os
parâmetros k e α do círculo inferior ou de tração, ou seja (Chen e Liu, 1990):
))sin((
)sin(φ
φα+
=332 (2.14)
))sin((
)cos(Ckφ
φ+
=33
6 (2.15)
onde:
C,φ Parâmetros de resistência do material
No caso 2D, para os problemas em estado plano de deformação, os
parâmetros k e α do critério de Drucker-Prager podem ser aproximados pelas
seguintes Equações (Chen e Liu, 1990):
)(tan
)tan(
φ
φα2129 +
= (2.16)
)(tan
Ckφ2129
3
+= (2.17)
2.5. Formulação Numérica da Análise Limite pelo MEF
O objetivo principal da Análise Limite é determinar a carga que leva uma
estrutura ao colapso (carga de colapso), onde a carga de colapso é determinada
como produto da carga inicial por um fator de colapso. As formulações existentes
na Análise Limite para o cálculo da carga de colapso são baseadas nos teoremas
de Análise Limite (inferior ou superior). A formulação pelo limite inferior
36
(formulação estática) deve satisfazer as condições de admissibilidade estática dos
campos de tensões; a formulação pelo limite superior (formulação cinemática)
deve satisfazer as condições de admissibilidade cinemática do campo de
velocidades; e, finalmente, a formulação mista tem a forma da formulação pelo
limite inferior, com variáveis em tensões, mas o campo de tensões não é
estaticamente admissível a não ser de forma aproximada.
A solução de problemas geotécnicos pode então ser feita através da
aplicação de qualquer uma das formulações estática, cinemática ou mista. Estas
formulações, podem também ser classificadas em formulação forte e formulação
fraca. A formulação é forte quando satisfaz explicitamente as condições dos
teoremas da Análise Limite, e é fraca quando as condições são satisfeitas através
do princípio dos trabalhos virtuais.
No presente trabalho é usada a formulação mista fraca com o uso do
princípio dos trabalhos virtuais para representar as equações de equilíbrio estático.
2.5.1. Condição de Equilíbrio
Na formulação mista fraca a condição de equilíbrio é garantida pelo
princípio dos trabalhos virtuais. Este princípio estabelece que, para qualquer
campo de velocidades cinematicamente admissível, o trabalho das forças internas
é igual ao trabalho das forças externas, como mostra a Equação 2.19.
dSdVdVS
T
VV
T ∫∫∫ += 00 TuFuσε αδαδδ &&& (2.19)
onde:
0F Força de volume inicial
0T Força de superfície inicial
α Fator (escalar) que multiplica as forças iniciais
u&δ Campo de velocidades virtuais
ε&δ Campo de velocidades de deformação virtuais
σ Campo de tensões
37
Pela teoria de elementos finitos, os campos são interpolados a partir de valores
nodais:
Tu
TT Nuu && δδ = (2.20)
TTT Buε)
&& ˆδδ = (2.21)
σNσ ˆσ= (2.22)
onde:
σ Campo de tensões
σ Campo de tensões nodais
u& Campo de velocidades
u& Campo de velocidades nodais
ε& Campo de velocidades de deformação
uN Matriz de interpolação de velocidades
σN Matriz de interpolação de tensões
B)
Matriz de compatibilidade cinemática obtida de uN
Substituindo-se as equações (2.20), (2.21) e (2.22) em (2.19), tem-se que:
∫ ∫∫ +=V S
Tu
TTu
T
V
TT dSdVdV 00ˆˆˆˆ TNuFNuσNBu &&
)& δαδαδ σ (2.23)
Fatorizando:
0ˆˆ00 =
⎥⎥⎦
⎤
⎢⎢⎣
⎡⎟⎟⎠
⎞⎜⎜⎝
⎛+−⎟⎟
⎠
⎞⎜⎜⎝
⎛∫ ∫∫V S
Tu
Tu
V
TT dSdVdV TNFNNBu ασδ σ
)& (2.24)
Como 0ˆ ≠Tu&δ , deslocamentos virtuais arbitrários, a Eq. 2.24 fica como:
0ˆ 00 =⎟⎟⎠
⎞⎜⎜⎝
⎛+−⎟⎟
⎠
⎞⎜⎜⎝
⎛∫ ∫∫V S
Tu
Tu
V
T dSdVdV TNFNσNB ασ
) (2.25)
38
Fazendo:
dVV
TσNBB ∫=
) (2.26)
∫ ∫+=V S
Tu
Tu dSdV 00 TNFNb (2.27)
A Equação 2.25 pode ser escrita como:
bσB α=ˆ (2.28)
Note-se que, na Equação 2.28, as variáveis do problema a serem
determinadas são as tensões nodais e o fator de colapso. Quando os elementos são
considerados de tensão constante, a matriz σN é igual à matriz identidade I , o
que conduz a:
dVV
T∫= BB)
(2.29)
2.5.2. Condições de Contorno
As condições de contorno precisam ser fornecidas para a solução do
problema. Para nós com campos de velocidades prescritos 0ˆ == uu && δδ , então,
deve-se eliminar de B e b (Equação 2.28) as linhas correspondentes às
velocidades prescritas.
2.5.3. Condição de Escoamento
Uma das condições para que o campo de tensões seja estaticamente
admissível (requisito do teorema de limite inferior) é que ele satisfaça um critério
de escoamento (Equação 2.30).
0)( ≤σg (2.30)
39
2.5.4. Problema de Análise Limite
Das Equações 2.28 e 2.30, a solução do problema pela Análise Limite (AL)
consiste em determinar o campo de tensões correspondentes a um fator α, que
maximiza a carga que a estrutura pode suportar sem violar a condição de
escoamento do material. O problema a ser resolvido consiste em um problema de
otimização como é apresentado a seguir e cuja solução será mostrada no Capítulo
3 deste trabalho.
max. αα =)(f (2.31)
0)(. ≤σgst (2.32)
bBσ α= (2.33)
2.5.5. Elementos Finitos Implementados
No presente trabalho foram implementados dois tipos de elementos finitos:
o elemento quadrilateral de 4 nós em 2D e o elemento hexaédrico de 8 nós em 3D.
Estes elementos são mostrados nas Figuras 2.7 (a) e (b). As funções de
interpolação uN utilizadas para interpolar as velocidades são as mesmas funções
de interpolação utilizadas para interpolar deslocamentos na formulação
convencional do MEF em deslocamentos.
40
(a) (b)(a) (b) Figura 2.7 – Elemento finito: (a) quadrilateral (2D), (b) hexaédrico (3D).
)1)(1(41 ssrrN iiiu ++= i=1,..,4 (2.34)
)1)(1)(1(81 ttssrrN iiiiu +++= i=1,..,8 (2.35)
onde:
uiN Funções de interpolação de velocidades
r,s,t Coordenadas paramétricas
ri,si,ti Coordenadas nos pontos nodais paramétricas
Para se interpolar os campos de tensões usa-se uma matriz identidade
(Equação 2.36) que corresponde a se adotar um campo de tensões constante
dentro do elemento.
IN =σ (2.36)
A escolha desses campos de tensões foi baseada em estudos feitos em
trabalhos anteriores (Gonzaga, 1997) (Farfan, 2000) e (Carrión, 2004).
3 SOLUÇÃO NUMÉRICA DA ANÁLISE LIMITE
O problema matemático resultante da aplicação de teoria da Análise Limite (AL)
associado ao Método dos Elementos Finitos (MEF) Equações (2.28) e (2.30),
pode ser representado pelas Equações 3.1-3.3. Estas equações configuram um
problema de otimização e que pode ser resolvido usando técnicas de otimização.
Max αα =)(f (3.1)
0)(.. ≤xgts (3.2)
0=− αbxB (3.3)
onde:
ℜ∈α Fator de colapso nℜ∈x Campo de tensões
mℜ∈)(xg Função de escoamento pxnℜ∈B Matriz de equilíbrio
pℜ∈b Vetor de cargas nodais iniciais
Na tentativa de resolver eficientemente o problema da Análise Limite
expresso pelas Equações (3.1), (3.2) e (3.3), foram pesquisados e testados vários
otimizadores comerciais existentes baseados nas técnicas de programação
matemática e otimizadores baseados em algoritmos genéticos. O uso de
otimizadores baseados em algoritmos genéticos para resolver o problema expresso
pelas Equações (3.1), (3.2) e (3.3) não foi bem sucedido e os otimizadores
comerciais baseados nas técnicas de programação matemática mostraram ser
ineficientes para problemas de grande escala como mostram os resultados de
testes apresentados na seção 3.1. Finalmente, foi implementado um otimizador
específico para o tipo de problema da análise limite baseado em técnicas de
programação matemática (PM).
42
3.1. Otimizadores Matemáticos Testados
Os otimizadores baseados em técnicas de programação matemática que
foram pesquisados e testados são LINGO (Lindo system inc., 1997), MINOS
(Murtagh e Saunders, 1998) e LANCELOT (Conn, Gould e Toin, 1992).
Os problemas de otimização podem ser classificados como de pequena,
média e de grande escala. Segundo Murtagh e Saunders (1998) os problemas são
considerados de “pequena escala” quando o número de restrições lineares mais o
número de restrições não lineares são menores que 100, de “média escala” quando
o número de restrições são aproximadamente 1000 ou 2000 e de “grande escala”
quando o número de restrições são maiores que 5000. O objetivo deste trabalho é
resolver problemas da Análise Limite de grande escala, para os quais os
otimizadores comerciais são testados.
Para testar a eficiência dos otimizadores indicados, um problema geotécnico
de estabilidade de talude em 2D é formulado. A Figura 3.1 mostra a geometria, as
condições de contorno e as propriedades do material do problema a ser analisado.
Figura 3.1 – Problema para teste de otimizadores.
Para se ter uma idéia comparativa do desempenho dos otimizadores, sete
análises foram feitas mantendo constantes os seguintes parâmetros: geometria do
problema, condições de contorno, propriedades do material (Figura 3.1) e os
mesmos recursos de hardware (Pentium IV com processador Intel de 3.07 GHz e
memória RAM de 1.4 Gb).
43
Para aumentar a complexidade do problema de otimização, a geometria do
problema (Figura 3.1) foi modelada com malhas de Elementos Finitos com
número de elementos variáveis, como mostrado pela Figura 3.2.
Figura 3.2 – Malhas: (a) 28, (b) 64, (c) 126, (d) 225, (e) 360, (f) 500 e (g) 750 elementos.
3.1.1. Otimizador Lingo
O programa LINGO desenvolvido pelo Lindo System Inc. (Lindo system
inc., 1997) é um otimizador para resolver problemas de otimização não lineares.
A grande vantagem deste otimizador é a sua facilidade de uso já que o
problema pode ser formulado facilmente no seu editor de texto ou ser importado a
partir de um arquivo de texto. O uso do programa é fácil porque tanto a função
objetivo, as restrições lineares e não lineares são escritas em forma de equações
44
explícitas. Uma outra vantagem deste programa é o pouco número de parâmetros
a serem controlados pelo usuário.
Na Tabela 3.1 apresentam-se os resultados obtidos a partir dos testes
realizados para o problema da Figura 3.1 com malhas da Figura 3.2, onde C é a
coesão, φ ângulo de atrito e γ peso especifico do material; n é o número de
variáveis do problema, m é número de restrições não lineares e p número de
restrições lineares; Mem é a quantidade de memória usada pelo otimizador, iter é
o número de iterações, t é o tempo requerido para otimizar e α é o fator de
colapso obtido.
Análise C φ γ Elem. Nós n m p m+p Mem.(MB) α iter. t(seg)1 25 5 17 28 40 85 28 56 84 0.200 3.76732 93 02 25 5 17 64 81 193 64 128 192 0.450 3.58688 245 23 25 5 17 126 150 379 126 252 378 0.875 3.52012 482 134 25 5 17 225 256 676 225 450 675 1.674 3.49650 857 1105 25 5 17 360 399 1081 360 720 1080 2.790 3.47713 1492 7756 25 5 17 500 546 1501 500 1000 1500 3.906 3.46887 3524 23057 25 5 17 750 806 2251 750 1500 2250 5.668 0.027280 43 2453
MalhaMaterial LingoProblema
Tabela 3.1 – Testes realizados com o programa Lingo.
Nas Figuras 3.3a, 3.3b, 3.3c e 3.3d ilustram-se os resultados obtidos com o
otimizador LINGO, onde, a Figura 3.3(b) indica a variação do número de
iterações requeridas pelo otimizador com o incremento de complexidade do
problema (número de elementos da malha). A Figura 3.3(c) indica o desempenho
do otimizador em termos de tempo em segundos com o incremento de
complexidade do problema. A Figura 3.3(a) indica a memória requerida pelo
otimizador com o incremento de complexidade do problema. A Figura 3.3(d)
indica o fator de colapso obtido com o refinamento da malha.
Os resultados indicam que o uso da memória cresce linearmente com a
complexidade do problema, enquanto que o número de iterações tem uma
variação quase linear até malhas com 35 elementos aproximadamente e logo
aumenta fortemente. O resultado mais importante é o desempenho do otimizador
em termos de tempo, a Figura 3.3c indica que o tempo de otimização cresce
exponencialmente com o número de elementos da malha. A Figura 3.3d mostra a
variação do fator de colapso com o refinamento da malha, este comportamento era
de se esperar, porque equações de equilíbrio obtidas por malhas menos refinadas
45
obtidas por elementos finitos tendem a fornecer modelos mais rígidos e
resistentes, superestimando as cargas de colapso. Isso explica porque o fator de
colapso se aproxima pelo limite superior ao fator de colapso real como acontece
na formulação pelo limite superior.
Uso da Memória - LINGO
0
1
1
2
2
3
3
4
4
5
0 100 200 300 400 500 600
Elementos da Malha
Série1
Mem
oria
(MB
)
Figura 3.3a – Variação da memória usada pelo otimizador LINGO.
Iterações - LINGO
0
500
1000
1500
2000
2500
3000
3500
4000
0 100 200 300 400 500 600
Elementos da Malha
Série1
Itera
ções
Figura 3.3b – Variação de número de iterações do otimizador LINGO.
46
Desempenho - LINGO
0
500
1000
1500
2000
2500
0 100 200 300 400 500 600
Elementos da Malha
Série1
Tem
po(s
eg)
Figura 3.3c – Desempenho do otimizador LINGO.
Fator de Colapso - LINGO
3.45
3.50
3.55
3.60
3.65
3.70
3.75
3.80
0 100 200 300 400 500 600
Elementos da Malha
Série1
Fato
r de
Col
apso
Figura 3.3d – Variação do fator de colapso obtido com LINGO.
Dos testes realizados com otimizador LINGO conclui-se que, para
problemas pequenos com malhas de elementos finitos com número de elementos
de ate 200, a eficiência deste otimizador é aceitável como pode ser observado
pelas análises 1, 2, 3 e 4 da Tabela 3.1, mas, para malhas com número de
elementos maiores de 200, o tempo de otimização cresce exponencialmente
(Figura 3.3c). Da Tabela 3.1, pode-se observar que para malha com 750 elementos
este otimizador converge a uma solução não ótima, portanto, é de se esperar que o
número de iterações e o tempo empregado não sejam reais e por isso a análise 7
não é usada para comparações.
47
3.1.2. Otimizador Minos
MINOS é um programa de otimização para resolver problemas de
otimização lineares e não-lineares, desenvolvido na universidade de Stanford por
Bruce A. Murtagh e Michael A. Saunders (Murtagh e Saunders, 1998).
O programa MINOS apresenta uma maior dificuldade para seu uso em
comparação ao programa LINGO. Entre as dificuldades que apresenta MINOS
estão a quantidade de parâmetros, 94 no total, que o usuário precisa controlar, uso
de arquivo com formato MPS (Mathematical Programming System) para a
entrada de dados das restrições lineares e um arquivo com formato SPECS
(Specification) para especificar a configuração dos parâmetros.
O teste com este otimizador foi feito para o problema apresentado pela
Figura 3.1, com malhas de Elementos Finitos apresentadas na Figura 3.2 e nas
mesmas condições que os testes realizados com o otimizador LINGO.
Os resultados dos testes realizados com o otimizador MINOS são
apresentados na Tabela 3.2 e nas Figuras 3.4a, 3.4b, 3.4c e 3.4d são ilustrados
graficamente estes resultados.
Análise C φ γ Elem. Nós n m p m+p Mem(MB) α iter. t(seg)1 25 5 17 28 40 85 28 56 84 253 3.767329 225 12 25 5 17 64 81 193 64 128 192 253 3.586883 599 43 25 5 17 126 150 379 126 252 378 253 3.566153 1602 64 25 5 17 225 256 676 225 450 675 253 3.497155 3556 165 25 5 17 360 399 1081 360 720 1080 253 3.478615 6561 726 25 5 17 500 546 1501 500 1000 1500 253 3.470797 10182 2027 25 5 17 750 806 2251 750 1500 2250 295 3.453981 18067 671
MalhaMaterial MinosProblema
Tabela 3.2 – Testes realizadas com o otimizador Minos.
48
Uso da Memória - MINOS
250255260265270275
280285290295300
0 100 200 300 400 500 600 700 800
Elementos da Malha
Série1
Mem
ória
(MB
)
Figura 3.4a – Variação da memória usada pelo otimizador MINOS.
Número de Iterações - MINOS
02000400060008000
10000
1200014000160001800020000
0 100 200 300 400 500 600 700 800
Elementos da Malha
Série1
Itera
ções
Figura 3.4b – variação de número de iterações do otimizador MINOS.
49
Desempenho - MINOS
0
100
200
300
400
500
600
700
800
0 100 200 300 400 500 600 700 800
Elementos da Malha
Série1
Tem
po(s
eg)
Figura 3.4c – Desempenho do otimizador MINOS.
Fator de Colapso - MINOS
3.40
3.45
3.50
3.55
3.60
3.65
3.70
3.75
3.80
0 100 200 300 400 500 600 700 800
Elementos da Malha
Série1
Fato
r de
Col
apso
Figura 3.4d – Variação do fator de colapso obtido com MINOS.
Os resultados dos testes realizados mostram que este otimizador faz uma
alocação prévia de 253 MB de memória como apresentado nas análises 1, 2, 3, 4,
5 e 6 da Tabela 3.2 e somente quando a aplicação precisa de mais memória, o
programa faz uma realocação de memória como mostrado na análise 7 da Tabela
3.2.
Este otimizador apresenta um melhor desempenho que o otimizador
LINGO, mas para malhas com mais 300 elementos o tempo de desempenho
também cresce quase exponencialmente (Figura 3.4c).
50
O teste também mostra a mesma variação do fator de colapso com o
incremento de número de elementos da malha obtidos com o otimizador LINGO
(Figura 3.4d).
3.1.3. Otimizador Lancelot
O otimizador LANCELOT (Linear And Nonlinear Constrained Extended
Lagrangian Otimization Techniques), segundo o manual do usuário (Conn, Gould
e Toin, 1992) é um programa para resolver problemas de otimização não linear de
grande escala.
Este programa foi desenvolvido pelos seguintes institutos de pesquisa:
Facultés Universitaires Notre-Dame de la Paix, Namur – Bélgica; Harwell
Laboratory, UK; IBM, USA; Rutherford Appleton Laboratory, UK e a University
of Waterloo, Canada.
Este programa tem muitos parâmetros a serem controlados pelo usuário e
para a especificação de dados utiliza arquivos com formato SIF (Standard Input
Format). Este programa apesar de ser multi-plataforma não apresenta suporte ou
não tem uma versão para o sistema operacional Windows. Os sistemas
operacionais para os quais este programa apresenta suporte são:
CRAY Unicos
DEC OSF1
DEC Ultrix
DEC VMS
HP HP-UX
IBM AIX
IBM CMS
IBM DOS (Salford and Waterloo Fortran)
Linux g77
Silicon Graphics IRIX
Sun SunOS and Solaris
51
No processo a pesquisa e adequação do código de LANCELOT para rodar
no sistema operacional Windows, encontrou-se um trabalho de testes de eficiência
e comparação realizados entre o programa LANCELOT e MINOS. O trabalho é
denominado “A Numerical Comparison Between the LANCELOT and MINOS
Package for Large-Scale Nonlinear Optimization” (Bongartz et all., 1997).
Na Tabela 3.3 apresentam-se os resultados de testes realizadas por Bongartz
e outros (Bongartz et all., 1997), onde LP indica problemas lineares, QP
problemas quadráticos, BC problemas sem restrição ou com restrições de
contorno simplesmente, LC problemas com restrições lineares e NC problemas
com restrições não lineares. A primeira coluna indica número total de problemas
testados e as colunas LANCELOT e MINOS indicam o tempo total acumulado
em segundos empregados pelos otimizadores.
Tabela 3.3 – Comparação da eficiência entre Lancelot e Minos (Bongartz et all, 1997).
Para os objetivos do presente trabalho de resolver eficientemente o
problema da Análise Limite com restrições lineares e não lineares, os resultados
mais importantes da Tabela 3.3 são as duas últimas linhas, onde se mostra que
para problemas com restrições lineares (LC) o programa MINOS tem uma melhor
eficiência com um tempo acumulado de 219 segundos contra 1526 segundos do
programa LANCELOT e para problemas com restrições não lineares (NC) o
programa MINOS também tem uma melhor eficiência com um tempo acumulado
de 1093 segundos contra 15754 segundos de LANCELOT.
Pelos resultados apresentados no trabalho citado era de se esperar que não
seria obtida nenhuma melhora em termos da eficiência (tempo de otimização) com
o uso do programa LANCELOT e se decidiu não seguir com a pesquisa sobre o
uso deste otimizador.
52
3.1.4. Comparação de Desempenho de Otimizadores Testados
Na Figura 3.5, apresenta-se uma comparação do desempenho dos
otimizadores pesquisados e testados. A Figura 3.5a mostra que o otimizador
LINGO faz um melhor uso da memória do computador, mas o resultado mais
importante é que MINOS tem um melhor desempenho em termos de tempo
empregado para otimizar o problema (Figura 3.5c). Ambos os otimizadores
mostraram a mesma variação do fator de colapso com o número de elementos da
malha (Figura 3.5d).
Uso da Memória
0
50
100
150
200
250
300
350
0 100 200 300 400 500 600 700 800
Elementos da Malha
LingoMinos
Mem
ória
(MB
)
Figura 3.5a – Comparação de uso da memória entre Lingo e Minos
Número de Iterações
02000400060008000
100001200014000160001800020000
0 100 200 300 400 500 600 700 800
Elementos da Malha
LingoMinos
Itera
ções
Figura 3.5b – Comparação de número de Iterações entre Lingo e Minos
53
Desempenho
0
500
1000
1500
2000
2500
0 100 200 300 400 500 600 700 800
Elementos da Malha
Lingo
Minos
Tem
po(s
eg)
Figura 3.5c – Comparação de desempenho entre Lingo e Minos
Fator de Colapso
3.40
3.45
3.50
3.55
3.60
3.65
3.70
3.75
3.80
0 100 200 300 400 500 600 700 800
Elementos da Malha
LingoMinos
Fato
r de
Col
apso
Figura 3.5d – Comparação de variação de fator de colapso, obtidos com Lingo e Minos.
Da pesquisa e dos testes realizados com os otimizadores comerciais
existentes conclui-se que o otimizador LINGO pode ser usado eficientemente para
problemas da Análise Limite de pequena escala e o otimizadores MINOS poderia
ser usado eficientemente para problemas da Análise Limite de até média escala.
A Análise Limite de problemas em 3D tem uma complexidade maior se
comparada com problemas em 2D porque o número de variáveis e o número de
restrições é muito maior. Portanto, o uso destes otimizadores pesquisados é quase
inviável. Em geral para problemas de Análise Limite de grande escala o uso
destes otimizadores não é recomendável.
54
3.2. Otimizador Implementado
Após a tentativa de usar otimizadores comerciais como Lingo e Minos e
sendo o desempenho alcançado não satisfatório para resolver problemas
geotécnicos pela Análise Limite com malhas suficientemente refinadas ou
problemas de grande escala, optou-se neste trabalho por desenvolver um
otimizador exclusivo para o tipo de problema da Análise Limite. Assim, nesta
pesquisa, implementou-se o algoritmo de Pontos Interiores. Sendo o objetivo da
tese a aplicabilidade prática da Programação Matemática para resolver
eficientemente o problema da análise limite, neste trabalho são expostos apenas
alguns conceitos que permitem o entendimento do algoritmo implementado.
3.2.1. Condições de Otimalidade
As condições de otimalidade, conhecidas também como as condições de
Karush, Kuhn e Tucker (Condições KKT), são o resultado teórico mais importante
no campo da programação não linear. Estas condições devem de ser satisfeitas
para a solução ótima de qualquer problema. Estas condições são as bases para o
desenvolvimento de muitos algoritmos computacionais e proporciona um critério
de parada para muitos outros (Castillo et all, 2002).
Para )(xgS ∇= , as condições de otimalidade para o problema de
otimização expresso pelas equações (3.1), (3.2) e (3.3), são escritas como segue:
0BµSλ =+ tt (3.4)
1−=bµ t (3.5)
0bxB =− α (3.6)
0)( =xgλ t (3.7)
0≥jλ (3.8)
0)( ≤xg (3.9)
55
onde: mℜ∈λ e pℜ∈µ são os multiplicadores de Lagrange associados
respectivamente às restrições (3.3) e (3.2). As condições (3.6) e (3.9) são
conhecidas como condições de factibilidade primal, a condição (3.7) é conhecida
como condição de complementaridade e a condição (3.8) indica a não
negatividade do multiplicador de Lagrange e é conhecida como condição de
factibilidade dual.
3.2.2. Algoritmo de Pontos Interiores
Este algoritmo foi proposto por Herskovits (1986 e 1989) (Farfan, 2000).
Ele já foi usado em vários trabalhos para resolver problemas da Análise Limite
como em Borges (1991), Zouain e outros (1993), Lyamin e Sloan (1997), e Farfan
(2000).
Os algoritmos são chamados de Pontos Interiores por que todo ponto gerado
deve de estar na região viável, ou seja, para cada iteração deve ser preservada a
factibilidade das restrições de desigualdade 0)( ≤xg e 0≥λ . A condição 0≥λ é
garantida pelo critério de atualização utilizado ao final de cada iteração. Portanto,
deve-se garantir a admissibilidade plástica 0)( ≤xg .
Na literatura pesquisada encontraram-se dois algoritmos de pontos
interiores; o primeiro, conhecido como Newton e Técnica de Deflexão e o
segundo conhecido como Newton e Técnica de relaxação-contração. Neste
trabalho, propõe-se uma outra técnica, denominada Vetorial.
Como discutido na seção 3.2.4 os três algoritmos se diferenciam na técnica
utilizada para garantir a factibilidade das desigualdades. Sendo o objetivo deste
trabalho implementar um otimizador eficiente, os três algoritmos foram
implementados e testados.
3.2.3. Inicialização
A inicialização das variáveis primárias (Equações 3.10 e 3.11) é feita
considerando que o ponto inicial deve satisfazer as equações de equilíbrio e a
admissibilidade plástica do material. Para satisfazer a condição de otimalidade
56
(Equação 3.8) a variável dual λ deve ser positiva, portanto, esta variável é
iniciada conforme sugerido por Borges (1991).
00 =x (3.10)
00 =α (3.11)
)(1
00 xλ
g−= (3.12)
3.2.4. Direção de Busca
Para o cálculo da direção de busca, o algoritmo de Newton-Rapson é
aplicado nas equações de otimalidade (3.4-3.7) e para )]([ 2 kk xgλH ∇= , ki
ki λxgG /)(= (matriz diagonal), kkk
x xxd −= +1 e kkkd ααα −= +1 são obtidas as
equações incrementais:
0µBλSdH =++ ++ 11 kkkx (3.13)
11 −=+kt µb (3.14)
0=− kkx dαbdB (3.15)
0λGdS =+ +1kkx (3.16)
As Equações (3.13)-(3.16) formam um sistema de quatro blocos de
equações com quatro grupos de incógnitas onde as variáveis ou incógnitas a serem
determinadas são kxd , 1+kλ , 1+kµ e kdα . A solução deste sistema de equações é
discutida na seção 3.2.8
Para o cálculo da direção de busca encontraram-se na literatura duas
técnicas e uma terceira técnica é proposta neste trabalho.
A primeira técnica implica o cálculo de ângulo de deflexão e o sistema de
equações é resolvido duas vezes (Zouain e outros, 1993). Na segunda técnica,
fatores de relaxação-contração são usados para determinar a direção viável e o
sistema de equações é resolvido somente uma vez (Borges, 1991). Uma terceira
técnica vetorial é proposta neste trabalho, onde a direção viável é determinada
mediante álgebra vetorial a partir do cálculo do ângulo de deflexão e o sistema de
equações é resolvido uma única vez.
57
3.2.4.1. Técnica de deflexão
Nesta técnica, um ângulo de deflexão θ é calculado a partir da solução do
sistema de equações (3.13-3.16). O ângulo de deflexão é calculado de tal forma a
garantir uma direção viável para kdα e kxd como mostrada pela Equação 3.17,
onde ]1,0[∈β é um parâmetro prefixado (Zouain et al,1993). Uma demonstração
completa para o cálculo de θ é apresentada por Borges (1991) e Zouain e outros
(1993).
∑ ⎟⎠⎞
⎜⎝⎛−
=+
ki
ki
kd
λλ 1
)1( αβθ (3.17)
O sub-índice i indica que o somatório é calculado para todos os elementos
da malha.
Na Figura 3.6 o ângulo de deflexão θ é ilustrado graficamente para uma
restrição )(xg , onde a direção 0d é uma direção não viável calculada a partir da
solução do sistema de equações 3.13-3.17.
Figura 3.6 – Técnica de deflexão.
Após o cálculo do ângulo de deflexão θ , a direção de busca viável 1d é
calculada. Para isso o sistema de equações (3.13-3.16) é resolvido novamente mas
desta vez considerando o ângulo de deflexão como mostrado pelas Equações
(3.18-3.21).
0µBλSdH =++ ++ 11 kkkx (3.18)
11 −=+kt µb (3.19)
58
0=− kkx dαbdB (3.20)
eλGdS θ−=+ +1kkx (3.21)
onde e é um vetor com todos os componentes unitários.
3.2.4.2. Técnica de relaxação-contração
Neste algoritmo, a factibilidade das tensões é garantida propondo-se uma
contração uniforme para as tensões e o fator de colapso. O procedimento consiste
em incrementar as tensões e o fator de colapso a partir da solução das Equações
(3.13-3.16), relaxada de um fator r . Em seguida calcula-se um fator de contração
c para garantir a factibilidade das tensões. A contração do fator de colapso pelo
mesmo fator c garante que a condição de equilíbrio (3.15) seja preservada ao
final da iteração, enquanto o fator de relaxação r é calculado para garantir a
condição de factibilidade na iteração.
Nesta técnica o novo ponto é calculado como mostrado na Equação 3.22
sem calcular direção de busca viável nem o comprimento de passo.
)(1 dxx rc kk +=+ (3.22)
Para ilustrar esta técnica apresenta-se a Figura 3.7, onde a direção 0d é
calculada a partir da solução das Equações (3.13-3.16) e kx são as variáveis da
iteração anterior. Esta técnica consiste em que para um fator de relaxação r o
vetor 0dAM r= é calculado, a seguir o vetor AMOAOM += é calculado e
finalmente o vetor OMON c= é calculado. Este processo é realizado para
diferentes valores do fator de relaxação r até que se cumpra a seguinte condição
de acréscimo kk αα >+1 . O fator de contração c neste processo é determinado
pela busca linear na direção do vetor OM .
59
Figura 3.7 – Técnica de Relaxação-Contração.
3.2.4.3. Técnica Vetorial - Proposto
Pelo fato de que, resolver duas vezes o sistemas de equações como proposto
pela primeira técnica (de deflexão) é computacionalmente custoso e isto afeta
diretamente o desempenho do otimizador implementado, neste trabalho propõe-se
uma técnica vetorial para o cálculo da direção de busca viável. Esta técnica é
descrita a seguir.
Da solução do sistema de Equações (3.13-3.16) a direção tangente não
viável kxdd =0 é determinada como mostra a Figura 3.8, a seguir um ângulo de
deflexão é calculado como expresso pela Equação 3.17 e finalmente a direção
viável 1d é calculada vetorialmente mediante álgebra vetorial, como a soma de
dois vetores (Figura 3.8), onde um vetor unitário na direção negativa do gradiente
é determinado pela Equação 3.23, a seguir a norma euclidiana do vetor oposto ao
ângulo de deflexão é calculada pela Equação 3.24 e logo a direção de busca viável
1d é obtida como a soma de dois vetores conforme Equação 3.25.
g(x)
=0
Região
viáve
l
g(x)<
0
g(x)
>0
1+kx
kx
0d
1d
)(xg∇−
tû
θ
g(x)
=0
Região
viáve
l
g(x)<
0
g(x)
>0
1+kx
kx
0d
1d
)(xg∇−
tû
θ
Figura 3.8– Técnica Vetorial - Proposta.
60
)()(ˆ
xgxgu
∇∇−
= (3.23)
)tan(0 θd=t (3.24)
udd ˆ01 t+= (3.25)
onde )(xg∇ , é o vetor gradiente com n elementos formados pelos gradientes das
restrições )(xig de cada elemento, u é um vetor unitário no sentido oposto ao
vetor gradiente e t é o módulo do vetor oposto ao ângulo de deflexão θ .
3.2.5. Comprimento de Passo
O comprimento de passo indicado na Figura 3.9 é determinado somente
para a primeira e a terceira técnica.
Determinada a direção de busca viável 1d na etapa anterior, é necessário
determinar o comprimento de passo. As restrições não lineares neste trabalho são
apresentadas pelos critérios de escoamento (Equações 2.8 e 2.9). Por sugestão de
Herskovits (2008), estas equações foram transformadas na sua forma quadrática
equivalente e substituindo na Equação 3.26 é obtida uma equação quadrática em
função de s como expresso na Equação 3.27, onde ),( 1 φdfa = ,
),,,( 1 φCfc k dx= e ),,,( fk Cfq γφx= . Finalmente, o comprimento de passo s
é obtido resolvendo a equação quadrática 3.27.
0)(.)( 1 =−+ kf
k gsg xdx γ (3.26)
02 =++ qscsa (3.27)
onde fγ é um parâmetro (Borges, 1991), C e φ são as propriedades do material.
61
Figura 3.9 – Comprimento de passo s .
3.2.6. Atualização das Variáveis
Calculada a direção viável e o comprimento de passo nas etapas anteriores
as variáveis x e α são atualizadas conforme Equações 3.28 e 3.29.
Graficamente, este processo de atualização é ilustrado na Figura 3.10, onde se vê
que o 1+kx é atualizado como a soma dos vetores kx e 1ds .
kx
kk sdxx +=+1 (3.28)
kkk ds ααα +=+1 (3.29)
Figura 3.10 – Atualização da variável x .
Pela condição de Karush-Kuhn-Tucker (Equação 3.8) a variável dual λ tem
que ser sempre positiva. A atualização desta variável é feita da seguinte forma:
62
),max(1∞
+ = λγλλ λki
ki (3.30)
onde },...,,max{ 21km
kk λλλλ =∞ e λγ é um parâmetro (Borges, 1991).
3.2.7. Teste de Desempenho dos Algoritmos
Para ilustrar o desempenho dos algoritmos, um talude 2D como indicado na
Figura 3.11 é analisado. As análises são feitas para as mesmas condições, ou seja,
considerando-se as mesmas propriedades do material coesão (C=25 kN/m²),
ângulo de atrito (φ = 5°) e peso específico (γ = 17 kN/m³); a mesma geometria
com o mesmo número de elementos da malha; o mesmo critério de escoamento,
Mohr-Coulomb; e usando o mesmo computador (Pentium IV com processador de
3.07 GHz e memória RAM de 1.4 Gb).
Na Tabela 3.4, apresentam-se as comparações de fator de colapso, número
de iterações, tempo de otimização e o mecanismo de ruptura obtido com cada um
dos algoritmos.
Figura 3.11 Geometria do problema para teste de algoritmos (malha de 25 elementos).
63
Algoritmo α Iter t(seg) Superficie de Falha
Deflexão 3.717250 13 1.040
Relax.-Cont. 3.717627 18 2.250
Vetorial 3.668298 13 0.687
Tabela 3.4 – Desempenho dos algoritmos implementados
3.2.8. Manipulação de sistemas lineares a serem resolvidos
Para resolver os sistemas de Equações (3.13-3.16) e (3.18-3.21), na
literatura pesquisada foram encontradas duas formas de reduzir estes sistemas
para logo serem resolvidos por algum método direto ou indireto, onde o sistema a
ser resolvido tem a mesma ordem que o número de restrições lineares do
problema. A primeira, é uma manipulação matricial global como mostrado por
Lyamin & Sloan (1997) ou por Farfan (2002). A segunda é uma manipulação
matricial por elementos como mostrado por Borges (1991) ou Zouain e outros
(1993). Neste trabalho, é proposta uma terceira forma de resolver estes sistemas
mediante uma solução direta sem manipulação prévia.
Para 0λλ =k , 0xx =k e 0λλ =k , as equações (3.13), (3.14), (3.15) e (3.16)
podem ser re-escritas como as Equações (3.31), (3.32), (3.33) e (3.34)
representado um sistema de 4 blocos de equações com 4 grupos variáveis ( kxd ,
1+kλ , 1+kµ e kαd ).
0µBλSHd =++ ++ 11 kkTkx (3.31)
01 =+ +kkx λGdS (3.32)
64
0=− kkx dαbdB (3.33)
11 −=+kt µb (3.34)
Neste trabalho foram implementadas e testadas as três formas de resolver o
sistema de Equações 3.31-3.34 descritas acima, tendo como objetivo a busca de
um melhor desempenho do otimizador implementado. As vantagens e
desvantagens de cada uma das formas de resolver são descritas a seguir e
finalmente é apresentada uma comparação de desempenho testado.
3.2.8.1. Manipulação matricial global
Nesta forma de manipulação as matrizes H , S , G e B são calculadas
para todo o sistema, onde nxnℜ∈H , nxmℜ∈S , mxmℜ∈G e nxpℜ∈B ; n é o
número de variáveis do problema, m é o número de restrições não lineares e p é
o número de restrições lineares. O vetor kxd pode ser explicitado na Equação 3.31
gerando a Equação 3.35. Substituindo-se em seguida kxd da Equação 3.35 na
Equação 3.32 e considerando que 1−= HSQ e GSQW −= T a Equação 3.36 é
obtida.
)( 111 ++− +−= kkTkx µBλSHd (3.35)
111 +−+ −= kTk uQBWλ (3.36)
Para QWQHD 11 −− −= T e substituindo a Equação 3.36 em 3.35, a Equação
3.37 é obtida. Substituindo a Equação 3.37 na Equação 3.33, e para TBDBK = e 11 )( −+= kk dαuv , a Equação 3.38 é obtida.
1+−= kTkx uDBd (3.37)
bKv −= (3.38)
A matriz K pode ser chamado de uma pseudo matriz de rigidez já que ela
relaciona forças aplicadas e velocidades.
A Equação 3.38 é um sistema de equações lineares de ordem igual ao
número de restrições lineares ( p ) e pode ser resolvido por qualquer dos métodos
discutidos na seção 3.2.9.
65
Conhecido o vetor v e substituindo kk dαvu =+1 na Equação 3.34, a variável
kdα é determinada e, logo, 1+ku é conhecido também. Substituindo 1+ku nas
Equações 3.37 e 3.36 as variáveis kxd e 1+kλ são calculados.
O processo de solução do sistema de Equações 3.18, 3.19, 3.20 e 3.21
considerando o ângulo de deflexão θ são similares e pode ser encontrado em
Lyamin e Sloan (1997) ou Farfan (2000).
Observa-se que o uso deste processo implica realizar operações como a
obtenção da inversa de matriz de ordem n (número de variáveis do problema),
assim como o produto de matrizes e produto de matrizes por vetores. Sendo estas
operações computacionalmente caras é de se esperar que esta implementação seja
ineficiente como ilustrado com os testes realizados.
3.2.8.2. Manipulação matricial por elementos
Uma segunda forma de manipulação na solução de sistema de Equações
(3.31), (3.32), (3.33) e (3.34) é apresentada por Borges (1991) e por Zouain e
outros em (1993). Nesta forma, a matriz K da Equação 3.38 é obtida pela
montagem de matrizes de rigidez elástico-plástica de cada elemento finito.
Nesta forma de manipulação as variáveis H , S , Q , W , epD e eK são
calculadas para cada elemento finito como indicado nas equações a seguir:
)]([ 2 xgλH ∇= (3.39)
)(xgS ∇= (3.40)
kλxgG /)(= (3.41)
1−= HSQ (3.42)
GSQW −= T (3.43)
QWQHD 11 −− −= Tep (3.44)
BDBK epTe = (3.45)
onde, para os elementos finitos implementados, quadrilateral em 2D e hexaédrico
em 3D, H é uma matriz de 3x3 para problemas em 2D e de 6x6 para problemas
em 3D, S e Q são vetores de 3 elementos para problemas em 2D e de 6
66
elementos para problemas em 3D, G e W são matrizes de 1x1, epD é a matriz
conhecida como módulo elástico-plástico de 3x3 para problemas em 2D e de 6x6
para problemas em 3D, B é a matriz de equilíbrio para cada elemento finito
obtido pelo MEF com dimensões de 8x3 para problemas 2D e de 24x6 para
problemas em 3D e eK é a matriz de rigidez elástico-plástica de 8x8 para
problemas em 2D e de 24x24 para problemas em 3D.
Da Equação 3.31, a variável kxd para cada elemento pode ser expressa como
mostra a Equação 3.46. Substituindo (3.46) em (3.32), a Equação 3.47 é obtida.
)( 111 ++− +−= kkTkx µBλSHd (3.46)
111 +−+ −= kTk uQBWλ (3.47)
Substituindo a Equação 3.47 em 3.46 e para módulo elástico-plástico epD
expresso pela Equação 3.44 a Equação 3.48 é obtida. Substituindo a Equação 3.48
na Equação 3.33, para epK expresso pela Equação 3.45 e para kekαdvu =+1 a
Equação 3.49 é obtida.
1+−= kTepkx uBDd (3.48)
eee bvK = (3.49)
A partir da matriz de rigidez elástico-plástica eK de cada elemento (Equação
3.49), a matriz de rigidez elástico-plástico K para todo o sistema é montada. Esta
montagem é similar à montagem de matriz de rigidez global da análise
convencional por Elementos Finitos. Uma vez montada a matriz K e sendo b
vetor dos carregamentos nodais montado também a partir dos carregamentos de
cada elemento, a Equação 3.49 é escrita para todo o sistema como indicado em
Equação 3.50.
bKv = (3.50)
A Equação 3.50 é similar à Equação 3.38 e pode ser resolvida por algum
dos métodos discutidos na seção 3.2.9. Determinado o vetor v as variáveis kdα e
1+ku são determinadas como expressas pelas equações 3.51 e 3.52.
67
1).( −= vbd kα (3.51)
vu kk dα=+1 (3.52)
Sendo o vetor 1+ku um vetor global para todo o sistema, este vetor é
decomposto facilmente para cada um dos elementos e logo os vetores kxd e 1+kλ
para cada elemento são calculados pelas Equações 3.48 e 3.47. Este mesmo
processo também é seguido para resolver as Equações de 3.31 a 3.34, mas
considerando o ângulo de deflexão θ . Uma abordagem completa pode ser
encontrada em Borges (1991) e Zouain e outros (993).
Esta forma de manipulação para resolver o sistema de Equações (3.31-3.34)
é mais eficiente que o primeiro porque faz um melhor uso da memória do
computador e os cálculos como a inversa da matriz, o produto de matrizes e o
produto de matrizes por vetores são feitos no nível de cada elemento e não para
todo o sistema.
3.2.8.3. Solução direta sem manipulação - proposta
Como uma alternativa às duas formas anteriores de resolver o sistema de
Equações (3.31-3.34), neste trabalho propõe-se resolver diretamente este sistema
de equações como mostra a Equação 3.53.
O sistema de Equações (3.18-3.21) ou (3.31-3.34) considerando o ângulo de
deflexão θ pode ser resolvido também como indicado na Equação 3.54.
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
=
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢
⎣
⎡
−− +
+
10
1
1
000
µλd
b00b00B
00GS0BSH
k
k
k
kx
t
tt
dα
(3.53)
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧−
=
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢
⎣
⎡
−− +
+
10
1
1
0e
0
µλd
b00b00B
00GS0BSH
θ
αk
k
k
kx
t
tt
d
(3.54)
68
Esta forma de resolver os sistemas de Equações (3.31-3.34) tem vantagens e
desvantagens. A vantagem é que não é necessário realizar nenhum tipo de cálculo
ou manipulação prévia à solução do sistema. A desvantagem é que o sistema de
equações a ser resolvido é de ordem de )1( +++ pmn , onde ( n +1) é o número
de variáveis, m é o número de restrições não lineares de desigualdade e p é o
número de restrições lineares de igualdade; isso indicaria que é necessário realizar
mais esforço computacional em termos de FLOPS (FLoating-point OPerationS) e
também uso de mais recursos do computador como a memória para armazenar a
matriz de coeficientes.
Sendo a matriz dos coeficientes uma matriz esparsa, as desvantagens
indicadas no parágrafo anterior como número de FLOPS e requerimento de mais
memória para armazenar a matriz de coeficientes são resolvidas neste trabalho
mediante o uso de uma técnica de tratamento de matrizes esparsas como vetores,
como discutido na seção 3.3.1.
3.2.8.4. Teste de Desempenho das Manipulações
Neste trabalho, as três formas de manipular e resolver o sistema de
Equações (3.31-3.34) foram implementadas e testadas. A eficiência das formas de
manipular e resolver os sistemas de equações é medida em termos de quantidade
de memória requerida e de tempo empregado na manipulação e solução do
sistema.
Os testes foram feitos utilizando cada uma das manipulações para o
problema apresentado na Figura 3.1 com malhas apresentadas na Figura 3.2. Nos
testes, para facilitar a identificação de cada método, a manipulação matricial
global é simplesmente denominada como “Matriz”, a manipulação matricial por
elementos é denominada como “Elementos” e a solução direta sem manipulação é
denominada como “Direta”.
a) Teste de manipulação matricial global
A Tabela 3.5 apresenta os resultados dos testes realizados com a
implementação mediante a manipulação matricial global. Nas Figuras 3.12 e 3.13
estes resultados são apresentados graficamente para uma melhor interpretação.
69
Como se pode verificar nas Figuras 3.12 e 3.13 tanto a quantidade de
memória requerida como o tempo empregado para resolver o sistema cresce
exponencialmente neste método. Devido ao fato de que estes valores têm um
crescimento exponencial os testes somente foram feitos para malhas de 28, 64 e
126 elementos.
Da Figura 3.13 pode-se observar também que entre 80 a 90 % de tempo
total é empregado na manipulação do sistema e de 10 a 20 % de tempo é
empregado na solução do sistema resultante.
Mem.Elem. Nós n m p m+p (MB) Manip. Resolv. total
28 40 85 28 56 84 0.54 0.125 0.015 0.14064 81 193 64 128 192 1.50 0.750 0.250 1.000126 150 379 126 252 378 47.00 6.844 1.172 8.016225 256 676 225 450 675360 399 1081 360 720 1080500 546 1501 500 1000 1500750 806 2251 750 1500 2250
Matriz
Tempo(seg)MalhaManipulação Problema
Tabela 3.5 – Teste da manipulação matricial global.
Uso da Memória
05
101520253035404550
0 100 200 300 400 500 600 700 800
Elementos da Malha
Matriz
Mem
ória
(MB
)
Figura 3.12 – Memória requerida pela manipulação matricial global.
70
Desempenho
0
1
2
3
4
5
6
7
8
9
0 100 200 300 400 500 600 700 800
Elementos da Malha
Manip.
Resolv.
Total
Tem
po(s
eg)
Figura 3.13 – Tempo empregado na manipulação e solução do sistema.
b) Teste de manipulação matricial por elementos
A Tabela 3.6 apresenta os resultados dos testes realizados com a
implementação mediante a manipulação matricial por elementos. Nas Figuras 3.14
e 3.15 estes resultados são apresentados graficamente para uma melhor
interpretação.
Como se pode observar nas Figuras 3.14 e 3.15, esta forma de resolver o
sistema de equações tem um melhor desempenho que a técnica anterior.
A Figura 3.14 mostra que o uso de memória tem um crescimento
exponencial suave, mas a Figura 3.15 mostra um crescimento exponencial em
relação ao tempo requerido para resolver o sistema.
A Figura 3.15 mostra também que 90 % do tempo é empregado na
manipulação, sendo o restante do tempo empregado na solução do sistema
resultante.
Mem.Elem. Nós n m p m+p (MB) Manip. Resolv. total
28 40 85 28 56 84 0.45 0.000 0.015 0.01564 81 193 64 128 192 0.81 0.078 0.031 0.109
126 150 379 126 252 378 2.00 0.547 0.140 0.687225 256 676 225 450 675 7.00 3.044 0.438 3.482360 399 1081 360 720 1080 16.00 13.828 1.125 14.953500 546 1501 500 1000 1500 29.00 50.281 4.812 55.093750 806 2251 750 1500 2250
Elementos
Tempo(seg)MalhaManipulação Problema
Tabela 3.6 – Teste da manipulação matricial por elementos.
71
Uso da Memória
0
5
10
15
20
25
30
35
0 100 200 300 400 500 600 700 800
Elementos da Malha
Elementos
Mem
ória
(MB
)
Figura 3.14 – Memória requerida para a solução do sistema.
Desempenho
0
10
20
30
40
50
60
0 100 200 300 400 500 600 700 800
Elementos da Malha
Manip.
Resolv.
Total
Tem
po(s
eg)
Figura 3.15 - Tempo empregado na manipulação e solução do sistema.
c) Teste de solução direta sem manipulação - proposta
Os resultados dos testes realizados com esta forma de resolver os sistemas
de equações são apresentados na Tabela 3.7 e nas Figuras 3.16 e 3.17 estes
resultados são apresentados graficamente.
Como esta forma de resolver o sistema de equações não tem manipulação, a
coluna Manip da Tabela 3.7 tem todos os valores iguais a zero, sendo portanto o
tempo total igual ao tempo empregado pelo resolvedor.
72
A Figura 3.16 indica que o requerimento de memória tem um crescimento
linear com o número de elementos da malha. A Figura 3.17 mostra que o
desempenho em termos de tempo requerido para resolver o sistema tem uma
tendência quase linear com pequenas oscilações.
Estes resultados indicam que esta técnica tem muito melhor desempenho,
tanto em relação ao uso da memória quanto em relação ao tempo requerido para
resolver o sistema.
Mem.Elem. Nós n m p m+p (MB) Manip. Resolv. total
28 40 85 28 56 84 0.45 0 0.0930 0.09364 81 193 64 128 192 0.60 0 0.1250 0.125
126 150 379 126 252 378 0.90 0 0.2340 0.234225 256 676 225 450 675 1.40 0 0.3750 0.375360 399 1081 360 720 1080 2.00 0 0.6870 0.687500 546 1501 500 1000 1500 2.70 0 0.9370 0.937750 806 2251 750 1500 2250 4.00 0 1.5620 1.562
Direta
Tempo(seg)Malha ProblemaManipulação
Tabela 3.7 – Teste da solução direta sem manipulação.
Uso da Memória
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
0 100 200 300 400 500 600 700 800
Elementos da Malha
Direta
Mem
ória
(MB
)
Figura 3.16 - Memória requerida para a solução do sistema.
73
Desempenho
0.0
0.2
0.40.6
0.8
1.0
1.21.4
1.6
1.8
0 100 200 300 400 500 600 700 800
Elementos da Malha
Manip.
Resolv.
Total
Tem
po(s
eg)
Figura 3.17 – Tempo empregado na solução do sistema.
d) Comparação das técnicas de solução
Nesta seção apresenta-se uma comparação dos testes realizados com as três
técnicas implementadas para resolver sistemas de equações. A Tabela 3.8
apresenta um resumo dos testes realizados com cada uma das técnicas e as Figuras
3.18 e 3.19 apresentam as comparações do desempenho em termos de uso da
memória e de tempo total requerido para manipular e resolver o sistema de
equações.
A Figura 3.18 mostra que a solução direta proposta neste trabalho é muito
mais eficiente em relação ao uso da memória do computador. A Figura 3.19
mostra também que a solução direta proposta neste trabalho é muito mais eficiente
em comparação às duas técnicas encontradas na literatura pesquisada. O resultado
mais importante mostrado pela Figura 3.19 é que o desempenho (em tempo) com
a solução direta proposta, cresce quase linearmente com o refinamento da malha,
este fato influencia diretamente no desempenho do otimizador implementado.
74
Mem.Elem. Nós n m p m+p (MB) Manip. Resolv. total
28 40 85 28 56 84 0.54 0.125 0.015 0.14064 81 193 64 128 192 1.5 0.750 0.250 1.000126 150 379 126 252 378 47.00 6.844 1.172 8.016225 256 676 225 450 675360 399 1081 360 720 1080500 546 1501 500 1000 1500750 806 2251 750 1500 225028 40 85 28 56 84 0.45 0.000 0.015 0.01564 81 193 64 128 192 0.81 0.078 0.031 0.109126 150 379 126 252 378 2.00 0.547 0.140 0.687225 256 676 225 450 675 7.00 3.044 0.438 3.482360 399 1081 360 720 1080 16.00 13.828 1.125 14.953500 546 1501 500 1000 1500 29.00 50.281 4.812 55.093750 806 2251 750 1500 225028 40 85 28 56 84 0.45 0 0.093 0.09364 81 193 64 128 192 0.60 0 0.125 0.125126 150 379 126 252 378 0.90 0 0.234 0.234225 256 676 225 450 675 1.40 0 0.375 0.375360 399 1081 360 720 1080 2.00 0 0.687 0.687500 546 1501 500 1000 1500 2.70 0 0.937 0.937750 806 2251 750 1500 2250 4.00 0 1.562 1.562
Elementos
Direta
Matriz
Tempo(seg)MalhaManipulação Problema
Tabela 3.8 – Comparação de resultados com as três técnicas de solução.
Uso da Memória
05
101520253035404550
0 100 200 300 400 500 600 700 800
Elementos da Malha
Matriz
Elementos
Direta
Mem
ória
(MB
)
Figura 3.18 – Memória usada pelas técnicas.
75
Desempenho
0
2
4
6
8
10
12
14
16
0 100 200 300 400 500 600 700 800
Elementos da Malha
Matriz
Elementos
Direta
Tem
po(s
eg)
Figura 3.19 – Tempo empregado pelas técnicas.
3.2.9. Resolvedores Implementados
A eficiência do otimizador implementado depende muito da solução
eficiente de sistema de equações lineares, por isso neste trabalho foram
implementados e testados diferentes algoritmos diretos e indiretos que permitem
resolver sistemas lineares da forma bAx = .
Os primeiros testes do otimizador implementado foram feitos para resolver
problemas pequenos com malha de 8 e 25 elementos usando o algoritmo de Gauss
Jordan (GJ) (Hoffman, 1992 e Tsao, 1989). Nos testes para aplicações com
malhas refinadas, o otimizador implementado usando Gauss-Jordan mostrou-se
ineficiente, por isso, foram implementados e testados outros métodos de solução
de sistemas lineares.
Entre os métodos diretos, o algoritmo de fatoração tLL. de Cholesky (Ames,
2000 e Alkire, 2002) e fatoração LU de Doolittle (Gómez e Burguest, 2004)
foram implementados e testados. Entre os métodos indiretos o método Quase
Newton BFGS(Broyden-Fletcher-Goldfarb-Shanno) QN-BFGS (Bathe, 1996) e o
método dos Gradientes Conjugados (CG) (Sandoval, 2006), foram implementados
e testados também.
Uma das vantagens dos métodos indiretos é que são adequados para aplicar
técnicas como tratamento da matriz esparsa e pré-condicionadores. Por isso,
76
foram implementados também o algoritmo dos Gradientes Conjugados com
tratamento da matriz esparsa (CSR-CG) e o método Quase Newton BFGS com
tratamento da matriz esparsa (CSR-QN-BFGS).
Dos testes realizados com os métodos implementados, o método dos
Gradientes Conjugados com tratamento da matriz esparsa apresentou melhor
desempenho, portanto neste trabalho somente é descrito o algoritmo dos
Gradientes Conjugados. Os métodos diretos e o método BFGS foram
implementados com base nas referências indicadas e pelo fato desses métodos
serem amplamente conhecidos e difundidos na literatura, eles não foram
apresentados no presente trabalho.
3.2.9.1. Método dos gradientes conjugados (CG)
O método dos Gradientes Conjugados, introduzido pela primeira vez por
(Hestenes e Stifel, 1952), é uma técnica de otimização e tornou-se um dos mais
populares métodos para a solução de sistemas lineares da forma bAx = , onde o
problema é resolvido como um problema de minimização de uma função objetiva
quadrática )(xf (Equação 3.55). A condição de primeira ordem para o mínimo de
uma função )(xf impõe que o gradiente da função objetivo bAxxf −=∇ )( seja
igual a zero, como mostra a Equação (3.56). Este método minimiza a função
)(xf , onde o valor de x que minimiza a função é a solução do sistema bAx = .
xbAxxxf tt −=21)( (3.55)
0)( =−=∇ bAxxf (3.56)
O algoritmo de Gradiente Conjugado pré-condicionado (Figura 3.20),
implementado no presente trabalho é a apresentado no trabalho de Sandoval em
2006 (Sandoval, 2006), onde M é a matriz pré-condicionadora.
77
Figura 3.20 – Algoritmo de Gradiente Conjugado pré-condicionado (Sandoval, 2006).
78
3.2.9.2. Teste de desempenho de resolvedores implementados
Os testes de eficiência dos métodos implementados foram feitos para os
problemas mostrados pela Figura 3.21. Na Tabela 3.9 são apresentados os
resultados dos testes dos métodos com melhor desempenho. Estes resultados são
mostrados graficamente na Figura 3.22. Da Tabela 3.9 pode-se observar que o
método Quase Newton BFGS mostra-se ineficiente. Por isso não foram realizados
mais testes e não foram apresentados resultados com esse algoritmo na Figura
3.22.
(b)
(d)
(f)
(a)
(c)
(e)
(b)
(d)
(f)
(b)
(d)
(f)
(a)
(c)
(e)
(a)
(c)
(e) Figura 3.21 – Malhas: (a) 8, (b) 25, (c) 64, (d) 150, (e) 400 e (f) 676 elementos.
79
A Figura 3.22 mostra que o método baseado na decomposição LU tem um
melhor desempenho para sistemas com número de variáveis N menores que 2400
(aproximadamente), mas para N maiores que 2400 o método de eliminação de
Gauss-Jordan (GJ) tem melhor eficiência.
Técnica de tratamento de matriz esparsa foi aplicada aos métodos indiretos,
para o qual o método de Gradiente Conjugado (CG) apresenta um ganho de
eficiência como discutido na seção 3.3.1. Pré-condicionadores também foram
implementados e aplicados no método de Gradiente Conjugado para melhorar a
eficiência como descrito na seção 3.3.2.
Malha N GJ LU QN_BFGS CG
a 47 0.000 0.000 0.0 0.0b 151 0.063 0.032 1.0 0.0c 385 0.829 0.547 5.0 5.0d 901 10.219 9.656 55.0 9.0e 2381 284.672 185.984 994.0f 6460 4854.000 6492.000 5958.0
Tempo (seg)
Tabela 3.9 – Comparação do desempenho dos métodos implementados.
Desempenho dos Métodos
0
1000
2000
3000
4000
5000
6000
7000
0 1000 2000 3000 4000 5000 6000 7000
N (No de variaveis do sistema)
GJ
LU
CG
T(seg)
Figura 3.22 – Comparação do desempenho dos métodos implementados.
80
3.2.10. Resolvedor SAMG Testado
Na procura por encontrar um método eficiente para resolver os sistemas
lineares na implementação do otimizador, neste trabalho tentou-se usar o
resolvedor comercial SAMG(Algebraic Multigrid Methods for System)
desenvolvido pelo Institute of Algorithms and Scientific Computation (Klaus &
Tonja, 2005). Este programa usa o Método Multigrid para resolver sistema
lineares.
A tentativa de resolver o problema expresso pelas Equações 3.53 e 3.54
usando este programa não foi bem sucedida.
3.3. Melhora do desempenho
Como o desempenho do otimizador implementado depende do desempenho
dos algoritmos que servem para resolver o sistema de equações lineares, neste
trabalho foram pesquisadas, implementadas e testadas técnicas como tratamento
da matriz esparsa e uso de pré-condicionadores. Estas técnicas são descritas a
seguir.
3.3.1. Tratamento de Matriz Esparsa
Uma matriz esparsa é aquela que apresenta muitos elementos iguais a zero,
e elementos diferentes de zero podem ser armazenados em uma estrutura de dados
especial ou em vetores (Tsal, 1998).
Existem várias técnicas ou formatos para tratamento de uma matriz esparsa
com vetores, neste trabalho usou-se o formato CSR (Compressed Sparse Row)
para armazenar a matriz de coeficientes em vetores.
3.3.1.1. Formato CSR(Compressed Sparse Row)
Neste formato, uma matriz esparsa é representada por 3 vetores como é
mostrado na Figura 3.23a (SMAILBEGOVIC et all, 2006), onde AN é um vetor
com os valores diferentes de zero da matriz A , AJ é um vetor com os índices das
81
colunas da matriz A e AI são os limites dos índices de colunas do vetor AN
para cada fila de A .
Neste trabalho, a matriz de coeficientes das Equações (3.53) e (3.54) são
armazenados em vetores no formato CSR. A escolha de este formato foi porque é
adequado para realizar o produto de uma matriz por um vetor e sua
implementação é simples demais como mostrado pelo código na Figura 3.23b.
Figura 3.23a – Armazenamento da matriz esparsa (SMAILBEGOVIC et all, 2006).
Figura 3.23b – Produto de uma matriz esparsa A por um vetor d .
3.3.1.2. Teste de desempenho de CG com tratamento da matriz esparsa
A Tabela 3.10 apresenta os resultados dos testes feitos com os métodos de
Gradiente conjugado CG e o método Quase-Newton BFGS com tratamento da
matriz esparsa pelo formato CSR. Observa-se que com o método Quase Newton
BFGS obteve-se uma melhora na eficiência de 50% aproximadamente e com o
método de Gradiente Conjugado CG a melhora de eficiência foi ainda maior.
Na Figura 3.24 são apresentados graficamente os resultados dos testes, onde
se pode visualizar o ganho obtido pelo método de Gradiente Conjugado com
tratamento da matriz esparsa CSR-CG. É de se esperar que este ganho influencie
diretamente na eficiência do otimizador implementado.
82
Malha N GJ LU QN_BFGS CSR+QN_BFGS CG CSR_CG
a 47 0.000 0.000 0.0 0.00 0.0 0.000b 151 0.063 0.032 1.0 0.00 0.0 0.000c 385 0.829 0.547 5.0 3.00 5.0 0.600d 901 10.219 9.656 55.0 27.00 9.0 1.400e 2381 284.672 185.984 994.0 12.000f 6460 4854.000 6492.000 5958.0 359.000
Tempo (seg)
Tabela 3.10 – Desempenho dos métodos com o tratamento da matriz esparsa.
Desempenho de CG com CSR
0
1000
2000
3000
4000
5000
6000
7000
0 1000 2000 3000 4000 5000 6000 7000
N (No de variaveis do sistema)
GJ
LU
CG
CSR_CG
T(seg)
Figura 3.24 – Desempenho do método CG com tratamento da matriz esparsa CSR.
3.3.2. Precondicionamento
A idéia do precondicionamento é transformar um sistema linear bAx = em
outro equivalente com condições espectrais mais favoráveis, onde o número de
iterações requeridas para a convergência é reduzido (Cervantes & Mejía, 2004).
Pré-condicionadores são técnicas utilizadas para acelerar a convergência da
solução de sistemas lineares (Equação 3.57) pelo método de Gradiente Conjugado
(CG) e, portanto aumenta o desempenho deste método em termos de tempo.
bAx = (3.57)
83
Neste trabalho foram implementados e testados três pré-condicionadores
encontrados na literatura e o uso de pré-condicionadores mistos proposto neste
trabalho.
Os pré-condicionadores implementados e testados são Escala Diagonal(DS)
(Pini e Gambolati, 1990), Escala Simétrica (SS) (Jennings e Malik, 1978) e
Fatoração Incompleta de Cholesky (ICF) (Meijerink e Van der Vorst, 1977). Os
seguintes pré-condicionadores mistos são propostos neste trabalho DS+SS e
ICF+SS. É importante mencionar que os pré-condicionadores mistos não
aparecem na literatura e foi uma tentativa deste trabalho de encontrar um pré-
condicionador eficiente.
Campos (1999) indica que, “o precondicionamento implica em alterar a
matriz de coeficientes do sistema, fazendo com que os autovalores desta matriz
sejam mais próximos e, logo, o sistema mais estável, reduzindo o número de
iterações necessárias para a solução do mesmo.
3.3.2.1. Escala Diagonal (DS)
A escala diagonal é um pré-condicionador muito simples, fácil de ser obtido
e implementado. A matriz pré-condicionadora M apresentada no algoritmo da
Figura 3.20 é uma matriz diagonal cujos elementos são formados pela diagonal da
matriz de coeficientes, como mostrado na Equação 3.58 (Pini e Gambolati, 1990).
)(AdiagM = (3.58)
3.3.2.2. Escala Simétrica (SS)
Este pré-condicionador foi apresentado por Jennings e Malik (1978) e
consiste em escalonar o sistema linear apresentado pela Equação 3.57 como
mostra a Equação 3.59. A matriz pré-condicionadora D é uma matriz diagonal
cujos elementos são obtidos como indicado na Equação 3.60.
DbxDADD =−1 (3.59)
84
5.0)( −= iiii ad (3.60)
Da Equação 3.59 para DADB = , xDy 1−= e Dbc = a Equação 3.52 é
transformada para a Equação 3.61. A Equação 3.61 é resolvida pelo método do
Gradiente Conjugado (Figura 3.20) com a matriz IM = .
cBy = (3.61)
Conhecido o vetor y , o vetor x é facilmente determinado pela relação
Dyx = .
3.3.2.3. Fatoração Incompleta de Cholesky (ICF)
A fatoração incompleta de Cholesky foi proposta por Meijerink e Van der
Vorst (1977). Este pré-condicionador basea-se na decomposição de Cholesky para
resolver sistemas lineares (Equação 3.57), onde a matriz de coeficientes A é
decomposta como o indicado na Equação 3.62.
LLA T= (3.62)
onde L é uma matriz triangular inferior. Esta decomposição na sua forma
original, requer do cálculo da raiz quadrada e dos elementos da diagonal, para
evitar esse cálculo, Meijerink e Van der Vorst (1977) apresentam o algoritmo de
Cholesky modificado como mostra a Equação 3.63.
DLLA T= (3.63)
onde D é uma matriz diagonal e as matrizes L e D são obtidas como mostrado
nas expressões:
iiii
i
k kkikjkjiii
LD
NijDLLAL
1
,/)1(
=
=−= ∑ −
(3.64)
85
No processo de decomposição, a matriz A perde o padrão de distribuição de
zeros, ou seja, elementos nulos podem tornar-se não nulos após a decomposição.
Para que a fatoração seja incompleta deve-se manter a mesma distribuição de
zeros da matriz original. Isto é conseguido fazendo-se com que os elementos nulos
na matriz original permaneçam nulos na matriz fatorizada.
Benzi & Tuma em 2001 indicam que uma fatoração incompleta pode falhar
para uma matriz geral de tipo esparsa positiva e definida (SPD) devido à
ocorrência de pivôs não positivos (Benzi & Tuma, 2001).
Na literatura pesquisada encontraram-se trabalhos indicando que o pré-
condicionador ICF tem o melhor desempenho, entre eles Pinheiro (1998), Campos
(1999).
3.3.2.4. Pré-condicionadores mistos - proposto
Neste trabalho propõe-se o uso de pré-condicionadores mistos para melhorar
o desempenho do método do gradiente conjugado na solução de sistemas lineares.
Na tentativa de encontrar um pré-condicionador eficiente para os tipos de
problema a serem resolvidos neste trabalho, foram implementados e testados dois
tipos de pré-condicionadores mistos. A escala diagonal misturada com escala
simétrica (DS+SS) e a escala diagonal misturada com fatoração incompleta de
Cholesky (DS+ICF).
Os pré-condicionadores mistos são fáceis de ser implementados. Para o pré-
condiconador DS+SS, primeiro é aplicado o pré-condicionador SS e logo o DS.
Similarmente para o pré-condicionador DS+ICF, primeiro aplica-se o ICF e logo
o DS.
Para o tipo de problema deste trabalho, o tipo de pré-condicionador que
mostrou o melhor desempenho é o pré-condicionador misto proposto DS+SS
como indicam os resultados dos testes realizados (Figura 3.25).
3.3.2.5. Teste de Desempenho de Pré-condicionadores
A Tabela 3.11 apresenta os resultados obtidos com os pré-condicionadores
implementados. Na Figura 3.25 estes resultados são ilustrados graficamente, onde
86
pode-se observar que o método do gradiente conjugado com pré-condicionador
misto DS+SS teve melhor desempenho para nosso tipo de problema.
Malha N CSR+CG DS+CSR+CG SS+CSR+CG ICF +CSR+CG DS+SS+SCR+CG ICF+SS+SCR+CG
a 47 0.0 0.000 0.000 0.000 0.000 0.000b 151 0.0 0.016 0.016 0.046 0.015 0.047c 385 0.6 0.063 0.063 0.437 0.063 0.375d 901 1.4 0.188 0.109 0.297 0.157 0.281e 2381 12.0 4.781 3.109 1.719 2.172 1.407f 6460 359.0 101.607 104.297 444.734 82.266 409.437
Temp(seg)
Tabela 3.11 – Teste de pré-condicionadores implementados.
Desempenho dos Precondicionadores
0
50
100
150
200
250
300
350
400
450
500
0 1000 2000 3000 4000 5000 6000 7000
N (No de variaveis do sistema)
CSR+CG
DS+CSR+CG
SS+CSR+CG
ICF +CSR+CG
DS+SS+SCR+CG
ICF+SS+SCR+CG
T(seg)
Figura 3.25 – Desempenho de CG com os pré-condicionadores implementados.
87
3.4. Teste de desempenho do Otimizador Implementado
Com o objetivo de ter uma idéia comparativa do desempenho do otimizador
implementado, o mesmo problema em 2D da Figura 3.1 com malhas da Figura 3.2
foi analisado usando o otimizador implementado e os resultados são comparados
com os otimizadores Lingo e Minos. Outro problema em 3D com resultados da
análise já conhecidos (obtidos usando o otimizador MINOS 5.5) é analisado e a
solução é obtida com o otimizador implementado.
3.4.1. Teste com problema em 2D
O problema da Figura 3.1 com malhas da Figura 3.2 que com o qual os
programas Lingo e Minos foram testados, foi analisado utilizando o otimizador
implementado.
Os resultados obtidos com o programa Geolima versão 2.0 cujo otimizador
foi implementado neste trabalho são apresentados na Tabela 3.12. Para fim de
comparação, estes resultados são apresentados graficamente nas Figuras 3.26,
3.27, 3.28 e 3.29, onde se observa que o programa Geolima com otimizador
implementado é mais eficiente tanto no uso de memória quanto no tempo de
otimização.
Elem. Nós Mem α iter. t(seg) Mem α iter. t(seg) Mem α iter. t(seg)28 40 0.20 3.76732 93 0 253 3.767329 225 1 0.45 3.766621 23 2.42264 81 0.45 3.58688 245 2 253 3.586883 599 4 0.60 3.584825 24 5.484
126 150 0.88 3.56612 482 13 253 3.566153 1602 6 0.90 3.560147 31 11.92225 256 1.67 3.49650 857 110 253 3.497155 3556 16 1.40 3.492299 39 27.48360 399 2.79 3.47713 1492 775 253 3.478615 6561 72 2.00 3.471841 27 27.97500 546 3.91 3.46887 3524 2305 253 3.470797 10182 202 2.70 3.463164 28 62.67750 806 5.67 0.027280 43 2453 295 3.453981 18067 671 4.00 3.447660 22 92.72
GeolimaMinosMalha Lingo
Tabela 3.12 – Resultados do teste em 2D do otimizador implementado.
88
Uso da Memória
0
50
100
150
200
250
300
350
0 100 200 300 400 500 600 700 800
Elem entos da Malha
Lingo
Minos
Geolima
Mem
ória
(MB)
Figura 3.26 – Comparação de uso da memória pelos otimizadores.
Número de Iterações
02000400060008000
100001200014000160001800020000
0 100 200 300 400 500 600 700 800
Elementos da Malha
Lingo
Minos
GeolimaItera
çõe
s
Figura 3.27 – Comparação de número de iterações.
89
Desempenho
0
500
1000
1500
2000
2500
0 100 200 300 400 500 600 700 800
Elementos da Malha
Lingo
Minos
Geolima
Tem
po(s
eg)
Figura 3.28 – Comparação do desempenho dos otimizadores.
Fator de Colapso
3.40
3.45
3.50
3.55
3.60
3.65
3.70
3.75
3.80
0 100 200 300 400 500 600 700 800
Elementos da Malha
Lingo
Minos
Geolima
Fat
or d
e C
olap
so
Figura 3.29 – Variação de fator de colapso.
90
3.4.2. Teste com problema em 3D
Com a finalidade de validar o otimizador implementado, um problema 3D
de frente de escavação de túneis (Figura 3.30) é analisado pelo programa
GEOLIMA 2.0 com seu próprio otimizador implementado. Este problema é
interessante, porque se tem resultados da simulação física do comportamento de
frente de escavação, feitas mediante ensaios de laboratório em modelos 3D em
escala reduzida. A simulação física foi feita no Laboratório da Mitsubishi Heavy
Industries Ltda., Japão (Sterpi et al., 1996). Este problema foi usado também para
validar os resultados da Análise Limite com GEOLIMA 1.0 implementada na
dissertação de mestrado e que usava o otimizador Minos 5.5.
Neste trabalho, a análise deste problema é feita com GEOLIMA 2.0 e os
resultados são comparados com os obtidos com o otimizador MINOS.
A análise foi feita mantendo todas as condições, ou seja para a mesma
geometria (Figura 3.30), a mesma malha (Figura 3.31) (676 elementos com 945
nós), as mesmas condições de contorno, as mesmas propriedades do material
coesão (C=50 kN/m²), ângulo de atrito (φ = 5°) e peso especifico (γ = 19.5
kN/m³); o mesmo critério de escoamento Drucker-Prager (parâmetros de
aproximação do círculo superior); com um carregamento inicial de oγ =1.0 kN/m³
e usando o mesmo computador (Pentium IV com dois processadores de 3.07 GHz
e memória RAM de 1.4 Gb).
O problema a ser resolvido pelo otimizador tem um total de 4057 variáveis,
1727 restrições lineares e 676 restrições não lineares. A ordem do sistema de
equações lineares a ser resolvido em cada iteração é de 6460.
O tempo requerido pelo Minos foi de 16 horas 43 minutos e o fator de
colapso obtido foi de 25.012=α .
O tempo requerido pelo otimizador implementado foi de 17 minutos com 49
seg. e o fator de colapso obtido foi de 24.6685=α . A Figura 3.32 apresenta o
mecanismo de ruptura do problema obtido pela Análise Limite com GEOLIMA
2.0 e é muito similar à obtida com Minos e também à obtida ao modelo físico
(Figura 3.33).
91
Figura 3.30 – Geometria da estrutura a ser analisada.
Figura 3.31 – Malha de elementos finitos (676 elementos e 945 nós).
92
Figura 3.32 – Mecanismo de colapso da estrutura obtido pelo programa GEOLIMA 2.0.
Figura 3.33 – Mecanismo de ruptura obtido a partir de ensaios em modelo físico em
escala reduzida (Sterpi,1996)
4 ANÁLISE DE CONFIABILIDADE COM ANÁLISE LIMITE
A avaliação da segurança das estruturas geotécnicas tem sido sempre um
dos objetivos da Engenharia Geotécnica. A forma convencional de quantificar a
segurança de uma estrutura foi e é feita ainda pela análise determinística, baseado
no conceito de fator de segurança. No entanto, o cálculo de fator de segurança
pela análise determinística pode ser muito conservador, pelo fato de que no seu
cálculo não são levadas em consideração as incertezas das variáveis como
carregamentos e propriedades do material.
Por outra parte, a confiabilidade é uma medida do nível de segurança de
uma estrutura que leva em consideração as incertezas das variáveis envolvidas
considerando-as como variáveis aleatórias.
Para materiais fabricados artificialmente como aço ou concreto é
perfeitamente possível estimar os parâmetros estatísticos de suas propriedades, a
partir de uma base de dados de experiências anteriores. No entanto, para materiais
de origem natural como o solo e rocha, cuja formação obedece a diferentes
processos físicos e químicos em diferentes condições de pressão e temperatura,
não é possível ter uma base de dados confiável, já que em cada lugar as
propriedades do material são completamente diferentes ou aleatórias.
Pelo fato de que as propriedades de solo são obtidas mediante ensaios de
campo ou de laboratório, incertezas de diferentes origens estão presentes nas
propriedades dos materiais assim obtidas o que torna a realização de uma análise
de confiabilidade necessária na análise de estruturas geotécnicas. Por outra parte,
conhecer os parâmetros estatísticos das propriedades dos materiais requer realizar
vários ou muitos ensaios de campo ou de laboratório, o qual inviabiliza o uso da
análise de confiabilidade em estruturas de pequeno porte, pelo custo econômico
que implica realizar estes ensaios. No entanto, para estruturas de grande porte e
importância, onde, uma possível falha ou colapso possa causar grandes perdas
econômicas, de vidas humanas, danos ecológicos e sociais, realizar a análise de
confiabilidade é perfeitamente possível, necessária e viável.
94
A confiabilidade é o complemento da probabilidade de falha e é avaliada,
então, a partir do cálculo da probabilidade de falha. Para problemas com função
de falha linear e dependente de poucas variáveis aleatórias com distribuição
normal, a probabilidade de falha pode ser calculada pelo método de integração
direta, como mostrado no trabalho de mestrado (Carrión, 2004). Mas, na prática, a
função de falha pode não ser linear e as variáveis aleatórias podem não ter uma
distribuição normal, como mostrado por Baecher & Christian em 2003, onde a
partir de dados de laboratório, verifica que o ângulo de atrito φ para alguns tipos
de solo tem uma distribuição Beta (Baecher & Christian, 2003). Sendo assim, são
necessários métodos alternativos que permitam o cálculo de confiabilidade para
uma função de falha qualquer dependente de variáveis aleatórias com distribuição
qualquer.
Para a análise de confiabilidade considerando variáveis aleatórias com
distribuição qualquer, algumas técnicas já foram desenvolvidas, tais como o
método FORM (First Order Reliability Method) e o método SORM (Second
Order Reliability Method) (Melchers, 2002). Estes métodos já foram amplamente
pesquisados, testados em diferentes trabalhos como (Kiureghian, 1994), (Wang e
Grandhi, 1994), (Val et all, 1995), (Imai e Frangopol, 1999), (Lee, 2000), (Pereira,
2007), (Lopes, 2007) e (Almeida, 2008).
O objetivo deste capítulo é descrever os conceitos fundamentais de
confiabilidade, apresentar um resumo dos diferentes métodos de cálculo,
apresentar um resumo do processo de cálculo pelo método FORM usado neste
trabalho e finalmente a Análise de Confiabilidade com a Análise Limite é
ilustrada com dois exemplos de aplicação. Um primeiro exemplo é aplicado para
um problema de talude 2D e o segundo exemplo para um problema de talude
confinado 3D. Para complementar o entendimento deste capítulo, no Apêndice A,
são apresentados os conceitos básicos da estatística.
O método FORM usado no presente trabalho, foi implementado, testado e
aplicado em várias pesquisas do Departamento de Engenharia Civil da PUC-Rio,
entre eles (Pereira, 2007), (Lopes, 2007) e (Almeida, 2008). No presente trabalho
usa-se este algoritmo simplesmente como usuário para a função de falha dada pela
Equação 4.4 e a gradiente da função de falha dada pelas equações (4.39-4.42).
95
4.1. Conceitos Fundamentais da Análise de Confiabilidade
Sendo o objetivo do presente trabalho a aplicabilidade prática da
confiabilidade em estruturas geotécnicas, nesta seção são apresentados somente
alguns conceitos de interesse para o presente trabalho.
4.1.1. Incertezas
A incerteza é a falta de conhecimento a priori de forma exata do resultado
de uma medição ou previsão da segurança no caso de uma estrutura. É dizer
aplica-se para prever ou quantificar a influência da não exatidão das medidas
físicas realizadas (como na medição das propriedades dos materiais) e na previsão
de eventos (como a segurança ou falha de uma estrutura). Este termo é utilizado
em um vasto número de campos, como a filosofia, estatística, economia, ciência,
engenharia, etc.
Na engenharia geotécnica, pelo fato que as propriedades de solo são obtidas
mediante ensaios de campo ou de laboratório, as incertezas nas medições podem
ser devido a diferentes fatores como erros humanos, má calibração dos
equipamentos de laboratório, alteração dos estados das amostras, etc. Estas
incertezas impossibilitam que uma estrutura apresente uma segurança absoluta ou
que a medida de algum nível de segurança seja calculada de forma absoluta ou
determinística.
Para levar em conta estas incertezas, as propriedades do material são
consideradas como variáveis aleatórias (ver seção A.1), e a análise de
confiabilidade pode ser realizada para conhecer o nível de segurança da estrutura,
dependente destas variáveis.
4.1.2. Função de Falha
A falha de uma estrutura é um estado que significa que a estrutura atingiu
condições indesejáveis, podendo ocasionar colapso total ou parcial (estado limite
último) ou então, interrupção do seu uso normal (estado limite de serviço). A
96
função matemática que representa este estado é conhecida como a função de falha.
Esta função, é também conhecida na literatura como a função de estado limite, de
performance ou de margem de segurança.
Para um conjunto de variáveis aleatórias iX com ni ,...2,1= expressa pelo
vetor X (Equação 4.1). Pode-se dizer que, a função de falha é uma função
matemática dependente das variáveis aleatórias (Equação 4.2). Esta função, define
a região segura e a região de falha (não segura) de uma estrutura em função das
variáveis aleatórias envolvidas com a estrutura (Figura 4.1)
TnXXX },...,,{ 21=X (4.1)
),...,,()( 21 nXXXFF =X (4.2)
Na Figura 4.1, este conceito é ilustrado graficamente para uma função de
falha )(XF dependente de duas variáveis aleatórias 1X e 2X . Esta função é
construída de tal forma que 0)( <XF representa a região de falha, 0)( >XF
representa a região segura, e 0)( =XF representa a superfície de falha, ou seja
uma superfície de separação entre os estados de falha e de segurança da estrutura.
1X
2X0)( <XF
0)( >XFRegião segura
0)( =XF
Região de falha
1X
2X0)( <XF
0)( >XFRegião segura
0)( =XF
Região de falha
Figura 4.1 – Função de falha.
97
Para a análise de confiabilidade de estruturas geotécnicas, o vetor das
variáveis aleatórias está em função dos parâmetros dos materiais e dos
carregamentos atuantes sobre a estrutura como expresso pela seguinte equação.
T
nmmmCCC },...,,,,...,,,,...,,,,...,,{ 21212121 ωωωγγγϕϕϕ=X (4.3)
onde m é o número de materiais da estrutura, n é o número de carregamentos
atuantes; iC (coesão), iϕ (ângulo de atrito) e iγ (peso específico) são as
propriedades do material para mi ,...,2,1= ; e jω são os carregamentos externos
para nj ,...,2,1= .
A função de falha usada para a Análise de Confiabilidade com a Análise
Limite é expressa pela seguinte equação:
cbwCaF jiii −−= ),(),()( ωγϕαX (4.4)
onde α é o fator de colapso determinado pela Análise Limite que depende dos
parâmetros de resistência iC e iϕ dos materiais, w representa o carregamento
aleatório que pode estar em função de cargas de gravidade iγ ou cargas externas
iω , a é uma constante que representa o carregamento inicial na Análise Limite,
b é um fator constante que multiplica aos carregamentos aleatórios, e c é uma
constante que representa aos carregamentos não aleatórios.
4.1.3. Função Densidade de Probabilidade Conjunta
Quando se trabalha com mais de uma variável aleatória, então, a função
densidade de probabilidade conjunta ),...,,( 21 nxxxp é usada para descrever as
variáveis aleatórias. Assim, a função densidade de probabilidade conjunta das
variáveis aleatórias nXXX ,...,, 21 , deve satisfazer as seguintes condições:
0),...,,( 21 ≥nxxxp (4.5)
98
1...),...,,(... 2121 =∫ ∫ ∫+∞
∞−
+∞
∞−
+∞
∞−nn dxdxdxxxxp (4.6)
nnn
b
a
b
a
b
ann bXabXadxdxdxxxxp
n
n
≤≤≤≤=∫ ∫ ∫ ,...,Pr[...),...,,(... 1112121
1
1
2
2
] (4.7)
onde ],...,,Pr[ 222111 nnn bXabXabXa ≤≤≤≤≤≤ é probabilidade que as
variáveis iX sejam maiores que ia e menores que ib para ni ,...,2,1= ; onde ia e
ib são constantes.
A função densidade de probabilidade conjunta de várias variáveis aleatórias,
fornece uma completa informação de probabilidade sobre a função densidade
probabilidade de quaisquer das variáveis aleatórias. Assim, a função densidade de
probabilidade (PDF) de cada uma das variáveis aleatórias é chamada de função
densidade de probabilidade marginal.
Duas ou mais variáveis aleatórias são estatisticamente independentes
quando a função densidade de probabilidade conjunta é igual ao produto das
funções densidade de probabilidade marginal de cada uma das variáveis aleatórias
(Equação 4.8).
)()...()(),...,,( 2121 nn xpxpxpxxxp = (4.8)
Para o caso de duas variáveis aleatórias a função densidade de probabilidade
conjunta e as funções densidade de probabilidade marginal são ilustradas na
Figura 4.2 (Melchers, 2002).
99
1x
2x
),( 21 xxp
0=F
0<F0>F
)( 2xp
)( 1xp
),( 21 xxp
1X
2X
1x
2x
),( 21 xxp
0=F
0<F0>F
)( 2xp
)( 1xp
),( 21 xxp
1X
2X
Figura 4.2 – Função densidade de probabilidade conjunta (Melchers, 2002).
4.1.4. Probabilidade de Falha
De uma forma geral pode-se dizer que a probabilidade de falha é uma
integral n-dimensional da função densidade probabilidade conjunta, como
expressa pela seguinte equação:
∫∫∫∫<
=0
2121 ...),...,,(...F
nnf dxdxdxxxxpP (4.9)
onde ),...,,( 21 nxxxp é a função densidade de probabilidade conjunta da variável
aleatória X e 0<F representa a região de falha.
100
Quando os parâmetros estatísticos ou a função densidade de probabilidade
da função de falha são conhecidos ou podem ser determinados, a probabilidade de
falha é como ilustrada pela Figura 4.3, e calculada pela Equação 4.10.
f
)( fpF
F
fPÁrea =
f
)( fpF
F
fPÁrea = fPÁrea =
Figura 4.3 – Probabilidade de falha.
∫∞−
=0
)( dffpPf (4.10)
onde fP é a probabilidade de falha e )( fp é a função densidade de probabilidade
da função de falha F .
A probabilidade de falha também pode ser avaliada a partir do cálculo do
índice de confiabilidade β (Seção 4.1.6). Desde logo, existem métodos que
avaliam a probabilidade de falha em função de índice do confiabilidade, como
mostra a Seção (4.2).
4.1.5. Confiabilidade
A confiabilidade é uma metodologia científica aplicada para conhecer a
estimativa da segurança associada a uma estrutura, de forma a assegurar que esta
cumpra sua função sem falhar durante sua vida útil. Esta metodologia, leva em
101
consideração as incertezas associadas às variáveis envolvidas com a estrutura,
considerando estas como variáveis aleatórias. Portanto, a confiabilidade é a
probabilidade de uma estrutura desempenhar sua função, sem falhar por um
determinado período de tempo e dentro das condições de uso especificado. Deve-
se ter em conta que, para uma estrutura operar com um nível especificado (alvo)
de confiabilidade, é necessário que ele tenha sido adequadamente projetada,
construída e operada.
Matematicamente, a confiabilidade de uma estrutura é definida como o
complemento da probabilidade de falha, como expressa pela seguinte equação:
fPC −= 1 (4.11)
onde C é a confiabilidade e fP é a probabilidade de falha (seção 4.1.4).
Da Equação 4.11, pode-se observar que o cálculo de confiabilidade depende
diretamente do cálculo da probabilidade de falha (Seção 4.1.4). Pelo que, o
cálculo de probabilidade de falha torna-se fundamental no cálculo de
confiabilidade.
4.1.6. Índice de Confiabilidade
Pelas dificuldades apresentadas no cálculo de probabilidade de falha pela
integração n-dimensional (Equação 4.9), alguns métodos fazem o cálculo de
probabilidade de falha a partir da avaliação de índice de confiabilidade. O índice
de confiabilidade é definido como a relação entre o valor esperado e o desvio
padrão da função de falha (Equação 4.12). Esta relação, é definida apenas para
funções de falha lineares e para variáveis com distribuição normal e não considera
outro tipo de distribuição, pelo que métodos que calculam a probabilidade de falha
a partir de índice de confiabilidade considerando outros tipos de distribuições,
precisam transformar para uma distribuição normal equivalente (seção 4.1.8).
102
))(()(XFVar
FE ><=
Xβ (4.12)
onde, β é o índice de confiabilidade, >< )(XFE é o valor esperado da função
de falha e ))(( XFVar é a variância da função de falha.
A interpretação gráfica do índice de confiabilidade β é apresentada na
Figura 4.4, onde se observa que este índice é um fator que multiplicado ao desvio
padrão, mede a distancia entre a origem e o valor médio da função de falha.
Fβσ
)0( >Fseguro)0( <Ffalha
F f
)( fpF
o
AreaPf =
Fβσ
)0( >Fseguro)0( <Ffalha
F f
)( fpF
o
AreaPf = AreaPf =
Figura 4.4 – Função densidade de probabilidade.
Da Figura 4.4, pode-se observar também que existe uma relação entre a
probabilidade de falha fP e o índice de confiabilidade β . Esta relação é
apresentada na Tabela 4.1, onde, podem-se observar os valores de probabilidade
de falha para diferentes valores de índice de confiabilidade.
Em métodos baseados no índice de confiabilidade, a probabilidade de falha
pode ser estimada a partir da Tabela 4.9, ou então, pode ser calculada como
expressa pela seguinte equação.
)( β−Φ=fP (4.13)
onde, Φ é função cumulativa normal padrão (Equação 4.16).
103
β Pf
5.200 1.00E-074.750 1.00E-064.270 1.00E-053.720 1.00E-043.090 1.00E-032.320 1.00E-021.280 1.00E-010.841 2.00E-010.524 3.00E-010.253 4.00E-010.000 5.00E-01-0.254 6.00E-01-0.525 7.00E-01-0.842 8.00E-01-1.286 9.00E-01-1.645 9.50E-01-2.327 9.90E-01
Tabela 4.1 – Índice de confiabilidade e probabilidade de falha.
4.1.7. Espaço Reduzido
Espaço reduzido é um espaço onde as variáveis aleatórias com média
diferente de zero )0( ≠X e desvio padrão diferente de um ( )1≠Xσ , são
transformados para um espaço equivalente onde a média é igual a zero )0( =Y e
o desvio padrão é igual a um )1( =Yσ . Esta transformação é feita pela seguinte
equação.
X
XXYσ−
= (4.14)
onde Y é a variável aleatória no espaço reduzido e X é a variável aleatória no
espaço original com média X e desvio padrão Xσ .
A Figuras 4.5, mostra a transformação de uma variável aleatória normal de
espaço original ( 0≠X e 1≠Xσ ) para espaço reduzido ( 0=Y e 1=Yσ ). A
distribuição normal no espaço reduzido é conhecida como distribuição normal
104
padrão, onde, sua função densidade de probabilidade e sua função distribuição de
probabilidade são dadas pelas equações que seguem:
2
2
21)(
y
ey−
=π
φ (4.15)
dzeyy z
∫∞−
−=Φ 2
2
21)(π
(4.16)
A Figura 4.6, ilustra a transformação de duas variáveis aleatórias com
distribuição quaisquer ( 0≠X e 1≠Xσ ) para o espaço reduzido com
distribuição normal padrão ( 0=Y e 1=Yσ ).
)(xp
xX
baX
x ==
σ
0Espaço original
)(yp
y0=Y
10==
y
Yσ
Espaço reduzido (normal padrão)
)(xTy =
)(xp
xX
baX
x ==
σ
0Espaço original
)(yp
y0=Y
10==
y
Yσ
Espaço reduzido (normal padrão)
)(xTy = )(xTy =
Figura 4.5 – Espaço original e espaço reduzido para uma variável.
105
1x
2x
Espaço original
1y
2y
Espaço reduzido (normal padrão)
)(xy T=
1x
2x
Espaço original
1y
2y
Espaço reduzido (normal padrão)
)(xy T= )(xy T=
Figura 4.6 – Espaço original e reduzido para duas variáveis.
4.1.8. Distribuição Normal Equivalente
Uma variável aleatória X com distribuição diferente da normal, pode ser
transformada para uma distribuição normal equivalente no ponto de pesquisa *x
como é ilustrada pelas Figuras 4.7 e 4.8.
*x x
)(xp
normalPDF
)( *xpXdePDF
*x x
)(xp
normalPDF
)( *xpXdePDF
Figura 4.7 – Funções densidade de probabilidade PDF.
106
*x x
)(xP
0.1normalCDF
XdeCDF)( *xP
*x x
)(xP
0.1normalCDF
XdeCDF)( *xP
Figura 4.8 – Funções distribuição de probabilidade CDF.
Obter uma distribuição normal equivalente implica obter a média e o desvio
padrão desta distribuição. Estas grandezas, são calculadas igualando-se as funções
densidade de probabilidade PDF e de distribuição CDF da distribuição real de X
e da distribuição normal equivalente no ponto de pesquisa *x , como expressa
pelas seguintes equações.
)(1 **
xpx
Nx
Nx
Nx
=⎟⎟⎠
⎞⎜⎜⎝
⎛ −σµ
φσ
(4.17)
)( **
xPx
Nx
Nx =⎟⎟⎠
⎞⎜⎜⎝
⎛ −Φ
σµ
(4.18)
onde, (.)φ e (.)Φ correspondem às funções PDF e CDF da distribuição normal
padrão respectivamente, (.)p e (.)P correspondem, respectivamente, às funções
PDF e CDF da distribuição não normal da variável X e Nxσ e N
xµ são,
respectivamente, a média e desvio padrão da normal equivalente no ponto *x .
As variáveis Nxσ e N
xµ , são calculadas através da resolução das Equações
4.17 e 4.18 e são expressas como:
[ ]{ }
)()(
*
*1
xpxPN
x
−Φ=φσ (4.19)
107
[ ])( *1* xPx Nx
Nx
−Φ−= σµ (4.20)
onde, (.)1−Φ corresponde a inversa da distribuição cumulativa normal padrão.
4.2. Métodos de Cálculo
Segundo a Equação 4.13, a confiabilidade é avaliada a partir do cálculo da
probabilidade de falha. Para o cálculo da probabilidade de falha existem diferentes
métodos desenvolvidos que a seguir são descritos brevemente.
Para os casos onde é conhecida a função densidade de probabilidade
conjunta(seção 4.1.3) ou a função densidade de probabilidade marginal de cada
variável aleatória (variáveis estatisticamente independentes), a probabilidade de
falha pode ser determinada pela avaliação da integral n -dimensional (Equação
4.9). Este método é conhecido como método de integração direta e pode ser
aplicado para problemas com poucas variáveis aleatórias e com a mesma
distribuição.
Sagrilo em 2003, indica que a avaliação da Equação 4.9, não é muito
simples, uma vez que ela envolve a avaliação de uma integral n -dimensional num
domínio complexo, onde n é o número de variáveis aleatórias. Mesmo com
desenvolvimento de técnicas modernas de integração numérica e com
computadores cada vez mais eficientes, na prática, a avaliação da Equação 4.9 por
integração, tem se restringido a variáveis com 5 ou 6 variáveis aleatórias no
máximo(Sagrilo, 2003).
Método de Montecarlo, é um método de simulação que também pode ser
usado para o cálculo da probabilidade de falha, como mostrado na dissertação de
mestrado (Carrion, 2004). A característica fundamental deste método é a geração
de dados amostrais para as variáveis aleatórias num experimento computacional.
Conjuntos de valores aleatórios são gerados para as variáveis aleatórias e a
resposta da estrutura é então avaliada para cada um dos dados aleatórios gerados.
O grande inconveniente deste método é a grande quantidade de avaliação da
resposta da estrutura, o que torna este método computacionalmente caro, quando a
108
avaliação da resposta da estrutura, tem custo computacional alto, como é o caso da
Análise Limite Numérica pelo MEF.
A probabilidade de falha também pode ser avaliada pelo método Estatístico
Linear como mostrado na dissertação de mestrado (Carrión, 2004). Este método
baseia-se simplesmente na média, desvio padrão e coeficiente de correlação das
variáveis aleatórias, onde os parâmetros estatísticos da função de falha (Equações
4.22 e 4.23) são aproximados a partir da aproximação da função de falha pela
serie de Taylor de primeira ordem (Equação 4.21).
...)()()(1
+−∂∂
+= ∑=
n
iii
i
XXXFXFXF (4.21)
)()( XFXFEF >≈<= (4.22)
TXF XXFVar )(.).()(2 FSF ∇∇≈=σ (4.23)
onde, F e Fσ são a média e o desvio padrão da função de falha )(XF , )(XF∇
é a gradiente da função de falha avaliado na média e XS é a matriz covariância
das variáveis aleatórias (Equação A.14).
Calculado os parâmetros estatísticos da função de falha, neste método, a
probabilidade de falha é determinada pela Equação 4.12.
Método MVFOSM (Mean Value First Order Second Moment). Este
método também é baseado na aproximação de primeira ordem por série de Taylor
da função de falha na media das variáveis aleatórias (Equação 4.21) (Cornell,
1969). É dizer, os parâmetros estatísticos são calculados do mesmo modo que para
o método Estatístico Linear, mas a probabilidade de falha é calculada a partir da
avaliação do índice de confiabilidade β , onde, este índice, é determinado pela
seguinte equação:
F
Fσ
β = (4.24)
109
Logo, a probabilidade de falha pode ser aproximada a partir da Tabela 4.1
ou calculada como )( β−Φ=fP , onde Φ é a função cumulativa normal padrão.
Ditlevesen (1973) e Veziano (1974) indicam que, o método MVFOSM dá
uma idéia tosca do nível de confiabilidade, isso pode não ser aceito como uma
aproximação válida devido à falta de formulação invariante e sua falha para
incorporar a informação da distribuição das variáveis aleatórias (Ditlevesen, 1973
e Veziano, 1974). Lee em 2000, indica que o índice de confiabilidade definido
pela Equação (4.24) falha ao não ser constante para uma função de performance
formulada de forma diferente mas mecanicamente equivalente(Lee, 2000).
O problema descrito no parágrafo anterior, foi superado pelo método
AFOSM (advanced first-order second moment), transformando as variáveis
aleatórias de espaço original para espaço reduzido (Hasofer and Lind, 1974).
O método AFOSM proposto por Hasofer e Lind é aplicável para variáveis
aleatórias com distribuição normal. Neste método, o índice de confiabilidade para
o cálculo de probabilidade de falha é definido como a mínima distância da origem
à superfície de falha no espaço reduzido, como ilustrado na Figura 4.9, para uma
função de falha não linear. Como primeiro passo, a variável em coordenadas reais
deveria de ser transformada para coordenadas reduzidas (seção 4.1.7). Logo a
função de falha no sistema de coordenadas originais é transformada também para
sistema de coordenadas reduzidas como:
),...,,()( 21 nYYYFYF = (4.25)
onde, iY são as variáveis aleatórias e )(YF a função de falha, ambos no espaço
reduzido.
110
1Y
2Y
β
0
*y
0)( <YFRegião não segura
0)( >YF
Região segura
0)( =YFFunção de falha
1Y
2Y
β
0
*y
0)( <YFRegião não segura
0)( >YF
Região segura
0)( =YFFunção de falha
0)( =YFFunção de falha
Figura 4.9 – Variáveis em coordenadas reduzidas: função de falha não linear.
Na Figura 4.9, o ponto *y que representa a distância mínima da origem à
função de falha é nomeado como ponto de projeto ou ponto mais provável de
falha MPP(most probable point). Para a função de falha não linear, o cálculo da
mínima distância torna-se um problema de otimização, e é formulado como segue:
0)(...
=
=
YFtsMinimize t yyβ (4.26)
O problema de otimização (Equação 4.26) foi resolvido por Shinozuka
(1983), usando o método de multiplicadores de Lagrange, de onde, a mínima
distância no espaço reduzido que representa o índice de confiabilidade é avaliada
pela seguinte equação.
∑
∑
=
=
⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
−=n
i i
i
n
ii
YF
YFy
1
2
1
*
β (4.27)
onde as derivadas ⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
iYF são avaliadas no ponto de projeto, ),...,,( **
2*1 nyyy .
111
O índice de confiabilidade de Hasofer-Lind para o cálculo de probabilidade
de falha é avaliado de forma exata, somente, quando todas as variáveis aleatórias
têm uma distribuição normal.
A grande desvantagem dos métodos anteriormente descritos é que eles não
consideram a informação referente ao tipo de distribuição das variáveis aleatórias
e que a avaliação da probabilidade de falha fP pode ser exata só para variáveis
aleatórias com distribuição normal.
Métodos baseados no ponto de projeto foram desenvolvidos como
alternativa aos métodos anteriormente descritos, entre estes métodos tem-se o
método FORM (first-order reliability method) e o método SORM(second-order
reliability method). A idéia do método SORM é a mesma que do método FORM,
a diferença entre estes dois métodos está na aproximação feita para a superfície de
falha no espaço reduzido. No método FORM a função de falha é aproximada no
ponto de projeto *y como superfície de falha linear pela série de Taylor de até
primeira ordem (Equação 4.21). No caso do método SORM, a função de falha é
aproximada como superfície quadrática pela série de Taylor de até segunda
ordem. Vários trabalhos ou pesquisas já foram realizados usando estes métodos,
entre eles Lee (2000), Lopez (2007), Pereira (2007) e Almeida (2008).
Lee(2000) indica que apesar do SORM poder prover melhores resultados
que FORM no mesmo problema, a necessidade de usar aproximação pelo método
SORM parece não ser significante. Chega a esta conclusão principalmente da
comparação do esforço computacional requerido e da melhora nos resultados
obtidos. O método SORM requer das derivadas parciais de segunda ordem da
função de falha que pode ser difícil de avaliar para estruturas complexas e requer
tempos de cálculo não necessários. Por outra parte, não é evidentemente óbvio
que os resultados podem ser mais próximos que os obtidos pelo método FORM.
Portanto, FORM pode ser um método mais prático que SORM para problemas de
engenharia, considerando a simplicidade conceitual e eficiência(Lee, 2000).
Segundo Almeida (2008) indica que o método FORM propicia, na maioria
dos problemas, uma precisão satisfatória com um tempo da análise computacional
reduzido quando comparado a outros métodos, o que justifica sua larga utilização
nas diversas análises de confiabilidade (Almeida, 2008).
112
Pelo indicado nos parágrafos anteriores, neste trabalho tenta-se aplicar o
método FORM, para avaliar a confiabilidade de Estruturas Geotécnicas.
4.2.1. Método FORM (First-Order Reliability Method)
O método de confiabilidade de primeira ordem FORM, é também conhecido
na literatura como a aproximação de Rackwitz-Fissler (Rackwitz e Fiessler,
1978). A grande vantagem deste método é que considera o tipo de distribuição das
variáveis aleatórias.
Rackwitz e Fiessler (1978) sugerem que, quando o problema envolve
variáveis aleatórias com distribuição não normal, eles podem ser resolvidos pela
transformação das variáveis não normais em variáveis normais equivalentes
(Seção 4.1.8).
Este método, é usado extensamente em muitos trabalhos de pesquisa. No
entanto, o cálculo só pode dar resultados aproximados para o caso de função de
falha com alta não linearidade. Não obstante, esta aproximação tem sido
considerada como um método efetivo na Análise de Confiabilidade por causa de
sua simplicidade e versatilidade.
4.2.1.1. Transformação de Variáveis
Para problemas de confiabilidade, onde, a função de falha depende de
variáveis aleatórias X com distribuição não normal e são estatisticamente
dependentes ou correlacionadas, existem métodos para transformar estas variáveis
em variáveis aleatórias Y com distribuições normais estatisticamente
independentes ou não correlacionadas.
As variáveis aleatórias estatisticamente independentes (não correlacionados)
com distribuição diferente da normal, podem ser transformadas para distribuição
normal equivalente no ponto de pesquisa como mostrado na Seção 4.1.8, e quando
as variáveis são estatisticamente dependentes (correlacionados) também é possível
usar a mesma transformação para obter as normais equivalentes, mas neste caso,
os coeficientes de correlação entre as variáveis originais devem também ser
113
corrigidos para coeficientes de correlação equivalentes. A correção dos
coeficientes de correlação é apresentada na Tabela 4.2. Uma referência completa
sobre a correção de coeficiente de correlação pode ser encontrada em Kiureghian
e Liu (Kiureghian e Liu, 1986) ou Lopes (Lopes, 2007).
Tipo III (min.)(Weibull)
Tipo II(máximos)
Tipo I(mínimo)
Tipo I (Max)(Gumbel)
Uniforme
Rayleigh
Lognormal
Normal
Normal
Var. ( j )Var. ( i )Coeficiente de Correlação Equivalente
Distribuição
Tipo III (min.)(Weibull)
Tipo II(máximos)
Tipo I(mínimo)
Tipo I (Max)(Gumbel)
Uniforme
Rayleigh
Lognormal
Normal
Normal
Var. ( j )Var. ( i )Coeficiente de Correlação Equivalente
Distribuição
ijEij ρρ =
ij
i
iEij ρ
δδρ
)1ln( 2+=
ijEij ρρ 014.1=
ijEij ρρ 023.1=
ijEij ρρ 031.1=
ijEij ρρ 031.1=
ijjjEij ρδδρ )364.0238.0030.1( 2++=
ijjjEij ρδδρ )328.0195.0031.1( 2++=
Tabela 4.2 – Coeficiente de correlação equivalente.
Na Tabela 4.2, jδ é o coeficiente de variação (ver seção A.4.6) e pode-se
observar que o coeficiente de correlação equivalente Eijρ não depende do ponto
*x onde a transformação está sendo realizada.
Para transformar variáveis normais correlacionadas em variáveis normais
estatisticamente independentes (não correlacionadas), o método amplamente
usado na análise de confiabilidade é o método conhecido como a Transformação
de Nataf (Kiureghian and Liu, 1986).
A transformação de Nataf indica que para um vetor X de variáveis
aleatórias normais correlacionadas entre si, um outro vetor Y de variáveis
114
aleatórias normais padrão estatisticamente independentes pode ser obtido pela
seguinte equação de transformação:
)( mXJY −= (4.28)
11 −−=∂∂
= σLJXY (4.29)
onde, m é o vetor com as médias das normais equivalentes e J é o Jacobiano da
transformação, L é uma matriz triangular inferior obtida pela fatoração de
Choleski (Equações 4.30-4.32) da matriz coeficiente de correlação ρ (Equação
A.16) e σ é a matriz do desvio padrão (Equação A.10), das normais equivalentes.
TLLρ .= (4.30)
⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢
⎣
⎡
=
nnnn LLL
LLL
L
MOMM
21
2221
11
00000
L (4.31)
11
11
,...,10.1
1
1
2
1
1
11
11
>−=
<<⎟⎟⎠
⎞⎜⎜⎝
⎛−=
===
∑
∑−
=
−
=
iLL
ikLLL
L
niLL
i
jijii
k
jkjijik
kkik
ii
ρ
ρ
(4.32)
onde ikρ é o coeficiente de correlação entre as variáveis iX e kX .
115
4.2.1.2. Pesquisa de Ponto de Projeto
O objetivo no cálculo de probabilidade de falha pelo método FORM é
calcular o índice de confiabilidade β . Segundo a Figura 4.10, o índice de
confiabilidade é a distância entre a origem no espaço reduzido e o ponto mais
provável de falha ( *y ), sendo que, o ponto mais provável se encontra na menor
distância da origem do espaço reduzido à superfície de falha. Este problema é
formulado como um problema de minimização da distância entre o origem e a
superfície de falha com uma restrição. Este problema é formulado então como
segue:
0)(..||min=YFts
Y (4.33)
O problema expresso pela Equação 4.33, pode ser resolvido por qualquer
algoritmo de otimização, mas o algoritmo mais usado na análise de confiabilidade
é aquele desenvolvido por Hasofer and Lind (1974) e melhorado por Rackwitz e
Fiessler (1978). Este algoritmo, comumente conhecido como HL-RF, apresenta a
seguinte expressão de recorrência.
[ ] TkkkTkk
k YFYFVYFYF
Y )()()(|)(|
12
1 ∇−∇∇
=+ (4.34)
onde )( kYF∇ é a gradiente da função de falha no espaço reduzido e )( kYF é o
valor da função de falha, ambos avaliados no ponto kY .
As seguintes relações são necessárias no processo de cálculo pelo método
HL-RF:
11 −−= σLJ (4.35)
)()( XY FF = (4.36)
)( mXJY −= (4.37)
116
)()()( 1 XJY FF T ∇=∇ − (4.38)
onde )(XF∇ é a gradiente da função de falha no espaço original avaliado no
ponto X .
Para a Análise de Confiabilidade com Análise Limite, onde as variáveis
aleatórias são expressas pela Equação 4.3 e a função de falha pela Equação 4.4, a
gradiente da função de falha )(XF∇ é calculada por diferenças finitas como
segue:
[ ]i
iiiii
i CCCCa
CF
ηφαφηα ),(),( −+
=∂∂ (4.39)
[ ]i
iiiii
i
CCaFηφ
φαηφφαφ
),(),( −+=
∂∂ (4.40)
bF
i
−=∂∂γ
(4.41)
bF
j
−=∂∂ϖ
(4.42)
onde, η é o fator de perturbação.
Segundo a literatura o fator de perturbação pode estar entre 212 1010 −− ≤≤η . No presente trabalho usou-se 410−=η .
1x
Espaço original -
2x0)( <XF
0)( =XF
0)( >XF
x
1y
Espaço reduzido -
2y
β
0)( <YF0)( >YF
0)( =YF
FORM
y
1x
Espaço original -
2x0)( <XF
0)( =XF
0)( >XF
x
1x
Espaço original -
2x0)( <XF
0)( =XF
0)( >XF
x
1y
Espaço reduzido -
2y
β
0)( <YF0)( >YF
0)( =YF
FORM
y
1y
Espaço reduzido -
2y
β
0)( <YF0)( >YF
0)( =YF
FORM
y Figura 4.10 – Espaço original e reduzido para duas variáveis.
117
4.2.1.3. Processo de Cálculo
O processo de cálculo da Análise de Confiabilidade pelo método FORM é
resumido pelo seguinte fluxograma.
Inicio
ãodistribuiçX iijXi i,,,, δρσ
Dados
ijEij ψρρ =
0=k iki XX =
[ ]{ })(
)(*
*1
i
iNX xp
xPi
−Φ=φσ
[ ])( *1*i
NXi
NX xPX
i
−Φ−= σµ
m σ ρLLL =← T
)()( XFYF =
)()( 1 XFJYF ∇=∇ −
)( kkk mXJY −=
[ ] TkkkTk
k
k YFYFYYFYF
)()(.)()(
12
1 ∇−∇∇
=+Y
|||||||| 1
k
kk
YYYerro −
=+
tolerro >
|| 1+= kYβ
)(1 β−Φ= −fP
fPC −=1
Fim
1+= kk
1−= LJ 1−σ
sim
não
Calcular
Inicializar
Calcular
Montar
Calcular
Calcular
Incremento
Calcular
Inicio
ãodistribuiçX iijXi i,,,, δρσ
Dados
ijEij ψρρ =
0=k iki XX =
[ ]{ })(
)(*
*1
i
iNX xp
xPi
−Φ=φσ
[ ])( *1*i
NXi
NX xPX
i
−Φ−= σµ
m σ ρLLL =← T
)()( XFYF =
)()( 1 XFJYF ∇=∇ −
)( kkk mXJY −=
[ ] TkkkTk
k
k YFYFYYFYF
)()(.)()(
12
1 ∇−∇∇
=+Y
|||||||| 1
k
kk
YYYerro −
=+
tolerro >
|| 1+= kYβ
)(1 β−Φ= −fP
fPC −=1
Fim
1+= kk
1−= LJ 1−σ1−= LJ 1−σ
sim
não
Calcular
Inicializar
Calcular
Montar
Calcular
Calcular
Incremento
Calcular
Figura 4.11 – Fluxograma da Análise de Confiabilidade pelo método FORM.
118
4.3. Exemplos de Aplicação
Com a finalidade de ilustrar a aplicabilidade da Análise de Confiabilidade
com a Análise Limite, dois exemplos de aplicação são apresentados a seguir. No
primeiro e segundo exemplo são apresentados os resultados do cálculo de
confiabilidade.
Para o cálculo da função de falha (Equação 4.4) o programa GEOLIMA é
chamado desde o algoritmo do método FORM e para o cálculo da gradiente da
função de falha por diferenças finitas (Equações 4.39 e 4.40) o programa
GEOLIMA é chamado também mais duas vezes. É dizer, para estruturas formadas
somente por um tipo de material, para cada iteração do método FORM o
programa GEOLIMA é chamado 3 vezes e para estruturas formados por dois
materiais diferentes o programa GEOLIMA seria chamado cinco vezes.
4.3.1. Talude 2D
Este exemplo tenta mostrar a aplicabilidade da Análise de Confiabilidade
com a Análise Limite em uma aplicação 2D. A malha do problema a ser analisada
é apresentada na Figura 4.12. Os dados do problema considerados como variáveis
aleatórias são: Coesão média de 151 =X kN/m2, ângulo de atrito médio
5.222 =X ° e peso especifico médio 163 =X kN/m3; considerou-se coeficiente
de variação de 5% para as três variáveis e o coeficiente de correlação entre a
coesão e ângulo de atrito de 5%. O tipo de distribuição considerado para as
variáveis são Normal para coesão, Lognormal para ângulo de atrito e Gumbel
Tipo I para o peso específico.
Como se pode observar nos resultados da Análise de Confiabilidade, a
convergência pelo método FORM alcançada na iteração 208. O índice de
confiabilidade calculado pelo método é 1.575045=β , a probabilidade de falha
calculada é 5.76%=fP e a confiabilidade é %26.94=C .
O ponto mais provável de falha MPP determinado pelo método FORM é
15.0421 =X , 24.3122 =X e 15.8843 =X . A visualização gráfica das zonas de
119
plastificação, assim como a superfície de falha, para o ponto mais provável de
falha são apresentadas nas Figuras 4.13 e 4.14.
Esta análise foi feita apenas com uma malha de 25 elementos, que gera um
problema de otimização de pequena escala, com 75 restrições no total. O tempo
requerido pelo programa GEOLIMA para resolver o problema em cada chamada
foi de 1 seg, fazendo um total de 3 seg para cada iteração do método FORM. O
tempo total requerido ate alcançar a convergência foi de 624 seg.
Figura 4.12 – Malha de elementos finitos 2D (25 elementos com 36 nós).
RELIABILITY ANALISYS REPORT
Method : FORM
Project : These
Structure: Slop 2D
User : mcp
DATA
| Var. | Mean | Desv.Std. |Distribution | Xo |
| X1 | 15.0 | 0.75 | Normal | 15.0 |
| X2 | 22.5 | 1.12 | LogNormal | 22.5 |
| X3 | 16.0 | 0.80 | Gumbel Tipo I | 16.0 |
120
Correlation Matrix
1.00 0.05 0.00
1.00 0.00
1.00
Iterations = 208
Most Probable Point
X1 = 15.042
X2 = 24.312
X3 = 15.884
Failure Function: F(X) = 4.203146e-02
Reliability Index: beta = 1.575045
Failure Probability: Pf = PHI(-beta ) = 5.76%
Reliability: C = 94.24%
Figura 4.13 – Zonas de plastificação no MPP(Most Probable Point).
Figura 4.14 – Superfície de falha no MPP.
121
Para se ter uma idéia da importância dos valores iniciais das variáveis, uma
nova análise é feita, mas desta vez os valores iniciais das variáveis já não são as
médias como na análise anterior, os valores iniciais correspondem a os valores da
iteração 190 da primeira análise ( 041285.151 =oX , 291812.242 =oX ,
883648.153 =oX ). Para este caso, o número de iterações necessárias para a
convergência foi de 18 e os resultados foram as mesmas da primeira análise. Isso
indica a importância do ponto de início para acelerar a convergência da análise, o
qual indica que a estimativa de ponto de início diferente das médias, deve ser
pesquisada em futuros trabalhos.
RELIABILITY ANALISYS REPORT
Method : FORM
Project : These
Structure: Slop 2D
User : mcp
DATA
| Var. | Mean | Desv.Std. |Distribution | Xo |
| X1 | 15.0 | 0.75 | Normal | 15.041285 |
| X2 | 22.5 | 1.12 | LogNormal | 24.291812 |
| X3 | 16.0 | 0.80 | Gumbel Tipo I | 15.883648 |
Correlation Matrix
1.00 0.05 0.00
1.00 0.00
1.00
Iterations = 18
Most Probable Point
X1 = 15.042
X2 = 24.312
X3 = 15.884
Failure Function: F(X) = 4.203146e-02
Reliability Index: beta = 1.575045
Failure Probability: Pf = PHI(-beta ) = 5.76%
Reliability: C = 94.24%
122
4.3.2. Talude Confinado 3D
Este segundo exemplo mostra a aplicabilidade da Análise de Confiabilidade
com a Análise Limite em um problema 3D. A malha do problema de talude
confinado a ser analisada é apresentada na Figura 4.15. Os dados do problema
considerados como variáveis aleatórias são: Coesão média de 251 =X kN/m2,
ângulo de atrito médio 52 =X ° e peso especifico médio 193 =X kN/m3;
considerou-se coeficiente de variação de 5% para as três variáveis e coeficiente de
correlação para coesão e ângulo de atrito de 5%. O tipo de distribuição
considerado para as variáveis são Normal para coesão, Lognormal para ângulo de
atrito e Gumbel Tipo I para o peso específico.
O método FORM alcança a convergência na iteração 241. O índice de
confiabilidade calculado pelo método é 1.7801=β , a probabilidade de falha
calculada é .75%3=fP e a confiabilidade é %25.96=C .
O ponto mais provável de falha MPP determinado pelo método FORM é
25.1101 =X , 5.4582 =X e 18.8443 =X . As zonas de plastificação e a
superfície de falha para o ponto mais provável de falha são apresentadas nas
Figuras 4.16 e 4.17.
Esta análise foi feita apenas com uma malha de 27 elementos, que gera um
problema um problema de otimização de pequena escala, com 81 restrições em
total. O tempo requerido pelo programa GEOLIMA para resolver o problema em
cada chamada foi de 3 seg, fazendo um total de 9 seg para cada iteração do
método FORM. O tempo total requerido ate alcançar a convergência foi de 2169
seg.
123
Figura 4.15 – Malha de elementos finitos 3D (27 elementos com 64 nós).
Figura 4.16 – Zonas de plastificação no MPP.
124
Figura 4.17 – Superfície de falha no MPP.
Dos exemplos apresentados no presente Capítulo, pode-se concluir que a
convergência pelo método FORM é alcançada em número de iterações maior que
200, o qual implica um custo computacional muito grande, porque para cada
iteração é necessário fazer três Análises Limite. Este fato inviabiliza o uso da
Análise de Confiabilidade pelo método FORM com a Análise Limite para
problemas com malhas muito refinadas.
O método requer de mais esforço computacional, que o requerido pelo
método de Monte Carlo. Pelo método de Monte Carlo é necessário avaliar entre
200 a 400 vezes a resposta da estrutura pela Análise Limite. Pelo método FORM
seria necessário avaliar mais de 600 vezes a resposta da estrutura, isso porque para
cada iteração é necessário avaliar 3 vezes a resposta da estrutura, uma para o
cálculo da função de falha (Equação 4.4) e duas para o cálculo da gradiente da
função de falha(Equações 4.39 e 4.40).
Em geral, para estruturas formadas por mais de um tipo de material o
número de avaliações necessárias da resposta da estrutura seria duas vezes o
número de materiais mais um, para cada iteração do método FORM; para
estruturas formados por cinco materiais, como é o caso da aplicação 4 do seguinte
Capítulo, seriam necessárias 11 avaliações da resposta da estrutura pela Análise
Limite, para cada iteração do método FORM.
6 CONCLUSÕES E SUGESTÕES
6.1. Conclusões
Neste trabalho foi realizado o desenvolvimento e a implementação de
algoritmos em conjunto com o programa GEOLIMA (Carrión, 2004) para realizar
a Análise Limite de problemas de gerande porte. Otimizadores comerciais como
LINGO, MINOS foram também utilizadas na pesquisa para efeitos de comparação
com os resultados obtidos com o otimizador desenvolvido.
Da pesquisa e dos testes realizados com os otimizadores comerciais,
conclui-se que o otimizador LINGO pode ser usado somente para resolver
problemas da Análise Limite de pequena escala e o otimizador MINOS pode ser
usado para resolver problemas de até média escala.
O resultados de testes realizados com os programas LINGO e MINOS
indicam que o tempo de otimização aumenta exponencialmente com o incremento
de número de elementos da malha, tornando inviável seu uso para a solução de
problemas de grande escala.
Dos testes realizados conclui-se que o fator de colapso diminui com o
refinamento da malha (Figuras 3.29). Esta diminuição é não linear e bastante
pronunciada até malhas com 150 elementos (aproximadamente) e logo a
diminuição tem uma tendência quase linear suave. Pelo fato que a formulação
implementada no presente trabalho é a formulação mista fraca, onde o campo de
tensões não é estaticamente admissível (a não ser de forma aproximada) e a
condição de equilíbrio é garantida pelo principio dos trabalhos virtuais, é de se
esperar que o fator de colapso se aproxime pelo limite superior ao fator de colapso
real com o refinamento da malha.
Pelo fato de os otimizadores comerciais pesquisados não apresentarem bom
desempenho para resolver problemas da Análise Limite de grande escala, neste
trabalho foi implementado um otimizador especifico para problemas da Análise
Limite, baseado no método de pontos interiores. O método de pontos interiores
162
pode ser aplicado para resolver problemas de otimização em geral, mas o
otimizador implementado neste trabalho é específico para resolver problemas da
Análise Limite.
Com o objetivo de melhorar a eficiência do otimizador implementado,
foram implementadas e testadas diferentes alternativas ou técnicas descritas na
literatura pesquisada e algumas alternativas são propostas neste trabalho para a
melhora da eficiência do otimizador implementado.
Na literatura pesquisada, foram encontrados dois algoritmos de pontos
interiores, o primeiro baseado em cálculo de ângulo de deflexão, e que precisa
resolver duas vezes o sistema de equações lineares para cada iteração e o segundo,
baseado no cálculo de fatores de relaxação e contração. Este segundo algoritmo
precisa resolver uma só vez o sistema de equações lineares mas se torna
ineficiente por que a convergência é alcançada em maior número de iterações.
Neste trabalho foi proposto e implementado um terceiro algoritmo, onde o sistema
de equações lineares é resolvido uma só vez e a direção de busca viável é
calculada por álgebra vetorial. Testes realizados com os três algoritmos mostram
que o algoritmo proposto fornece valores de fator de colapso ligeiramente
menores, mas que tem um ganho em desempenho (tempo de otimização).
Um ponto crítico encontrado na implementação do otimizador foi a solução
de sistemas lineares em cada iteração. Na literatura pesquisada foram encontradas
duas formas de manipular e resolver o sistema de equações resultantes. A primeira
é uma manipulação matricial global e a segunda uma manipulação matricial por
elementos. Neste trabalho propõe-se a solução direta sem manipulação de sistema
de equações lineares. Dos testes realizados (Tabela 3.8 e Figuras 3.18 e 3.19)
conclui-se que a solução direta com tratamento da matriz esparsa proposta neste
trabalho é muito mais eficiente tanto em uso de memória quanto em desempenho
em termos de tempo.
Para a solução de sistemas lineares foram implementados diferentes
métodos diretos e indiretos. Dentre os métodos testados, o dos gradientes
conjugados (CG) com tratamento da matriz esparsa e uso de pré-condicionador
misto (DS+SS) foi o mais eficiente.
Para melhorar a eficiência na solução de sistemas lineares, foram
implementados e testados os seguintes pré-condicionadores encontrados na
literatura, escala diagonal (DS), escala simétrica (SS) e fatoração incompleta de
163
Cholesky (ICF). Neste trabalho propõe-se também o uso de pré-condicionadores
mistos DS+SS e DS+FIC. Dos testes realizados conclui-se que o pré-
condicionador misto proposto DS+SS teve melhor desempenho.
Da comparação dos testes realizados em problema 2D, onde o problema é
resolvido usando otimizadores Lingo, Minos e Geolima com otimizador próprio,
conclui-se que o otimizador implementado faz melhor uso da memória como
também tem melhor desempenho em termos de tempo e o mais importante e que o
diferencia dos outros otimizadores é que o desempenho de GEOLIMA cresce
quase linearmente com o refinamento da malha.
Com o objetivo de validar e comparar a eficiência do otimizador
implementado, um problema 3D de frente de escavação de túneis é analisado. Este
problema tem seu comportamento conhecido a partir de uma simulação com um
modelo físico 3D em escala reduzida feito no Laboratório da Mitsubishi em Japão
e também da análise com o otimizador Minos apresentado em Carrión (2004). O
fator de colapso obtido com o otimizador Mino foi 012.25=α e com o
otimizador implementado foi 6685.24=α . O tempo empregado por Minos foi de
16 horas 43 min. e com o otimizador implementado foi de 17 min 49 seg. O
mecanismo de colapso obtido foi semelhante aquele obtido com Minos e muito
similar aquele obtido no modelo físico.
Propriedades de solo são obtidas mediante ensaios de campo ou de
laboratório, incertezas de diferente origem estão presentes nas propriedades dos
materiais assim obtidas. Assim, torna-se relevante a incorporação de estudos de
confiabilidade na análise de estruturas geotécnicas.
Existem vários métodos para avaliar a confiabilidade, mas a grande maioria
é desenvolvida para funções de falha lineares e considerando simplesmente uma
distribuição normal para as variáveis aleatórias. Neste trabalho, tentou-se aplicar
o método FORM(First Order Reliability Method) cuja formulação considera uma
função de falha qualquer e distribuições quaisquer das variáveis aleatórias, mas
resultados obtidos indicam que o método alcança a convergência em um número
muito alto de iterações (maior que 200), o qual resulta em custo computacional
maior ainda que o método de Monte Carlo, inviabilizando o uso deste método par
a cálculo de confiabilidade com a Análise Limite.
164
Dos resultados da primeira aplicação (talude infinito homogêneo), conclui-
se que o fator de colapso, o número de iterações e o tempo empregado pelo
programa em resolver o problema variam com o tipo de material da estrutura.
Além disso, zonas de plastifição, a localização e forma da superfície de ruptura e
o mecanismo de colapso variam com tipo de material.
Dos resultados da segunda aplicação (talude infinito heterogêneo), conclui-
se que o fator de colapso, número de iterações e o tempo requerido para a análise
variam com o refinamento da malha. Sendo o refinamento da malha importante
para identificar as zonas de plastifição e a superfície de ruptura.
Dos resultados da terceira análise (talude com percolação), pode-se concluir
que a presença de percolação de água no talude influencia na redução de fator de
colapso em um 23%, assim como o programa precisou de maior número de
iterações para resolver o problema e por conseguinte maior tempo de análise
também. E quanto às zonas de plastificação e superfície de falha, pode-se observar
das figuras 5.12(a,b) e 5.13(a,b) que a presença da percolação faz que tanto as
zonas de plastificação como a superfície de falha sejam mais superficiais.
Dos resultados da Análise Limite da barragem de terra(quarta aplicação),
conclui-se que, no caso da análise 1, a estrutura entraria em colapso somente no
caso em que o carregamento ( gravitacional ) fosse 2.46 vezes maior. No caso da
análise 2, a estrutura entraria em colapso se a pressão da água for maior em 5.24
vezes, o que indica que a barragem projetada é bastante conservadora em relação
à segurança.
Da quinta aplicação (talude confinado 3D) se pode concluir que, com o
programa implementado pode-se resolver problemas com malhas suficientemente
refinadas que geram problemas de otimização de grande escala, o qual já
possibilita a aplicabilidade da Análise Limite na solução de problemas reais.
A sexta aplicação se refere à análise de um depósito de rejeito em 3D.
Sendo o fator de colapso obtido para esta estrutura próximo a 1, não se pode
afirmar que a estabilidade do depósito está garantida, porque, efeito de algum
outro fenômeno como movimento sísmico na direção WE poderia desestabilizar o
depósito de rejeito. Assim, é recomendável diminuir a inclinação de talude com a
finalidade de garantir sua estabilidade para fins de abandono definitivo.
165
6.2. Sugestões para Futuras Pesquisas
A implementação da Análise Limite Numérica mediante processamento
paralelo devem de ser considerados nas futuras pesquisas, para a análise de
problemas de maior porte.
Refinamento da malha implica em custos computacionais na Análise
Limite, nas futuras pesquisas devem-se implementar o refinamento automático da
malha (autoadaptativo) nas zonas de maior plastificação.
Para garantir a estabilidade de estruturas geotécnicas são usadas reforços
(como geosintéticos e ancoragens). Em trabalhos futuros deve-se considerar estes
reforços na formulação de Análise Limite, assim como o uso de elementos de
contato para a represetação de meios fraturados.
Para modelar problemas com geometria complexa, em trabalhos futuros
devem-se incorporar outros tipos de elementos finitos como triangulares em 2D e
tetraédricos em 3D, assim como elementos de interface.
Exemplos de aplicação da Análise de Confiabilidade associada à Análise
Limite mostram que o algoritmo de otimização HL-RF (Hasofer and Lind –
Rackwitz and Fiessler) usado pelo método FORM, tem uma convergência pobre.
Assim, devem-se considerar a implementação de algoritmos mais efetivos para a
Análise de Confiabilidade pelo método FORM.
Um dos fatores para a uma convergência pobre do algorimo HL-LR pode
ser o cálculo aproximado da gradiente da função de falha por diferencias finitas,
portanto, em pesquisas futuras devem-se considerar o uso de superfícies de
resposta, geradas das soluções do problema pela Análise Limite. Outras técnicas
como o uso de ponto de partida diferente das médias das variáveis devem de ser
pesquisado, como o objetivo de acelerar a convergência.
7 REFERÊNCIAS BIBLIOGRÁFICAS
ALMEIDA A. F. Projeto Ótimo Baseado em Confiabilidade de Pórticos Planos de Concreto armado. Tese de Doutorado, Departamento de Engenharia Civil – PUC-Rio, Rio de Janeiro – RJ – Brasil. 2008. ALKIRE B. Cholesky Factorization of Aumented Positive Definite Matrices. Electrical Engineering Department, University of California Los Angeles UCLA, 2002. AMES, W. F. “Numerical Methods” The Engineering Handbook. Ed. Richard C. Drof. Boca Raton CRC Press LLC, 2000. BAECHER G. B. and CHRISTIAN J. T. Reliability and Statistics in Geotechnical Engineering. John Wiley & Sons Ltd. England, 2003. BATHE K. J. Finete Element Procedures. Prentice-Hall, Inc. Printed in the United States of America, 1996. BENZI M. e TUMA M. A Robust Incomplete Factorization Preconditioner for Positive Definite Matrices. Numerical Linear Algebra with Applications. Copyright @ John Wiley & Sons, Ltd., 2001. BONGARTZ I., CONN A. R., GOULD N. I. M., SAUNDERS M. A. A Numerical Comparison Between the LANCELOT and MINOS Package for Large-scale Nonlinear Optimization. Report 97/13. Canadá: ProMIRA Software Inc., 1997. BORGES L. Formulação e Solução para Análise Limite com Superfície de Escoamento Não Linear. Tese de Doutorado. Departamento de Engenharia Mecânica – Pontifícia Universidade Católica de Rio de Janeiro PUC-Rio, Rio de Janeiro - Brasil, 1991. CAMPOS, J. L. E. Análise Numérica de Transporte de Contaminantes em Meios Porosos com Reações Químicas. Tese de Doutorado, Departamento de Engenharia Civil – PUC-Rio, Rio de Janeiro – RJ – Brasil. 1999. CARRION M. Análise Limite Tridimensional Determinística e Não Determinística. Dissertação de Mestrado. Departamento de Engenharia Civil – Pontifícia Universidade Católica de Rio de Janeiro PUC-Rio, Rio de Janeiro, 2004.
167
CARRION M, VARGAS E., VAZ E., FARFAN A. Análise Limite Numérica em 3D de Frente de Escavação de Túneis. INFOGEO 5to Simpósio Brasileiro de Aplicações de Informática em Geotecnia. Belo Horizonte-Brasil, 2005. CASTILLO E., CORNEJO A., PEDREGAL P., GARCIA R. E ALGUACIL N. Formulación y Resolución de Modelos de Programación Matemática en Ingeniería y Ciencia. España, 2002. CERVANTES G. & MEJIA C. Precondicionamiento de Métodos Iterativos. Revista de la Academia Colombiana de Ciências. Colombia, 2004. CHEN W. & LIU X. Limit Analysis in Soil Mechanics. Amsterdam-Oxford-New York: Elsevier Scientific Publishing Company, 1990. CIMNE. GID 9.0 User Manual. Barcelona: International Center for Numerical Methods in Engineering. 2008. CIMNE. GID 9.0 Reference Manual. Barcelona: International Center for Numerical Methods in Engineering, 2008. CISMID – Centro Peruano Japonés de Investigación Sísmica y Mitigación de Desastres. Estudio Definitivo para el Abandono de los Depósitos de Relave de Yauliyacu, Bellavista e Antuquito. Vol. I, II e III. Lima-Perú, 1998. CISMID - Centro Peruano Japonés de Investigación Sísmica y Mitigación de Desastres. Estudio de los Aspectos Geotécnicos y Análisis Detallado del Estado Actual de los Terraplenes del Pulmón de Regulación Horaria Catalinayocc. Lima-Perú, 2002. CONN A. R., GOULD N. I. M. AND TOINT P. L. LANCELOT A Fortran Package for Large-Scale Nonlinear Optimization (Release A). Springer Series in Computational Mathematics. United State of America, 1992. CORNELL C. A. Probabilidty-based Structural Code. Journal of ACI, 66(12), pp. 974-985, 1969. DRUCKER D. C.; GREENBERG H. J. E PRAGER W. Extended Limit Design Theorems for Continuous Media. Q. Appli. Math., vol 9. 1952. DITLEVSEN, O. Structural Reability and the Invariance Problem. Research Report, bf22, Solid Mechanics Div., University of Waterloo, Otario, Canada, 1973. FARFÁN A. Aplicações da Análise Limite a Problemas Geotécnicos Modelados como Meios Contínuos Convencionais e Meios de Cosserat. Tese de Doutorado. Rio de Janeiro: PUC-Rio, 2000. GOMEZ J e BURGUEST J. Curso de Cálculo Numérico. Departamento de Física Atômica Molecular y Nuclear. Universidad de Valencia, Facultat de Física. España, 2004.
168
GONZAGA L. Estudo Numérico de Problemas de Estabilidade em Materiais Geotécnicos Através da Análise Limite. Tese de Doutorado. Rio de Janeiro: PUC-Rio, 2000. GVOZDEV A. A. The Determination of the Collapse Load for Statically Indeterminate System Undergoing Plastic Deformations. Int. J. Mech. Sci., vol 1. (Tradução do original de 1938), 1960. HASOFER A. M. and LIND N. C. Exact and Invariant Second-moment code format. Journal of Engineering Mechanics, ASCE 100 (EM1). PP. 111-121, 1974. HERSKOVITS J. Comunicação Persoal. 2008. HESTENES M. R. and STIFEL E. Methods of Conjugate Gradient for Solving Linear Systems. Journal of Research of National Bureau of Standars. 49, 409-436, 1952. HOFFMAN J. D. Numerical Methods for Engineers and Scientists. Second Edition, Revised and Expanded. Department of Mechanical Engineering, Purdue University. New York, 1992. IMAI K. and FRANGOPOL D. Geometrically Nonlinear Finite Element Reliability Análysis of Structural Systems. PERGAMON - Computers & Structures. Elsevier Science Ltd., 1999. JENNINGS A. and MALIK G. M. The Solution of Sparce Linear Equations by the Conjugate Gradient Method. International Journal for Numerical Methods in Engineering, 12:141-158, 1978. KIUREGHIAN A. D. and LIU P. L. Structural Reliability Under Incomplete Probability Information. Journal of Engineering Mechanics (ASCE), Vol. 110, No 3, 1986. KIUREGHIAN A. Structural Reliability Methods for Seismic Safety Assessment. Department of Civil Engineering, University of California, Berkeley, USA, 1994. STUBEN K. and CLEES T.. SAMG USER’s Manual Release 22c. Fraunhofer Institute SCAI Scholoss Birlinghoven. Germany, 2005. STUBEN K. SAMG Data Structure and File Format Specification. Fraunhofer Institute SCAI Scholoss Birlinghoven. Germany, 2005. LANCELLOTTA R. Technical University of Turin. Department of Structural Engineering. Geotechnical Engineering. Bologna: Zanichelli Editores S.P.A., 1995.
169
LEE S. Y. Static and Dinamic Reliabilidty Analysis off Frama and Shear Walls Structural Systems. Doctor of Philosophy These. Department of Civil Engineering and Engineering Mechanics of the Universitty of Arizona, 2000. LINDO SYSTEMS INC. Lingo Manual for Release 5.3. Chicago: Lindo Systems Inc., 1997. LOPES, M. T. de A. Análise de Confiabilidade de Estruturas Aplicada ao Projeto de Reforço à Força Cortante de Vigas em Concreto Armado com Compósito de Fibras de Carbono. Tese de Doutorado. Pontifícia Universidade Católica de Rio de Janeiro PUC-Rio, 2007.
LYAMIN A. V. & SLOAN S. W. A comparison of linear and nonlinear programming formulation for lower bound limit analysis. Em Pietruszczak & Pande (eds), Numerical Models in Geomechanics: pp. 367-373. Rotterdam: Balkema. 1997. MEIJERINK J. A. and VAN DER VORST H. A. An Interecative Solution Method for Linear Systems of which the Coefficient Matrix is a Symetric m-matrix. Mathematics of Computation, 31: 148-155. 1977. MELCHERS, R. E. Structural Reliability Analysis and Prediction. 437p. John Wiley & Sons, New York. 2002.
MURTAGH B. A. AND SAUNDERS M. A. Minos 5.5 User’s Guide. California: Stanford University, 1998. PEREIRA A. Projeto Ótimo Baseado em Confiabilidade: Aplicação a Treliças espaciais. Tese de Doutorado, Departamento de Engenharia Civil – Pontifícia Universidade Católica de Rio de Janeiro PUC-Rio, Rio de Janeiro, 2007. PINHEIRO C. Aplicação de um Algoritmo Genético no Estudo das Perdas de Controle de Tensão em Sistemas Elétricos de Potencia. Dissertação de Mestrado. Universidade Federal de Minas Gerais, 1998. PINI G. and GAMBOLATI G. Is Simple Diagonal Scaling the Best Preconditioner for Conjugate Gradient on Supercomputers?. Advanced in Water Resources, 12(3):147-153, 1990. SAGRILO L. Apostila do Curso de Confiabilidade Estrutural. UFRJ, Rio de Janeiro, Brasil, 2003. SANDOVAL M. L. Métodos Iterativos Eficientes Para Problemas de Convección-Difusión Transitorios. Tesis Doctoral, Departamento de Matemática Aplicada III UPC. Barcelona, 2006. SMAILBEGOVIC F., GAYDADJIEV G. and VASSILIANDIS S. Sparse Matrix Storage Format. Coputer Engineering Laboratory, Electrical Engineering, Mathematics and Computer Sciense. Mekelweg 4, 2628CD Delf, 2006.
170
SHINOZUKA M. Basic Analysis of Structural Safety. Journal of Structural Engineering Division, ASCE, 109(3), pp. 721-740, 1983. RACKWITZ R. and FIESSLER B. Structural Reliability Under Combined Random Load Sequences. Computer and Structures, Vol. 9, 1978. TSAL R. J. Numerical Modeling of Multi-Fan Air Distribution System for Research Laboratories. Netsal & Associates. Fountain Valley CA 92708, 1998. TSAO N. On Equivalence of Gaussian Elimination and Gauss-Jordan Reduction in Solving Linear Equations. NASA Technical Memorandum 101466. Wayne State University and Institute for Computational Mechanics in Propulsion Lewis Research Center. Michigan and Ohio, 1989. VAL D., BLJUGER F. and YANKELEVSKY D. Optimization Problem Solution in Reliability Analysis of Reinforced Concrete Structures. National Building Research Institute, Faculty of Civil Engineering, Technion – Israel, 1995. VENEZIANO D. Contribution to Second Moment Reliability Theory. Research Report, Department of Civil Engineering, MIT, 1974. WANG L., and GRANDHI R. Safety Index Calculation Using Intervening Variables for Structural Reliability Analysis. Department of Mechanical and Materials Engineering, Wright State University, Dayton – U.S.A., 1994. ZOUAIN N., HERSKOVITS J., BORGES L. and FEIJO R. An Interactive Algorithm for Limit Analysis with Nonlinear Yield Functions. Int. J. Solids Structure Vol. 30. Great Britain, 1993.
A CONCEITOS ESTATÍSTICOS
A.1. Variáveis Determinísticas e Aleatórias
Se uma medição ou um experimento é feito repetidamente, mantendo todas
as condições possíveis, e os resultados obtidos são os mesmos, então a variável é
conhecida como determinística. Todavia, se os resultados obtidos nas medições
variam ou são diferentes, então as variáveis são conhecidas como aleatórias.
Na Engenharia Geotécnica, os parâmetros dos materiais geológicos são
obtidos mediante ensaios de laboratório ou de campo. Sendo os valores obtidos
diferentes para cada ensaio, estes parâmetros são de natureza aleatória.
A.2. Espaço Amostral, Evento e Valor Observado
O espaço amostral é o conjunto formado por todos os possíveis e diferentes
resultados observados de uma variável aleatória. Um evento é um subconjunto
formado por uma ou mais observações do espaço amostral. O valor numérico de
cada medição de uma variável aleatória é conhecido como valor observado e é
parte de seu espaço amostral.
A.3. Medidas de Tendência Central
Na Estatística existem diferentes medidas de tendência central, como a
mediana, a moda e diferentes tipos de médias. Esta pesquisa refere-se só à média
aritmética, por ser de interesse para os objetivos deste trabalho.
172
A.3.1. Momento estatístico de ordem m
Matematicamente o momento estatístico de ordem m é definido como a
somatória da potência m dos valores observados do espaço amostral dividida pelo
número de valores observados da variável aleatória.
∑=
=n
j
mi
jmi
s xn
X1
1 (A.1)
onde:
ij x j-ésimo valor observado da variável aleatória Xi
m Potência à qual estão elevados os valores da variável
n Número de valores observados da variável aleatória mi
s X Momento estatístico de ordem m da amostra
A.3.2. Média aritmética
A média aritmética é um momento estatístico de primeira ordem (m=1) e é
talvez a medida de tendência central mais usada.
∑=
=n
ji
ji
s xn
X1
1 (A.2)
onde:
ij x j-ésimo valor observado da variável aleatória Xi
n Número de valores observados da variável aleatória
is X Média da amostra da variável aleatória Xi
173
A.3.3. Características da média aritmética
• A unidade de medida da média aritmética é a mesma que a unidade dos
valores observados da variável aleatória.
• É um valor único.
• É fácil de compreender e aplicar.
• Utiliza todos os valores da variável em seu cálculo.
• É afetada pela mudança de algum valor da variável.
A.4. Medidas de Dispersão da Variável Aleatória
A.4.1. Desvio
O desvio é a diferença entre o valor observado de uma variável aleatória e a
média aritmética da amostra.
XxDesvio i−= (A.3)
A medida dos desvios das variáveis aleatórias tem duas características muito
importantes:
• A somatória de todos os desvios dos valores da variável aleatória é sempre
igual a zero.
01
=−∑=
]Xx[ is
n
ji
j (A.4)
• A somatória dos quadrados dos desvios dos valores da variável aleatória é
sempre um valor mínimo.
mínimo]Xx[ is
n
ji
j =−∑=
2
1
(A.5)
174
A.4.2. Momento estatístico central de ordem m
É definido como a somatória da mma potência dos desvios dos valores
observados da variável aleatória dividida pelo número de valores observados.
mi
sn
ji
jmi
s ]Xx[n
X −= ∑=1
1 (A.6)
onde: mi
s X Momento estatístico central de ordem m
ij x j-ésimo valor observado da variável aleatória Xi
is X Média da amostra da variável aleatória Xi
m Ordem a variável aleatória
n Número de valores observados da variável aleatória
A.4.3. Desvio absoluto médio
É um momento estatístico central de ordem um ( m=1 )
]Xx[n
X is
n
ji
ji
s −= ∑=1
1 (A.7)
onde:
is X Desvio absoluto médio
A.4.4. Variância
A variância da mostra é o momento estatístico central de ordem dois (m=2)
2
1
1 ]Xx[n
)x(Var is
n
ji
ji
s −= ∑=
(A.8)
175
Características da variância:
• A variância é sempre um valor positivo.
• Se todos os valores de uma variável forem iguais então a variância é zero.
A.4.5. Desvio padrão
Uma desvantagem da variância é sua unidade de medida, o quadrado da
unidade de medida dos valores da variável. Como essa unidade não explica nada
sobre as características dos valores da amostra, é definido o desvio padrão, que
mantém a unidade de medida da variável.
O desvio padrão da amostra é definido como a raiz quadrada da variância da
amostra.
2
1
1 ]Xx[n i
sn
ji
jxi
−= ∑=
σ (A.9)
onde:
ixσ Desvio padrão da variável aleatória Xi
Características do desvio padrão:
• O desvio padrão é sempre um número positivo.
• Se os valores de uma variável forem iguais, o desvio padrão é zero.
• O desvio padrão depende da soma dos quadrados dos desvios da variável;
portanto, quanto menor o desvio padrão, mais os valores da variável
aproximam-se da sua média.
• Se o desvio de um valor da variável for menor que o desvio padrão da
amostra, então esse valor estará mais próximo da média do que outro valor
com desvio maior.
• Quanto mais os valores da variável se afastarem de sua média, maior serão os
desvios e, conseqüentemente, maior será o desvio padrão e mais aberta será a
distribuição de freqüências da variável.
176
• Se duas variáveis tiverem a mesma média e desvios padrões diferentes, a
distribuição da variável com maior desvio padrão será mais aberta que a de
variável com menor desvio padrão.
• Para n variáveis aleatórias a matriz do desvio padrão é expressa pela seguinte
equação.
⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢
⎣
⎡
=
nx
x
x
σ
σσ
L
MOMM
00
000000
2
1
σ (A.10)
A.4.6. Coeficiente de variação
A comparação da dispersão de duas ou mais distribuições pelo simples
confronto de seus desvios padrões nem sempre é suficiente, pois as amostras
podem ter unidades diferentes ou, tendo a mesma unidade, seus valores de média
podem estar bastante afastados.
O coeficiente de variação é uma medida relativa de dispersão que permite a
comparação de distribuições, pois é definido como a relação do desvio padrão por
unidade de média. Na comparação de duas variáveis, a variável que tiver menor
coeficiente de variação tem menor dispersão, variabilidade ou incerteza.
i
s
xx X
i
i
σδ = (A.11)
onde:
ixδ Coeficiente de variação da variável aleatória Xi
177
A.5. Medidas de Correlação de Variáveis Aleatórias
As medidas estatísticas que medem a tendência e a relação linear entre as
variáveis aleatórias são a covariância e o coeficiente de correlação.
A.5.1. Covariância
A covariância resume em um único número a tendência e a relação linear
entre duas variáveis e é definida como a média dos produtos dos desvios das
variáveis.
]Xx[]Xx[n
)X,X(Cov sjn
j
sj22
11121
1−−= ∑
=
(A.12)
onde:
Cov(X1,X2) Covariância das variáveis aleatórias X1 e X2
1xj , 2xj Valores observados das variáveis aleatórias X1 e X2
1X , 2X Média dos valores observados das variáveis aleatórias X1 e X2
n Número de valores observados das variáveis aleatórias X1 e X2
Características da covariância:
• As duas variáveis devem ter o mesmo número de valores observados.
• O valor da covariância pode ser qualquer valor real, pois ela pode ser negativa,
positiva ou zero.
• A unidade de medida é o resultado do produto das unidades dos valores das
variáveis aleatórias.
• A covariância de uma variável e ela mesma é a própria variância da variável.
• A permutação das variáveis não altera o resultado da covariância, se os
mesmos pares de valores forem mantidos: )X,X(Cov)X,X(Cov 1221 = .
178
• Se as variáveis X1 e X2 forem estatisticamente independentes, então a
covariância destas variáveis será igual a zero.
• Se o resultado da covariância das variáveis X1 e X2 for igual a zero, não se
pode afirmar que as duas variáveis sejam estatisticamente independentes.
• Para a, b, c e d constantes e X, Y e Z variáveis aleatórias, então sempre se
cumpre que:
)Y,Z(Cov)Y,X(Cov)Y,ZX(Cov)Y,X(acCov)dcY,baX(Cov
)Y,X(aCov)Y,aX(Cov)Y,X(Cov)Y,X(Cov
)a,X(Cov
+=+=++
=−=−
= 0
(A.13)
• Para n variáveis aleatórias dependentes entre si, a matriz covariância XS de
estas variáveis é expressa pela seguinte equação:
⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢
⎣
⎡
=
)(),(),(
)(),()(
21
221
1
nmm
X
XVarXXCovXXCov
XVarXXCovSimétricoXVar
L
MOMMS (A.14)
A.5.2. Coeficiente de correlação
O coeficiente de correlação permite a comparação de duas variáveis, por ser
independente das unidades de medida.
O coeficiente de correlação de duas variáveis aleatórias é definido como a
covariância das amostras dividida pelo desvio padrão de cada amostra
kj
kjxx
kjXX
XXCovσσ
ρ),(
, = (A.15)
onde:
kj XX ,ρ Coeficiente de correlação das variáveis aleatórias Xj e Xk
179
O coeficiente de correlação entre duas variáveis aleatórias varia no intervalo
de 11 +<≤− ijρ . Uma interpretação gráfica deste conceito é apresentada na
Figura A.1, e o grado de dependência entre variáveis aleatórias é classificada
segundo o valor de coeficiente de correlação como mostrada na Tabela A.1.
Por questão prática, a notação formal de coeficiente de correlação entre duas
variáveis ji XX ,ρ é escrita em alguns casos simplesmente como ijρ .
ix ixix
ix
jx
jx jx jx
0=ijρ 0=ijρ 0=ijρ
ix
jx
1+=ijρix
jx
1−=ijρ 10 << ijρ
ix ixix
ix
jx
jx jx jx
0=ijρ 0=ijρ 0=ijρ 0=ijρ 0=ijρ 0=ijρ
ix
jx
1+=ijρix
jx
1+=ijρix
jx
1−=ijρ ix
jx
1−=ijρ 10 << ijρ 10 << ijρ
Figura A.1 – Interpretação gráfica de coeficiente de correlação.
Intervalo de ρ Grau de dependência
0.0 - 0.3 Baixo
0.3 - 0.5 Médio
0.5 - 0.7 Importante
0.7 - 0.9 Forte
0.9 - 1.0 Muito Forte
Tabela A.1 – Grão de dependência das variáveis aleatórias.
180
Características do coeficiente de correlação:
• Os valores do coeficiente de correlação estão limitados entre -1 e +1.
• O coeficiente de correlação é um valor único para o espaço amostral.
• O coeficiente de correlação de uma variável e ela mesma é igual a um 1, =xxρ
• A permutação de variáveis não altera o resultado de coeficiente de correlação
se os mesmos pares de valores forem mantidos: xyyx ,, ρρ = .
• Se as variáveis X e Y forem estatisticamente independentes, então o
coeficiente de correlação destas variáveis será igual a zero.
• Se o resultado do coeficiente de correlação das variáveis X e Y for igual a zero,
não se pode afirmar que as duas variáveis sejam estatisticamente
independentes.
• Para n variáveis aleatórias dependentes entre si, a matriz coeficiente de
correlação ρ é expressa pela seguinte equação:
⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢
⎣
⎡
=
nnnn XXXXXX
XXXX
XX Simétrico
,,,
,,
,
21
2212
11
ρρρ
ρρρ
L
MOMMρ (A.16)
A.6. Caracterização de Variáveis Aleatórias
A caracterização de variáveis aleatórias é feita por intermédio de funções
que usam as distribuições de probabilidade sobre o espaço amostral. Estas funções
são basicamente duas, a função densidade de probabilidade e a função distribuição
de probabilidade.
A.6.1. Função Densidade de Probabilidade (PDF)
181
A função densidade de probabilidade é uma função matemática contínua
que permite descrever as variáveis aleatórias. Esta função satisfaz as seguintes
condições:
• A função densidade de probabilidade é sempre maior ou igual a zero.
0≥)x(p (A.17)
• A área abaixo da função densidade de probabilidade é igual a um.
∫+∞
∞−
=1dx)x(p (A.18)
• A probabilidade de que um valor da variável aleatória observada esteja entre
dois valores numéricos xa e xb é denotada por ]xXxPr[ ba ≤≤ e é expressa
como:
∫=≤≤b
a
x
xba dx)x(p]xXxPr[ (A.19)
Qualquer função matemática que satisfaça as 3 condições anteriores pode
ser considerada uma função densidade de probabilidade (PDF).
Existem várias funções que são usadas freqüentemente como função
densidade de probabilidade, como exemplo, a função densidade de probabilidade
normal é apresentada nesta seção e na Tabela A.1, são apresentados as funções
mais usadas.
A função matemática da função densidade probabilidade normal é expressa
pela Equação (A.20) e sua representação gráfica é mostrada na Figura A.2.
⎩⎨⎧
−∞>>∞≤≤∞−
=−
−
bcx
cxp e c
bx
,021)(
2)(21
π (A.20)
onde:
bX = Média
182
c=σ Desvio padrão
)(xp
xa
)(ap
AreaaP =)(
b
)(xp
xa
)(ap
AreaaP =)(
b Figura A.2 – Função densidade de probabilidade normal.
A.6.2. Função Distribuição de Probabilidade
A função distribuição de probabilidade é conhecida também na literatura
como função de distribuição cumulativa CDF (Cumulative Distribution Function),
esta função é definida como a integral da função densidade de probabilidade(PDF)
e indica a probabilidade de que uma variável aleatória X tenha um valor menor
ou igual que a (Equação A.21).
∫∞−
=≤=a
dxxpaXaP )(]Pr[)( (A.21)
onde:
)(aP Função distribuição de probabilidade
Pr[X ≤ a] Probabilidade de que a variável X seja menor ou igual que a
p(x) Função densidade de probabilidade
183
Uma função distribuição de probabilidades deve de satisfazer as seguintes
condições:
0.0)( =−∞P (A.22)
1)(0 ≤≤ XP (A.23)
1)( =∞P (A.24)
Para o caso da função densidade de probabilidade normal da Figura A.2, sua
função distribuição de probabilidade é expressa pela Equação A.25, e sua
representação gráfica é mostrada pela Figura A.3.
∫∞−
⎟⎠⎞
⎜⎝⎛ −
−=
xc
bz
dzec
xP2
21
21)(π
(A.25)
x
)(xP
0.1
a
)(aP
x
)(xP
0.1
a
)(aP
Figura A.3 – Função distribuição de probabilidade normal.
Na Tabela A.2, apresenta-se um resumo das funções densidade
probabilidade (PDF) e as funções distribuição de probabilidade (CDF) mais
usadas.
Tipo III (min.)(Weibull)
Tipo II(máximos)
Tipo I(mínimo)
Tipo I (Max)(Gumbel)
Uniforme
Rayleigh
Exponencial
Lognormal
Normal
E(X)P(X)p(x)
Desvio PadrãoMédiaCDFPDFTipo de
Distribuição
Tipo III (min.)(Weibull)
Tipo II(máximos)
Tipo I(mínimo)
Tipo I (Max)(Gumbel)
Uniforme
Rayleigh
Exponencial
Lognormal
Normal
E(X)P(X)p(x)
Desvio PadrãoMédiaCDFPDFTipo de
Distribuição
⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛ −
−2
21exp
21
σµ
πσx
⎟⎠⎞
⎜⎝⎛ −
Φσµx
⎟⎠⎞
⎜⎝⎛ + 2
21exp ξλ
)( XVar
µ σ
⎟⎟
⎠
⎞
⎜⎜
⎝
⎛⎟⎟⎠
⎞⎜⎜⎝
⎛ −−
2)ln(
21exp
21
ξλ
πξx
x ⎟⎟⎠
⎞⎜⎜⎝
⎛ −Φ
ξλ)ln( x
1)exp()( 2 −ξXE
)exp( xλλ − )exp(1 xλ−−λ1
λ1
⎟⎟
⎠
⎞
⎜⎜
⎝
⎛⎟⎟⎠
⎞⎜⎜⎝
⎛−
2
2 21exp
RR
xxσσ ⎟
⎟
⎠
⎞
⎜⎜
⎝
⎛⎟⎟⎠
⎞⎜⎜⎝
⎛−−
2
21exp1
R
xσ Rσπ
⎟⎟⎠
⎞⎜⎜⎝
⎛
2 Rσπ⎟⎟⎠
⎞⎜⎜⎝
⎛−
22
ab −1
abax
−−
2ba +
12ab −
( )( ))(exp)(exp uxux −−−−− ααα ( )( ))(expexp ux −−− αα5772.0
+uα
π6
( )( ))(exp)(exp uxux −−− ααα ( )( ))(expexp1 ux −−− αα5772.0
−uα
π6
⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛−⎟
⎠⎞
⎜⎝⎛
+ kk
xv
xv
vk exp
1
⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛−
k
xvexp ⎟
⎠⎞
⎜⎝⎛ −Γ
kv 11
21
1121 2⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛ −Γ−⎟
⎠⎞
⎜⎝⎛ −Γ
kkv
⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛−⎟
⎠⎞
⎜⎝⎛
+ kk
vx
vx
vk exp
1
⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛−−
k
vxexp1 ⎟
⎠⎞
⎜⎝⎛ +Γ
kv 11
21
1121 2⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠⎞
⎜⎝⎛ +Γ−⎟
⎠⎞
⎜⎝⎛ +Γ
kkv
Tabela A.2 – Funções densidade probabilidade (PDF) e de distribuição de probabilidade (CDF) ( Lopes, 2007).
Na Tabela A.2, os valores de k e a função Gamma )(kΓ nas distribuições
tipo II e tipo III, são avaliadas pelas seguintes equações:
09.1−= ik δ (A.26)
∫∞
−−=Γ0
1)( dxxek kx (A.27)
onde iδ , é o coeficiente variação (seção A.4.6)
A.6.3. Coeficiente de Inclinação de uma Distribuição
O coeficiente de inclinação de uma amostra indica a inclinação de sua
distribuição de freqüências. Matematicamente, é definido como a relação entre o
momento estatístico central de terceira ordem dividido pelo cubo do desvio
padrão.
3
3
ix
iXInclinaçãode.Coef
σ= (A.28)
Interpretação:
• Se o coeficiente de inclinação for igual a zero, a distribuição de freqüências
será simétrica.
• Se o coeficiente de inclinação for negativo, a distribuição de freqüências terá
inclinação esquerda ou negativa.
• Se o coeficiente de inclinação for positivo, a distribuição de freqüências terá
inclinação direita ou positiva.
186
A.6.4. Coeficiente de Curtose
O coeficiente de curtose de uma amostra permite comparar a distribuição de
freqüências da amostra com a distribuição normal. Este coeficiente é definido
como o momento estatístico central de quarta ordem, dividido pela quarta
potência do desvio padrão.
4x
4i
i
XCurtosede.Coef
σ= (A.29)
Se duas distribuições de freqüências têm a mesma dispersão e inclinação,
isso não será suficiente para supor que as duas distribuições têm a mesma forma.
Interpretação:
• Se o coeficiente de curtose for igual a zero, então a distribuição de freqüências
será a própria distribuição normal.
• Se o coeficiente de curtose for negativo, então a distribuição será achatada,
plana.
• Se o coeficiente de curtose for positivo, então a distribuição de freqüências
será concentrada ao redor da média, distribuição com pico.
A.7. Esperança Matemática de uma Variável Aleatória
A esperança matemática é uma operação que consiste na integração
ponderada de uma função de variável aleatória.
∫+∞
∞−
>=< dx)x(p)x(g)X(gE (A.30)
onde:
g(X) Função da variável aleatória
187
p(x) Função de ponderação
E<g(X)> Esperança matemática da função de variável aleatória
Esta operação pode ser visualizada na Figura A.4, que mostra que a
esperança de uma função de variável aleatória é igual a área abaixo da curva
definida por g(x) vezes p(x).
Casos especiais da esperança matemática:
• A esperança da variável aleatória é sua média aritmética.
Xdx)x(xpXE =>=< ∫+∞
∞−
(A.31)
• A esperança matemática de uma variável aleatória multiplicada por uma
constante é a média multiplicada pela constante.
XCXCE >=< (A.32)
• A esperança do quadrado da variável aleatória é igual ao quadrado da média
da amostra.
222 Xdx)x(pxXEa
=>=< ∫+
∞−
(A.33)
• A variância da variável aleatória é definida como:
2222 )X(xEdx)x(p)Xx()XX(E)X(Vara
−><=−>=−<= ∫+
∞−
(A.34)
• O momento estatístico central de ordem m é definido como:
∫+∞
∞−
−>=−< dx)x(p)Xx()XX(E mm (A.35)
188
(a)
(b)
(c)
(a)
(b)
(c) Figura A.4 – (a) Função de variável aleatória; (b) Função de ponderação; (c) Esperança
matemática da função de variável aleatória.
A esperança matemática para uma função geral g dependente de n variáveis
aleatórias X1, X2,...Xn e sua correspondente função densidade de probabilidade
p(x1,x2,...,xn) é definida como:
∫ ∫+∞
∞−
+∞
∞−
>=< nnnn dx...dx)x,...,x,x(p)x,...,x,x(g...)X,...,X,X(gE 1212121 (A.36)
Casos especiais mais usados:
• Covariância:
∫ ∫+∞
∞−
+∞
∞−
−−>=−−<= kjkjkkjjkkjjkj dxdx)x,x(p)Xx)(Xx(.)XX)(XX(E)X,X(Cov (A.37)
• Coeficiente de correlação:
>−<>−<
>−−<=
22 )()(
))((
kkjj
kkjjjk
XXEXXE
XXXXEρ (A.38)
Recommended