Elementos Finitos na mecânica estrutural

O objetivo do método dos elementos finitos na mecânica estrutural é encontrar uma função aproximada que satisfaça o sistema de equações diferenciais da teoria da elasticidade, obedecendo às condições de contorno do problema específico, com um grau de precisão satisfatório.

Por exemplo: uma barra de aço sujeita a esforços, e fixada em uma das extremidades. Para geometrias simples podemos aplicar as fórmulas de resistência de materiais. Mas com seções variáveis e/ou detalhes como furos ou chavetas, isso não é possível. Sabemos pela teoria da elasticidade que a relação entre tensões e deformações resulta em um sistema de equações diferenciais parciais. A função vetorial u(x,y,z) é a solução procurada e representa o deslocamento de cada ponto do sólido em relação à condição não deformada. Conhecendo-se u, é possível, através de suas derivadas, determinar as deformações e tensões em cada ponto da barra.

Um dos métodos de resíduos ponderados, normalmente o de Galerkin, é usado para determinar uma solução aproximada. Entretanto, não se procede na forma usual desse método, arbitrando-se um polinômio e determinando seus coeficientes por integração ao longo de todo o volume da barra.

Em vez disso, a barra é dividida em subdomínios, formando uma malha, contendo elementos volumétricos (os elementos finitos) e nós na intercessão das linhas da malha (pode haver mais nós além desses em elementos mais complexos).

A função solução aproximada para os deslocamentos, válida para cada elemento, é
$(u_{x}(x,y,z),u_{y}(x,y,z),u_{z}(x,y,z))$=$(\sum _{i=1}^{n_{e}}{(u_{xi}*N_{i})},\sum _{i=1}^{n_{e}}{(u_{yi}*N_{i})},\sum _{i=1}^{n_{e}}{(u_{zi}*N_{i})})$,
onde $n_{e}$ é o número de nós do elemento.

Cada uma das funções $N_{i}(x,y,z)$ é definida de maneira a assumir o valor 1 nas coordenadas correspondentes ao nó $=i$, e 0 para os demais nós do elemento. São também iguais a zero em qualquer ponto da malha fora do elemento. $u_{xi},u_{yi},u_{zi}$ são os coeficientes a determinar, conforme o método de Galerkin.

Como pode ser observado, do modo como as funções $N_{i}(x,y,z)$ são construídas, os deslocamentos de cada nó em cada direção coordenada coincidem com um dos coeficiente a determinar. Isso porque, quando $x_{k},y_{k},z_{k}$ são as coordenadas do nó$=k$, todas as funções presentes nos somatórios, exceto $N_{k}(x,y,z)$ são zero, e $N_{k}(x,y,z)=1$, portanto u$(x_{k},y_{k},z_{k})=(u_{xk},u_{yk},u_{zk})$.

Montagem do Sistema de equações lineares

Equações relativas aos elementos

Os coeficientes $u_{xi},u_{yi},u_{zi}$ são o resultado da solução de um sistema de equações lineares, como descrito nos métodos de resíduos ponderados, aplicando-se as condições de contorno do problema específico.

Para o mesmo exemplo acima, de uma barra de aço sujeita a determinados esforços, as condições de contorno essenciais são os nós em que não há deslocamento algum, ou são fixos em algumas direções coordenadas. Ou ainda quando esse deslocamento é diferente de zero mas pré-definido.

As condições de contorno de Neumann, no caso de problemas elásticos em pequenas deformações, resultam de que as deformações são derivadas parciais dos deslocamentos, ex: $\epsilon _{xx}={\frac {\partial {u_{x}}}{\partial x}}$, $\epsilon _{yy}={\frac {\partial {u_{y}}}{\partial y}}$ e $\epsilon _{zz}={\frac {\partial {u_{z}}}{\partial z}}$. E que existe uma relação linear entre as tensões e deformações, ex: $\sigma _{xx}=(\lambda +2*\mu )*\epsilon _{xx}+\lambda *\epsilon _{yy}+\lambda *\epsilon _{zz}$ para material isotrópico, onde $\lambda$ e $\mu$ são as constantes de Lamé do material.

Assim, as condições de contorno de Neumann dependem do conhecimento das tensões aplicadas em regiões da superfície de alguns dos elementos.

