Upload
buikhue
View
213
Download
0
Embed Size (px)
Citation preview
INSTITUTO DE ENGENHARIA DE ESTRUTURAS,
TERRITÓRIO E CONSTRUÇÃO
Placa-Leg – Programa de Aplicação de um Modelo Híbrido-Misto
de Tensão à Análise Elastoplástica de Placas
Luís Mendes; Luís M.S.S. Castro
– Julho de 2006 –
Relatório ICIST
DTC 08/2006
I C I S T Av. Rovisco Pais, 1049-001 Lisboa Tel: 218 418 244
i
ÍNDICE
1 INTRODUÇÃO ............................................................................................................. 1
2 FORMULAÇÃO DE PLACAS .................................................................................... 1 2.1 Introdução.............................................................................................................. 1 2.2 Definição do problema .......................................................................................... 2 2.3 Relações de equilíbrio ........................................................................................... 3 2.4 Relações de compatibilidade ................................................................................. 4 2.5 Relações constitutivas ........................................................................................... 4
2.5.1 Fase elástica...................................................................................................... 5 2.5.2 Fase plástica...................................................................................................... 6
2.6 Integração das relações constitutivas..................................................................... 9
3 MODELO HÍBRIDO-MISTO DE TENSÃO ........................................................... 11 3.1 Introdução............................................................................................................ 11 3.2 Definição da aproximação ................................................................................... 11 3.3 Relações de equilíbrio ......................................................................................... 11 3.4 Relações de compatibilidade ............................................................................... 12 3.5 Relações constitutivas ......................................................................................... 13
3.5.1 Fase elástica.................................................................................................... 13 3.5.2 Fase plástica.................................................................................................... 13
3.6 Sistema governativo ............................................................................................ 15 3.6.1 Modelo elástico .............................................................................................. 15 3.6.2 Modelo elastoplástico..................................................................................... 15
3.7 Análise incremental ............................................................................................. 16 3.7.1 Processo iterativo – Método de Newton-Raphson ......................................... 18 3.7.2 Organização do sistema governativo.............................................................. 20
4 IMPLEMENTAÇÃO DO MODELO ........................................................................ 23 4.1 Definição da aproximação ................................................................................... 23
4.1.1 Funções de aproximação ................................................................................ 23 4.1.2 Comentários.................................................................................................... 24
4.2 Controlo da cedência ........................................................................................... 25 4.2.1 Controlo em pontos de colocação................................................................... 25 4.2.2 Controlo em células plásticas ......................................................................... 26
4.3 Critérios de cedência ........................................................................................... 27 4.3.1 Critério de Von Mises .................................................................................... 28 4.3.2 Critério de Drucker-Prager ............................................................................. 29
4.4 Cálculo dos operadores estruturais ...................................................................... 29 4.4.1 Transformação de coordenadas ...................................................................... 30 4.4.2 Operador de flexibilidade ............................................................................... 32 4.4.3 Operador de compatibilidade no domínio ...................................................... 34 4.4.4 Operador de compatibilidade na fronteira ...................................................... 38 4.4.5 Operador N∗ ................................................................................................... 41 4.4.6 Operador M∗ .................................................................................................. 42 4.4.7 Vector das forças de massa............................................................................. 43 4.4.8 Vector das forças na fronteira......................................................................... 45
4.5 Cálculo dos campos de tensões e deslocamentos ................................................ 47 4.6 Algoritmo do modelo elastoplástico.................................................................... 47
ii
5 PROGRAMA DE CÁLCULO AUTOMÁTICO ...................................................... 52 5.1 Introdução............................................................................................................ 52 5.2 Estrutura do ficheiro de dados ............................................................................. 53 5.3 Funcionamento do programa............................................................................... 54
5.3.1 Configuração .................................................................................................. 54 5.3.2 Executar o programa....................................................................................... 56
6 EXEMPLOS DE APLICAÇÃO................................................................................. 62
7 REFERÊNCIAS BIBLIOGRÁFICAS ...................................................................... 75
ANEXOS ANEXO A – POLINÓMIOS DE LEGENDRE
ANEXO B – APROXIMAÇÃO DOS INCREMENTOS DOS PARÂMETROS
PLÁSTICOS
ANEXO C – CRITÉRIOS DE CEDÊNCIA
iii
ÍNDICE DE FIGURAS
Figura 2.1: Elemento de placa............................................................................................................ 2
Figura 3.1: Estrutura em consola exemplificativa............................................................................ 21
Figura 4.1: Técnicas de controlo da cedência. ................................................................................. 26
Figura 4.2: Representação de uma célula plástica e dos pontos de integração, nos diversos
referenciais considerados. ................................................................................................................ 27
Figura 4.3: Representação gráfica dos critérios de cedência [Zienkiewicz et al. 1991]................... 29
Figura 4.4: Transformação de coordenadas para um elemento genérico trapezoidal de 4 nós. ....... 30
Figura 4.5: Funções de aproximação da geometria do elemento. .................................................... 30
Figura 4.6: Ordenamento dos termos de um bloco do operador de flexibilidade F. ........................ 34
Figura 4.7: Ordenamento dos termos de um bloco do operador vA . ............................................... 38
Figura 4.8: Elemento mestre utilizado no programa. ....................................................................... 40
Figura 4.9: Ordenamento dos termos de um bloco do operador Aγ . ............................................... 41
Figura 4.10: Representação esquemática do método de Newton-Raphson e do método de Newton-
Raphson modificado......................................................................................................................... 43
Figura 4.11: Ordenamento dos termos do vector das forças de massa............................................. 45
Figura 4.12: Ordenamento dos termos do vector das forças das forças de fronteira........................ 46
Figura 4.13: Algoritmo do incremento de carga. ............................................................................. 50
Figura 4.14: Algoritmo do processo iterativo. ................................................................................. 51
Figura 5.1: Menu principal do programa.......................................................................................... 54
Figura 5.2: Interface do módulo de configuração. ........................................................................... 56
Figura 5.3: Interface do módulo de pré-processamento. .................................................................. 57
Figura 5.4: Interface do módulo de visualização do problema em análise. ..................................... 58
Figura 5.5: Interface do módulo de processamento. ........................................................................ 59
Figura 5.6: Interface do módulo de pós-processamento................................................................... 61
Figura 6.1: Exemplo 1 – Definição das características da placa quadrada em consola. .................. 62
Figura 6.2: Exemplo 2 – Definição das características da viga bi-encastrada. ................................ 62
iv
Figura 6.3: Exemplo 1 – Visualização global-elastoplástico dos resultados, =0.215λ . .................. 64
Figura 6.4: Exemplo 1 – Visualização global-elastoplástico dos resultados, =0.45186λ ................ 65
Figura 6.5: Exemplo 1 – Visualização deslocamentos (uv/ug), =0.45186λ . .................................... 66
Figura 6.6: Exemplo 1 – Visualização das células plásticas e dos pontos de integração, =0.45186λ .
.......................................................................................................................................................... 66
Figura 6.7: Exemplo 1 – Visualização dos campos de tensões, ao longo do corte GHI, =0.45186λ .
.......................................................................................................................................................... 67
Figura 6.8: Exemplo 1 – Visualização dos campos de deslocamentos uv, ao longo do corte GHI,
=0.45186λ . ....................................................................................................................................... 67
Figura 6.9: Exemplo 1 – Visualização dos campos de deslocamentos uγ, ao longo do corte
GHI, =0.45186λ . ............................................................................................................................... 67
Figura 6.10: Exemplo 2 – Visualizações de diversas grandezas estruturais, =0.12λ . ..................... 70
Figura 6.11: Exemplo 2 – Visualização dos campos de tensões e deslocamentos ao longo do corte
CD, =0.12λ . ..................................................................................................................................... 71
Figura 6.12: Exemplo 2 – Visualizações de diversas grandezas estruturais, =0.17567λ . ................ 72
Figura 6.13: Exemplo 2 – Visualização dos campos de tensões e deslocamentos ao longo do corte
CD, =0.17567λ ................................................................................................................................. 73
v
ÍNDICE DE TABELAS
Tabela 5.1: Listagem das rotinas de cálculo do programa. .............................................................. 52
Tabela 5.2: Listagem dos ficheiros adicionais. ................................................................................ 53
Tabela 5.3: Listagem dos ficheiros de configuração das saídas gráficas. ........................................ 53
Tabela 5.4: Formato do ficheiro de dados........................................................................................ 55
Tabela 6.1: Exemplo 1 – Ficheiro de dados. .................................................................................... 63
Tabela 6.2: Exemplo 2 – Ficheiro de dados. .................................................................................... 68
Tabela 6.3: Exemplo 2 – Ficheiro de dados (cont.).......................................................................... 69
1
1 INTRODUÇÃO
O presente relatório tem por objectivo servir como documento de suporte ao programa de
cálculo denominado Placa-Leg. O programa foi desenvolvido no âmbito da dissertação de
Mestrado em Engenharia de Estruturas [Mendes 2002] e contou com o apoio do Núcleo de
Análise de Estruturas, do Instituto de Engenharia de Estruturas, Território e Construção
do Instituto Superior Técnico (ICIST).
O programa é caracterizado por recorrer ao modelo de elementos finitos híbrido misto de
tensão para efectuar a análise elastoplástica de placas (estados planos de tensão e de
deformação) e por utilizar polinómios de Legendre para efectuar a aproximação das
grandezas estruturais. Neste programa encontram-se implementados o critério de cedência
de Von Mises [Von Mises 1913] e o critério de Drucker-Prager [Drucker et al. 1952]. O
controlo da cedência pode ser efectuado através da utilização da técnica dos pontos de
colocação, da definição de células plásticas ou de ambas as técnicas em simultâneo.
Para além desta introdução, o relatório encontra-se organizado em 6 secções. Nas primeiras
secções apresentam-se os aspectos gerais da formulação de placas e uma descrição sucinta
do modelo híbrido misto de tensão (HMT). Posteriormente, apresentam-se as
características da aproximação e define-se o cálculo dos operadores estruturais. Por último,
apresentam-se alguns aspectos práticos relacionados com a utilização do programa e ainda
dois exemplos de aplicação. Este documento termina com a apresentação das referências
bibliográficas e dos anexos.
2 FORMULAÇÃO DE PLACAS
2.1 Introdução
Os elementos de placa são utilizados para modelar situações nas quais o comportamento
estrutural se pode considerar como um estado plano de tensão (EPT) ou como um estado
plano de deformação (EPD).
Esta secção tem como objectivo definir de forma sucinta as variáveis e as relações
fundamentais para os elementos de placa. Textos mais aprofundados sobre o assunto
podem ser encontrados em diversas publicações [Reddy 1993; Zienkiewicz et al. 1991].
2
2.2 Definição do problema
Considere-se o elemento de placa representado na Figura 2.1, constituído por um material
homogéneo e isotrópico, com um domínio V, delimitado por uma fronteira Γ , e definido
num referencial cartesiano.
A fronteira pode ser dividida na fronteira estática, σΓ , onde são conhecidas as tensões
aplicadas, e na fronteira cinemática, uΓ , onde são impostos deslocamentos.
Figura 2.1: Elemento de placa.
O elemento de placa encontra-se sujeito a um conjunto de forças de massa e de tensões
aplicadas na fronteira. Definem-se os vectores das forças de massa b e das tensões na
fronteira t através das suas componentes cartesianas:
x
y
bb
b⎡ ⎤
= ⎢ ⎥⎣ ⎦
, (2.1)
x
y
tt
t⎡ ⎤
= ⎢ ⎥⎣ ⎦
. (2.2)
O campo de tensões em qualquer ponto fica definido pelo vector onde se reúnem as
componentes independentes do tensor das tensões:
xx
yy
xy
σσ σ
σ
⎡ ⎤⎢ ⎥= ⎢ ⎥⎢ ⎥⎣ ⎦
. (2.3)
Analogamente, as componentes do tensor das deformações podem ser agrupadas no vector:
xx
yy
xy
εε ε
γ
⎡ ⎤⎢ ⎥= ⎢ ⎥⎢ ⎥⎣ ⎦
. (2.4)
3
O campo de deslocamentos do elemento de placa fica definido em cada ponto através das
componentes do vector dos deslocamentos:
x
y
uu
u⎡ ⎤
= ⎢ ⎥⎣ ⎦
. (2.5)
2.3 Relações de equilíbrio
A condição de equilíbrio no domínio pode ser obtida impondo o equilíbrio de um volume
infinitesimal, obtendo-se [Arantes e Oliveira 1969]:
, 0ij i jbσ + = . (2.6)
Tratando-se de um problema de elasticidade plana, as equações de equilíbrio reduzem-se a:
, ,
, ,
00
xx x xy y x
xy x yy y y
bb
σ σσ σ
+ + =⎧⎨ + + =⎩
. (2.7)
As relações anteriores podem ser escritas na forma matricial, através de:
D 0bσ + = , (2.8)
onde o operador diferencial de equilíbrio D fica definido por:
( ) ( )0
( ) ( )0
x yD
y x
∂ ∂⎡ ⎤⎢ ⎥∂ ∂⎢ ⎥=
∂ ∂⎢ ⎥⎢ ⎥∂ ∂⎣ ⎦
. (2.9)
A condição de equilíbrio na fronteira estática impõe que o campo de tensões esteja em
equilíbrio com as tensões exteriores aplicadas. Recorrendo à fórmula de Cauchy, obtém-se
[Arantes e Oliveira 1969]:
ij i jn tγσ = . (2.10)
No caso das placas, as equações anteriores reduzem-se a:
x xx y xy x
y yy x xy y
n n t
n n t
σ σ
σ σ
+ =⎧⎪⎨ + =⎪⎩
, (2.11)
ou na forma matricial:
N tσ = , (2.12)
onde N representa a matriz que reúne as componentes da normal exterior unitária à
fronteira, vindo definida por:
4
0
0x y
y x
n nN
n n⎡ ⎤
= ⎢ ⎥⎣ ⎦
. (2.13)
2.4 Relações de compatibilidade
Na hipótese dos pequenos deslocamentos e das pequenas deformações, as condições de
compatibilidade no domínio podem ser expressas na forma [Timoshenko et al. 1982]:
( ), ,12ij i j j iu uε = + . (2.14)
Para problemas de elasticidade plana, as relações de compatibilidade definidas em (2.14)
podem ser expressas da seguinte forma:
,
,
, ,
xx x x
yy y y
xy x y y x
uu
u u
ε
ε
γ
⎧ =⎪
=⎨⎪ = +⎩
, (2.15)
ou na forma matricial:
D uε ∗= , (2.16)
onde o operador diferencial de compatibilidade D∗ vem definido por:
( ) 0
( )0
( ) ( )
x
Dy
y x
∗
⎡ ⎤∂⎢ ⎥∂⎢ ⎥
⎢ ⎥∂= ⎢ ⎥∂⎢ ⎥⎢ ⎥∂ ∂⎢ ⎥∂ ∂⎣ ⎦
. (2.17)
A condição de fronteira cinemática é definida por:
x
y
uu
uγ
⎡ ⎤= ⎢ ⎥⎣ ⎦
, (2.18)
onde uγ representa o campo de deslocamentos imposto.
2.5 Relações constitutivas
As relações constitutivas permitem relacionar as tensões (campo estático) com as
deformações (campo cinemático).
Considerando materiais com comportamento elastoplástico, as relações constitutivas
apresentam um primeiro troço elástico, onde a deformação total tε é elástica e portanto
5
reversível eε , e um segundo troço elastoplástico, caracterizado pelo aparecimento de
deformações plásticas irreversíveis, pε . Pode então escrever-se:
Troço elástico: t eε = ε , (2.19)
Troço elastoplástico: t e pε = ε + ε . (2.20)
A decomposição das deformações totais numa parcela elástica e numa outra plástica é uma
das características da definição cinemática da plasticidade [Almeida 1991]. Este tipo de
formulação ajusta-se bem ao modelo de elementos finitos utilizado, o qual utiliza uma
formulação de flexibilidade para as relações constitutivas elásticas.
2.5.1 Fase elástica
Na fase elástica, a relação constitutiva do material é do tipo linear, podendo ser definida no
formato de flexibilidade (2.21), ou no formato de rigidez (2.22) [Arantes e Oliveira 1969]:
θε = f σ + ε , (2.21)
pσ = k ε + σ . (2.22)
Nas equações anteriores, θε e pσ representam respectivamente os campos de deformações
e de tensões instalados antes de se iniciar o carregamento, podendo ser utilizados para
permitir a consideração de deformações térmicas e tensões residuais.
As matrizes f e k reúnem os parâmetros elásticos do material. Tendo em conta as hipóteses
de base consideradas, estas matrizes são simétricas e positivas definidas, podendo ser
escritas na forma:
( )
( )
( )
1 -ν -ν 0 0 0
-ν 1 -ν 0 0 0
-ν -ν 1 0 0 0
0 0 0 2 1 + ν 0 0
0 0 0 0 2 1 + ν 0
0 0 0 0 0 2 1 + ν
1f =E
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
, (2.23)
(1 - ν) ν ν 0 0 0
ν (1 - ν) ν 0 0 0
ν ν (1 - ν) 0 0 0
(1 - 2 ν)E 0 0 0 0 02(1 + ν)(1 - 2ν)
(1 - 2 ν)0 0 0 0 0
2(1 - 2 ν)
0 0 0 0 02
k =
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
, (2.24)
6
onde E representa o módulo de elasticidade e ν o coeficiente de Poisson.
2.5.2 Fase plástica
De seguida apresenta-se a formulação adoptada para o comportamento de material em fase
plástica. Informações mais detalhadas sobre este tema podem ser encontradas em diversos
trabalhos publicados [Chen 1994; Kachanov 1974; Zienkiewicz et al. 1991].
Para caracterizar o comportamento do material em fase plástica é necessário definir:
1. A Condição de Cedência, que permite definir o limite a partir do qual se iniciam as
deformações plásticas irreversíveis.
2. A Lei de Endurecimento, que permite definir a evolução da superfície de cedência ao
longo do processo de deformação.
3. A Lei de Escoamento, que permite definir a direcção e o sentido do incremento da
deformação plástica.
4. A Condição de Consistência, que permite definir os incrementos dos potenciais
plásticos e da deformação plástica, para que não se viole a condição de cedência e a
lei de endurecimento.
5. As Condições de Complementaridade, que permitem relacionar a evolução do campo
estático e do campo cinemático plástico.
Hipóteses base:
No modelo implementado considerou-se uma lei de escoamento associada, o que significa
que se utilizou a mesma função para definir a cedência e os potenciais plásticos. Esta lei é
frequentemente adoptada nos trabalhos do domínio da plasticidade [Castro 1996]. No
entanto, alguns autores consideram-na inadequada para simular o comportamento de
materiais não metálicos, como o betão e alguns solos [Chen 1994].
Considerou-se também como válido o Postulado de Drucker, também conhecido como
critério de estabilidade do material. Este critério estabelece que “para que um material seja
estável, o trabalho realizado pelas deformações plásticas, durante um ciclo fechado de
tensões, é sempre não negativo”. Informações mais detalhadas podem ser encontradas em
[Mendes 2002].
7
Condição de Cedência:
A condição de cedência permite definir o limite a partir do qual se iniciam as deformações
plásticas. Matematicamente, materializa-se na denominada função de cedência, f:
( )ˆ ˆ 0f σ = σ - σ , (2.25)
onde σ̂ representa a tensão de comparação associada à condição de cedência e 0σ
representa um parâmetro, ou um conjunto de parâmetros, que caracterizam a envolvente de
cedência do material (ex. tensão de cedência à tracção ou ao corte, ângulo de atrito interno,
coesão), usualmente determinados por via experimental. A função de cedência toma
valores negativos quando o material se encontra em regime elástico e valor nulo quando o
material se encontra em cedência. Valores positivos representam um domínio inacessível
para o material.
A condição de cedência fica definida impondo que a função de cedência tome valores não
positivos:
0f ≤ , (2.26)
Lei de Endurecimento:
O endurecimento não foi implementado nesta versão do programa. Informação sobre como
incorporar o endurecimento no modelo pode ser encontrada em outros trabalhos dos
autores [Castro 1996; Mendes 2002].
Lei de Escoamento:
A lei de escoamento define a direcção e o sentido do incremento de deformação plástica,
não fornecendo no entanto qualquer indicação sobre a sua intensidade.
A lei de escoamento define-se como:
pij *
ij
gdε = dεσ∂∂
, (2.27)
pelo que a direcção do incremento das deformações plásticas é normal à superfície
associada ao potencial plástico g, definida por ( )0g = . Por esta razão, a lei de escoamento
é muitas vezes denominada por condição de normalidade.
8
A variável escalar *ε determina a intensidade do incremento. Devido à obrigatoriedade do
incremento da deformação plástica ser sempre não negativo, também *ε terá de tomar
valores não negativos, definindo-se desta a forma a condição de escoamento plástico:
*dε 0≥ . (2.28)
Para assegurar ainda a verificação do postulado de Drucker, o sentido do incremento das
deformações plásticas deve ser tal que assegure um valor não-negativo para o trabalho.
Conclui-se então que os incrementos da deformação plástica devem ser normais e
exteriores à superfície de cedência [Chen 1994; Mendes 2002].
Condição de Consistência:
Durante um processo de deformação plástica, independentemente de haver endurecimento,
o estado de tensão tem de permanecer “sobre” a superfície de cedência, o que se traduz
por:
( ), 0df σ κ = , (2.29)
ou seja:
ijij
f fdf = dσ + dκ= 0σ κ∂ ∂∂ ∂
. (2.30)
Analisando a expressão anterior, é possível estabelecer a seguinte conclusão para o caso de
não haver endurecimento 0fκ∂
=∂
:
0ijij
f dσσ∂
= ⇒∂
o vector incremento de tensão é tangente à superfície de cedência.
Condições de Complementaridade:
Definem-se as condições de complementaridade para relacionar a evolução do campo
estático f , com o campo cinemático plástico, pε , tendo em consideração as situações de
carregamento possíveis. Estas condições podem-se definir através de [Castro 1996]:
pf dε = 0 , (2.31)
pdf dε = 0 . (2.32)
Informações mais detalhadas podem ser encontradas em [Castro 1996; Mendes 2002].
9
Nesta fase aproveita-se para simplificar a notação utilizada, tirando partido do facto de se
ter adoptado uma lei de escoamento associada. Define-se então de forma indiscriminada
através de φ∗ , a função de cedência e a função dos potenciais plásticos:
= f = gφ∗ . (2.33)
Outra notação muito utilizada na literatura da especialidade, diz respeito à representação
através de n∗ do vector das normais à superfície definida pelo valor nulo da função dos
potenciais plásticos:
ij
gn =σ∗∂∂
, (2.34)
2.6 Integração das relações constitutivas
As relações constitutivas na fase plástica definidas anteriormente encontram-se expressas
em termos das suas variações infinitesimais. Para efectuarmos uma análise elastoplástica
incremental, necessitamos de as definir em termos de incrementos finitos. Para isso
admite-se que o carregamento é aplicado e a plasticidade se desenvolve num intervalo de
tempo fictício suficientemente pequeno , , 1,j jτ τ∗ ∗ +⎡ ⎤⎣ ⎦ , denominado intervalo de tempo
convencional de plasticidade [Castro 1996].
Em cada intervalo, o incremento das deformações plásticas é definido por:
*, j+1
*, jp p *
τ
τΔε = ε dτ∫ . (2.35)
Salienta-se que, tendo em conta as hipóteses adoptadas, nenhuma das relações definidas é
dependente do tempo. Assim sendo, o termo pε está na realidade associado às suas
variações no incremento de carga ( )p p *dε = ε dτ e não a uma variação no tempo
propriamente dita.
Introduzindo a lei de escoamento (2.27), a qual estabelece a forma como se processam as
variações infinitesimais das deformações plásticas, obtém-se:
*, j+1
*, jp * * *
τ
τΔε = n ε dτ∫ . (2.36)
O modo como se efectua a integração indicada em (2.36) é de extrema importância, tendo
influência não só nos resultados obtidos, como na própria estruturação do modelo
elastoplástico. Nos trabalhos de Crisfeld [Crisfield 1991] e Zienkiewicz et al. [Zienkiewicz
10
et al. 1991] encontra-se abordado com maior profundidade este assunto. Informações mais
detalhadas podem também ser encontradas em [Mendes 2002].
Várias estratégias podem ser utilizadas para o cálculo de (2.36), podendo ser classificados
em métodos explícitos e implícitos [Zienkiewicz et al. 1991]. Os métodos implícitos
utilizam grandezas que são apenas determinadas durante ou após o incremento de carga.
Assim sendo, os dados base que permitem definir e iniciar o incremento de carga, são
produto do próprio incremento de carga, levando a que estes métodos necessitem estar
associados a algoritmos iterativos, tipo Newton-Raphson.
O método implícito mais utilizado é o método implícito de Euler, na designação anglo-
saxónica backward-Euler [Crisfield 1991]. Este método consiste em corrigir, após cada
incremento de carga, os pontos situados fora da superfície de cedência, através de
alterações aos campos de tensões e da introdução de deformações plásticas proporcionais à
normal da superfície de cedência no final do incremento. Gera-se então um processo
iterativo que termina quando a correcção final conduz a um estado de tensão que satisfaça
a condição de cedência, ou quando a diferença seja tão pequena quanto se queira.
Aplicando o método implícito de Euler para efectuar a integração (2.36), obtém-se:
p *, j+1 *Δε = n Δε , (2.37)
ficando as deformações plásticas no final do passo de carga, definidas por:
p, j+1 p, j p p, j *, j+1 *ε = ε +Δε = ε + n Δε . (2.38)
A condição de escoamento plástico (2.39), de cedência (2.40), de consistência (2.41) e de
complementaridade (2.42) e (2.43) escritas na forma incremental, são definidas por [Castro
1996]:
*Δε 0≥ , (2.39)
* *+Δ 0φ φ ≤ , (2.40)
t* *, j+1Δ = Δn σφ , (2.41)
* *Δε = 0φ , (2.42)
* *Δ Δε = 0φ . (2.43)
11
3 MODELO HÍBRIDO-MISTO DE TENSÃO
3.1 Introdução
Esta secção tem como objectivo introduzir de uma forma sucinta o modelo de elementos
finitos utilizado para a resolução do problema genérico formulado na secção anterior.
Informação mais detalhada sobre o modelo híbrido-misto de tensão pode ser encontrada
nas publicações [Castro 1996; Freitas et al. 1999].
O modelo apresentado não engloba os operadores que permitem considerar o
endurecimento. Noutros trabalhos dos autores apresenta-se a formulação do modelo que
engloba o endurecimento [Castro 1996].
3.2 Definição da aproximação
Neste modelo são aproximados simultaneamente e de forma independente os campos de
tensões (3.1) e de deslocamentos no domínio (3.2) de cada elemento. São também
aproximados os incrementos dos parâmetros plásticos (3.3) no domínio e o campo de
deslocamentos na fronteira estática (3.4) [Castro 1996].
pσ = S X +σ em Vv , (3.1)
v v v pu =U q +u em V , (3.2)
* * *Δε = P Δe em V , (3.3)
γ γ σu =U q em Γ . (3.4)
As grandezas duais associadas às variáveis discretas introduzidas são as deformações
generalizadas e , as forças de massa generalizadas vQ , os incrementos dos potenciais
plásticos generalizados *ΔΦ , e as forças na fronteira generalizadas Qγ , definidas por
[Castro 1996]:
tve S dVε= ∫ , (3.5)
tv vQ U b dV= ∫ , (3.6)
t* * *ΔΦ = P Δφ dV∫ , (3.7)
tQ U t dγ γ γ σ= Γ∫ . (3.8)
3.3 Relações de equilíbrio
A condição de equilíbrio no domínio é dada por [Castro 1996]:
12
tv v pA X Q Q= − − , (3.9)
onde:
(D )tv v vA S U dV= ∫ , (3.10)
tv vQ U b dV= ∫ , (3.11)
Dtp v pQ U dVσ= ∫ . (3.12)
A condição de equilíbrio na fronteira estática pode-se definir através de [Castro 1996]:
tpA X Q Qγ γ γ= + − , (3.13)
onde:
(N )tvA S U dγ γ σ= Γ∫ , (3.14)
tQ U t dγ γ γ σ= Γ∫ , (3.15)
Ntp pQ U dγ γ σσ= Γ∫ . (3.16)
A expressão (3.13) corresponde à imposição ponderada da condição de fronteira estática.
Pode ser encarada também como resultando da imposição ponderada do equilíbrio entre
elementos adjacentes, definindo-se neste contexto a fronteira estática como abrangendo
também as fronteiras inter-elementares.
3.4 Relações de compatibilidade
A condição de compatibilidade no domínio é definida por [Castro 1996]:
A Av v ppe q q e eγ γ γ= − + + − , (3.17)
onde:
(N )tv ue S u dγ γ= Γ∫ , (3.18)
(D )tpp v pe S u dV= ∫ . (3.19)
A condição de fronteira cinemática é imposta localmente. A continuidade dos
deslocamentos numa fronteira comum a dois elementos é igualmente imposta localmente,
quando se garante que estes partilham a aproximação para o campo de deslocamentos
definida ao longo da fronteira comum.
13
3.5 Relações constitutivas
Conforme foi visto anteriormente, no domínio elastoplástico o campo de deformações
totais recebe a contribuição de uma parcela elástica reversível e de uma parcela plástica
irreversível (2.20). O mesmo sucede com as variáveis discretas do modelo de elementos
finitos:
t e pe e e= + , (3.20)
Através do mesmo critério energético utilizado para definir (3.5), obtém-se a definição de
cada uma das parcelas de deformação [Castro 1996]:
te v ee S dVε= ∫ , (3.21)
tp v pe S dVε= ∫ . (3.22)
3.5.1 Fase elástica
As relações constitutivas na fase elástica (2.21) são impostas na forma dos resíduos
pesados, utilizando-se para efectuar a ponderação as funções de aproximação do campo de
esforços no domínio [Castro 1996]:
( ) 0tv eS f dVθε σ ε− − =∫ . (3.23)
Introduzindo a aproximação (3.1) e a definição (3.21):
t t te v v v p ve S f S X dV S f dV S e dVθσ= + +∫ ∫ ∫ , (3.24)
que pode ser escrita na forma:
e pee F X e eθ= + + , (3.25)
onde:
F tv vS f S dV= ∫ , (3.26)
tpe v pe S f dVσ= ∫ , (3.27)
tve S dVθ θε= ∫ . (3.28)
3.5.2 Fase plástica
A lei de escoamento (2.37) é imposta ponderadamente, utilizando as funções de
aproximação do campo de esforços no domínio [Castro 1996]:
( )p , 1 0tv jS n dVε ε∗ + ∗Δ − Δ =∫ . (3.29)
14
Introduzindo a aproximação (3.3), obtém-se a igualdade:
, 1 0t tv p v jS dV S n P e dVε ∗ + ∗ ∗Δ − Δ =∫ ∫ , (3.30)
que pode ser escrita na forma:
pe N e∗ ∗Δ = Δ , (3.31)
onde:
, 1tv jN S n P dV∗ ∗ + ∗= ∫ . (3.32)
A condição de escoamento plástico (2.39) define-se no modelo de elementos finitos,
através da substituição da variável contínua pela discreta correspondente [Castro 1996]:
0e∗Δ ≥ . (3.33)
A condição de cedência (2.48) é imposta na forma de resíduos pesados, utilizando-se as
funções de aproximação dos incrementos dos parâmetros plásticos para efectuar a
ponderação [Castro 1996]:
( ) 0tP dVφ φ∗ ∗ ∗+ Δ ≤∫ . (3.34)
Tendo em conta a definição (3.7), obtemos:
0∗ ∗Φ + ΔΦ ≤ . (3.35)
A condição de consistência (2.41) é imposta na forma de resíduos pesados, utilizando as
funções de aproximação dos incrementos dos parâmetros plásticos para efectuar a
ponderação [Castro 1996]:
( ), 1 0t tjP n dVφ σ∗ ∗ ∗ +Δ − Δ =∫ . (3.36)
Tendo em conta a definição dos potenciais plásticos generalizados (3.7), introduzindo a
aproximação do incremento do campo de tensões (3.1) e a aproximação dos incrementos
dos parâmetros plásticos (3.3), obtém-se:
, 1t t
j vP n S X dV∗ ∗ ∗ +ΔΦ = Δ∫ , (3.37)
que pode ser escrito na forma:
tN X∗ ∗ΔΦ = Δ , (3.38)
As condições de complementaridade (2.42) e (2.43), são impostas em média [Castro 1996]
através de:
0dVφ ε∗ ∗Δ =∫ , (3.39)
15
0dVφ ε∗ ∗Δ Δ =∫ . (3.40)
Introduzindo a aproximação dos incrementos dos parâmetros plásticos (3.3) obtém-se:
0P e dVφ∗ ∗ ∗Δ =∫ , (3.41)
0P e dVφ∗ ∗ ∗Δ Δ =∫ , (3.42)
que conduz a:
0t e∗ ∗Φ Δ = , (3.43)
0t e∗ ∗ΔΦ Δ = . (3.44)
3.6 Sistema governativo
3.6.1 Modelo elástico
A condição de compatibilidade no domínio (3.17) e as relações constitutivas (3.25), estão
ambas escritas em termos das deformações generalizadas, podendo ser agrupadas numa
única equação [Castro 1996]:
F A Av v pe ppX q q e e e eγ γ γ θ+ − = + − − − . (3.45)
Juntando as condições de equilíbrio no domínio (3.9) e na fronteira (3.13), obtém-se o
sistema governativo elementar:
F A A
A 0 0A 0 0
v pe pptv v v pt
p
X e e e eq Q Qq Q Q
γ γ θ
γ γ γ γ
⎡ ⎤ ⎡ ⎤⎡ ⎤− − − −⎢ ⎥ ⎢ ⎥⎢ ⎥ = − −⎢ ⎥ ⎢ ⎥⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥− − +⎣ ⎦ ⎣ ⎦⎣ ⎦
(3.46)
Salienta-se que devido à preservação da dualidade estática-cinemática e da reciprocidade
nas relações constitutivas no modelo discreto, a matriz do sistema governativo elástico é
simétrica.
3.6.2 Modelo elastoplástico
O modelo é definido no instante final de um passo de carga , , 1,j jτ τ∗ ∗ +⎡ ⎤⎣ ⎦ , a que
designaremos de forma simplificada por instante ( )1j + . Posteriormente, será redefinido
de forma a facilitar a implementação da análise incremental e iterativa.
O sistema governativo elastoplástico pode ser obtido do sistema elástico, adicionando as
condições do modelo elastoplástico. A primeira diferença em relação ao modelo elástico
16
reside no facto das deformações generalizadas totais te , agora englobarem também uma
parcela de deformações plásticas pe .
Desta forma, no final do passo de carga a expressão que engloba a imposição ponderada
das relações constitutivas, fica definida por:
1 , 1 , 1 , 1t j pe j j p je F X e e eθ+ + + += + + + . (3.47)
Estas relações podem-se juntar novamente com as relações de compatibilidade (3.17),
obtendo-se:
1 , 1 , 1 , 1 , 1 , 1 , 1 , 1F A Aj v v j j j pp j j pe j p jX q q e e e e eγ γ γ θ+ + + + + + + ++ − = + − − − − . (3.48)
As relações de equilíbrio (3.9) e (3.13), são introduzidas no sistema através de:
1 , 1 , 1tv j v j p jA X Q Q+ + += − − , (3.49)
1 , 1 , 1t
j j p jA X Q Qγ γ γ+ + +− = − − . (3.50)
É ainda necessário que se verifique a condição de cedência (3.35), escrita no final do
incremento de carga.
Matematicamente, o sistema governativo no instante ( )1j + , fica traduzido pelo seguinte
conjunto de equações e inequações:
1 , 1 , 1 , 1 , 1 , 1
, 1 , 1 , 1
, 1 , 1 , 1
*, 1 * *, 1
- - - - -0 0 - -
- 0 0 -
0
v j j pp j j pe j p jtv v j v j p jt
j j p j
tj j
F A A X e e e e eA q Q QA q Q Q
P dV
γ γ θ
γ γ γ γ
φ
+ + + + + +
+ + +
+ + +
+ +
⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥+⎣ ⎦ ⎣ ⎦⎣ ⎦
Φ = ≤∫
. (3.51)
3.7 Análise incremental
Conforme foi referido, pretende-se efectuar uma análise onde o carregamento é aplicado
em sucessivos patamares de carga. Salienta-se que por carregamento se entende não só as
tensões ou cargas impostas à estrutura, mas também os deslocamentos ou deformações.
A análise incremental e iterativa é necessária para acompanhar o aparecimento dos
fenómenos plásticos e devido à utilização do método implícito de Euler na integração das
relações constitutivas. Neste método é adoptado um valor fixo da normal à superfície de
cedência durante todo o patamar de carregamento , 1jn∗ + , levando a que seja necessário
definir o carregamento a aplicar em incrementos suficientemente pequenos para que a
qualidade da solução não seja afectada.
17
O sistema definido em (3.51) para o instante 1j + , pode ser escrito em função de cada uma
das grandezas definidas no instante inicial do passo de carga j , acrescidas de um
incremento associado à nova etapa de carregamento.
As deformações plásticas generalizadas na forma incremental, são definidas por:
, 1 ,p j p j pe e e+ = + Δ . (3.52)
Introduzindo a lei de escoamento (3.31), obtém-se:
, 1 ,p j p je e N e+ ∗ ∗= + Δ . (3.53)
A imposição ponderada das relações de constitutivas no final do passo de carga, fica então
definida por:
1 , 1 , 1 ,t j pe j j p je F X e e e N eθ+ + + ∗ ∗= + + + + Δ . (3.54)
Estas relações podem combinar-se novamente com as relações de compatibilidade (3.17),
obtendo-se:
1 , 1 , 1 , 1 , 1 , 1 , 1 ,j v v j j j pp j j pe j p jF X A q A q N e e e e e eγ γ γ θ+ + + ∗ ∗ + + + ++ − + Δ = + − − − − . (3.55)
A condição de cedência fica definida por:
, 0j∗ ∗Φ + ΔΦ ≤ . (3.56)
O processo de carregamento deve ser tal que se verifiquem também as condições de
complementaridade (3.43) e (3.44), e a condição de escoamento plástico (3.33), definidas
respectivamente por:
( ), 0; 0
0
t
j
t
ee
e∗ ∗ ∗
∗
∗ ∗
⎧ Φ + ΔΦ Δ =⎪ Δ ≥⎨⎪ΔΦ Δ =⎩
. (3.57)
O sistema governativo elastoplástico na forma incremental fica então definido pelo
seguinte conjunto de equações e inequações:
( )
* , , , , ,,
, ,,
, ,*
*, * **, *
* *
- - - - - - - -0 0 0 - - - -
- 0 0 0 - -
00; ;
0
jv j pp j pp j pe j pe p j
v j vtv v j v p j p
jtj p j p
t
jj
t
X XF A A N e e e e e e e e e
q qA Q Q Q Q
q qA Q Q Q Q
e
e
e
γ γ γ θ θ
γ γγ γ γ γ γ
+Δ⎡ ⎤⎡ ⎤ ⎡ ⎤+Δ Δ Δ Δ⎢ ⎥+Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ = Δ Δ⎢ ⎥ ⎢ ⎥⎢ ⎥+Δ⎢ ⎥ ⎢ ⎥Δ + + Δ⎢ ⎥ ⎣ ⎦⎣ ⎦ Δ⎣ ⎦
⎧ Φ + ΔΦ Δ =⎪Φ + ΔΦ ≤ ⎨⎪ΔΦ Δ =⎩
* 0eΔ ≥
,(3.58)
18
3.7.1 Processo iterativo – Método de Newton-Raphson
Conforme já foi referido, na integração das relações constitutivas é necessário incorporar
um método iterativo. Foi implementado um dos algoritmos iterativos mais utilizados, o
método de Newton-Raphson. Neste método, a estimativa do valor de uma variável x que
satisfaça determinada equação ( ) 0f x = , pode ser obtida em cada iteração i através de:
( )( )
1i
i ii
f xx x
f x+ = −
′. (3.59)
Para a análise elastoplástica efectuada, aplicou-se este conceito muito simples generalizado
a um sistema de equações de várias variáveis. Para evitar inverter a matriz que reúne as
primeiras derivadas, é conveniente escrever o sistema de equações da seguinte forma:
( ) { } ( ){ }1i i i if x x x f x+⎡ ⎤′ − = −⎣ ⎦ , (3.60)
O algoritmo do método de Newton-Raphson resume-se a efectuar recursivamente
estimativas lineares definidas recorrendo às primeiras derivadas das relações fundamentais
do modelo, de forma a eliminar os erros (restos ou resíduos) inerentes ao facto das
estimativas serem lineares e obter desta forma a solução para o incremento de carga.
Conforme já foi referido, os restos indicam a diferença entre a solução obtida no final de
cada iteração, e aquela que respeita as condições fundamentais do problema. O seu cálculo
consiste na utilização das equações do sistema governativo presentes em (3.58), definidas
na forma da procura de raízes. Desta forma, no final da iteração i os restos são calculados
através de:
i1 , , , , 1
, 1 , 1 , 1
t2 v , 1 , 1
t3 , 1 , 1
4 1
=it it it
it
it
N N Ni i i i i
j v v j v j p j ji i i
pp j j pe j
Ni i
j v j p ji
Ni i
j j p ji
i ij
R F X X A q q A q q N e e e
e e e
R A X X Q Q
R A X X Q Q
R
γ γ γ γ
θ
γ γ γ
∗ ∗ +
+ + +
+ +
+ +
+
+ Δ + + Δ − + Δ + Δ + − +
+ + +
= + Δ + +
= − + Δ + −
=
⎧ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠ ⎝ ⎠
⎛ ⎞⎜ ⎟⎝ ⎠⎛ ⎞⎜ ⎟⎝ ⎠
Φ
∑ ∑ ∑
∑
∑
⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩
. (3.61)
Conforme foi referido, a verificação da condição de escoamento plástico é garantida à
partida e as condições de complementaridade referem-se apenas à forma como o modelo
19
elastoplástico deve evoluir, pelo que não é necessária a sua inclusão tanto no cálculo dos
restos, como no sistema de equações do método de Newton-Raphson.
O sistema de equações do método de Newton-Raphson é definido de forma a eliminar ou a
minimizar os restos. É escrito na forma apresentada em (3.60) para os incrementos das
variáveis do sistema governativo. É no entanto necessário incluir um conjunto de equações
que permitam definir a evolução dos potenciais plásticos, e não apenas calcular os seus
valores no final de cada iteração, como na expressão (3.61). Para este efeito recorre-se à
condição de consistência (3.38), que é escrita apenas em função dos incrementos das
variáveis. Desta forma obtemos a seguinte definição alternativa para o cálculo no final de
cada iteração de 4R :
4 , , 1i i t
j j jR N X∗ ∗ ∗ ∗ += Φ + ΔΦ = Φ + Δ . (3.62)
Introduzindo o conjunto de equações anteriores na definição dos restos, obtém-se um
sistema de equações determinado e simétrico.
Recordando a noção de diferencial de uma função de várias variáveis indicada no anexo I,
define-se o sistema de equações a adoptar em cada iteração, por:
( )
1 1 1 1
2 2 2 2 1
2
33 3 3 3
4
4 4 4 4
i i i i
v
i i i i i i
i ivi ii i i i
iv
i i i i
v
R R R RX q q e
R R R R X RX q q e q R
q RR R R RX q q e e R
R R R RX q q e
γ
γ
γ
γ
∗
∗
∗ ∗
∗
⎡ ⎤∂ ∂ ∂ ∂⎢ ⎥∂ ∂ ∂ ∂Δ⎢ ⎥⎢ ⎥ ⎡ ⎤∂ ∂ ∂ ∂ ⎡ ⎤Δ⎢ ⎥ ⎢ ⎥ ⎢ ⎥∂ ∂ ∂ ∂Δ Δ⎢ ⎥ ⎢ ⎥ ⎢ ⎥= −⎢ ⎥ ⎢ ⎥ ⎢ ⎥Δ∂ ∂ ∂ ∂⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥∂ ∂ ∂ ∂Δ Δ Δ ⎢ ⎥⎢ ⎥ ⎣ ⎦⎣ ⎦⎢ ⎥∂ ∂ ∂ ∂⎢ ⎥
⎢ ⎥∂ ∂ ∂ ∂Δ⎣ ⎦
. (3.63)
Na obtenção das derivadas indicadas nas equações anteriores, o cálculo menos evidente
corresponde à derivação do resto 1iR em relação aos esforços generalizados X [Bussamra
1999].
[ ] [ ]1R F X N eX X X ∗ ∗
∂ ∂ ∂= + Δ
∂ ∂ ∂, (3.64)
ou seja:
1R N eFX X
∗ ∗∂ ∂ Δ= +
∂ ∂. (3.65)
20
Desenvolvendo o segundo termo da equação anterior e omitindo por simplicidade todos os
índices, que indicariam tratar-se da iteração ( )i , obtém-se:
( )( )
( )( )
( )( ) ( )
tt t vv , 1 v , 1 1
SS SN j j j
dVn P e dV n dVeX X X X
φ εε σ∗
∗∗ + ∗ ∗ ∗ + ∗∗ ∗ +
⎛ ⎞∂∂ Δ⎜ ⎟⎜ ⎟∂ Δ ∂ Δ∂ Δ ∂⎝ ⎠= = =
∂ ∂ ∂ ∂
∫∫ ∫ . (3.66)
Tendo em conta:
( )( )
( )( )
( )( )
( )( )
i
vi i i iS
X X
σ
σ σ
∂∂ ∂ ∂= =
∂ ∂ ∂ ∂, (3.67)
podemos escrever (3.66), na forma [Bussamra 1999]:
( )( )
N= M
i ii
i
e
X∗ ∗
∗
∂ Δ
∂, (3.68)
onde:
2
** *2
1
i tv v
j
M S S dVφ εσ +
∂= Δ
∂∫ . (3.69)
As restantes derivadas são de cálculo imediato, pelo que a equação (3.63) pode ser escrita
na forma:
( ) ( )
* * 1
2
3
* 4*
-0 0 0
-- 0 0 0
0 0 0
i i i iv
t i iv vt i i
t i ii
F M A A N X RA q RA q R
e RN
γ
γ γ
⎡ ⎤+ ⎡ ⎤Δ ⎡ ⎤⎢ ⎥ ⎢ ⎥ ⎢ ⎥Δ⎢ ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ ⎢ ⎥Δ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ Δ Δ⎢ ⎥ ⎢ ⎥⎣ ⎦⎣ ⎦⎣ ⎦
. (3.70)
Usualmente, a matriz da equação anterior é denominada por Hessiana iH . A expressão
(3.70) representa o sistema de equações, a aplicar em cada iteração i, para um dado
incremento de carga 1j + .
3.7.2 Organização do sistema governativo
Para ilustrar a organização do sistema governativo global, apresenta-se a sua estrutura base
para o caso da placa em consola representada na Figura 3.1.
21
Figura 3.1: Estrutura em consola exemplificativa.
Os operadores e os vectores ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )i,m1 2 3 4, N , , , ,i i i, j i i j m
v γF ,A ,A R R R R∗ , que aparecem no
sistema governativo global, encontram-se associados aos (elementos, fronteiras, modos de
cedência activos) definidos respectivamente pelos índices ( ), ,i j m .
O sistema governativo apresentado em (3.71), representa a forma genérica do sistema de
equações a adoptar na resolução do problema em regime elástico (3.46):
( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
( )
( )
( )
( )
( ) ( )
( )
( )
( )
1 1 1,1 1,4 1,6v
2 2 2,2 2,4 2,5 2,7v
1v
2v
1,1
2,2
1,4 2,4
2,5
1,6
2,7
. . . . .
. . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
t
t
t
t
t t
t
t
t
F A -A -A -A
F A -A -A -A -A
A
A
-A
-A
-A -A
-A
-A
-A
γ γ γ
γ γ γ γ
γ
γ
γ γ
γ
γ
γ
⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢
⎣
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
1
2
1
1
1
2
4
5
6
7
v
v
X
X
q
q
q
q
q
q
q
q
γ
γ
γ
γ
γ
γ
⎤ ⎡ ⎤⎥ ⎢ ⎥⎥ ⎢ ⎥⎥ ⎢ ⎥⎥ ⎢ ⎥⎥ ⎢ ⎥⎥ ⎢ ⎥⎥ ⎢ ⎥⎥ ⎢ ⎥ =⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥
⎦ ⎣ ⎦
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( )
( ) ( )
( ) ( )
( ) ( )
( ) ( ) ( )
( ) ( )
( ) ( )
( ) ( )
1 1 1 1
2 2 2 2
1 1
2 2
1 1
2 2
4 4 4
5 5
6 6
7 7
pp pe
pp pe
v p
v p
p
p
p p
p
p
p
e e e ee e e e
Q QQ QQ QQ Q
Q Q Q
Q QQ QQ Q
γ θ
γ θ
γ γ
γ γ
γ γ γ
γ γ
γ γ
γ γ
− − −
− − −
− −
− −
− +=
− +
− + +
− +
− +
− +
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
. (3.71)
22
Na expressão (3.72), representa-se a forma genérica do sistema de equações do algoritmo
de Newton-Raphson, a adoptar na resolução do problema em regime elastoplástico (3.70).
Admite-se que se encontram activos dois modos de plastificação, um pertencente ao
elemento 1 e outro ao elemento 2. Omite-se por simplicidade os índices que indicariam
tratar-se da iteração i.
( ) ( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
( )
( )
( )
( )
( ) ( )
( )
( )
1 1 1 1,1 1,4 1,6 1,1* *
2 2 2 2,2 2,4 2,5 2,7 2,2* *
1
2
1,1
2,2
1,4 2,4
2,5
1,6
. . - . - . - . .
. . . - - - . - .
. . . . . . . . . . .
. . . . . . . . . . .
- . . . . . . . . . . .
. - . . . . . . . . . .
- - . . . . . . . . . .
. - . . . . . . . . . .
- . . . .
v
v
tv
tv
t
t
t t
t
t
F M A A A A N
F M A A A A A N
A
A
A
A
A A
A
A
γ γ γ
γ γ γ γ
γ
γ
γ γ
γ
γ
+
+
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )( )( )( )1
*
2*
1
2
1
1
1
2
4
5
6
2,7 7
1,1*
2,2*
. . . . . . .
. - . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
v
v
t
t
t
e
e
X
X
q
q
q
q
q
q
q
A q
N
N
γ
γ
γ
γ
γ
γ γ
Δ
Δ
Δ
Δ
⎡ ⎤ ⎡ ⎤Δ⎢ ⎥ ⎢ ⎥
Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥
Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥Δ⎢ ⎥ ⎢ ⎥
Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥
Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥Δ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎣ ⎦⎣ ⎦
=
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
11
21
12
22
13
23
43
53
63
73
14
24
R
R
R
R
R
R
R
R
R
R
R
R
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
= . (3.72)
Salienta-se a definição elementar dos operadores F e vA , enquanto que γA é definido para
cada deslocamento não impedido ao longo dos troços da fronteira estática. É neste
operador que se faz sentir a ligação entre os elementos, através da partilha da aproximação
23
para o campo de deslocamentos ao longo da fronteira comum. É o caso da fronteira 4 (7ª
coluna da matriz do sistema) onde se regista a participação dos dois elementos ( ) ( )1,4 2,4γ γA , A .
O operador N∗ e o vector e∗Δ , encontram-se definidos no domínio onde é efectuado o
controlo de cedência e englobam apenas os modos de plastificação que se encontram
activos.
4 IMPLEMENTAÇÃO DO MODELO
4.1 Definição da aproximação
4.1.1 Funções de aproximação
As funções de aproximação utilizadas são do tipo polinomial. A construção da
aproximação pode ser obtida, no caso de grandezas unidimensionais, agrupando os
polinómios por ordem crescente de grau, conforme se indica em (3.73). No caso de
aproximações em domínios bidimensionais, a construção da aproximação é resultante do
produto tensorial de polinómios unidimensionais (3.74), podendo ser agrupados com a
sequência indicada em (3.75):
( ) ( ) ( )0 1 nP x P x P x⎡ ⎤⎣ ⎦ , (3.73)
( )
( )
( ) ( )[ ] ( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
0 0 0 0 0 1 0
0 1
P x P y P y P x P y P x P y P x P y
P x P x P y P x P y P x P y
α α
α α α α α
=
⎡ ⎤ ⎡ ⎤⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦
, (3.74)
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )0 0 0 0P x P y P x P y P x P y P x P yα α α α⎡ ⎤⎣ ⎦ . (3.75)
Quanto à organização das matrizes que reúnem as funções de aproximação , ,v vS U Uγ ,
optou-se por agrupar primeiro todas as funções associadas a cada grandeza aproximada.
Estas matrizes são definidas com um número de linhas igual ao número de grandezas
aproximadas. Desta forma, a matriz (3.76) representa, por exemplo, o caso onde são duas
as grandezas aproximadas:
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )0 0
0 0
0 00 0
P x P y P x P yP x P y P x P y
α α
α α
⎡ ⎤⎢ ⎥⎣ ⎦
. (3.76)
Outras formas de organização são possíveis, influenciando no entanto a distribuição dos
termos no sistema governativo. Procurou-se garantir que a organização das matrizes
24
assegure uma distribuição eficiente dos coeficientes não-nulos, nomeadamente através da
optimização da proximidade à diagonal principal, conduzindo à minimização das
operações a efectuar quando se utilizam técnicas apropriadas para a resolução eficaz do
sistema de equações.
Utilizando aproximações com polinómios completos, a dimensão de qualquer uma das
matrizes anteriores é dada para o caso de uma aproximação unidimensional por (3.77), e
por (3.78) para o caso bidimensional:
( ), 1 maxn n g+⎡ ⎤⎣ ⎦ , (3.77)
( )2, 1 maxn n g⎡ ⎤+⎣ ⎦ , (3.78)
onde n representa o número de grandezas aproximadas e maxg o grau máximo das funções
consideradas nessa aproximação.
Utilizam-se para definir todas as aproximações séries completas de polinómios de
Legendre, os quais correspondem a soluções da equação diferencial de Legendre [Spiegel
et al. 1990]:
( ) ( ) ( )2(1 ) 2 ( 1) 0P x P n n Pγ γ γ γ′′ ′− − + + = . (3.79)
Estas funções polinomiais podem ser geradas recursivamente pela fórmula de Bonnet:
1 1( 1) ( ) (2 1) ( ) ( ) 0n n nn P n x P n Pγ γ γ+ −+ − + + = . (3.80)
Salienta-se que neste tipo de formulação se perde o conceito de elementos isoparamétricos,
pois utilizam-se funções diferentes para aproximar a geometria dos elementos e as
grandezas estruturais.
No Anexo A abordam-se com maior detalhe os aspectos relacionados com os polinómios
de Legendre.
4.1.2 Comentários
A escolha deste tipo de funções deve-se a duas razões principais: em primeiro lugar são
funções onde é possível explorar a ortogonalidade, o que conduz à geração de sistemas
governativos altamente esparsos, e consequentemente, assegura um menor esforço
computacional quando se implementam rotinas apropriadas ao tratamento e resolução deste
tipo de sistemas. Em segundo lugar, esta série completa de funções ortogonais está
associada a uma maior estabilidade numérica na definição das aproximações, o que
permite efectuar refinamentos tipo-p nos quais se pode aumentar de forma significativa o
25
grau máximo dos polinómios utilizados. Esta é uma vantagem importante se se tiverem em
consideração os problemas numéricos que geralmente ocorrem quando nos elementos
finitos de deslocamento se incrementa o grau das funções de interpolação de Lagrange. Por
outro lado, a exploração das propriedades dos polinómios ortogonais de Legendre
possibilita, na maioria das situações, a definição de expressões analíticas para o cálculo dos
operadores estruturais. Torna-se assim dispensável o recurso a quaisquer esquemas de
integração numérica, o que permite não só optimizar a velocidade de cálculo como também
maximizar a precisão numérica dos cálculos envolvidos.
4.2 Controlo da cedência
Neste trabalho foram utilizadas duas técnicas para o controlo da cedência: o controlo em
pontos de colocação e o controlo em células plásticas. Para que se verifique localmente a
condição de escoamento plástico (2.39) as funções que constituem a base de aproximação
P∗ são escolhidas de forma a apresentarem no seu domínio sempre valores não negativos.
Conforme já foi referido, o controlo da cedência é efectuado num referencial ( ),γ ψ ,
denominado por referencial local de plasticidade. É neste referencial que se definem as
grandezas associadas ao comportamento elastoplástico, como os potenciais plásticos ∗Φ ,
os incrementos dos parâmetros plásticos e∗Δ e o operador N∗ . Esta segunda mudança de
coordenadas foi efectuada por duas razões: em primeiro lugar se tivermos as células
plásticas definidas em intervalos [ ]1,1− a implementação dos métodos de integração
numérica (quadratura de Gauss ou de Lobatto) é facilitada e directa (ver Figura 4.2). A
segunda razão está relacionada com questões de coerência e de uniformidade nos
procedimentos de cálculo dos operadores estruturais.
4.2.1 Controlo em pontos de colocação
No método baseado na utilização de pontos de colocação, a plasticidade é controlada de
uma forma pontual. Foram utilizadas funções de Dirac ( )ˆ ,δ γ ψ [Bussamra 1999], como
funções de aproximação dos incrementos dos parâmetros plásticos. Cada uma destas
funções é definida num ponto, de forma a que ao efectuar-se o integral num dado domínio,
da função de Dirac multiplicada por outra função, obtém-se como resultado o valor nesse
ponto da segunda função. Assim sendo, para cada modo de cedência m, o potencial
26
plástico generalizado ∗Φ está apenas associado a um ponto ( )0 0,γ ψ do domínio (Figura
4.1a), tomando o valor da função dos potenciais plásticos nesse ponto,
( )0 0,m φ γ ψ∗ ∗=Φ . (3.81)
Considerou-se a possibilidade dos pontos de colocação se localizarem em ( )n x n pontos
de Gauss e de Lobatto e num conjunto de pontos equidistantes.
Figura 4.1: Técnicas de controlo da cedência.
4.2.2 Controlo em células plásticas
No controlo em células plásticas, foram implementados três conjuntos de funções de
aproximação para os incrementos dos parâmetros plásticos. O primeiro consiste na
consideração de uma única função constante de valor unitário, denominada adiante pelo
tipo CEL_1. O segundo conjunto foi definido utilizando quatro funções do mesmo tipo das
utilizadas na aproximação da geometria do elemento (Figura 4.1b), denominada por
CEL_4. Por último, considerou-se também um conjunto de 9 funções obtidas do produto
nas duas direcções de polinómios de segundo grau, CEL_9. Todas estas funções
encontram-se definidas no Anexo B.
O valor do potencial plástico ∗Φ associado ao modo m, pode ser calculado através de:
( )( ) ( ) ( )1 1
1 1
, , , 'tm mP J J d dγ ψ φ γ ψ γ ψ ψ γ∗ ∗ ∗
− −
Φ = ∫ ∫ , (3.82)
onde ( ),J γ ψ representa o Jacobiano de transformação de coordenadas entre o referencial
global e o local, enquanto que 'J representa o Jacobiano de transformação de
coordenadas entre o referencial local e local de plasticidade (ver Figura 4.2). Da análise da
Figura 4.2, facilmente se conclui que o Jacobiano não depende das coordenadas ( ),γ ψ ,
a)
b)
27
pois toma um valor constante igual à relação entre as áreas do elemento no referencial
local (2x2=4) e a área da célula plástica no referencial local de plasticidade.
Figura 4.2: Representação de uma célula plástica e dos pontos de integração, nos diversos referenciais
considerados.
A expressão (3.82) foi calculada numericamente, pois não é geralmente possível definir
matematicamente a forma como variam os potenciais plásticos no domínio da célula
plástica. Para efectuar a integração, utilizou-se o método da Quadratura de Gauss ou o
método da Quadratura de Lobatto, obtendo-se:
( ) ( ) ( )1 1
, , ' , 'n n
k kP W J Jγ ψ φ γ ψ γ ψ∗ ∗ ∗⎡ ⎤Φ = ⎣ ⎦∑∑ , (3.83)
onde n representa o número de pontos em cada direcção e 'W representa o peso associado
a cada ponto para a quadratura considerada.
4.3 Critérios de cedência
O comportamento dos materiais pode ser agrupado em dois grandes grupos: materiais
isotrópicos e materiais anisotrópicos. Os materiais isotrópicos podem ainda ser
classificados em materiais insensíveis à pressão hidrostática, caso da maioria dos metais, e
materiais sensíveis à pressão hidrostática, como os solos e o betão.
Diversos modelos foram desenvolvidos para caracterizar a cedência deste tipo de materiais.
Entre os mais simples encontram-se o critério de Tresca [Tresca 1864], o critério de von
Mises [Von Mises 1913], para os materiais insensíveis à pressão hidrostática, o critério de
Mohr-Coulomb [Coulomb 1873] e o critério de Drucker-Prager [Drucker et al. 1952], para
os materiais sensíveis à pressão hidrostática.
( ),J γ ψ 'J
28
Um dos objectivos do trabalho que originou o programa Placa-Leg consistia em incorporar
critérios de cedência que considerem a resistência combinada de atrito e coesão. Assim
sendo, para além do critério de Von Mises foi implementado o critério de Drucker-Prager.
O critério de Drucker-Prager pode ser calibrado para se aproximar do critério de Mohr-
Coulomb. Este critério é na maioria das situações preferível, pois tem uma definição
matemática mais atraente, permitindo contornar os problemas das arestas nas superfícies de
cedência.
Em seguida apresenta-se a definição genérica dos critérios de cedência incorporados no
programa Placa-Leg, deixando para o Anexo C a sua definição matemática mais
pormenorizada.
4.3.1 Critério de Von Mises
Este critério, inicialmente proposto por Von Mises em 1913 [Von Mises 1913], considera
que a cedência ocorre quando a tensão tangencial octaédrica, (C.6), atinge o valor de
23
kτ , onde kτ representa a tensão de cedência num estado de corte puro, ou seja num
estado de tensão onde o primeiro invariante do tensor das tensões seja nulo, 1 0I = .
Matematicamente, este critério fica definido por:
22 23 3
J kτ= , (3.84)
ou:
22J kτ= . (3.85)
No espaço das tensões principais, este critério fica definido por um prisma com uma elipse
como base, e com o eixo 1 2 3σ σ σ= = como geratriz (ver Figura 4.3).
O critério de Von Mises considera que a cedência do material depende apenas do segundo
invariante do tensor das tensões deviatóricas, 2J , não estando relacionada com a pressão
hidrostática instalada, 1I .
Calibrando o parâmetro kτ de modo a que a tensão de cedência num ensaio de tracção puro
1 2 3, 0, 0tσ σ σ σ= = = tome o valor cedσ , obtém-se:
3
cedkτσ
= . (3.86)
29
Introduzindo a relação anterior em (3.85), obtém-se a definição mais usual do critério de
Von Mises [Chen 1994; Zienkiewicz et al. 1991]:
( )2 ced 2,σ 3 0cedf J J σ= − = . (3.87)
Figura 4.3: Representação gráfica dos critérios de cedência [Zienkiewicz et al. 1991].
4.3.2 Critério de Drucker-Prager
O critério de Drucker-Prager foi inicialmente formulado em 1952 [Drucker et al. 1952],
resultando da modificação do critério de Von Mises, de forma a incluir o efeito da pressão
hidrostática, (ver Figura 4.3):
( )1 2 1 2, 0f I J I J kα= + − = , (3.88)
onde α e k são constantes que caracterizam o material, determinados de forma a que a
superfície cónica associada ao critério de Drucker-Prager, representada no espaço das
tensões principais, circunscreva a pirâmide hexagonal do critério de Mohr-Coulomb.
( )( )( )
( )( )( )
2 6;
3 3 3 3sen c cos
ksen senφ φ
αφ φ
= =− −
. (3.89)
No espaço das tensões principais, este critério fica definido por um prisma em forma de
cone com uma elipse como base, e com o eixo 1 2 3σ σ σ= = como geratriz.
4.4 Cálculo dos operadores estruturais
Nesta secção definem-se as expressões associadas ao cálculo dos operadores estruturais.
Salientam-se, neste contexto, os trabalhos de Pereira et al. [Pereira et al. 2000] e de Castro
[Castro 1996], que serviram como ponto de partida à implementação que aqui se reporta.
30
A aproximação é definida pelas variáveis , ,s uv ugg g g , que representam respectivamente o
grau máximo da aproximação dos campos de tensões, de deslocamentos no domínio e de
deslocamentos na fronteira.
4.4.1 Transformação de coordenadas
Os elementos utilizados são do tipo trapezoidal, com a geometria definida por 4 nós. Para
permitir a sistematização do cálculo dos operadores estruturais, efectua-se uma
transformação de coordenadas do referencial global ( ),x y para um elemento mestre
definido no referencial local ( ),ξ η . Foi implementado um elemento mestre definido em
[ ]1,1− (ver Figura 4.4).
Figura 4.4: Transformação de coordenadas para um elemento genérico trapezoidal de 4 nós.
Para a definição da geometria do elemento, considera-se uma aproximação do tipo bilinear
(ver Figura 4.5), definida por:
1
2
3
4
1 4(1 )(1 )
1 4(1 )(1 )
1 4(1 )(1 )
1 4(1 )(1 )
ψ ξ η
ψ ξ η
ψ ξ η
ψ ξ η
⎧ = − −⎪
= + −⎪⎨
= + +⎪⎪ = − +⎩
. (5.1)
Figura 4.5: Funções de aproximação da geometria do elemento.
1ψ 2ψ 3ψ 4ψ
31
A geometria do elemento nos eixos globais ( )1 2,x x pode ser definida na forma [Pereira
1993]:
[ ] [ ] [ ] [ ]1 2 3 41 2 3 4 , {1, 2}k k k k kx x x x x para kψ ψ ψ ψ= + + + = , (5.2)
onde [ ]ikx representa a coordenada k, definida no sistema de eixos global, do nó i do
elemento.
Salienta-se que, devido à forma como foi definida a aproximação da geometria do
elemento, a orientação dos lados locais é necessariamente definida no sentido anti-horário.
Introduzindo as aproximações definidas em (5.1) na expressão (5.2), obtém-se:
( )[ ] [ ] [ ] [ ]( ) [ ] [ ] [ ] [ ]( )
[ ] [ ] [ ] [ ]( ) [ ] [ ] [ ] [ ]( )
1 2 3 4 1 2 3 4
1 2 3 4 1 2 3 4
1 1, 4 4
1 14 4
k k k k k k k kk
k k k k k k k k
x x x x x x x xx
x x x x x x x x
ξξ η
η ξη
+ + + + − + + − +=+
− − + + + − + −. (5.3)
As equações anteriores podem ser escritas no formato [Pereira et al. 2000]:
( ) 0, , {1, 2}k k k k kx x para kξ η α η β ξ γ ξη= + + + = , (5.4)
onde:
[ ] [ ] [ ] [ ]( )[ ] [ ] [ ] [ ]( )[ ] [ ] [ ] [ ]( )
[ ] [ ] [ ] [ ]( )
1 2 3 40
1 2 3 4
1 2 3 4
1 2 3 4
1 4
1 4
1 4
1 4
k k k k k
k k k k k
k k k k k
k k k k
x x x x x
x x x x
x x x x
x x x xκ
α
β
γ
⎧ = + + +⎪⎪ = − − + +⎪⎨
= − + + −⎪⎪⎪ = − + −⎩
. (5.5)
A relação entre as variações infinitesimais das grandezas definidas nos sistemas de
coordenadas global e local, pode ser definida através da seguinte igualdade:
{ }, 1, 2i ii
x xdx d d iξ ηξ η∂ ∂
= + =∂ ∂
, (5.6)
ou na forma matricial:
1 1
1
2 2 2
x xdx ddx x x d
ξξ ηη
ξ η
∂ ∂⎡ ⎤⎢ ⎥∂ ∂⎡ ⎤ ⎡ ⎤⎢ ⎥=⎢ ⎥ ⎢ ⎥∂ ∂⎢ ⎥ ⎣ ⎦⎣ ⎦⎢ ⎥∂ ∂⎣ ⎦
. (5.7)
A matriz presente na equação anterior é usualmente denominada por matriz Jacobiana da
transformação de coordenadas. Para os elementos trapezoidais utilizados é definida por:
32
( ) 1 1 1 1
2 2 2 2
,Jβ γ η α γ ξ
ξ ηβ γ η α γ ξ
+ +⎡ ⎤= ⎢ ⎥+ +⎣ ⎦
. (5.8)
Define-se o Jacobiano da transformação como sendo o determinante da matriz Jacobiana,
o qual pode ser obtido através de:
1 1 2 2 1 1 2 2( , ) det[ ( , )] ( ) ( ) ( )( )J Jξ η ξ η β γ η α γ ξ α γ ξ β γ η= = + + − + + . (5.9)
Efectuando o desenvolvimento dos produtos e agrupando os termos, obtém-se:
1 2 1 2 1 2 1 2 1 2 2 1 1 2 1( , ) ( ) ( ) ( ) ( )J ξ η β α α β γ α α γ η β γ β γ ξ γ γ γ γ ξη2= − + − + − + − . (5.10)
Observando que o último termo se anula, é possível escrever a expressão anterior na forma
condensada [Pereira et al. 2000]:
0( , )J J J Jξ ηξ η ξ η= + + , (5.11)
onde:
0 1 2 1 2J β α α β= − , (5.12)
1 2 1 2Jη γ α α γ= − , (5.13)
1 2 2 1Jξ β γ β γ= − . (5.14)
4.4.2 Operador de flexibilidade
De acordo com o modelo de elementos finitos, o operador de flexibilidade é calculado para
cada elemento finito, e, através da expressão:
( )tev vF S f S dV= ∫ .
Efectuando a transformação de coordenadas para o elemento mestre definido no referencial
local ( ),ξ η , a definição do operador irá englobar o Jacobiano da transformação:
( )( ) ( ) ( )1 1
1 1
, , ,te
v vF S f S J d dξ η ξ η ξ η η ξ− −
= ∫ ∫ . (5.15)
Desta forma, o cálculo de um termo genérico do operador de flexibilidade, associado às
funções de aproximação definidas pelos índices i, j, m, n e ao termo klf da matriz que
reúne os parâmetros elásticos, pode ser efectuado a partir de:
( ) ( ) ( ) ( ) ( ) ( )1 1
-1 -1
, , , , , = ,ei j kl m nF i j k l m n P ξ P η f P ξ P η J ξ η d dη ξ∫ ∫ . (5.16)
33
Particularizando a expressão anterior para elementos trapezoidais, nos quais o Jacobiano da
transformação de coordenadas é definido em (5.11), um elemento genérico do operador de
flexibilidade pode ser calculado por:
( ) ( ) ( ) ( ) ( ) ( )1 1
01 1
, , , , ,ei j kl m nF i j k l m n P P f P P J J J d dξ ηξ η ξ η ξ η η ξ
− −
= + +∫ ∫ . (5.17)
A expressão anterior pode ser subdividida em três parcelas, relacionadas com os termos
que definem o Jacobiano. Tendo em consideração as expressões apresentadas no Anexo A,
obtém-se:
Termo 1:
( ) ( ) ( ) ( ) ( ) ( )1 1 1 1
0 01 1 1 1
,kl j n i m klf J P P d P P d f J se i m j nη η η ξ ξ ξ− − − −
= = =⎛ ⎞
∧⎜ ⎟⎝ ⎠∫ ∫ ∫ ∫ . (5.18)
Termo 2:
( ) ( ) ( ) ( )
( ) ( )
( ) ( )
1 1 1 1
1 1 1 1
, 12
, 12
kl j n i m
kli m
kli m
f J P P d P P d
iJ f se i m j n
mJ f se i m j n
ξ
ξ
ξ
η η η ξ ξ ξ ξ
λ λ
λ λ
− − − −
=
= + ∧ =
== − ∧ =
⎛ ⎞⎜ ⎟⎝ ⎠
⎧⎪⎪⎨⎪⎪⎩
∫ ∫ ∫ ∫
. (5.19)
Termo 3:
( ) ( ) ( ) ( )
( ) ( )
( ) ( )
1 1 1 1
1 1 1 1
, 12
, 12
kl j n i m
klj n
klj n
f J P P d P P d
jJ f se j n i m
nJ f se j n i m
η
ξ
ξ
η η η η ξ ξ ξ
λ λ
λ λ
− − − −
=
= + ∧ =
== − ∧ =
⎛ ⎞⎜ ⎟⎝ ⎠
⎧⎪⎪⎨⎪⎪⎩
∫ ∫ ∫ ∫
. (5.20)
Para cada termo não nulo de klf gera-se um bloco de termos quadrado, com dimensão
( )21 sg+ , no qual os termos anteriores se definem de acordo com o exemplo apresentado
na Figura 4.6, para 3sg = .
34
i 0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3m n / j 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 30 00 10 20 31 01 1 - TERMO 11 2 - TERMO 21 3 - TERMO 32 02 12 22 33 03 13 23 3
Figura 4.6: Ordenamento dos termos de um bloco do operador de flexibilidade F.
O operador de flexibilidade eF da estrutura é traduzido por uma matriz quadrada com a
dimensão:
º
2
13 (1 )
n elemis
ig
=
⎡ ⎤+⎣ ⎦∑ . (5.21)
4.4.3 Operador de compatibilidade no domínio
Conforme definido anteriormente, o operador de compatibilidade no domínio para o
elemento finito e, é definido pela expressão:
( )e tv v vA D S U dV= ∫ .
Para efectuar a transformação de coordenadas para o elemento mestre, adopta-se um
procedimento semelhante ao anterior. No entanto, é necessário calcular a relação entre o
operador de equilíbrio definido no sistema de coordenadas globais ( ),x y e o mesmo
operador definido nas coordenadas locais ( ),ξ η .
Recorrendo à regra da derivação da função composta (5.22), podemos relacionar para cada
elemento finito as derivadas parciais associadas às coordenadas globais, com as derivadas
associadas às coordenadas locais:
( ) ( ) ( )ξ ηγ ξ γ η γ
∂ ∂ ∂ ∂ ∂= +
∂ ∂ ∂ ∂ ∂, (5.22)
definindo-se o novo operador de equilíbrio ( ),D ξ η , da seguinte forma:
35
( ) ( ) ( ) ( )0( , )
( ) ( ) ( ) ( )0
x x y yD
y y x x
ξ η ξ ηξ η ξ η
ξ ηξ η ξ η
ξ η ξ η
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂⎡ ⎤+ +⎢ ⎥∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂⎢ ⎥=∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂⎢ ⎥+ +⎢ ⎥∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂⎣ ⎦
. (5.23)
Tendo em atenção que os termos , , ,x y x yξ ξ η η⎛ ⎞∂ ∂ ∂ ∂
⎜ ⎟∂ ∂ ∂ ∂⎝ ⎠ da expressão anterior podem ser
calculados com base nas expressões associadas à inversa da transformação de coordenadas
utilizada na definição da geometria do elemento, é possível relacionar as derivadas nos
dois sistemas de coordenadas, através de:
1
( )( )
( ) ( )x J
y
ξ
η
−
∂⎡ ⎤∂⎡ ⎤⎢ ⎥⎢ ⎥ ∂∂ ⎢ ⎥⎢ ⎥ =
∂ ∂⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥∂ ∂⎣ ⎦ ⎣ ⎦
, (5.24)
onde:
( )
, ,1
, ,
1,
y yJ
x xJη ξ
η ξξ η−
−⎡ ⎤= ⎢ ⎥−⎣ ⎦
. (5.25)
Desta forma, de acordo com as expressões apresentadas em (5.4) para os elementos
trapezoidais, o operador diferencial de equilíbrio pode ser definido através de:
( ) ( ) ( )1, ,,
eD x y DJ
ξ ηξ η
= , (5.26)
onde:
( ) 2 2 1 1 2 2 1 1
1 1 2 2 1 1 2 2
0 0( ) ( ),0 0
Dα γ ξ α γ ξ β γ η β γ η
ξ ηα γ ξ α γ ξ β γ η β γ ηξ η
+ − − + − −
− − + − − +
⎡ ⎤ ⎡ ⎤∂ ∂= −⎢ ⎥ ⎢ ⎥∂ ∂⎣ ⎦ ⎣ ⎦
.(5.27)
Num formato mais compacto, podemos definir:
( ) ( ) ( )( ) ( )
1 2
2 1
D , 0 D ,,
0 D , D ,D
ξ η ξ ηξ η
ξ η ξ η⎡ ⎤
= ⎢ ⎥⎣ ⎦
, (5.28)
onde:
( ) ( ) ( )
( ) ( ) ( )
1 2 2 2 2
2 1 1 1 1
( ) ( ),
( ) ( ),
D
D
ξ η α γ ξ β γ ηξ η
ξ η α γ ξ β γ ηξ η
+ +
− − − −
∂ ∂⎧ = −⎪ ∂ ∂⎪⎨ ∂ ∂⎪ = −⎪ ∂ ∂⎩
. (5.29)
36
O operador de compatibilidade no domínio pode ser então definido no referencial local
( ),ξ η , através de:
( )
( ) ( ) ( ) ( )1 1
1 1
1 , , , ,,
t
ev v vA D S U J d d
Jξ η ξ η ξ η ξ η η ξ
ξ η− −
⎛ ⎞⎜ ⎟=⎜ ⎟⎝ ⎠
∫ ∫ . (5.30)
Eliminando os termos referentes aos Jacobianos, o cálculo de um termo genérico do
operador de compatibilidade, associado às funções de aproximação definidas pelos índices i ,j, m, n, e ao termo do operador de equilíbrio kD definido em (5.29), pode ser efectuado
através de:
( ) ( ) ( ) ( )( ) ( ) ( )1 1
1 1
, , , , ,te
v k i j m nA i j k m n D P P P P d dξ η ξ η ξ η η ξ− −
= ∫ ∫ . (5.31)
Os termos associados a 1D , resultam do cálculo de:
( ) ( ) ( )( ) ( ) ( ) ( ) ( )( ) ( )
( )( ) ( ) ( ) ( )( )( ) ( ) ( ) ( )
1 1 1 1
2 21 1 1 1
1 1 1 1
2 21 1 1 1
i ij n m j n m
j jn i m n i m
P PP P d P d P P d P d
P PP d P P d P d P P d
ξ ξα η η η ξ ξ γ η η η ξ ξ ξ
ξ ξ
η ηβ η η ξ ξ ξ γ η η η ξ ξ ξ
η η
− − − −
− − − −
∂ ∂+ −
∂ ∂
∂ ∂−
∂ ∂
∫ ∫ ∫ ∫
∫ ∫ ∫ ∫. (5.32)
Tendo em consideração as expressões apresentadas no Anexo A, obtém-se:
( ) ( ) ( )( ) ( )
( ) ( ) ( ) ( )
1 1
21 1
22 , 10, . .
ij n m
i m se j n m i m i for par
PP P d P d
Termoc c
ξα η η η ξ ξ
ξ
α λ λ− −
= ∧ < ∧ +
∂=
∂
→⎧⎪⎨⎪⎩
∫ ∫, (5.33)
( ) ( ) ( )( ) ( )
( ) ( ) ( )( ) ( ) ( ) ( )
1 1
21 1
2
2
, 2
2 , 30, . .
ij n m
i m
se j n i m
se j n m i m i for par
PP P d P d
m Termo
Termoc c
ξγ η η η ξ ξ ξ
ξ
γ
γ λ λ
− −
= ∧ =
= ∧ < ∧ +
∂=
∂
→⎧⎪
→⎨⎪⎩
∫ ∫, (5.34)
( )( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
1 1
21 1
2 í2 , 4
0, . .
jn i m
j n se i m n j n j for mpar
PP d P P d
Termo
c c
ηβ η η ξ ξ ξ
η
β λ λ− −
= ∧ < ∧ +
∂− =
∂
− →⎧⎪⎨⎪⎩
∫ ∫, (5.35)
37
( ) ( )( ) ( ) ( )
( ) ( ) ( )( ) ( ) ( ) ( )
1 1
21 1
2
2
, 5
2 , 6
0, . .
j ni m
j n
se i m n j
se i m n j n j for par
P Pd P P d
n Termo
Termo
c c
η η ηγ η ξ ξ ξ
η
γ
γ λ λ
− −
= ∧ =
= ∧ < ∧ +
∂− =
∂
− →⎧⎪− →⎨⎪⎩
∫ ∫. (5.36)
De forma idêntica para os termos associados a 2D :
( ) ( ) ( )( ) ( ) ( ) ( ) ( )( ) ( )
( )( ) ( ) ( ) ( )( )( ) ( ) ( ) ( )
1 1 1 1
1 11 1 1 1
1 1 1 1
1 11 1 1 1
i ij n m j n m
j jn i m n i m
P PP P d P d P P d P d
P PP d P P d P d P P d
ξ ξα η η η ξ ξ γ η η η ξ ξ ξ
ξ ξ
η ηβ η η ξ ξ ξ γ η η η ξ ξ ξ
η η
− − − −
− − − −
∂ ∂−
∂ ∂
∂ ∂
∂ ∂
− +
+
∫ ∫ ∫ ∫
∫ ∫ ∫ ∫, (5.37)
obtendo-se:
( ) ( ) ( )( ) ( )
( ) ( ) ( ) ( )
1 1
11 1
12 , 10, . .
ij n m
i m se j n m i m i for par
PP P d P d
Termoc c
ξα η η η ξ ξ
ξ
α λ λ− −
= ∧ < ∧ +
∂− =
∂
− →⎧⎪⎨⎪⎩
∫ ∫, (5.38)
( ) ( ) ( )( ) ( )
( ) ( ) ( )( ) ( ) ( ) ( )
1 1
11 1
1
1
, 2
2 , 30, . .
ij n m
i m
se j n i m
se j n m i m i for par
PP P d P d
m Termo
Termoc c
ξγ η η η ξ ξ ξ
ξ
γ
γ λ λ
− −
= ∧ =
= ∧ < ∧ +
∂− =
∂
− →⎧⎪− →⎨⎪⎩
∫ ∫, (5.39)
( )( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
1 1
11 1
1 í2 , 4
0, . .
jn i m
j n se i m n j n j for mpar
PP d P P d
Termo
c c
ηβ η η ξ ξ ξ
η
β λ λ− −
= ∧ < ∧ +
∂=
∂
→⎧⎪⎨⎪⎩
∫ ∫, (5.40)
( )( ) ( ) ( ) ( )
( ) ( ) ( )( ) ( ) ( ) ( )
1 1
11 1
1
1
, 5
2 , 6
0, . .
jn i m
j n
se i m n j
se i m n j n j for par
PP d P P d
n Termo
Termo
c c
ηγ η η η ξ ξ ξ
η
γ
γ λ λ
− −
= ∧ =
= ∧ < ∧ +
∂=
∂
→⎧⎪
→⎨⎪⎩
∫ ∫. (5.41)
Para cada termo kD , gera-se um bloco de termos com dimensão ( ) ( )2 21 1s uvg x g+ + , com
a posição definida de acordo com o exemplo apresentado na Figura 4.7, para
3, 2s uvg g= = .
38
m 0 0 0 1 1 1 2 2 2i j / n 0 1 2 0 1 2 0 1 20 00 10 20 31 0 - TERMO 2 + TERMO 51 1 - TERMO 41 2 - TERMO 31 3 - TERMO 12 0 - TERMO 62 12 22 33 03 13 23 3
Figura 4.7: Ordenamento dos termos de um bloco do operador vA .
O operador de compatibilidade vA é traduzido por uma matriz com a dimensão:
º2
1
º2
1
º : 3 (1 )
º : 2 (1 )
n elemis
i
n elemiuv
i
N de linhas g
N de colunas g
=
=
⎡ ⎤+⎣ ⎦
⎡ ⎤+⎣ ⎦
∑
∑. (5.42)
4.4.4 Operador de compatibilidade na fronteira
De acordo com a secção relativa ao modelo de elementos finitos, o operador de
compatibilidade na fronteira associado ao elemento finito e, e à fronteira estática f, é
calculado através de:
( ), (N )e f tvA S U dγ γ σ= Γ∫ .
Para o cálculo deste operador, necessitamos da relação entre as componentes das normais
exteriores definidas no sistema de coordenadas globais ( ),x y e definidas em coordenadas
locais ( ),ξ η .
Prova-se que [Pereira 1993]:
( )( ) ( ) ( ) ( )1,
, ,,
x
y
nnJ J
nnξ
η
ξ ηξ η ξ η
ξ η−⎡ ⎤ ⎡ ⎤
=⎢ ⎥ ⎢ ⎥⎣ ⎦⎣ ⎦
. (5.43)
Atendendo a que:
( ) ( ) ( )1 , ,
, ,
1,,
y yJ
x xJη ξ
η ξ
ξ ηξ η
− −⎡ ⎤= ⎢ ⎥−⎣ ⎦
, (5.44)
a expressão (5.43) pode ser rescrita na forma:
39
( )( )
, ,
, ,
,,
x
y
y y nnx x nnη ξ ξ
η ξ η
ξ ηξ η
−⎡ ⎤ ⎡ ⎤ ⎡ ⎤=⎢ ⎥ ⎢ ⎥ ⎢ ⎥−⎣ ⎦ ⎣ ⎦⎣ ⎦
. (5.45)
De acordo com as expressões apresentadas para os elementos trapezoidais (5.4), obtém-se:
( )( )
2 2 2 2
1 1 1 1
,,
x
y
nnnnξ
η
ξ η α γ ξ β γ ηξ η α γ ξ β γ η
⎡ ⎤ + − − ⎡ ⎤⎡ ⎤=⎢ ⎥ ⎢ ⎥⎢ ⎥− − +⎣ ⎦ ⎣ ⎦⎣ ⎦
. (5.46)
Como as normais exteriores obtidas pela expressão (5.46) não são unitárias, é necessário
dividi-las pela sua norma, obtida através de:
( ) ( )2 2( , ) , ,x yn n nξ η ξ η ξ η= + , (5.47)
ficando o operador de compatibilidade na fronteira definido por:
( )( ) ( )
( ) ( ) ( )t
1
γ v γ2 2-1 x y
N ξ,ηA = S ξ,η U γ J γ dγ
n ξ,η +n ξ,η
⎛ ⎞⎜ ⎟⎜ ⎟⎝ ⎠∫ , (5.48)
onde:
( ) ( )
( ) ( ), 0 ,
0 , ,x y
y x
n nN
n nξ η ξ η
ξ η ξ η⎡ ⎤
= ⎢ ⎥⎣ ⎦
. (5.49)
Tendo em conta que o Jacobiano da transformação de coordenadas, pode ser obtido de:
( ) ( ) ( )2 2, ,x yJ n nγ ξ η ξ η= + , (5.50)
o operador de compatibilidade pode ser definido por:
( ) ( ) ( )( ) ( )1
,
1
, ,te f
vA N S U dγ γξ η ξ η γ γ−
= ∫ . (5.51)
Desta forma, um termo genérico do operador, associado às funções de aproximação
definidas pelos índices i, j, m, e à normal exterior ( ),kn ξ η , pode ser calculado através de:
( ) ( ) ( ) ( ) ( )( ) ( )1
,
1
, , , ,te f
k i j mA k i j m n P P P dγ ξ η ξ η γ γ−
= ∫ . (5.52)
De acordo com a expressão anterior, um termo genérico do operador de compatibilidade na
fronteira de Aγ , pode ser calculado através de:
( ) ( ) ( ) ( )( ) ( )1
,k
1
,te f
i j iA n P P P dγ ξ η ξ η γ γ−
= ∫ . (5.53)
Particularizando para cada um dos lados do elemento mestre definido Figura 4.8, obtem-se:
40
Figura 4.8: Elemento mestre utilizado no programa.
Lado I - ( )γ ξ= , ( )1η = − e ( ) ( ); 0 ; 1n nξ η = − :
( )
( )
I
I
2 2 2 2 2 2
1 1 1 1 1 1
01
x
y
nn
α γ ξ β γ β γα γ ξ β γ β γ
⎡ ⎤ + − + −⎡ ⎤ ⎡ ⎤⎡ ⎤= =⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥− − − − +−⎢ ⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦⎣ ⎦
, (5.54)
( )I 2 2 1 1
1 1 2 2
00
Nβ γ β γ
β γ β γ− − +⎡ ⎤
= ⎢ ⎥− + −⎣ ⎦, (5.55)
( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ), I I
1
1
1 1 ,e It
i j m jA N P P P d P N se i mγ ξ ξ ξ−
= − = − =∫ . (5.56)
Lado II - ( )γ η= , ( )1ξ = e ( ) ( ); 1 ; 0n nξ η = :
( )
( )
II
II
2 2 2 2 2 2
1 1 1 1 1 1
10
x
y
nn
α γ β γ η α γα γ β γ η α γ
⎡ ⎤ + − − +⎡ ⎤ ⎡ ⎤⎡ ⎤= =⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥− − + − −⎢ ⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦⎣ ⎦
, (5.57)
( )II 2 2 1 1
1 1 2 2
00
Nα γ α γ
α γ α γ+ − −⎡ ⎤
= ⎢ ⎥− − +⎣ ⎦, (5.58)
( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ), II II
1
1
1 1 ,e IIt
i j m iA N P P P d P N se j mγ η η η−
= = =∫ . (5.59)
Lado III - ( )γ ξ= , ( )1η = e ( ) ( ); 0 ; 1n nξ η = :
( )
( )
III
III
2 2 2 2 2 2
1 1 1 1 1 1
01
x
y
nn
α γ ξ β γ β γα γ ξ β γ β γ
⎡ ⎤ + − − − −⎡ ⎤ ⎡ ⎤⎡ ⎤= =⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥− − + +⎢ ⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦⎣ ⎦
, (5.60)
( )III 2 2 1 1
1 1 2 2
0N
0β γ β γ
β γ β γ− − +⎡ ⎤
= ⎢ ⎥+ − −⎣ ⎦, (5.61)
( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ), III III
1
1
N 1 1 N ,e IIIt
i j m jA P P P d P se i mγ ξ ξ ξ−
= = =∫ . (5.62)
Lado IV - ( )γ η= , ( )1ξ = − e ( ) ( ); 1 ; 0n nξ η = − :
41
( )
( )
IV
IV
2 2 2 2 2 2
1 1 1 1 1 1
10
x
y
nn
α γ β γ η α γα γ β γ η α γ
⎡ ⎤ − − − − +−⎡ ⎤ ⎡ ⎤⎡ ⎤= =⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥− + + −⎢ ⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦⎣ ⎦
, (5.63)
( )IV 2 2 1 1
1 1 2 2
00
Nα γ α γ
α γ α γ− + −⎡ ⎤
= ⎢ ⎥− − +⎣ ⎦, (5.64)
( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ), IV IV
1
1
1 1 ,e IVt
i j m iA N P P P d P N se j mγ η η η−
= − = =∫ . (5.65)
Para cada termo não nulo de N, gera-se um bloco de termos com dimensão
( ) ( )21 1s ug x g γ+ + , com a posição definida de acordo com o exemplo apresentado no
Figura 4.9, para ( )3, 2s ug g γ= = .
O operador de compatibilidade Aγ é traduzido por uma matriz com a dimensão:
º2
1
º . .
1
º : 3 (1 )
º : (1 )
n elemis
i
n g liu
i
N de linhas g
N de colunas g γ
=
=
⎡ ⎤+⎣ ⎦
⎡ ⎤+⎣ ⎦
∑
∑. (5.66)
i j / m 0 1 2 i j / m 0 1 20 0 0 00 1 0 10 2 0 20 3 0 31 0 1 01 1 1 11 2 1 2 - LADOS I e III1 3 1 32 0 2 0 - LADOS II e IV2 1 2 12 2 2 22 3 2 33 0 3 03 1 3 13 2 3 23 3 3 3
Figura 4.9: Ordenamento dos termos de um bloco do operador Aγ .
4.4.5 Operador N∗
Confome foi visto anteriormente, o operador N∗ é calculado para o elemento e, e para o
modo de cedência m, através de:
( ) ( ), te m e m mvN S n P dV∗ ∗ ∗= ∫ .
As expressões para o cálculo das normais à superfície dos potenciais plásticos n∗ ,
encontram-se definidas no Anexo C, para as diversas condições de cedência consideradas.
42
Salienta-se o facto de o operador N∗ não se encontrar definido no elemento, mas sim nos
pontos de colocação ou nas células plásticas. Apenas é necessário calcular este operador
para os modos de cedência activos. No caso do controlo da cedência ser efectuado com
base em células plásticas, efectua-se a transformação de coordenadas para o elemento
mestre definido no referencial local de plasticidade ( ),γ ψ .
O termo genérico associado às funções de aproximação definidas pelos índices ( ),i j , ao
termo k do vector das normais exteriores, ao elemento e e ao modo de cedência m, é
calculado através de:
( ) ( ) ( ) ( ) ( ) ( ) ( )1 1
,,
1 1
N , , , , , 'e m e e mi j ki j k P P n P J J d dγ ψ γ ψ γ ψ γ ψ ψ γ∗ ∗ ∗
− −
= ∫ ∫ . (5.67)
Mais uma vez, devido a não ser possível definir matematicamente a forma como varia na
célula plástica o vector das normais exteriores à superfície de cedência n∗ , as integrações
definidas em (5.67) foram efectuados numericamente, recorrendo-se para tal aos métodos
de quadratura de Gauss ou de Lobatto, através do somatório duplo:
( ) ( ) ( ) ( ) ( ) ( ),1 1
N , , , , , ' , 'n n
m e e mi j ki j k m P P n P W J Jγ ψ γ ψ γ ψ γ ψ∗ ∗ ∗⎡ ⎤= ⎣ ⎦∑∑ , (5.68)
onde n representa o número de pontos em cada direcção e 'W representa o peso associado
a cada ponto para a quadratura considerada.
No caso do controlo da cedência ser efectuado com base em pontos de colocação, a cada
modo de cedência está apenas associado a um ponto do domínio ( )0 0,γ ψ . Assim sendo,
um termo genérico associado às funções de aproximação definidas pelos índices ( ),i j , ao
termo k do vector das normais exteriores, e ao modo m, é calculado através de:
( ) ( ) ( ) ( ) ( )0 0 , 0 0N , , , , ,m e e mi j ki j k m P P n Pγ ψ γ ψ γ ψ∗ ∗ ∗= . (0.69)
4.4.6 Operador M∗
De acordo com a secção relativa ao modelo de elementos finitos, o operador M∗ é
calculado para o elemento finito e, através da expressão:
( )2
e ev v2
t eM S S dVφ ε∗∗ ∗
∂= Δ
∂σ∫ .
43
Este operador traduz o efeito das deformações plásticas que se desenvolvem durante o
processo iterativo no operador de flexibilidade generalizado. Caso não fosse tido em conta,
as relações constitutivas permaneceriam inalteradas, independentemente das deformações
plásticas que se fossem desenvolvendo durante todo o processo iterativo.
A sua implementação tem dois inconvenientes para o modelo numérico: em primeiro lugar
irá destruir a esparsidade do operador F , visto não ser possível estabelecer nenhuma
relação de ortogonalidade, conduzindo a uma matriz cheia. Em segundo lugar não é
possível definir expressões analíticas para o seu cálculo, sendo necessário recorrer a
procedimentos de integração numérica.
No entanto, o cálculo e a implementação deste operador em cada iteração são dispensáveis.
O modelo assim definido pode ser visto como a utilização do método de Newton-Raphson
modificado (Figura 4.10b), tendo a vantagem de ser mais eficiente do ponto de vista
numérico (maior esparsidade) e bastante menos exigente computacionalmente, mesmo
conduzindo a um maior número de iterações.
Figura 4.10: Representação esquemática do método de Newton-Raphson e do método de Newton-Raphson
modificado.
No programa Placa-Leg não foi implementado o operador M∗ . O desenvolvimento das
expressões para o cálculo deste operador podem ser encontradas em [Mendes 2002].
4.4.7 Vector das forças de massa
O vector das forças de massa associado ao elemento e é definido por:
( )tev vQ U b dV= ∫ ,
a) b)
44
e pode ser calculado recorrendo à transformação de coordenadas para o referencial local
definida pela expressão (5.70):
( ) ( )1 1
1 1
,te
v vQ U b J d dξ η η ξ− −
= ∫ ∫ . (5.70)
O termo genérico associado às funções de aproximação definidas pelos índices i, j e à
componente kb do vector das forças de massa, pode ser calculado através de:
( ) ( ) ( ) ( )1 1
1 1
, , ,ev i j kQ i j k P P b J d dξ η ξ η η ξ
− −
= ∫ ∫ . (5.71)
Particularizando para os elementos trapezoidais, um elemento genérico do vector das
forças de massa pode ser calculado por:
( ) ( ) ( ) ( )1 1
01 1
, ,v i j kQ i j k P P b J J J d dξ ηξ η ξ η η ξ− −
= + +∫ ∫ . (5.72)
A expressão anterior pode ser dividida em três parcelas, relacionadas com os termos que
definem o Jacobiano, obtendo-se:
Termo 1:
( ) ( ) ( ) ( ) ( ) ( )
( )
1 1 1 10 0
0 00 01 1 1 1
0
0
, 0
0, . .
k i j k i j
k
P PJ b P P d d J b P d P d
J b se i j
c c
ξ ηξ η η ξ ξ ξ η η
λ λ
λ
− − − −
=
⎧ = =⎪= ⎨⎪⎩
∫ ∫ ∫ ∫, (5.73)
Termo 2:
( ) ( ) ( ) ( ) ( ) ( )
( )( ) ( )
1 1 1 10 0
0 01 1 1 1
20
, 1 020, . .
k i j k i j
k
i
P PJ b P P d d J b P d P d
J b mse i j
c c
ξ ξ
ξ
ξ ηξ η ξ η ξ ξ ξ ξ η η
λ λ
λ λ
− − − −
=
⎧= ∧ =⎪= ⎨
⎪⎩
∫ ∫ ∫ ∫, (5.74)
Termo 3:
( ) ( ) ( ) ( ) ( ) ( )
( )( ) ( )
1 1 1 10 0
0 01 1 1 1
20
, 0 12
0, . .
k i j k i j
k
j
P PJ b P P d d J b P d P d
J b nse i j
c c
η η
η
ξ ηξ η η η ξ ξ ξ η η η
λ λ
λ λ
− − − −
=
⎧= ∧ =⎪= ⎨
⎪⎩
∫ ∫ ∫ ∫. (5.75)
45
A posição dos termos anteriores no vector das forças massa associado a um elemento,
encontra-se definido no Figura 4.11, para 3sg = .
i j0 00 10 20 31 01 11 2 - TERMO 11 3 - TERMO 22 0 - TERMO 32 12 22 33 03 13 23 3
Figura 4.11: Ordenamento dos termos do vector das forças de massa.
O vector terá um tamanho total de:
º
2
12 (1 )
n elemiuv
ig
=
⎡ ⎤+⎣ ⎦∑ . (5.76)
4.4.8 Vector das forças na fronteira
O vector das forças na fronteira associado à fronteira estática f é definido através de:
( )tfQ U t dγ γ γ σ= Γ∫ ,
e pode ser calculado recorrendo à transformação de coordenadas para o referencial local,
obtendo-se a expressão:
( ) ( )1
1
tfQ U t J dγ γ γ γ γ−
= ∫ . (5.77)
Atendendo à definição (5.50), o termo genérico associado à função de aproximação
definida pelo índice i e à componente ,ktγ do vector das forças na fronteira, pode ser
calculado através de:
( ) ( ) ( ) ( )1
2 2,
1
, , ,fi k x yQ m k P t n n dγ γξ ξ η ξ η γ
−
= +∫ . (5.78)
A consideração de cargas uniformes definidas nas fronteiras dos elementos, pode ser
implementada através da particularização da expressão (5.78) para cada um dos lados do
elemento mestre, obtendo-se:
46
Lado I - ( )γ ξ= e ( ) ( )2 2 1 1; ;x yn n β γ β γ= − − + :
( ) ( ) ( ) ( )
( ) ( )
12 2 0
, 2 2 1 101
2 2, 2 2 1 1
0
, 0
0, . .
k i
k
Pt P d
tse i
c c
γ
γ
ξβ γ β γ ξ ξ
λ
β γ β γλ
−
− + − + =
⎧ − + − +⎪ == ⎨⎪⎩
∫, (5.79)
Lado II - ( )γ η= , ( ) ( )2 2 1 1; ;x yn n α γ α γ= + − − :
( ) ( ) ( ) ( )
( ) ( )
12 2 0
, 2 2 1 101
2 2, 2 2 1 1
0
, 0
0, . .
k i
k
Pt P d
tse i
c c
γ
γ
ηα γ α γ η η
λ
α γ α γλ
−
+ + − − =
⎧ + + − −⎪ == ⎨⎪⎩
∫, (5.80)
Lado III - ( )γ ξ= , ( ) ( )2 2 1 1; ;x yn n β γ β γ= − − + :
( ) ( ) ( ) ( )
( ) ( )
12 2 0
, 2 2 1 101
2 2, 2 2 1 1
0
, 0
0, . .
k i
k
Pt P d
tse i
c c
γ
γ
ξβ γ β γ ξ ξ
λ
β γ β γλ
−
− + − + =
⎧ − + − +⎪ == ⎨⎪⎩
∫, (5.81)
Lado IV - ( )γ η= , ( ) ( )2 2 1 1; ;x yn n α γ α γ= − + − :
( ) ( ) ( ) ( )
( ) ( )
12 2 0
, 2 2 1 101
2 2, 2 2 1 1
0
, 0
0, . .
k i
k
Pt P d
tse i
c c
γ
γ
ξα γ α γ ξ ξ
λ
α γ α γλ
−
− + + − =
⎧ − + + −⎪ == ⎨⎪⎩
∫. (5.82)
Para cada fronteira estática, o termo genérico calculado pelas expressões anteriores, insere-
se na posição do vector das forças na fronteira apresentada na Figura 4.12.
m0 - LADOS I,II,III e IV12
Figura 4.12: Ordenamento dos termos do vector das forças das forças de fronteira.
47
O vector terá um tamanho total de:
º . .
1(1 )
n g liu
ig γ
=
⎡ ⎤+⎣ ⎦∑ . (5.83)
4.5 Cálculo dos campos de tensões e deslocamentos
Os campos de tensões σ , no elemento finito e, são calculados através da aproximação
(3.1). Assim sendo, para um ponto ( ),ξ η do domínio do elemento mestre, obtém-se:
( ) ( ) ( )0 0
( , ) , ( , )s sg g
ei j p
i j
P P X i jσ ξ η ξ η σ ξ η= =
⎡ ⎤= +⎣ ⎦∑∑ , (5.84)
onde ( ),X i j representa o peso associado às funções de aproximação do campo de
tensões no domínio do elemento e, de índices i e j, e sg representa o grau máximo das
funções utilizadas na definição da aproximação.
Os campos de deslocamentos vu no domínio do elemento finito e, são calculados através
da aproximação (3.2). Desta forma, para um ponto ( ),ξ η do elemento mestre, obtém-se:
( ) ( ) ( )0 0
( , ) , ( , )uv uvg g
ev i j v p
i j
u P P q i j uξ η ξ η ξ η= =
⎡ ⎤= +⎣ ⎦∑∑ , (5.85)
onde ( ),vq i j representa o peso associado às funções de aproximação do campo de
deslocamentos no domínio do elemento e, de índices i, j, e uvg representa o grau máximo
das funções.
Os campos de deslocamentos nas fronteiras estáticas uγ , são calculados através da
aproximação (3.3). Desta forma, para um ponto γ da fronteira f do elemento mestre,
obtém-se:
( ) ( )0
( )ug
fv i
iu P q i
γ
γγ γ=
⎡ ⎤= ⎣ ⎦∑ , (5.86)
onde ( )q iγ representa o peso associado à função de aproximação do campo de
deslocamentos na fronteira f, de índice i, e ug γ representa o grau máximo das funções.
4.6 Algoritmo do modelo elastoplástico
Nesta secção é abordado as questões relacionadas com o funcionamento do algoritmo
utilizado na resolução do problema não-linear. É discutida a forma adoptada para a
48
definição do incremento de carga, o processo de inicialização e actualização de variáveis, o
critério de convergência adoptado para o processo iterativo e por fim são apresentados os
diagramas de fluxo do algoritmo desenvolvido.
Conforme já foi referido, é conveniente não ter incrementos de carga muito grandes, pois
caso contrário os erros consequentes da utilização do método ímplicito de Euler poderão
influenciar a qualidade da solução. Outra restrição à dimensão do passo de carga é a
necessidade de se obter convergência ou que o número de iterações necessárias não seja
muito elevado.
No programa Placa-Leg não se implementou nenhum algoritmo de estimativa do erro, pelo
que a definição do incremento de carga foi feita através da análise do esforço de
convergência. Assim sendo, implementou-se um algoritmo que permite definir
automaticamente o incremento de carga, tendo em consideração o número de iterações
necessárias à convergência em passos anteriores. Caso o número de iterações ultrapasse o
desejado, é também possível definir regras para a subdivisão do actual incremento de
carga.
No início do passo de carga 1j + , as variáveis são inicializadas com os valores do final do
passo de carga anterior, à excepção dos incrementos dos parâmetros plásticos que são
inicializados com valor nulo, pois referem-se exclusivamente ao passo de carga em
questão:
01
0, 1 ,
0, 1 ,
0, 1 0
j j
v j v j
j j
j
X X
q q
q q
e
γ γ
+
+
+
∗ +
⎧ =⎪⎪ =⎪⎨
=⎪⎪⎪Δ =⎩
. (5.87)
No final da iteração itN a solução final pode ser obtida de:
( )
1
,1
,1
1
it
it
it
it
Ni i
ji
Ni iv v j v
i
Ni i
ji
Ni i
i
X X X
q q q
q q q
e e
γ γ γ
=
=
=
∗ ∗=
⎧= + Δ⎪
⎪⎪⎪ = + Δ⎪⎨⎪ = + Δ⎪⎪⎪Δ = Δ Δ⎪⎩
∑
∑
∑
∑
. (5.88)
49
No programa Placa-Leg optou-se por considerar que foi atingida a convergência quando o
maior valor absoluto dos restos das variáveis é inferior a uma dada tolerância _Tol R , tão
pequena quanto se queira:
{ }1 2 3 4, , , _i i i iR R R R Tol R≤ . (5.89)
Outra tolerância que é necessária definir, podendo condicionar fortemente a convergência
do processo iterativo, diz respeito ao valor do zero numérico na determinação dos modos
de cedência que se encontram activos e inactivos. Neste trabalho utilizou-se a seguinte
expressão para determinar se o modo de cedência m se encontra inactivo:
_ , _ 0m Tol com TolΦ ≤ − Φ Φ ≥ , (5.90)
todos os modos que não verifiquem esta condição, são considerados activos.
Figura 4.14 procura-se ilustrar através de diagramas de fluxo, o algoritmo utilizado para a
implementação do modelo de elementos finitos. A primeira destas figuras refere-se à
implementação dos incrementos de carga, enquanto que a segunda se refere
exclusivamente ao algoritmo iterativo implementado.
50
Figura 4.13: Algoritmo do incremento de carga.
PROCESSO ITERATIVO
INICIALIZAÇÃO DE VARIÁVEIS
INCREMENTO DE CARGA
{ }0 ,0 ,0 0v pX q q eγ= = = =
ACTUALIZAÇÃO DAS DEFORMAÇÕES PLÁSTICAS
, 1 ,p j p je e N e+ ∗ ∗= + Δ
NÃO
FIM
SIM
1jλ λ +=
FIM DO CARREGAMENTO ?
( )1 , 1 , 1, , ,j v j jX q q eγ+ + + ∗Δ
51
Figura 4.14: Algoritmo do processo iterativo.
ITERAÇÕES
FIM
SIM
NÃO
SOLUÇÃO SATISFAZ O CRITÉRIO
DE CONVERGÊNCIA ?
MONTAGEM DO SISTEMA DE EQUAÇÕES
INICIALIZAÇÃO DE VARIÁVEIS
CÁLCULO DOS RESTOS
{ } { }i R inH Var =Δ −⎡ ⎤⎣ ⎦
{ }
0 0 01 , 1 , , j+1 , j
0,
, ,0 ;
j j v j v j
j
X X q q q qe N N
γ γ+ +
∗ ∗ ∗
= = =
Δ = =
( )1 2 3 4, , ,i i i iR R R R
CÁLCULO DO CAMPO DE ESFORÇOS
CÁLCULO DOS POTENCIAIS PLÁSTICOS
,i t iP dVφ∗ ∗ ∗=Φ ∫
CÁLCULO DOS OPERADORES
i t ivN S n P dV∗ ∗ ∗= ∫
MÉTODO DE NEWTON-RAPHSON
1. Resolução do sistema de equações.
2. Actualização do vector solução.
i iS Xv pσ σ= +
52
5 PROGRAMA DE CÁLCULO AUTOMÁTICO
5.1 Introdução
Nesta secção descreve-se o funcionamento do programa de cálculo Placa-Leg. O programa
foi desenvolvido na plataforma Matlab [MathWorks Inc. 1999] e é constituído pelas
rotinas listadas na Tabela 5.1 e pelos ficheiros auxiliares indicados na Tabela 5.2.
Para visualizar os resultados recorre-se a um conjunto de rotinas de visualização
denominadas por ToolVis, igualmente escritas na plataforma Matlab. Estas rotinas
recorrem a ficheiros de configuração (extensão TDF) para efectuar as visualizações (ver
Tabela 5.3).
O programa foi desenvolvido na versão 6 (R12) do Matlab [MathWorks Inc. 1999] tendo
sido testado com êxito nas versões 7.0 (R13) e 7.1 (R14). A utilização de versões
anteriores ou posteriores poderá implicar a necessidade de se efectuarem ajustes no código
que se inclui no CD em anexo (versão de Maio de 2006). Contactando os autores através
de correio electrónico ([email protected] ou [email protected]) é possível obter o URL
onde se encontra a última versão do programa disponível para descarregar.
Rotina Objectivo Tipo \PLACA.m Arranque do programa Executável \VARIAVEIS.m Criação de variáveis globais Executável \CFG\PLACA_CFG.m Rotina de configuração Executável \EST\PLACA_EST.m Visualização do problema Executável \PRE\PLACA_PRE.m Rotina de pré-processamento Executável \PRE\Calc_GLOB_LOC.m Mapeamento do elemento mestre Função \PRO\PLACA_PRO.m Rotina de processamento Executável \PRO\Calc_F.m Cálculo do operador F Função \PRO\Calc_Av.m Cálculo do operador Av Função \PRO\Calc_Ag.m Cálculo do operador Ag Função \PRO\Calc_Qv.m Cálculo do vector Qv Função \PRO\Calc_Qg.m Cálculo do vector Qg Função \PRO\Monta_sist_e.m Montagem da matriz do sistema governativo Função \PRO\Monta_TI_e.m Montagem dos termos independentes do s.g. Função \PRO\Res_SIST.m Resolução do sistema governativo Executável \PRO\Gere_ANL.m Gere algoritmo não linear Função \PRO\Calc_Iteracoes.m Cálculo das iterações Função \PRO\Efectua_REGISTO.m Efectua registo de variáveis Função \POS\PLACA_POS.m Rotina de pós-processamento Executável \POS\Calc_DESLOCAMENTOS_uv.m Cálculo de uv Função \POS\Calc_DESLOCAMENTOS_ug.m Cálculo de ug Função \POS\Calc_ENERGIA_DEFORMACAO.m Cálculo da energia deform. Executável \TOOLVIS.V2.0\TOOLVIS.m Rotina toolvis Função \TOOLVIS.V2.0\CORES.m Def. gradiente de cores Executável
Tabela 5.1: Listagem das rotinas de cálculo do programa.
53
Ficheiro Objectivo \PLACA.fig Interface gráfica (menu principal) \CONFIG.OPT Configuração do programa (ficheiros e directórios) \CFG\PLACA_CFG.fig Interface gráfica (configuração) \EST\PLACA_EST.fig Interface gráfica (estrutura) \EST\PLACA_EST.OPT Configuração do programa (estrutura) \PRE\PLACA_PRE.fig Interface gráfica (pré-processamento) \PRO\PLACA_PRO.fig Interface gráfica (processamento) \PRO\PLACA_PRO.OPT Configuração do programa (processamento) \PRO\Final.wav Ficheiro wav para assinalar o fim dos cálculos \POS\PLACA_POS.fig Interface gráfica (pós-processamento) \POS\PLACA_POS.OPT Configuração do programa (pós-processamento) \TOOLVIS.V2.0\CORES.fig Interface gráfica (gradiente de cores) \TOOLVIS.V2.0\CM\*.CM Definição dos mapas de cores
Tabela 5.2: Listagem dos ficheiros adicionais.
Ficheiro Objectivo ESTRUTURA.TDF Visualização do problema em análise(Interface gráf.) PLACA_VIS.TDF Visualização do problema em análise (executável) CELULAS_CEL.TDF Visualização das células plásticas CELULAS_PC.TDF Visualização dos pontos de colocação CELULAS_PC_CEL.TDF - CELULAS_PC_PI.TDF - CELULAS_PC_PI_CEL.TDF - CELULAS_PI.TDF Visualização dos pontos de integração CELULAS_PI_CEL.TDF - CORTE_s.TDF Visualização das tensões segundo um corte CORTE_uv.TDF Visualização dos deslocamentos uv segundo um corte CORTE_ug.TDF Visualização dos deslocamentos ug segundo um corte DEMO.TDF Ficheiro de exemplo DESLOCAMENTOS.TDF Visualização dos deslocamentos (standard) DESLOCAMENTOS_uv.TDF Visualização dos deslocamentos uv DESLOCAMENTOS_uv_2fig.TDF Visualização dos deslocamentos uv (2 figuras) DESLOCAMENTOS_ug.TDF Visualização dos deslocamentos ug GLOBAL-ELÁSTICO.TDF Visualização standard - problema elástico GLOBAL-ELASTOPLASTICO.TDF Visualização standard - problema elastoplástico GLOBAL-ELASTOPLASTICO_SEM_TENSOES.TDF
Visualização standard - problema elastoplástico (sem tensões)
REGISTOS_COM_MARKERS.TDF Visualização dos registos com markers REGISTOS_SEM_MARKERS.TDF Visualização dos registos sem markers TENSOES.TDF Visualização das tensões TENSOES_CELULAS.TDF Visualização das tensões nas células plásticas TENSOES_IN_CORES.TDF Visualização das tensões – entrada dos mapas de cores TENSOES_OUT_CORES.TDF Visualização das tensões – saída dos mapas de cores TENSOES_sxx.TDF Visualização das tensões, σxx TENSOES_syy.TDF Visualização das tensões, σyy TENSOES_sxy.TDF Visualização das tensões, σxy
Tabela 5.3: Listagem dos ficheiros de configuração das saídas gráficas.
5.2 Estrutura do ficheiro de dados
O ficheiro de dados é um ficheiro ASCII e recorre a blocos identificadores para estruturar
os dados, permitindo uma maior flexibilidade na ordenação da informação. Na Tabela 5.4
encontra-se exemplificada a estrutura do ficheiro de dados.
54
5.3 Funcionamento do programa
O programa Placa-Leg pode funcionar através da utilização da interface gráfica
apresentada Figura 5.1. No entanto, alguns comandos podem ser executados directamente
na command window do Matlab.
Figura 5.1: Menu principal do programa.
5.3.1 Configuração
Para funcionar, o programa necessita das rotinas indicadas na Tabela 5.1, dos ficheiros
adicionais indicados na Tabela 5.2, da configuração das visualizações indicada na Tabela
5.3, e ainda das rotinas de visualização Toolvis (versão 2.0 de Janeiro de 2006). Estas
últimas rotinas podem armazenar-se num directório à parte (recomendado), assim como os
ficheiros de dados.
Na primeira execução do programa deverá ser efectuada a configuração, através do
seguinte procedimento:
1. Direccionar a current directory do Matlab para o directório principal onde o
programa se encontra instalado;
55
[Titulo] CONSOLA QUADRADA [Nos] 4 0 0 1 0 1 1 0 1 [Lados] 4 1 2 0 0 2 2 20 2 3 0 0 2 2 20 3 4 0 0 2 2 20 4 1 -1 -1 2 2 20 [Elementos] 1 1 2 3 4 1 0.3 1 2 20 20 [Aproximacoes] 2 [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16] [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15] [Fmassa] 0 [Ffront] 1 3 2 -1 [tipo] 1 [Macro_Celulas] 1 0 0 1 0 1 1 0 1 15 2 2 1 2 0 1 [Registos] 13 1 1 0 0.5 2 1 0 0.5 3 1 0 0.5 4 1 0 1 4 1 0.25 0.75 5 1 0 1 5 1 0.25 0.75 7 1 1 0.5 8 1 1 0.5 9 1 1 0.5 13 1
Marca identificadora do bloco de dados “Titulo” Título do problema Marca identificadora do bloco de dados “Nos” Número de nós Coordenada x; Coordenada y do nó 1 Marca identificadora do bloco de dados “Lados” Número de lados Nó inicial; nó final; cond. apoio x, cond. apoio y, aprox. ug(x); aprox. ug(y), número de intervalos para pós-processamento Nota: Condição de apoio: 0–livre; -1–fixo. Não é possível
considerar apoios elásticos. Marca identificadora do bloco de dados “Elementos” Número de elementos Lado 1, lado 2, lado 3, lado 4, E, Poisson, aprox. σ, aprox. uv, número de intervalos σ; número de intervalos uv Marca identificadora do bloco de dados “Aproximacoes” Número de aproximações Graus da aproximação n Marca identificadora do bloco de dados “Fmassa” Número de forças de massa Elemento, direcção, valor Marca identificadora do bloco de dados “Ffront” Número de forças na fronteira Fronteira, tipo, valor Marca identificadora do bloco de dados “tipo” Tipo Nota: 1-EPT; 2-EPD Marca identificadora do bloco de dados “Macro_Celulas” Número de macro-células Marca identificadora do bloco de dados “Registos” Número de registos Tipo elemento posição_x posição_y Nota: Tipos: 1 – sxx; 2 – syy; 3 – sxy ; 4 – scomp (mcel); 5 –
fi(mcel) ; 6 – fi(modo) ; 7 – uv_x ; 8 – uv_y ; 9 – |uv| ; 10 – ug_x ; 11 – ug_y; 12 – |ug| ; 13 – U;
Nota: O tipo 13, refere-se ao registo da energia de deformação e não é necessário definir o elemento e a posição, uma vez que se aplica à estrutura na sua globalidade.
Número de passos entre registos
Tabela 5.4: Formato do ficheiro de dados.
56
2. Digitar na command window o comando “PLACA”, o que fará aparecer o menu
principal (Figura 5.1);
3. Clicar no botão “Configuração”, que fará aparecer a janela de configuração (Figura
5.2);
4. Clicar em “Procurar o directório principal” que deverá fazer aparecer na caixa de
texto o directório principal onde o programa se encontra instalado;
5. Clicar em “Guardar configuração” para gravar as novas definições no ficheiro
“CONFIG.OPT”;
6. Clicar em “Actualizar o Matlab Path”, para adicionar os directórios ao caminho de
busca do programa.
Os directórios adicionados por defeito ao caminho de busca, são os seguintes: “.\”;
“.\CFG”; “.\EST”; “.\PRE”; “.\POS”; “.\PRO”; “.\EXEMPLOS”; “.\TOOLVIS.v2.0”;
“.\TOOLVIS.v2.0\CM'” e “.\TOOLVIS.v2.0\TDF”, onde “.\” se refere ao directório
principal onde o programa se encontra instalado. Alterações podem ser efectuadas na sub-
rotina “BT_ACTUALIZA_MATLAB_PATH_Callback” do ficheiro “PLACA_CFG.m”.
Figura 5.2: Interface do módulo de configuração.
5.3.2 Executar o programa
Para efectuar os cálculos e visualizar os resultados da análise, deve-se seguir o seguinte
procedimento:
1. Pré-processamento;
2. Visualização do problema em análise (opcional);
3. Processamento;
4. Pós-processamento.
57
5.3.2.1 Pré-processamento
O pré-processamento pode ser efectuado acedendo ao respectivo módulo (Figura 5.3),
clicando em “Pré-Processamento” no menu principal ou digitando “PLACA_PRE”.
Posteriormente deve-se seguir o seguinte procedimento:
1. Clicar em “Clear all” e em “Variáveis”, para limpar o Workspace e para visualizar
as variáveis globais, respectivamente;
2. No grupo de opções “Ficheiros e directórios”, definir o “Directório de trabalho” e o
“Ficheiro de dados” (Nota: Clicando na lista de ficheiros à direita é possível buscar
o directório e o ficheiro pretendido);
3. Clicar em “Ler” para carregar o ficheiro de dados;
4. Escolher o tipo de cálculo como “Elástico” caso se pretenda efectuar o cálculo
elástico e “Elastoplástico” caso se pretenda efectuar o cálculo elastoplástico. Em
seguida clicar em “Executar”, o que iniciará os cálculos de pré-processamento;
5. Clicar em “Menu Principal”.
Figura 5.3: Interface do módulo de pré-processamento.
58
Eventuais alterações aos dados podem ser feitas utilizando os restantes grupos de opções
(“Gerais”; “Nós”; “Lado”, etc). Para que as alterações sejam gravadas no ficheiro de dados
basta clicar em “Gravar”. Estes menus podem igualmente ser utilizados para criar um novo
ficheiro de dados.
5.3.2.2 Visualizar o problema em análise
Depois de ter sido efectuado o pré-processamento, clicando em “Visualizar a estrutura” no
menu principal ou digitando “PLACA_EST”, é possível visualizar uma representação
gráfica do problema carregado em memória (Figura 5.4). Através dos menus, é possível
activar ou desactivar várias opções de visualização.
Figura 5.4: Interface do módulo de visualização do problema em análise.
5.3.2.3 Processamento
Depois de ter sido efectuado o pré-processamento, o processamento é efectuado acedendo
ao respectivo módulo (Figura 5.5) acedido através do menu principal, ou digitando
“PLACA_PRO” na command windows. Posteriormente, deve-se seguir o seguinte
procedimento:
59
1. Activar a opção “Calcular operadores elásticos” e “Calcular a energia de
deformação” (opcional) no grupo de opções “Gerais”;
2. Activar a opção “Elástico” ou “Elastoplástico” em “Tipo de cálculo”;
3. No caso de uma análise elastoplástica, definir os incrementos de carga a serem
efectuados na análise incremental, utilizando o botão “Adicionar” e uma das três
opções disponíveis;
4. Clicar em “Executar”, que iniciará os cálculos de processamento;
5. Clicar em “Menu Principal”.
Desactivando a opção “Calcular operadores elásticos”, é possível evitar recalcular os
operadores elásticos ( , ,vF A Aγ ), escolha especialmente útil em situações onde foi alterado
apenas o carregamento.
Figura 5.5: Interface do módulo de processamento.
60
O programa permite continuar um cálculo anteriormente efectuado através da opção
“Continuar” em “Sequência”. Caso o problema não esteja carregado em memória, através
do botão “Ler Passo”, lendo os ficheiros “Dados_Pre.mat”, “Elast.mat” e um dado passo
de carga “P_x.mat” é possível recuperar todas as variáveis e continuar o carregamento. A
opção “Reiniciar” em “Sequência” é utilizada para iniciar um novo cálculo.
Como o cálculo é efectuado incrementalmente, é necessário estipular os incrementos de
carga a considerar, utilizando o botão “Adicionar” e as três opções para definir os valores
numéricos dos incrementos. É também possível efectuar registos de diversas grandezas
estruturais (e.g. tensões, deslocamentos, etc.). Para tal, deverão ser utilizadas as opções
indicadas em “Registo” e que podem ficar guardadas nos ficheiros de dados gravando-o
posteriormente na interface de pré-processamento.
O programa inclui um procedimento de ajuste dos incrementos de carga através da análise
do esforço de convergência, que permite também definir um número máximo de iterações
a partir do qual o incremento de carga é divido em incrementos mais pequenos.
Outros parâmetros que influenciam fortemente a evolução e até a convergência do cálculo
são as tolerâncias. Na interface gráfica é possível definir as seguintes tolerâncias: Tol_zero
– representa o valor do zero numérico; Tol_zero_sol – representa o valor do zero numérico
para o vector solução (não utilizado); Tol_zero_iter – representa o valor da tolerância do
resíduo máximo para definir a convergência do processo iterativo, calculado através de
(5.91); Tol_zero_fi – representa o zero numérico na determinação dos modos de cedência
que se encontram activos e inactivos, efectuada recorrendo à expressão (5.92):
{ }1 2 3 4, , , _i i i iR R R R Tol R≤ , (5.91)
_ , _m Tol com TolΦ Φ 0Φ ≤ − ≥ . (5.92)
Recomenda-se a utilização dos seguintes valores: Tol_zero=1e-15; Tol_zero_sol=0;
Tol_zero_iter=1e-12; Tol_zero_fi=1e-13.
5.3.2.4 Pós-processamento
Depois do processamento, o pós-processamento é efectuado acedendo ao respectivo
módulo (Figura 5.6), através do menu principal ou digitando “PLACA_POS”.
Posteriormente, deve seguir-se o seguinte procedimento:
1. Escolher o tipo de visualização que se pretende efectuar (“Global”; “Tensões”;
“Deslocamentos”, etc.).
61
2. Definir as opções de visualização do grupo escolhido. A cada opção encontra-se
associado um ficheiro TDF.
3. Activar a opção “Calcular”
4. Activar a opção “Visualizar”
5. Clicar em “Executar”, que fará aparecer uma janela com a visualização.
É possível limitar os elementos e lados representados nas visualizações através dos grupos
de opções no topo da janela.
É necessário efectuar o cálculo das variáveis de pós-processamento apenas uma vez, razão
pela qual se deve desactivar a opção “Calcular” nas visualizações posteriores.
Digitar o comando “TOOLVIS GLOBAL-ELÁSTICO.TDF” na command window, é o
mesmo que executar a visualização “Global”, “Global elástico” com a opção calcular
desactivada.
No grupo de opções “Opções” é possível definir, entre outros parâmetros, o factor de
escala das deformadas.
Figura 5.6: Interface do módulo de pós-processamento.
62
6 EXEMPLOS DE APLICAÇÃO
Nesta secção apresentam-se dois exemplos de aplicação. No primeiro caso analisa-se uma
placa quadrada em consola, recorrendo a uma malha de 1 elemento (ver Figura 6.1) e a
células plásticas do tipo CEL_1. O ficheiro de dados encontra-se representado na Tabela
6.1: Exemplo 1 – Ficheiro de dados.Tabela 6.1 e os resultados obtidos nas Figuras 6.4 a
6.7. O segundo exemplo diz respeito à análise de uma viga bi-encastrada, discretizada em 5
elementos finitos (ver Figura 6.2). O ficheiro de dados encontra-se representado nas
Tabelas 6.2 e 6.3 e os resultados obtidos nas Figuras 7.9 a 7.14.
A análise destes e outros exemplos pode ser consultada em [Mendes 2002].
( )1.0
sv
uv
uγ
ced
EPT
a = 1.00
E = 1.00
ν= 0.30
g =16
g =15
g =15
von Mises σ =
Figura 6.1: Exemplo 1 – Definição das características da placa quadrada em consola.
( )1.0sv uv uγ cedEPT; E = 1.00; ν= 0.30; g =16; g =15; g = 15; von Mises σ =
Figura 6.2: Exemplo 2 – Definição das características da viga bi-encastrada.
63
[Titulo] CONSOLA QUADRADA [Nos] 4 0 0 1 0 1 1 0 1 [Lados] 4 1 2 0 0 2 2 20 2 3 0 0 2 2 20 3 4 0 0 2 2 20 4 1 -1 -1 2 2 20 [Elementos] 1 1 2 3 4 1 0.3 1 2 20 20 [Aproximacoes] 2 [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16] [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15] [Fmassa] 0 [Ffront] 1 3 2 -1 [tipo] 1 [Macro_Celulas] 1 0 0 1 0 1 1 0 1 15 2 2 1 2 0 1 [Registos] 23 1 1 0.5 0.5 2 1 0.5 0.5 3 1 0.5 0.5 1 1 0 0.5 2 1 0 0.5 3 1 0 0.5 4 1 0 1 4 1 0.25 0.75 5 1 0 1 5 1 0.25 0.75 7 1 1 0.5 8 1 1 0.5 9 1 1 0.5 7 1 0.25 0.75 8 1 0.25 0.75 9 1 0.25 0.75 7 1 0 1 8 1 0 1 9 1 0 1 10 2 1 0.5 11 2 1 0.5 12 2 1 0.5 13 1
Tabela 6.1: Exemplo 1 – Ficheiro de dados.
64
Figura 6.3: Exemplo 1 – Visualização global-elastoplástico dos resultados, =0.215λ .
65
Figura 6.4: Exemplo 1 – Visualização global-elastoplástico dos resultados, =0.45186λ .
66
Figura 6.5: Exemplo 1 – Visualização deslocamentos (uv/ug), =0.45186λ .
Figura 6.6: Exemplo 1 – Visualização das células plásticas e dos pontos de integração, =0.45186λ .
67
Figura 6.7: Exemplo 1 – Visualização dos campos de tensões, ao longo do corte GHI, =0.45186λ .
Figura 6.8: Exemplo 1 – Visualização dos campos de deslocamentos uv, ao longo do corte GHI, =0.45186λ .
Figura 6.9: Exemplo 1 – Visualização dos campos de deslocamentos uγ, ao longo do corte GHI, =0.45186λ .
Salienta-se que as abcissas dos gráficos da Figura 6.9 se encontram definidas segundo a
orientação dada na definição da fronteira, que neste caso corresponde a (-x).
68
[Titulo] VIGA ENCASTRADA [Nos] 12 0 0 3 0 6 0 7 0 10 0 13 0 0 1 3 1 6 1 7 1 10 1 13 1 [Lados] 16 1 2 0 0 2 2 20 2 3 0 0 2 2 20 3 4 0 0 2 2 20 4 5 0 0 2 2 20 5 6 0 0 2 2 20 7 8 0 0 2 2 20 8 9 0 0 2 2 20 9 10 0 0 2 2 20 10 11 0 0 2 2 20 11 12 0 0 2 2 20 1 7 -1 -1 2 2 20 2 8 0 0 2 2 20 3 9 0 0 2 2 20 4 10 0 0 2 2 20 5 11 0 0 2 2 20 6 12 -1 -1 2 2 20 [Elementos] 5 1 12 6 11 1 0.3 1 2 20 20 2 13 7 12 1 0.3 1 2 20 20 3 14 8 13 1 0.3 1 2 20 20 4 15 9 14 1 0.3 1 2 20 20 5 16 10 15 1 0.3 1 2 20 20 [Aproximacoes] 2 [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14] [0 1 2 3 4 5 6 7 8 9 10 11 12 13] [Fmassa] 0 [Ffront] 1 8 2 -1 [tipo] 1
Tabela 6.2: Exemplo 2 – Ficheiro de dados.
69
[Macro_Celulas] 13 0 0 1 0 1 1 0 1 8 2 2 1 3 0 1 1 0 2 0 2 1 1 1 6 2 2 1 3 0 1 2 0 3 0 3 1 2 1 3 2 2 1 3 0 1 3 0 4 0 4 1 3 1 3 2 2 1 3 0 1 4 0 5 0 5 1 4 1 3 2 2 1 3 0 1 5 0 6 0 6 1 5 1 6 2 2 1 3 0 1 6 0 7 0 7 1 6 1 8 2 2 1 3 0 1 7 0 8 0 8 1 7 1 6 2 2 1 3 0 1 8 0 9 0 9 1 8 1 3 2 2 1 3 0 1 9 0 10 0 10 1 9 1 3 2 2 1 3 0 1 10 0 11 0 11 1 10 1 3 2 2 1 3 0 1 11 0 12 0 12 1 11 1 6 2 2 1 3 0 1 12 0 13 0 13 1 12 1 8 2 2 1 3 0 1 [Registos] 41 1 1 0 1 2 1 0 1 3 1 0 1 4 1 0 1 5 1 0 1 7 1 0 1 8 1 0 1 9 1 0 1 1 1 0 0.5 2 1 0 0.5 3 1 0 0.5 4 1 0 0.5 5 1 0 0.5 7 1 0 0.5 8 1 0 0.5 9 1 0 0.5 1 3 6.5 0.5 2 3 6.5 0.5 3 3 6.5 0.5 4 7 6.5 0.5 5 7 6.5 0.5 7 3 6.5 0.5 8 3 6.5 0.5 9 3 6.5 0.5 1 3 6.5 1 2 3 6.5 1 3 3 6.5 1 4 7 6.5 1 5 7 6.5 1 7 3 6.5 1 8 3 6.5 1 9 3 6.5 1 1 3 6.5 0 2 3 6.5 0 3 3 6.5 0 4 7 6.5 0 5 7 6.5 0 7 3 6.5 0 8 3 6.5 0 9 3 6.5 0 13 1
Tabela 6.3: Exemplo 2 – Ficheiro de dados (cont.).
70
Figura 6.10: Exemplo 2 – Visualizações de diversas grandezas estruturais, =0.12λ .
71
Figura 6.11: Exemplo 2 – Visualização dos campos de tensões e deslocamentos ao longo do corte CD,
=0.12λ .
72
Figura 6.12: Exemplo 2 – Visualizações de diversas grandezas estruturais, =0.17567λ .
73
Figura 6.13: Exemplo 2 – Visualização dos campos de tensões e deslocamentos ao longo do corte CD,
=0.17567λ .
75
7 REFERÊNCIAS BIBLIOGRÁFICAS
Almeida, J. P. B. M. (1991) - "Modelos de Elementos Finitos para a Análise Elastoplástica", Doutoramento em Engenharia Civil, Universidade Técnica de Lisboa, Instituto Superior Técnico, Lisboa.
Arantes e Oliveira, E. R. (1969) - "Resistência dos Materiais - Livro 2: Elementos da Teoria da Elasticidade", AEIST, Lisboa.
Bussamra, F. (1999) - "Elementos Finitos Híbridos-Trefftz: Um Modelo Elastoplástico Tridimensional ", Doutoramento em Engenharia, Escola Politécnica da Universidade de São Paulo, São Paulo.
Castro, L. M. S. S. (1996) - "Wavelets e Séries de Walsh em Elementos Finitos", Doutoramento em Engenharia Civil, Universidade Técnica de Lisboa, Instituto Superior Técnico, Lisboa.
Chen, W. (1994) - "Constitutive Equations for Engineering Materials, Volume 2: Plasticity and Modeling", Elsevier,
Coulomb, C. A. (1873) - "Sur une application des regles de maximis et minimis a quelques problemes de statique relatifs a l'architecture", Acad. R. Sci. Mem. Math. Phys. par divers svants, Vol. 7, pp. 343-382.
Crisfield, M. A. (1991) - "Non-Linear Finite Element Analysis of Solids and Structures, Volume 1", John Wiley & Sons Ltd, ----.
Drucker, D. C. e W. Prager (1952) - "Soil Mechanics and Plastic Analysis of Limit Design", Quarterly of Applied Mathematics, Vol. 10, pp. 157-165.
Freitas, J. A. T., J. P. B. M. Almeida, et al. (1999) - "Non-conventional Formulations for the Finite Element Method", Computational Mechanics - Springer, Volume 23(Number 5-6), pp. 488-501.
Kachanov, L. M. (1974) - "Fundamentals of Theory of Plasticity", Mir Publishers, Moscow.
Lode, W. (1926) - "Versuche ueber den Einfluss der mitt leren Hauptspannung auí das Fliessen der Metalle Eisen Kupfer und Nickel", Zeitschrift fuer Physik, 36, 913-939.
MathWorks Inc., T. (1999) - "MATLAB - The Language Of Technical Computing.
Mendes, L. M. (2002) - "Modelos de Elementos Finitos Híbridos-Mistos de Tensão na Análise Elastoplástica de Estruturas Laminares Planas", Mestrado em Engenharia de Estruturas, Instituto Superior Técnico, Universidade Técnica de Lisboa, Lisboa.
Pereira, E. M. B. R. (1993) - "Elementos Finitos de Tensão, Aplicação à Análise Elástica de Estruturas", Doutoramento em Engenharia Civil, Universidade Técnica de Lisboa, Instituto Superior Técnico, Lisboa.
Pereira, E. M. B. R. e J. A. T. Freitas (2000) - "Numerical Implementation of a Hybrid-Mixed Finite Element Model for Reissner-Mindlin Plates", Computers & Structures – Elsevier Science Ltd(Number 74), pp. 323-334.
Reddy, J. N. (1993) - "An Introduction to The Finite Element Method", McGraw-Hill Inc., New York.
76
Spiegel , R. M. e L. Abellanas (1990) - "Fórmulas e Tabelas de Matemática Aplicada", McGraw-Hill,
Timoshenko, S. P. e S. Woinowsky-Kriefer (1982) - "Theory of Elasticity ", McGraw-Hill, Tóquio.
Tresca, H. (1864) - "Sur l'ecoulement des corps solides soumis a de fortes pression", Compt. Rend, Vol. 59, p.754.
Von Mises, R. (1913) - "Mechanik der festen Koerper in Plastisch deformabelm Zustand", Geottinger Nachr. Math. Phys., pp. 582-592.
Zienkiewicz, O. C. e R. L. Taylor (1991) - "The Finite Element Method - Volume 1 - Basic Formulation and Linear Problems", McGraw-Hill Inc., Londres.
Zienkiewicz, O. C. e R. L. Taylor (1991) - "The Finite Element Method - Volume 2 - Solid and Fluid Mechanics, Dynamics and Non-Linearity", London.
77
ANEXO A – POLINÓMIOS DE LEGENDRE
A.1 Introdução
As funções de Legendre de ordem (n), correspondem às soluções da equação diferencial de
Legendre:
2(1 ) 2 ( 1) 0x y x y n n y′′ ′− − + + = . (A.1)
Somos desta forma conduzidos à definição de funções do tipo polinomial, que podem ser
geradas pela fórmula de Rodriguez [Spiegel et al. 1990]:
( )21( ) 12 !
n n
n n n
dP x xn dx
= − , (A.2)
onde Pn(x) representa o polinómio de Legendre de ordem n.
A.2 Propriedades dos polinómios de Legendre
Verifica-se que polinómios de grau consecutivo são alternadamente funções pares e
ímpares.
( ) ( 1) ( )nn nP x P x− = − . (A.3)
Os integrais definidos no intervalo [-1,1] do produto de dois polinómios de Legendre,
podem ser calculados a partir das igualdades:
1
11
1
( ) ( ) 0 ,
2( ) ( ) ,2 1
m n
m n
P x P x dx para m n
P x P x dx para m nn
−
−
⎧= ≠⎪
⎪⎨⎪ = =⎪ +⎩
∫
∫. (A.4)
Assim sendo, conclui-se que os polinómios de Legendre são ortogonais no intervalo [-1,1] .
A.3 Fórmulas geradoras de polinómios de Legendre
Várias fórmulas podem ser utilizadas para gerar os polinómios de Legendre,
correspondendo a maioria a fórmulas de recorrência [Spiegel et al. 1990]:
( ) ( ) ( ) ( )1 1n n nP x x P x n P x+′ ′− = + , (A.5)
( ) ( ) ( )1n n nx P x P x n P x−′ ′− = , (A.6)
( ) ( ) ( ) ( )1 1 2 1n n nP x P x n P x+ −′ ′− = + , (A.7)
78
( ) ( ) ( ) ( )211 n n nx P x n x P x n P x−′− = − . (A.8)
Utilizou-se neste trabalho a fórmula recursiva de Bonnet:
1 1( 1) ( ) (2 1) ( ) ( ) 0n n nn P x n x P x n P x+ −+ − + + = , (A.9)
com ( )0 1P x = e ( )1P x x= .
Por forma a que se verifique a condição (A.10), a fórmula de Bonnet é escalada por nλ
(A.11), obtendo-se o formato alternativo apresentado em (A.12):
[ ]1
2
1
( ) 1nP x dx−
=∫ , (A.10)
2 12n
nλ += , (A.11)
1 11 1
1 2 1( ) ( ) ( )n n nn n n
n n nP x x P x P xλ λ λ+ −
+ −
+ += − , (A.12)
com ( )0 0P x λ= e ( )1 1P x xλ= :
Na Figura A-1 apresenta-se a representação gráfica dos polinómios de Legendre
unidimensionais de grau igual ou inferior a 14, e nas Figuras A-2 e A-3, os polinómios de
Legendre bidimensionais até ao terceiro grau, gerados através da expressão (A.12).
79
Figura A.1: Representação gráfica dos polinómios de Legendre unidimensionais.
80
Figura A.2: Representação gráfica dos polinómios de Legendre bidimensionais.
81
Figura A.3: Representação gráfica dos polinómios de Legendre bidimensionais (cont.).
82
A.4 Expressões para as integrações analíticas
Uma das principais vantagens associadas à utilização destes polinómios consiste em ser
possível dispensar os processos numéricos no cálculo das integrações, pois é sempre
possível determinar expressões analíticas para os operadores estruturais (em problemas
onde se assumem as hipóteses de linearidade física e geométrica), mesmo para os casos em
que não se consegue tirar directamente partido da ortogonalidade.
Apresentam-se de seguida as expressões indicadas por Pereira et al. [Pereira et al. 2000]:
1
10ij 1
1
( ) ( ) 1 ,
( ) ( ) 0 , . .
i j
i j
P x P x dx se i j
P x P x dx c c
−
−
⎧= =⎪
⎪= ⎨⎪ =⎪⎩
∫
∫A ; (A.13)
1
1
11ij
1
1
1
( ) ( ) , 12
( ) ( ) , 12
( ) ( ) 0 , . .
i ji j
i ji j
i j
iP x P x x dx se i j
jP x P x x dx se i j
P x P x x dx c c
λ λ
λ λ
−
−
−
⎧= = +⎪
⎪⎪⎪= = = −⎨⎪⎪⎪ =⎪⎩
∫
∫
∫
A ; (A.14)
( )( )( )
( )( )
( )
1 3 22
21
12
12ij 1
2
1
12
1
4 6 1( ) ( ) ,2 2 1 2 3
1( ) ( ) , 2
2 2 1
1( ) ( ) , 2
2 2 1
( ) ( ) 0 , . .
i ji
i ji j
i ji j
i j
i iP x P x x dx se i ji i
j jP x P x x dx se i j
j
i iP x P x x dx se i j
i
P x P x x dx c c
λ
λ λ
λ λ
−
−
−
−
⎧ + −= =⎪ − +⎪
⎪ −⎪ = = −−⎪⎪= ⎨−⎪ = = +⎪ −
⎪⎪
=⎪⎪⎩
∫
∫
∫
∫
A ; (A.15)
83
( )( )( )( )
( )( )
( )( )
( )( )( )( )
13
1
213
21
213 3ij 2
1
13
1
1
1
1 2( ) ( ) , 3
2 2 1 2 3
3 2( ) ( ) , 1
2 4 9
3 2( ) ( ) , 1
2 4 9
1 2( ) ( ) , 3
2 2 1 2 3
( ) ( )
i ji i
i ji j
i ji j
i ji j
i j
j j jP x P x x dx se i j
j j
j jP x P x x dx se i j
j
i iP x P x x dx se i j
i
i i iP x P x x dx se i j
i i
P x P x x
λ λ
λ λ
λ λ
λ λ
−
−
−
−
−
− −= = −
− −
−= = −
−
−= = = +
−
− −= = +
− −
∫
∫
∫
∫
∫
A
3 0 , . .dx c c
⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪
=⎪⎪⎩
; (A.16)
1
10ij 1
1
( )( ) 2 ,
( )( ) 0 , . .
ji i j
ji
P xP x dx se i j i j é ímpar
xP x
P x dx c cx
λ λ−
−
∂⎧= < ∧ +⎪ ∂⎪= ⎨
∂⎪ =⎪ ∂⎩
∫
∫B ; (A.17)
1
11
1ij
11
1
( )( ) ,
( )( ) 2 ,
( )( ) 0 , . .
ji
ji i j
ji
P xP x x dx i se i j
xP x
P x x dx se i j i j é parx
P xP x x dx c c
x
λ λ
−
−
−
∂⎧= =⎪ ∂⎪
⎪ ∂⎪= = < ∧ +⎨ ∂⎪⎪ ∂⎪ =
∂⎪⎩
∫
∫
∫
B ; (A.18)
( )12
1
1 22
12ij 1
2
11
2
1
( ) 1( ) , 1
2
( ) 3 1( ) , 12
( )( ) 2 , 1
( )( ) 0 , . .
ji
i j
ji
i j
ji i j
ji
P x i iP x x dx se i j
x
P x j jP x x dx se i jx
P xP x x dx se i j i j é ímpar
xP x
P x x dx c cx
λ λ
λ λ
λ λ
−
−
−
−
∂⎧ −= = +⎪ ∂⎪
⎪ ∂ − −⎪ = = −∂⎪
= ⎨∂⎪ = < − ∧ +⎪ ∂⎪
⎪ ∂=⎪
∂⎩
∫
∫
∫
∫
B . (A.19)
85
ANEXO B – APROXIMAÇÃO DOS INCREMENTOS DOS PARÂMETROS PLÁSTICOS
B.1 Aproximação constante
A aproximação constante é definida por uma única função, com valor unitário em todo o
seu domínio de definição. Matematicamente é definida por:
( ), 1P γ ψ∗ = . (B.1)
B.2 Aproximação linear
A aproximação linear é definida pelas quatro funções (B.2), sempre não negativas. Cada
função toma o valor unitário num dos cantos da célula plástica (ver Figura B.1):
( ) ( )( )
( ) ( )( )
( ) ( )( )
( ) ( )( )
1
2
3
4
1, 1 1
41
, 1 141
, 1 141
, 1 14
P
P
P
P
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
∗
∗
∗
∗
= − −
= + −
= + +
= − +
⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩
. (B.2)
B.3 Aproximação Quadrática
A aproximação quadrática é definida por nove funções (ver Figura B.2), sempre não
negativas. Cada função toma valor unitário num dos cantos ou a meio dos lados das células
plásticas (Figura F-2):
86
( ) ( ) ( )
( ) ( )( )
( ) ( ) ( )
( ) ( ) ( )( ) ( )( )( ) ( ) ( )
( ) ( ) ( )
( ) ( )( )
( ) ( ) ( )
2 2
2 2
2 2
2 2
22
2 2
1
22 2
2 23
24 2
5
6
7
8
9
1, 1 1
161
, 1 141
, 1 1161
, 1 14
, 1 1
1, 1 1
41
, 1 1161
, 1 141
, 1 116
P
P
P
P
P
P
P
P
P
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
γ ψ γ ψ
∗
∗
∗
∗
∗
∗
∗
∗
∗
= − −
= − −
= + −
= − −
= − −
= + −
= − +
= − +
= + +
⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩
Figura B.1: Aproximação Linear.
∗1P
∗2P
∗3P
∗4P
87
Figura B.2: Aproximação Quadrática.
∗1P
∗2P
∗4P
∗5P ∗
6P
∗7P
∗8P ∗
10P
∗3P
89
ANEXO C – CRITÉRIOS DE CEDÊNCIA
Este anexo tem como objectivo definir as expressões relacionadas com as condições de
cedência implementadas no programa. Para além da definição da própria condição de
cedência, é necessário definir as suas derivadas, pois são necessárias ao cálculo dos
operadores:
tvN S P dVφ∗
∗ ∗
∂=
∂σ∫ ,
2
tv v2M S S e dVφ∗
∗ ∗
∂= Δ
∂σ∫ .
Recordam-se as seguintes definições, úteis para expressar matematicamente os critérios de
cedência:
Vector das componentes independentes do tensor das tensões:
txx yy zz xy yz xzσ σ σ σ σ σ σ⎡ ⎤= ⎣ ⎦ . (C.1)
Primeiro invariante do tensor das tensões:
1 xx yy zzI σ σ σ= + + . (C.2)
Tensões normal e tangencial médias:
( ) 1
3 3xx yy zz
mIσ σ σ
σ+ +
= = , (C.3)
225m Jτ = . (C.4)
Tensões normal e tangencial octaédricas:
1
3oct mIσ σ= = , (C.5)
223oct Jτ = . (C.6)
Tensor das tensões deviatóricas, e o seu segundo e terceiro invariante:
3
kk ijij ijs
σ δσ= − , (C.7)
( ) ( ) ( )2 2 2 2 2 22
1 12 6ij ij xx yy yy zz zz xx xy yx zxJ s s σ σ σ σ σ σ σ σ σ⎡ ⎤= = − + − + − + + +⎢ ⎥⎣ ⎦
, (C.8)
90
313 ij jk ki ijJ s s s det s= = . (C.9)
A magnitude da componente hidrostática e deviatórica:
1 3 33 oct m
Iξ σ σ= = = , (C.10)
22 3 5oct mJρ τ τ= = = . (C.11)
O ângulo de Lode [Lode 1926]:
( ) 33 22
3 33 , ,2 6 6
Jsen comJ
π πθ θ ⎡ ⎤= − ∈ −⎢ ⎥⎣ ⎦. (C.12)
É possível definir a seguinte expressão para o cálculo das tensões principais ( )1 2 3, ,σ σ σ ,
com base nos três valores de θ contidos no intervalo ,6 6π π⎡ ⎤−⎢ ⎥⎣ ⎦
,
( )1
2 1 2 3
3
23
2 2 ,3
43
m
m
m
sen
sen com
sen
θ πσ σσ σ θ σ σ σσ σ
θ π
⎡ ⎤⎛ ⎞+⎜ ⎟⎢ ⎥⎡ ⎤ ⎡ ⎤ ⎝ ⎠⎢ ⎥⎢ ⎥ ⎢ ⎥= + > >⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎛ ⎞⎣ ⎦ ⎣ ⎦ ⎢ ⎥+⎜ ⎟⎢ ⎥⎝ ⎠⎣ ⎦
. (C.13)
Um estado de tensão é definido, na sua forma mais simples, através de um ponto P
pertencente ao espaço das tensões principais ( )1 2 3, ,σ σ σ . Para materiais com
comportamento isotrópico, pode ser igualmente definido através da sua componente no
eixo hidrostático 1 2 3σ σ σ= = , por 1 3 33 oct m
Iξ σ σ= = = , pela sua componente no
plano das tensões deviatóricas, perpendicular a esse eixo, obtida de
22 3 5oct mJρ τ τ= = = , e ainda, através do ângulo de Lode θ , medido a partir da
direcção positiva do eixo da maior tensão principal 1σ .
Assim sendo, a definição de um estado de tensão pode ser inequivocamente definido,
recorrendo a qualquer uma das formas: ( )1 2 3, ,σ σ σ , ( )1 2 3, ,I J J , ( ), ,oct octσ τ θ e ( ), ,ξ ρ θ .
91
C.1 Critério de Von Mises
Define-se para o critério de Von Mises [Von Mises 1913]:
2
3
33
2 3222
xy
yz
xz
a
b
cJ
φ
σ
σ
σ
∗∂=
∂σ
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
, (C.14)
2
2 2 2
2
2 2 2
22
2 2 22
22
2 2 2
2
2 2
2
2
1 1 13 6 6
1 13 6
13 3
1
. 1
1
xy yz xz
xy yz xz
xy yz xz
xy xy yz xy xz
yz yz xz
xz
A A AA A B AC
J J J
B B BB B C
J J J
C C CC
J J JJ
J J J
SimJ J
J
σ σ σ
σ σ σ
σ σ σφ
σ σ σ σ σ
σ σ σ
σ
∗
− − − − − − − −
− − − − − −
− − − −∂
=∂σ
− − −
− −
−
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢⎢⎣ ⎦
⎥⎥
, (C.15)
onde:
2 2 2
, ,6 6 6
a b cA B C
J J J= = = , (C.16)
2 , 2 , 2xx yy zz xx yy zz xx yy zza b cσ σ σ σ σ σ σ σ σ= − − = − + − = − − + . (C.17)
92
C.2 Critério de Drucker-Prager
Define-se para o critério de Drucker-Prager [Drucker et al. 1952]:
2
311
31 10 2 30 20 2
2
xy
yz
xz
a
b
cJ
φα
σ
σ
σ
σ
∗∂= +
∂
⎡ ⎤⎢ ⎥
⎡ ⎤ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎢ ⎥
⎢ ⎥⎣ ⎦
, (C.18)
2
2 2 2
2
2 2 2
22
2 2 22
22
2 2 2
2
2 2
2
2
1 1 13 6 6
1 13 6
11 3
1
. 1
1
xy yz xz
xy yz xz
xy yz xz
xy xy yz xy xz
yz yz xz
xz
A A AA A B AC
J J J
B B BB B C
J J J
C C CC
J J JJ
J J J
SimJ J
J
σ σ σ
σ σ σ
σ σ σφσ σ σ σ σ σ
σ σ σ
σ
∗
− − − − − − − −
− − − − − −
− − − −∂
=∂
− − −
− −
−
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢⎢⎣ ⎦
⎥⎥
, (C.19)
onde:
2 2 2
, ,6 6 6
a b cA B C
J J J= = = , (C.20)
2 , 2 , 2xx yy zz xx yy zz xx yy zza b cσ σ σ σ σ σ σ σ σ= − − = − + − = − − + . (C.21)