Para modelar o sistema de equações lineares, no caso de nosso exemplo, o ponto de partida é o sistema de equações diferenciais do modelo elástico:

${\frac {\partial \sigma _{xx}}{\partial x}}+{\frac {\partial \sigma _{yx}}{\partial y}}+{\frac {\partial \sigma _{zx}}{\partial z}}+b_{x}=0$

${\frac {\partial \sigma _{xy}}{\partial x}}+{\frac {\partial \sigma _{yy}}{\partial y}}+{\frac {\partial \sigma _{zy}}{\partial z}}+b_{y}=0$

${\frac {\partial \sigma _{xz}}{\partial x}}+{\frac {\partial \sigma _{yz}}{\partial y}}+{\frac {\partial \sigma _{zz}}{\partial z}}+b_{z}=0$

Se o peso próprio do sólido for desprezível frente às tensões aplicadas (b = 0), e trabalhando apenas com a primeira das equações (para as demais o processo é análogo):

${\frac {\partial {((\lambda +2*\mu )*\epsilon _{xx}+\lambda *\epsilon _{yy}+\lambda *\epsilon _{zz}})}{\partial x}}+{\frac {\partial {(\mu *\epsilon _{yx}})}{\partial y}}+{\frac {\partial {(\mu *\epsilon _{zx}})}{\partial z}}=0$

Explicitando as deformações como derivadas dos deslocamentos:

$(\lambda +2*\mu )*{\frac {\partial ^{2}{u_{x}}}{\partial x^{2}}}+\lambda *{\frac {\partial ^{2}{u_{y}}}{\partial y*\partial x}}+\lambda *{\frac {\partial ^{2}{u_{z}}}{\partial z*\partial x}}+\mu *{\frac {\partial ^{2}{u_{y}}}{\partial x*\partial y}}$ + $\mu *{\frac {\partial ^{2}{u_{z}}}{\partial x*\partial z}}=0$

Agora, de acordo com o método de Galerkin, devemos substituir u pela função aproximada, e a equação diferencial resulta em um resíduo em vez de zero.

Esse procedimento é feito para cada um dos elementos:
$(\lambda +2*\mu )*\sum _{i=1}^{n_{e}}{(u_{xi}*{\frac {\partial ^{2}{N_{i}}}{\partial x^{2}}})}+\lambda *\sum _{i=1}^{n_{e}}{(u_{yi}*{\frac {\partial ^{2}{N_{i})}}{\partial y*\partial x}})}+\lambda *\sum _{i=1}^{n_{e}}{(u_{zi}*{\frac {\partial ^{2}{N_{i})}}{\partial z*\partial x}})}$ + $\mu *\sum _{i=1}^{n_{e}}{(u_{yi}*{\frac {\partial ^{2}{N_{i})}}{\partial x*\partial y}})}+\mu *\sum _{i=1}^{n_{e}}{(u_{zi}*{\frac {\partial ^{2}{N_{i})}}{\partial x*\partial z}})}=res$

Onde $n_{e}$ é o número de nós de cada elemento.
E finalmente, o resíduo deve ser multiplicado por cada uma das funções de interpolação, e a integral ao longo do volume igualada a zero:

$\int _{V}N_{1}*res\,dV=0$

$\int _{V}N_{2}*res\,dV=0$

.
.

$\int _{V}N_{n}*res\,dV=0$

Para inserir as condições de contorno correspondentes às tensões aplicadas, utiliza-se o método de integração por partes. Por exemplo, o primeiro termo da primeira equação, lembrando que $u_{xi}$ são constantes, pode ser expresso como:

$\sum _{i=1}^{ne}(u_{xi}*\int _{V}{(N_{1}*{\frac {\partial ^{2}{N_{i}}}{\partial x^{2}}})}\,dV)$

A expressão dentro da integral pode ser integrada por partes, resultando em:

$\sum _{i=1}^{ne}(u_{xi}*\int _{V}{\,d(N_{1}*{\frac {\partial {N_{i}}}{\partial x}})}-u_{xi}*\int _{V}{({\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}})}\,dV)$

Na primeira integral, o termo $\sum _{i=1}^{ne}(u_{xi}*{{\frac {\partial {N_{i}}}{\partial x}})}=\epsilon _{xx}$, ou seja a deformação xx do elemento. Isso porque $\epsilon _{xx}={\frac {\partial {u_{x}}}{\partial x}}$, e $u_{x}=\sum _{i=1}^{n_{e}}{(u_{xi}*N_{i})}$.
Aplicando-se o teorema da divergência:

$\sum _{i=1}^{ne}(\int _{S}{(n*N_{1}*\epsilon _{xx})}\,dS-u_{xi}*\int _{V}{({\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}})}\,dV)$, onde n é a normal à superfície.

Usando o mesmo processo para todas as equações, e aplicando as relações lineares entre tensões e deformações:

$(\lambda +2*\mu )*\sum _{i=1}^{ne}({u_{xi}*\int _{V}{\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}}}\ dV)+\lambda *\sum _{i=1}^{ne}({u_{yi}*\int _{V}{\frac {\partial {N_{1}}}{\partial y}}*{\frac {\partial {N_{i}}}{\partial y}}}\ dV)$ + $\lambda *\sum _{i=1}^{ne}({u_{zi}*\int _{V}{\frac {\partial {N_{1}}}{\partial z}}*{\frac {\partial {N_{i}}}{\partial z}}}\ dV)+\mu *\sum _{i=1}^{ne}({u_{yi}*\int _{V}{\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}}}\ dV)$ + $\mu *\sum _{i=1}^{ne}({u_{zi}*\int _{V}{\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}}}\ dV)=\sum _{i=1}^{ne}\int _{S}({n*N_{1}*\sigma _{xx})}\,dS$ + $\sum _{i=1}^{ne}\int _{S}{(n*N_{1}*\sigma _{yx})}\,dS+\sum _{i=1}^{ne}\int _{S}{(n*N_{1}*\sigma _{zx})}\,dS$

O lado direito da equação expressaria, se não houvesse o termo $N1$, a força total na direção x sofrida pelo elemento. A multiplicação pela função de interpolação $N1(x)$ aloca que porção dessa força atua no nó 1. Agrupando os termos com coeficientes comuns no lado esquerdo:
$\sum _{i=1}^{ne}{(u_{xi}*(\lambda +2*\mu )*\int _{V}{\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}}}\ dV)$+$\sum _{i=1}^{ne}({u_{yi}*(\lambda +\mu )*\int _{V}({\frac {\partial {N_{1}}}{\partial y}}*{\frac {\partial {N_{i}}}{\partial y}}+{\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}}})\ dV)$+$\sum _{i=1}^{ne}({u_{zi}*(\lambda +\mu )*\int _{V}({\frac {\partial {N_{1}}}{\partial z}}*{\frac {\partial {N_{i}}}{\partial z}}+{\frac {\partial {N_{1}}}{\partial x}}*{\frac {\partial {N_{i}}}{\partial x}}})\ dV)=F_{1x}$

Efetuando cálculos similares para todas as integrais dos resíduos geradas pelas 3 equações diferenciais, temos um sistema relativo ao elemento.

Equações relativas à malha

Repete-se o mesmo procedimento para todos os demais elementos, gerando $3*n_{e}*El$ equações, onde $El$ é o número de elementos da malha. O número de equações é maior que o número de incógnitas (os coeficientes a determinar: $u_{xi},u_{yi},u_{zi}$) por causa dos nós em comum aos elementos adjacentes.

Para gerar o sistema da malha, as equações com o mesmo lado direito (ex:$F_{1x}$ no exemplo acima) são somadas, gerando assim uma equação para cada nó e cada direção coordenada (x, y, z).

No lado esquerdo das equações, os termos que multiplicam os coeficientes, contendo uma combinação de constantes de Lamé e derivadas parciais das funções de interpolação, formam a matriz $K$ de rigidez global do sistema. Essa matriz é sempre simétrica e positiva definida.

Entretanto, como não foram ainda inseridas as condições de contorno essenciais (nós com deslocamento zerado ou pré fixado), o sistema está ainda indeterminado, o que se verifica pelo fato de que a soma de qualquer linha ou coluna da matriz ser zero. Fisicamente, isso significa que não posso determinar o deslocamento dos nós em função de tensões aplicadas, se o sólido não está fixo, e pode girar ou se mover livremente.

O procedimento para modificar o sistema, levando em conta essas condições de contorno é exemplificada abaixo para um sistema de 3 equações e 3 incógnitas:

$3x+2y-5z=12$

$2x+6y-8z=22$

$-5x-8y+13z=-34$

A soma de qualquer das linhas ou colunas do lado esquerdo das equações do sistema é zero. Para resolver o sistema é necessário saber o valor de uma das incógnitas. Se por exemplo y=2:

$3x+2*2-5z=12$

$2x+6*2-8z=22$

$-5x-8*2+13z=-34$

Passando as constantes para o lado direito e eliminando a segunda equação:

$3x-5z=12-2*2$

$-5x+13z=-34+8*2$

O que resulta em y=2 e z=-1

O procedimento portanto é: eliminar da matriz as colunas e linhas correspondentes ao coeficiente conhecido, e alterar correspondentemente o lado direito das equações não eliminadas.

Em forma compacta, o sistema a ser resolvido é: $K.X=F$, onde $K$ é a matriz de rigidez global, com as condições essenciais de contorno, $X$ é o vetor de coeficientes, e $F$ é o vetor de forças, também modificado para as mesmas condições de contorno.

Funções de interpolação e mapeamento

Como os elemento tem nós com coordenadas x,y,z específicas, as funções $N_{i}(x,y,z)$ são diferentes para cada um deles. Entretanto não há necessidade de calculá-las várias vezes. Usa-se o procedimento de definir as funções para um elemento de referência.

Para o caso de um elemento cúbico de 8 nós, esse elemento é centrado em (0,0,0), com o comprimento das arestas = 2, e com o plano das faces perpendicular aos eixos coordenados. As coordenadas, (denominadas locais) do nós são portanto:

Elemento de referência em coordenadas locais e mapeamento a um elemento distorcido em coordenadas reais

$1)(1,-1,-1)$

$2)(1,1,-1)$

$3)(-1,1,-1)$

$4)(-1,-1,-1)$

$5)(1,-1,1)$

$6)(1,1,1)$

$7)(-1,1,1)$

$8)(-1,-1,1)$

As funções $N_{i}(x_{L},y_{L},z_{L})$ devem ser definidas para assumir o valor 1 no nó correspondente e 0 nos demais:

$N1=1/8*(x_{L}+1)*(y_{L}-1)*(z_{L}-1)$

$N2=-1/8*(x_{L}+1)*(y_{L}+1)*(z_{L}-1)$

$N3=1/8*(x_{L}-1)*(y_{L}+1)*(z_{L}-1)$

$N4=-1/8*(x_{L}-1)*(y_{L}-1)*(z_{L}-1)$

$N5=-1/8*(x_{L}+1)*(y_{L}-1)*(z_{L}+1)$

$N6=1/8*(x_{L}+1)*(y_{L}+1)*(z_{L}+1)$

$N7=-1/8*(x_{L}-1)*(y_{L}+1)*(z_{L}+1)$

$N8=1/8*(x_{L}-1)*(y_{L}-1)*(z_{L}+1)$

Qualquer outro elemento cúbico de 8 nós da malha, com outras dimensões de aresta, e ainda que esteja distorcido com arestas de comprimentos diferentes e/ou ângulos diferentes de 90 graus pode ser mapeado pelo elemento de referência. Cada ponto de coordenadas $(x_{L},y_{L},z_{L})$ do elemento de referência corresponderá ao ponto $N_{i}(x,y,z)$ do elemento real da malha através das funções:

$x=\sum _{i=1}^{8}N_{i}(x_{L},y_{L},z_{L})*x_{i}$

$y=\sum _{i=1}^{8}N_{i}(x_{L},y_{L},z_{L})*y_{i}$

$z=\sum _{i=1}^{8}N_{i}(x_{L},y_{L},z_{L})*z_{i}$

onde $x_{i},y_{i},z_{i}$ são as coordenadas de cada nó do elemento da malha.

Para calcular a matriz K de rigidez dos elementos, é necessário saber os valores das derivadas parciais das funções $N_{i}(x,y,z)$ em relação a cada um dos eixos coordenados (derivadas reais). E depois disso, integrar cada termo da matriz ao longo do volume do elemento.

Cálculo das derivadas reais

Como as derivadas de $N_{i}(x_{L},y_{L},z_{L})$ (derivadas locais) são conhecidas, bastando aplicar as regras de derivação às funções definidas acima, podemos saber as derivadas reais usando a regra da cadeia para funções de várias variáveis:

${\frac {\partial N_{i}}{\partial x}}={\frac {\partial N_{i}}{\partial x_{L}}}*{\frac {\partial x_{L}}{\partial x}}+{\frac {\partial N_{i}}{\partial y_{L}}}*{\frac {\partial y_{L}}{\partial x}}+{\frac {\partial N_{i}}{\partial z_{L}}}*{\frac {\partial z_{L}}{\partial x}}$.

As expressões para ${\frac {\partial N_{i}}{\partial y}}$ e ${\frac {\partial N_{i}}{\partial z}}$ são similares.

Em forma compacta, o vetor das derivadas reais ${\frac {\partial N_{i}}{\partial X}}=J_{X_{L}->X}*{\frac {\partial N_{i}}{\partial X_{L}}}$, onde $J_{X_{L}->X}$ é a matriz jacobiana da transformação e ${\frac {\partial N_{i}}{\partial X_{L}}}$ é o vetor das derivadas locais.

A matriz jacobiana $J_{X_{L}->X}$é:

${\begin{Bmatrix}{\frac {\partial x_{L}}{\partial x}}&{\frac {\partial y_{L}}{\partial x}}&{\frac {\partial z_{L}}{\partial x}}\\{\frac {\partial x_{L}}{\partial y}}&{\frac {\partial y_{L}}{\partial y}}&{\frac {\partial z_{L}}{\partial y}}\\{\frac {\partial x_{L}}{\partial z}}&{\frac {\partial y_{L}}{\partial z}}&{\frac {\partial z_{L}}{\partial z}}\\\end{Bmatrix}}$

Entretanto, pelas funções de mapeamento de x, y e z, o que pode ser calculado são as derivadas das coordenadas reais em relação às locais:

${\frac {\partial x}{\partial x_{L}}}=\sum _{i=1}^{8}{\frac {\partial N_{i}(x_{L},y_{L},z_{L})}{\partial x_{L}}}*x_{i}$

${\frac {\partial y}{\partial y_{L}}}=\sum _{i=1}^{8}{\frac {\partial N_{i}(x_{L},y_{L},z_{L})}{\partial y_{L}}}*y_{i}$

${\frac {\partial z}{\partial z_{L}}}=\sum _{i=1}^{8}{\frac {\partial N_{i}(x_{L},y_{L},z_{L})}{\partial z_{L}}}*z_{i}$

e não seu inverso como exige a matriz jacobiana.
Por isso, calcula-se a matriz jacobiana $J_{X->X_{L}}$: ${\begin{Bmatrix}{\frac {\partial x}{\partial x_{L}}}&{\frac {\partial y}{\partial x_{L}}}&{\frac {\partial z}{\partial x_{L}}}\\{\frac {\partial x}{\partial y_{L}}}&{\frac {\partial y}{\partial y_{L}}}&{\frac {\partial z}{\partial y_{L}}}\\{\frac {\partial x}{\partial z_{L}}}&{\frac {\partial y}{\partial z_{L}}}&{\frac {\partial z}{\partial z_{L}}}\\\end{Bmatrix}}$

inverte-se essa matriz 3x3, e aplica-se a transformação: ${\frac {\partial N_{i}}{\partial X}}=J_{X->X_{L}}^{-1}*{\frac {\partial N_{i}}{\partial X_{L}}}$

Cálculo das integrais

Com os resultados da seção anterior, as derivadas reais são espressas em coordenadas locais, referentes ao elemento de referência. Mas elas devem ser integradas ao longo do volume do elemento real. Essa transformação de variáveis exige a substituição de $dV$ por $\det {J}*dV_{L}$, onde $\det {J}$ é o determinante da matriz jacobiana, e $dV_{L}$ é o elemento de integração de volume do cubo de referência.

Para codificação em computador, o cálculo analítico dessas integrais não é prático. Em vez disso usa-se a quadratura de Gauss, que transforma a integral num somatório. Os limites das fórmulas de Gauss (entre -1 e 1) são os limites do cubo de referência. Essa integral numérica é exata no caso de elementos cúbicos não distorcidos angularmente.