Análise da influência de microestruturas heterogêneas na resposta macromecânica do problema bidimensional de placas

Analysis of the heterogeneous microstructures influence on the macromechanical response of the two-dimensional problem of plates

RESUMO

Neste trabalho apresentam-se análises do problema bidimensional de placas compostas por materiais heterogêneos com acoplamento do MEC (Método dos Elementos de Contorno) e MEF (Método dos Elementos Finitos), através de uma abordagem multi-escala. O MEC é adotado para modelar o problema definido no macro-contínuo, que neste trabalho é dado pela análise não-linear do problema bidimensional de placas, enquanto o problema de equilíbrio na micro-escala (representada pelo Elemento de Volume Representativo - EVR) é resolvido pelo MEF. É importante notar, que na formulação não-linear do MEC é adotado o operador tangente consistente ao longo do processo iterativo, a fim de reduzir o esforço computacional, que é muito importante em uma análise multi-escala. Um EVR deve estar associado a cada ponto do macro-contínuo, onde se faz necessário conhecer as tensões e o tensor constitutivo a fim de solucionar o problema não-linear da placa. Para solucionar o problema de equilíbrio do EVR, devem-se impor ao mesmo, condições de contorno em termos de flutuações dos deslocamentos. Depois de resolver o problema de equilíbrio do EVR, a passagem do micro-contínuo para o macro-contínuo é feita adotando-se técnicas de homogeneização para os campos de tensões e do tensor constitutivo, que permitem calcular os respectivos valores homogeneizados para um ponto do macro-contínuo. Nos exemplos numéricos são definidos diferentes EVRs, os quais podem ter inclusões ou vazios no seu domínio, sendo cada fase modelada por modelos constitutivos próprios. Os resultados confirmam as potencialidades de aplicação da modelagem proposta.

Palavras-chave:
Modelagem multi-escala; homogeneização; elementos de contorno; problema bidimensional de placas

ABSTRACT

A multi-scale modelling for analyzing the stretching problem of plates composed of heterogeneous materials using a coupling BEM (Boundary Element Method) and FEM (Finite Element Method) and a multiscale approach is presented. The BEM is adopted to model the macro-continuum problem (represented by the non-linear analysis of the stretching problem) while the equilibrium problem at micro-scale (represented by the Representative Volume Element - RVE) is solved by a FEM formulation. It is important to note that in the non-linear analysis of the stretching problem the consistent tangent operator has been considered along the iterative procedure, what is very important to reduce the computational effort. One RVE has to be assigned to every point of the macro-continuum where the stresses and constitutive tangent operator are required to solve the non-linear plate problem. To solve the RVE equilibrium problem, boundary conditions in terms of displacement fluctuations have to be imposed. After solving the equilibrium problem of the RVE the micro-to-macro transition can be made by applying homogenization techniques for the fields of stress and constitutive tensor in order to compute their homogenized values for a macro-continuum point. In the numerical examples, to define the material microstructure, different RVEs composed of inclusions or voids are considered, each one modeled by properly constitutive models. The results confirm the potentialities of the proposed modeling.

Keywords:
Multi-scale modeling; homogenization; boundary elements; stretching plate problem

1. INTRODUÇÃO

Apesar dos materiais usados em Engenharia (metais, polímeros, compósitos, concretos e madeira) terem diferentes microestruturas, em nível macroscópico observam-se características similares do seu comportamento mecânico, como por exemplo: elasticidade, viscosidade, deformação plástica, ruptura frágil, ruptura dúctil, etc. Por causa dessas semelhanças, nas análises numéricas de estruturas usualmente são adotados modelos constitutivos baseados na mecânica do contínuo e na termodinâmica dos sólidos. No entanto, é importante observar que os processos de deformação e ruptura ocorrem em nível microscópico. Assim, muitos trabalhos foram desenvolvidos a fim de analisar o fenômeno dissipativo dos materiais na micro-escala, ver [¹–⁵]. Neste contexto, modelar o comportamento de estruturas compostas por materiais heterogêneos em diferentes escalas é muito importante para sua melhor representação, ver [⁶–¹⁰]. Em muitos casos, a análise não-linear convencional não consegue representar de forma precisa o comportamento de tais estruturas, ver [¹¹].

A formulação multi-escala apresentada neste trabalho utiliza-se do conceito de existência de um Elemento de Volume Representativo (EVR), o qual representa o comportamento de um ponto da macroestrutura procurando evidenciar os fenômenos dissipativos que ocorrem na microestrutura e, por conta de técnicas de homogeneização, leva informações dos tensores de tensão e constitutivo do micro para o macro-contínuo, ver [¹²–¹⁷]. Cada EVR pode possuir vazios e/ou inclusões com diferentes comportamentos mecânicos levando, portanto, a uma heterogeneidade da microestrutura do material.

Em termos gerais, numa análise multi-escala, após a aplicação de um estado de excitação no macro-contínuo, um estado de deformação é obtido, sendo o mesmo aplicado a um ponto do macro-contínuo. Essa deformação dita macroscópica, por sua vez, é aplicada ao ponto que é representado por um EVR. Com condições de contorno particulares, o EVR deforma-se e o seu problema de equilíbrio é resolvido obtendo-se os campos de deformação, tensão e deslocamentos no EVR, o que requer um processo iterativo de solução. A solução numérica apresenta resultados de tensão e relação constitutiva representativas do EVR através da aplicação de técnicas de homogeneização e médias volumétricas, permitindo passar informações atualizadas do micro para o macro-contínuo. Um novo estado de excitação é aplicado no macro-contínuo, continuando o processo incremental. O Método dos Elementos Finitos (MEF) é utilizado para modelar a microestrutra, enquanto que o Método dos Elementos de Contorno (MEC) é utilizado para modelar a macroescala, ver [¹²–¹³].

Ao se utilizar o MEC para a obtenção da solução da macro-escala, deve-se discretizar o contorno da placa em elementos nos quais os deslocamentos e esforços são aproximados. Além disso, na análise não-linear necessita-se também discretizar o domínio da placa em células, onde são aproximadas as forças normais iniciais ou inelásticas. Assim, para a se obter a solução em multi-escala, deve-se definir um EVR em cada nó de célula do macro.

Nos exemplos numéricos deste trabalho são considerados EVRs com inclusões elásticas ou vazios definidos no seu interior, adotando-se os critérios de von Mises ou Mohr Coulomb para governar o comportamento não-linear do material da matriz. A influência da distribuição e proporcionalidade dessas heterogeneidades da microestrutura na macroestrutura é foco do presente trabalho. Portanto, este trabalho é um desenvolvimento dos trabalhos anteriores dos autores voltado para o problema bidimensional de chapa composta por materiais heterogêneos e com a inclusão do modelo de Mohr-Coulomb com vistas a futuras aplicações em materiais compósitos cerâmicos, como o concreto. Além disso, procura-se evidenciar a importância de se utilizar uma modelagem como a aqui proposta, para simular o comportamento mecânico da microestrutura e seu reflexo na macroestrutura, o que não é razoável em se pensar quando modelos constitutivos fenomenológicos são empregados diretamente na macroestrutura, como ilustrado no primeiro exemplo deste trabalho.

O trabalho está dividido em 4 seções. Na primeira seção, uma breve introdução sobre o tema é apresentada. Na seção 2, a formulação do MEC para a análise da macroescala é abordada, assim como os conceitos da formulação multiescala utilizando MEF são apresentados. O algoritmo de solução do problema encerra a seção 2. Já os exemplos numéricos são discutidos na seção 3. Por fim, a seção 4 apresenta algumas conclusões e considerações finais.

2. MATERIAIS E MÉTODOS

2.1 FORMULAÇÃO DO MEC PARA ANÁLISE NÃO-LINEAR DO PROBLEMA BIDIMENSIONAL DE PLACAS

2.1.1 Relações Básicas

Para definir o problema considere uma placa de contorno externo Γ e domínio Ω, onde as variáveis são definidas segundo um sistema de coordenadas cartesianas onde x₁ e x₂ são as direções no plano da placa. Assume-se que a placa suporta apenas cargas distribuídas no plano da sua superfície média. Como o presente trabalho trata de análise não-linear, todas as variáveis do problema serão definidas em termos de taxas, ou seja, suas derivadas no tempo, sendo (ẋ) = dx/dt. Assim, as variáveis relacionadas ao problema bidimensional são:forças (ṗ_n e ṗ_s) e deslocamentos (${\dot{u}}_n$ e ${\dot{u}}_s$), onde (n, s) é o sistema de coordenadas locais, sendo n e s, respectivamente, as direções normal e tangencial ao contorno da placa. As equações básicas do problema bidimensional serão omitidas aqui, mas as mesmas podem ser encontradas em [¹⁸–²⁴].

As taxas das forças internas de membrana Ṅ_ij são obtidas integrando-se as tensões ${\dot{\sigma }}_{ij}$ ao longo da espessura t_h da placa, as quais, considerando-se estado plano de tensão, resultam em:

(1) $$\begin{array}{cc}{\dot{N}}_{ij}={\int }_{-t_h/2}^{t_h/2}{\dot{\sigma }}_{ij}dz=t{\dot{\sigma }}_{ij}& \text{i},\text{j}=1,2\end{array}$$

No contexto não-linear convencional, as tensões (${\dot{\sigma }}_{ij}$) são obtidas a partir do modelo constitutivo adotado. Porém, neste trabalho como é feita uma abordagem em multi-escala, essas tensões são calculadas após a solução do EVR, adotando-se uma técnica de homogeneização para se fazer a passagem da micro-escala para a macro-escala. Por outro lado, as deformações podem ser divididas em suas parcelas elásticas ${\dot{\epsilon }}_{ij}^e$ e plásticas ${\dot{\epsilon }}_{ij}^p$ como segue:

(2) $${\dot{\epsilon }}_{ij}={\dot{\epsilon }}_{ij}^e+{\dot{\epsilon }}_{ij}^p$$

Observe que, aplicando-se a lei de Hooke, as tensões ${\dot{\sigma }}_{ij}$ são relacionadas às deformações elásticas, sendo ${\dot{\epsilon }}_{ij}^e={\dot{\epsilon }}_{ij}-{\dot{\epsilon }}_{ij}^p$ e as forças internas de membrana de tentativa ${\dot{N}}_{ij}^e$ relacionadas à deformação total ${\dot{\epsilon }}_{ij}$. Portanto, as forças ${\dot{N}}_{ij}^e$ podem ser escritas em termos das derivadas dos deslocamentos totais como segue:

(3) $$\begin{array}{cc}{\dot{N}}_{\text{ij}}^{\text{e}}=\mu \left[{\dot{\text{u}}}_{\text{i},j}+{\dot{u}}_{\text{j},i}\right]+\frac{2\mu {\nu }^{\prime }}{1-2{\nu }^{\prime }}{\dot{u}}_{k,k}{\delta }_{\text{ij}}& \text{i},\text{j},\text{k}=1,2\end{array}$$

onde δ_ij é o delta de Kronecker; v′= v/(1 + v) para estado plano de tensão, sendo v o coeficiente de Poisson; μ é o modulo elástico de cisalhamento, que tem que ser multiplicada pela espessura t para estado plano de tensão.

Portanto, a taxa de forças internas de membrana inelásticas ${\dot{N}}_{ij}^p$ podem ser definidas como:

(4) $${\dot{N}}_{ij}^p={\dot{N}}_{ij}^e-{\dot{N}}_{ij}$$

2.1.2 Equações integrais e Algébricas

A equação integral do deslocamento para um ponto interno é obtida a partir do Teorema de Reciprocidade de Betti, que no caso do problema bidimensional de placas é dado por (ver mais detalhes em [²⁰]):

(5) $$\begin{array}{cc}{\int }_{\Omega }{\epsilon }_{ijk}^{\text{*}}{N}_{jk}d\text{Ω}={\int }_{\Omega }{N}_{ijk}^{\text{*}}{\dot{\epsilon }}_{jk}d\text{Ω}-{\int }_{\Omega }{\epsilon }_{ijk}^{\text{*}}{\dot{N}}_{ij}^{\left(p\right)}d\text{Ω}& \text{i},\text{j},\text{k}=1,2\end{array}$$

onde os termos com índice * são relacionados ao problema fundamental, sendo i a direção da carga fundamental.

Integrando-se a equação (5) por partes, se obtém a conhecida representação de deslocamento no plano da placa para um ponto interno:

(6) $$K_i{\dot{u}}_i=-{\int }_{\Gamma }\left({\dot{u}}_n{p}_{in}^{\text{*}}+{\dot{u}}_s{p}_{is}^{\text{*}}\right)d\text{Γ}+{\int }_{\Gamma }\left(u_{in}^{\text{*}}{\dot{p}}_n+u_{is}^{\text{*}}{\dot{p}}_s\right)d\text{Γ}+{\int }_{{\Omega }_{\text{b}}}\left(u_{in}^{\text{*}}{\dot{b}}_n+u_{is}^{\text{*}}{\dot{b}}_s\right)d\text{Ω}+{\int }_{\Omega }{\epsilon }_{ijk}^{\text{*}}{\dot{N}}_{jk}^{p}d\text{Ω}$$

onde $u_n^{*}$, $u_s^{*}$, $p_n^{*}$ e $p_s^{*}$ são valores fundamentais das forças e deslocamentos no plano da placa, Ω_b a região carregada da placa; o termo livre K_i(q) podem assumir os valores, respectivamente, para pontos internos e sobre o contorno não coincidente com um canto: K(q) = 1 e K(q) = 1/2.

Note que a solução não-linear é obtida a partir de um processo incremental de carga. Assim considerando-se Δt = t_n+1 – t_n um passo de tempo, relacionado a um incremento de carga, o problema consiste em achar a solução no passo de tempo t_n+1, sendo a solução já conhecida em t_n. Desse modo, nas equações anteriores os valores escritos em termos de taxas devem ser substituídos por seus respectivos valores em incrementos. A representação integral das derivadas dos deslocamentos Δu_i,ℓ é obtida derivando-se a equação (6), ou seja:

(7) $$\text{Δ}{u}_{i,\ell }=-{\int }_{\Gamma }\text{Δ}\left(u_n{p}_{in,\ell }^{\text{*}}+u_s{p}_{is,\ell }^{\text{*}}\right)d\text{Γ}+{\int }_{\Gamma }\left(u_{in,\ell }^{\text{*}}\text{Δ}{p}_n+u_{is,\ell }^{\text{*}}\text{Δ}{p}_s\right)d\text{Γ}+{\int }_{{\Omega }_b}\left(u_{in,\ell }^{\text{*}}\text{Δ}{b}_n+u_{is,\ell }^{\text{*}}\text{Δ}{b}_s\right)d{\text{Ω}}_{\text{b}}++\frac{\partial }{\partial x_1}{\int }_{\Omega }{\epsilon }_{ijk}^{\text{*}}\text{Δ}{N}_{jk}^{p}d\text{Ω}$$

Para transformar as equações integrais em algébricas, o contorno é discretizado em elementos com aproximação quadrática das variáveis, enquanto o domínio é discretizado em células com aproximação linear das forças de membrana inelásticas (ou plásticas). Ao longo do contorno têm-se quatro valores definidos: Δu_n, Δu_s, Δp_n e Δp_s, sendo dois desses valores dados como condição de contorno. Portanto, duas equações de deslocamentos serão escritas em cada nó de contorno. Note que o cálculo das integrais de domínio envolvendo o carregamento, assim como as integrais nas células são calculadas transformando-as em integrais de contorno. Além disso, a fim de melhorar a precisão numérica adota-se a técnica de sub-elementos na integração numérica das células e dos elementos de contorno (ver mais detalhes em [²¹]).

Para se obter a solução não-linear, deve-se ainda escrever três equações das forças de membrana elásticas de tentativa $\text{Δ}{N}_{jk}^{e\left(MEC\right)}$ nos nós das células. Essas equações são obtidas considerando-se a equação (3), onde as derivadas dos deslocamentos são definidas pela equação (7). Após escrever todas as equações necessárias, obtém-se um sistema de equações, que após aplicar as condições de contorno, pode ser escrito como (ver mais detalhes em [¹²]):

(8) $$\text{Δ}X=\text{Δ}L+R_N\text{Δ}{N}^P$$

onde o vetor ΔX contêm as incógnitas do contorno, ΔL representa a parte elástica dessas incógnitas, R_N representa as correções devido ao incremento das forças de membrana inelásticas ΔN^P.

Escrevendo-se as equações algébricas de $\text{Δ}{N}_{jk}^{e\left(MEC\right)}$ para todos os nós de células, obtém-se (ver mais detalhes em [¹²]):

(9) $$\text{Δ}{N}^{e\left(MEC\right)}=\text{Δ}K+S_N\text{Δ}{N}^P$$

onde ΔK é a solução elástica dada em termos de incremento de forças de membrana, S_N expressa o efeito do incremento de forças de membrana inelásticas ΔN^P.

Observe que ΔN^P, ΔN^e e ΔN definidos, respectivamente, nas equações (4), (3) e (1) são calculados localmente para um determinado ponto, isto é, são obtidos levando em conta apenas seu incremento de deformação Δε_ij e de tensões Δσ_ij (obtidas após a solução do EVR). Por outro lado, os incrementos elásticos de forças de membrana ΔN^e(MEC) definidos na equação (9) são calculados levando-se em conta os incrementos de forças de membrana inelásticas ΔN^P de todos os nós da placa. Neste contexto, pode-se também definir a equação algébrica do incremento de forças de membrana ΔN^(MEC) (ver mais detalhes em [¹²]):

(10) $$\text{Δ}{N}^{MEC}=C_N\text{Δ}\epsilon -\text{Δ}{K}_N-S_N\text{Δ}{N}^P+\text{Δ}N$$

onde os valores nodais ΔN^P são dados por:

(11) $$\text{Δ}{N}^P=C_N\text{Δ}\epsilon -\text{Δ}N$$

sendo Δε o incremento de deformações nodais da placa C_N a matriz elástica obtida a partir da equação (3).

2.1.3 Equação de equilíbrio e Operador Tangente Consistente

Para um incremento n, a equação de equilíbrio da placa é dada por:

(12) $$\text{Δ}{K}_n-\text{Δ}{N}_n^{MEC}=0$$

onde os vetores ΔK_n e $\text{Δ}{N}_n^{MEC}$ são definidos, respectivamente, nas equações (9) e (10).

Substituindo-se a equação (10) na equação (12), obtém-se a expressão final para a equação de equilíbrio ou equação de resíduos da placa:

(13) $$R_N\left(\text{Δ}{\epsilon }_n\right)=2\text{Δ}{K}_n-C_N\text{Δ}{\epsilon }_n+S_N\left(C_N\text{Δ}{\epsilon }_n-\text{Δ}{N}_n\right)-\text{Δ}{N}_n=0$$

Após aplicar à placa o incremento de deformações Δε_n e obter as tensões Δσ_ij para todos os nós de células, se a equação (13) não for satisfeita, haverá um resíduo R_N, isto é o incremento não será elástico. Nesse caso, a equação (13) será resolvida aplicando-se o método de Newton-Raphson, no qual se necessita de um processo iterativo para obter-se o valor do incremento Δε_n que satisfaz a equação de equilíbrio da placa. Seja uma iteração i onde o incremento de curvatura $\text{Δ}{\epsilon }_n^i$ é conhecido. O próximo incremento de tentativa $\text{Δ}{\epsilon }_n^{i+1}$ na iteração (i+1) é obtido adicionado as correções $\delta \text{Δ}{\epsilon }_n^{i+1}$, ou seja:

(14) $$\text{Δ}{\epsilon }_n^{i+1}=\text{Δ}{\epsilon }_n^i+\delta \text{Δ}{\epsilon }_n^{i+1}$$

onde as correções $\delta \text{Δ}{\epsilon }_n^{i+1}$ são calculadas linearizando-se a equação (13), sendo dadas por:

(15) $$\delta \text{Δ}{\epsilon }_n^{\text{i}+1}=-{\left[\frac{\partial R_N\left(\text{Δ}{\epsilon }_n^i\right)}{\partial \text{Δ}{\epsilon }_n^i}\right]}^{-1}{R}_N\left(\text{Δ}{\epsilon }_n^i\right)$$

onde $-\frac{\partial R_N(\Delta {\epsilon }_n^i)}{\partial \Delta {\epsilon }_n^i}$ o operador tangente consistente obtido derivando-se a equação (13), resultando em:

(16) $$-\frac{\partial R_N(\text{Δ}{\epsilon }_n^i)}{\partial \text{Δ}{\epsilon }_n^i}=S_N\left(C_{N(n)}^{ep(i)}-C_N\right)+C_N+C_{N(n)}^{ep(i)}$$

Na equação (16) $C_{N(n)}^{ep(i)}$ é a matriz que contém os tensores elasto-plásticos ${\left[C_n^{ep}\right]}_k^i$ para todos os nós de células. Para um nó k, ${\left[C_n^{ep}\right]}_k^i$ é dado por: ${\left[C_n^{ep}\right]}_k^i=\frac{\partial N_k^i}{\partial \text{Δ}{\epsilon }_k^i}=t{\left[C_{ep}\right]}_k^i$, sendo ${\left[C_{ep}\right]}_k^i=\frac{\partial \left({\sigma }_k^i\right)}{\partial {\epsilon }_k^i}$, isto, é, ${\left[C_{ep}\right]}_k^i$ é o tensor constitutivo que relaciona tensões e deformações, que na análise em multi-escala é obtido após solução do processo iterativo do EVR.

Para um incremento n, o processo numérico começa fazendo-se um passo elástico de tentativa, onde se calculam os incrementos de deformações para todos os nós de células. Então, esse incremento de deformações é imposto aos EVRs relacionados a cada nó de célula e o incremento de tensão Δσ_ij obtido após resolver o problema de equilíbrio do respectivo EVR discutido adiante. Calculadas as tensões em todos os nós de célula, as forças internas de membrana (equação (1)) podem ser calculadas e a equação (13) verificada. Se ela não for satisfeita, as correções $\delta \text{Δ}{\epsilon }_n^{i+1}$ (equação (15)) têm que ser calculadas continuando o processo na iteração i+1.

O critério de convergência adotado para o problema de equilíbrio da placa é dado por: $\frac{\sqrt{R_N^T{R}_N}}{\sqrt{{\text{ΔK}}_{\text{N}}^{\text{T}}\text{Δ}{K}_N}}\le tol$, onde tol é a tolerância adotada no processo iterativo.

2.2 MODELAGEM DA MICRO-ESCALA UTILIZANDO O MEF

A micro-escala é representada pelo EVR, cujo problema de equilíbrio é resolvido pelo Método dos Elementos Finitos. Essa formulação do MEF é discutida mais detalhadamente nos trabalhos [¹²–¹⁵], sendo neste trabalho abordados apenas seus conceitos principais.

2.2.1 Definição do EVR, Passagem do Micro para o Macro e Campo de deslocamentos no EVR

Seja o macrocontínuo representado na figura (1), que neste trabalho é definido pela placa, onde x é um ponto qualquer do macrocontínuo. Na análise em multi-escala, cada ponto x é representado pelo EVR (Elemento de Volume Representativo) (ver [¹²] – [¹⁵]), onde define-se V_μ como seu volume, Ω_μ como seu domínio e ∂Ω_μ como seu contorno, sendo y um ponto qualquer do mesmo. Neste trabalho, deve-se definir um EVR para cada ponto de Gauss definido ao longo da espessura e referente a um determinado nó de célula. Observe que a dimensão, discretização e definição dos nós no RVE são completamente independentes das dimensões ou discretização do problema estudado na macro-escala. É importante dizer que, na entrada de dados definem-se os nós e a discretização para apenas um EVR padrão, pois esses dados do EVR se mantêm constantes, ou seja, são independentes da posição que ele ocupe na placa.

Figura 1
Definição do macro-contínuo e da micro-estrutura

Note na figura (1) que o EVR pode ser composto por vazios (domínio ${\text{Ω}}_{\mu }^v$) e por partes sólidas (domínio ${\text{Ω}}_{\mu }^s$), sendo ${\text{Ω}}_{\mu }={\text{Ω}}_{\mu }^v\cup {\text{Ω}}_{\mu }^s$. Além disso, a parte sólida pode ser composta de várias fases, as quais podem ter módulo de elasticidade e coeficiente de Poisson diferentes uma da outra, além de poderem ser governadas por diferentes modelos constitutivos. Para simplificar, no que segue será considerado apenas o caso que os buracos não interceptam o contorno do EVR.

Note que se o EVR não possuir vazios, nem inclusões, isto é, for constituído por apenas uma fase, a resposta numérica da estrutura seria igual àquela obtida por uma análise não-linear convencional. Assim, a heterogeneidade do material é introduzida pela definição de inclusões ou vazios no interior da sua matriz.

Neste trabalho, o EVR também é considerado como meio contínuo e, portanto o conceito de tensão permanece válido em nível microscópico. Assim, a tensão microscópica pode ser escrita em termos de deformação como: σ_μ(y,t)=f_y(ε_μ(y,t)), onde f_y é o tensor constitutivo definido de acordo com o critério adotado ou é dado pela lei de Hooke se um comportamento elástico for adotado para a fase.

Assume-se que o tensor de deformação ε(x,t) assim como o tensor das tensões σ(x,t) referentes a um ponto x do macro-contínuo seja a média volumétrica de seus respectivos campos microscópicos (ε_μ=ε_μ(y,t) ou σ_μ=σ_μ(y,t)) do EVR associado a x, ou seja, aplicando-se uma técnica de homogeneização, para um instante qualquer t, tem-se:

(17) $$\epsilon (x,t)=\frac{1}{V_{\mu }}{\int }_{{\text{Ω}}_{\mu }}{\epsilon }_{\mu }(y,t)dV$$

(18) $$\sigma (x,t)=\frac{1}{V_{\mu }}{\int }_{{\text{Ω}}_{\mu }}{\sigma }_{\mu }(y,t)dV$$

onde os tensores ε = ε(x, t) e σ = σ(x, t) são designados, respectivamente, deformação e tensão homogeneizados.

Além disso, a deformação microscópica ε_μ pode ser escrita em termos do campo de deslocamento microscópico u_μ como segue:

(19) $${\epsilon }_{\mu }\left(y,t\right)={\nabla }^S{u}_{\mu }(y,t)$$

onde ∇^s é o operador gradiente simétrico.

Pelo processo de homogeneização, pode-se também definir o tensor constitutivo homogeneizado C^ep definido na equação (16) como:

(20) $$C^{ep}(x,t)=\frac{\partial \sigma (x,t)}{\partial \epsilon (x,t)}=\frac{\frac{1}{V_{\mu }}{{\displaystyle \int }}_{{\text{Ω}}_{\mu }}\partial {\sigma }_{\mu }(y,t)dV}{\partial \epsilon (x,t)}=\frac{\frac{1}{V_{\mu }}{{\displaystyle \int }}_{{\text{Ω}}_{\mu }}\partial f_y({\epsilon }_{\mu }(y,t))dV}{\partial \epsilon (x,t)}$$

Portanto, após alcançar convergência do processo iterativo definido no EVR, as tensões e o tensor constitutivo referente ao macro-contínuo podem ser obtidos a partir das equações (18) e (20).

No presente trabalho, o campo de deslocamentos u_μ do EVR é dividido da seguinte maneira:

(21) $$u_{\mu }(y,t)=\epsilon (x,t)y+{\tilde{u}}_{\mu }(y,t)$$

onde a parcela εy varia linearmente, sendo obtida multiplicando-se a deformação constante ε imposta pelo macro pelas coordenadas do ponto y do EVR; a parcela ${\tilde{u}}_{\mu }$ é denotada flutuação de deslocamentos e representa a variação de deformação no EVR, isto é, se a deformação no EVR é constante, tem-se ${\tilde{u}}_{\mu }$ nulo. Analogamente, a deformação microscópica pode ser escrita como:

(22) $${\epsilon }_{\mu }(y,t)=\epsilon (x,t)+{\tilde{\epsilon }}_{\mu }(y,t)$$

onde ε é a deformação homogênea imposta pelo macro-contínuo e ${\tilde{\epsilon }}_{\mu }$ é o campo das flutuações de deformações, sendo definida como: ${\tilde{\epsilon }}_{\mu }(y,t)={\nabla }^S{\tilde{u}}_{\mu }$.

A condição para se ter um campo de flutuações de deslocamentos cinematicamente admissível é que ${\tilde{u}}_{\mu }ϵ{\tilde{K}}_{\mu }^{*}$, sendo ${\tilde{K}}_{\mu }^{*}$ o espaço vetorial dos campos de flutuações de deslocamentos minimamente restringidos e cinematicamente admissíveis, definido como (ver mais detalhes em [¹⁵]):

(23) $${\tilde{K}}_{\mu }^{\text{*}}\equiv \left\{v/{\int }_{\partial {\text{Ω}}_{\mu }}v{\otimes }_{S}ndA=0\right\}$$

onde faz-se uso da seguinte expressão $u{\otimes }_{s}v=\frac{1}{2}\left(u\otimes v+v\otimes u\right)$, no cálculo das integrais, válida para quaisquer vetores u e v, sendo a expressão um produto tensorial resultando num tensor simétrico.

Portanto, escrevendo-se a equação (22) na forma de taxas, a taxa de deformação microscópica é dita cinematicamente admissível se:

(24) $$\begin{array}{cc}{\dot{\epsilon }}_{\mu }(y,t)={\nabla }^{\text{S}}{\dot{\tilde{u}}}_{\mu }=\dot{\epsilon }(x,t)+{\dot{\tilde{\epsilon }}}_{\mu }(y,t)& \forall {\dot{\tilde{u}}}_{\mu }\in V_{\mu }\end{array}$$

2.2.2 Equação de equilíbrio do EVR e condições de contorno

Sejam: σ_μ = σ_μ(y,t) a tensão em um ponto y do EVR, ou seja, a tensão microscópica; b = b(y,t) o campo de forças de volume que agem no EVR e t^e = t^e(y, t) o campo de forças externas de superfície que atuam no contorno ∂Ω_μ do EVR. O Princípio dos Trabalhos Virtuais estabelece que o EVR estará em equilíbrio se e somente se a seguinte equação variacional for satisfeita a cada instante t:

(25) $$\begin{array}{cc}{\int }_{{\Omega }_{\mu }^S}{\sigma }_{\mu }(y,t):{\nabla }^S\eta dV-{\int }_{{\Omega }_{\mu }^S}b(y,t).\eta dV+{\int }_{{\Omega }_{\mu }^v}{\sigma }_{\mu }(y,t):{\nabla }^S\eta dV-{\int }_{{\Omega }_{\mu }^v}b(y,t).\eta dV-{\int }_{\partial {\Omega }_{\mu }}{t}^e(y,t).\eta dA=0& \forall \eta \in V_{\mu }\end{array}$$

onde V_μ é o espaço de deslocamentos virtuais que satisfaz ${\text{V}}_{\text{μ}}={\tilde{K}}_{\mu }$ e, portanto, η é um campo arbitrário de deslocamentos virtuais.

Assumindo-se que as forças de volume que agem no buraco ou vazio são nulas, levando-se em conta o Princípio de macro-homogeneidade de Hill-Mandel (ver mais detalhes em [¹⁵]), considerando-se ainda as equações (19) e (22) e escrevendo-se a tensão em termos das deformações: σ_μ = f_y(ε_μ), onde f_y é o tensor constitutivo, a equação (25) pode ser reescrita da seguinte forma em termos de flutuação de deslocamentos:

(26) $$\int \begin{array}{cc}{}_{{\Omega }_{\mu }^S}{f}_y\left(\epsilon (x,t)+{\nabla }^S{\tilde{u}}_{\mu }(y,t)\right):{\nabla }^S\eta dV=0& \forall \eta \in V_{\mu }\end{array}$$

Para um campo η arbitrário, após a discretização do EVR em elementos, a seguinte equação de equilíbrio deve ser satisfeita para um passo de tempo Δt_n = t_n+1 – t_n e discretizaçãoh, sendo ${\tilde{u}}_{\mu (n+1)}={\tilde{u}}_{\mu (n)}+\Delta {\tilde{u}}_{\mu (n)}$:

(27) $$G^{n+1}={\int }_{{\Omega }_{\mu }^h}{B}^T{f}_y\left({\epsilon }^{n+1}+B{\tilde{u}}_{\mu }^{n+1}\right)dV=0$$

onde B é a matriz global que relaciona deslocamentos com deformações.

Assim, o problema de equilíbrio de EVR consiste em encontrar o campo de flutuações de deslocamentos ${\tilde{u}}_{\mu (n+1)}={\tilde{u}}_{\mu (n)}+\Delta {\tilde{u}}_{\mu (n)}$ tal que a cada instante t_n+1 a equação (27) seja satisfeita. Se o incremento n+1 é não-linear, a equação (27) é resolvida aplicando-se o método de Newton-Raphson, que consiste em encontrar as correções de flutuações $\delta {\tilde{u}}_{\mu }^{i+1}$ na iteraçãoi+1, tal que:

(28) $$F^i+K^i\delta {\tilde{u}}_{\mu }^{i+1}=0$$

sendoK a matriz de rigidez tangente e F o vetor das forças internas; no caso de se ter uma discretização com N_e elementos eles são dados por:

(29) $$F^i={\int }_{{\text{Ω}}_{\mu }^h}{B}^T{f}_y({\epsilon }_{n+1}+B{\tilde{u}}_{\mu }^i)dV={\displaystyle {\sum }_{e=1}^{Ne}{B}_e^T{\sigma }_{\mu }^{e(i)}{V}_e}$$

(30) $$K^i=\left[{\int }_{{\text{Ω}}_{\mu }^h}{B}^T{D}_{\mu }^{i}BdV\right]={\displaystyle {\sum }_{e=1}^{\mathrm{Ne}}{B}_e^T{D}_{\mu }^{e(i)}{B}_e{V}_e}$$

OndeV_e é o volume de um elemento e qualquer e $D_{\mu }^e$ é o tensor constitutivo tangente do elemento e, definido como:

(31) $$D_{\mu }^i=\left({\frac{d{f}_y}{d{\epsilon }_{\mu }}|}_{{\epsilon }_{\mu }={\epsilon }^{n+1}+B{\tilde{u}}_{\mu }^{i+1}}\right)$$

Após o cálculo das correções $\delta {\tilde{u}}_{\mu }^{i+1}$ pela equação (28), o próximo campo de flutuação de deslocamentos de tentativa a ser considerado na iteração i+1 referente à micro-estrutura, é dado por: ${\tilde{u}}_{\mu }^{i+1}={\tilde{u}}_{\mu }^i+\delta {\tilde{u}}_{\mu }^{i+1}$.

Finalmente, a formulação fica completa com a escolha do espaço V_μ, isto é, com a escolha das restrições cinemáticas a ser impostas ao EVR em termos de flutuação dos deslocamentos. Na formulação desenvolvida em [¹⁵], podem ser consideradas as seguintes condições de contorno no EVR: (i) modelo de Taylor, (ii) deslocamentos lineares, (iii) deslocamentos com variações periódicas, (iv) forças de superfície uniformes no contorno do EVR. Note que cada condição de contorno leva a uma resposta numérica diferente, definindo-se assim, diferentes modelos em multi-escala.

Note que após definir as condições de contorno, em termos de flutuação dos deslocamentos, a serem impostas ao EVR, o sistema dado pela equação (28) pode ser reduzido, pois alguns valores nodais de flutuação são prescritos nulos ou podem ser escritos em função de outros valores nodais. Portanto, a equação (28) pode ser escrita na seguinte forma reduzida:

(32) $$\delta {\tilde{u}}_R^{i+1}=-K_R^i{}^{-1}{F}_R^i$$

onde os vetores ${\left\{\delta {\tilde{u}}_R\right\}}^{i+1}$ e {F_R}ⁱ assim como a matriz [K_R]ⁱ são definidos de acordo com o modelo em multi-escala (ou condição de contorno). Para flutuações periódicas no contorno, elas serão definidos a seguir.

Neste trabalho será adotada apenas a condição de contorno que impões flutuações periódicas ao longo do contorno do EVR. Porém, outros tipos de condições de contorno são discutidos em [¹⁷].

O modelo com flutuações periódicas no contorno é adequado para descrever o comportamento de materiais que tem microestrutura periódica. No entanto, pode-se mostrar que se for utilizada uma discretização refinada, a resposta de qualquer material pode ser modelada por essa condição. Para definir o campo de flutuações neste modelo, considere a figura (2), onde são representados um EVR retangular e outro hexagonal.

Figura 2
Definição de EVRs para meio periódicos: célula retangular e hexagonal

Observe que nos EVRs definidos na figura (2) cada lado ${\text{Γ}}_{\text{i}}^{+}$ corresponde a um lado igual e oposto ${\text{Γ}}_{\text{i}}^{-}$, sendo $n_i^{+}$ a direção normal ao contorno ${\text{Γ}}_{\text{i}}^{+}$ e $n_i^{-}$ a direção normal ao contorno ${\text{Γ}}_{\text{i}}^{-}$, com $n_i^{+}=-n_i^{-}$. Assim, para cada ponto y⁺ pertencente ao contorno ${\text{Γ}}_{\text{i}}^{+}$ existe um ponto correspondente y⁻ do contorno ${\text{Γ}}_{\text{i}}^{-}$.

Neste modelo, adota-se que as flutuações do par de pontos y⁺ e y⁻ são iguais, ou seja:

(33) $$\begin{array}{cc}{\tilde{u}}_{\mu }\left(y^{+},t\right)={\tilde{u}}_{\mu }(y^{-},t)& \forall \left\{y^{+},y^{-}\right\}ϵ\partial {\text{Ω}}_{\mu }\end{array}$$

Além disso, são prescritas flutuações de deslocamentos nulas nos cantos. Portanto, neste modelo o espaço V_μ é adotado como:

(34) $$V_{\mu }=\left\{{\tilde{u}}_{\mu }\in {\tilde{K}}_{\mu }^{*}/{\tilde{u}}_{\mu }(y^{+},t)={\tilde{u}}_{\mu }(y^{-},t)\forall \left\{y^{+},y^{-}\right\}\in \partial {\text{Ω}}_{\mu }\right\}$$

A fim de satisfazer o princípio de Hill-Mandel as forças t^e devem ser anti-periódicas no contorno, isto é: t^e(y⁺,t) = −t^e(y⁻,t) e as forças de volume b(y,t) devem ser nulas no domínio do EVR. Neste modelo o sistema definido em (28) pode ser escrito da seguinte forma:

(35) $${\left\{\begin{array}{c}{F}_p\\ F_m\\ F_i\end{array}\right\}}^i+{\left[\begin{array}{ccc}{K}_{pp}& K_{pm}& K_{pi}\\ K_{mp}& K_{mm}& K_{mi}\\ K_{ip}& K_{im}& K_{ii}\end{array}\right]}^i{\left\{\begin{array}{c}\delta {\tilde{u}}_p\\ \delta {\tilde{u}}_m\\ \delta {\tilde{u}}_i\end{array}\right\}}_{\mu }^{i+1}=0$$

onde os índices p, m e I definem, respectivamente, pontos y⁺, y⁻ e internos.

Considerando que $\delta {\tilde{u}}_m=\delta {\tilde{u}}_p$, o sistema (35) pode ser reduzido e representado da seguinte forma:

(36) $${\left\{\begin{array}{c}{F}_p+F_m\\ F_i\end{array}\right\}}^i+{\left[\begin{array}{cc}{K}_{pp}+K_{pm}+K_{mp}+K_{mm}& K_{pi}+K_{mi}\\ K_{ip}+K_{im}& K_{ii}\end{array}\right]}^i{\left\{\begin{array}{c}\delta {\tilde{u}}_p\\ \delta {\tilde{u}}_i\end{array}\right\}}_{\mu }^{i+1}=0$$

isto é, nesse caso apenas as flutuações de deslocamentos dos pontos p e I são incógnitas no sistema de equações necessário resolver no processo iterativo para alcançar o equilíbrio do EVR.

Note que a partir da equação (36) podem ser obtidos os vetores e matrizes definidos na equação (32).

2.2.3 Tensões e Tensor Constitutivo Homogeneizados

A tensão homogeneizada é calculada pela equação (18) e, considerando-se que o EVR é composto de vazios e parte sólida (${\text{Ω}}_{\mu }={\text{Ω}}_{\mu }^v\cup {\text{Ω}}_{\mu }^s$), resulta em:

(37) $$\sigma =\sigma \left(x,t\right)=\frac{1}{V_{\mu }}{\int }_{{\text{Ω}}_{\mu }^S}\sigma \mu \left(y,t\right)dV+\frac{1}{V_{\mu }}{\int }_{{\text{Ω}}_{\mu }^v}\sigma \_\mu \left(y,t\right)dV$$

Considerando-se o teorema de Green e discretizando o EVR em elementos finitos, a equação (37) pode ser escrita da seguinte forma num passo de tempo Δt_n = t_n+1 – t_n relacionado a um determinado incremento de carga:

(38) $${\sigma }_{n+1}=\frac{1}{V_{\mu }}\left[{\int }_{\partial {\Omega }_{\mu }^h}{t}_{n+1}^e{\otimes }_{s}ydA-{\int }_{{\Omega }_{\mu }^{S(h)}}{b}_{n+1}{\otimes }_{s}ydV\right]$$

onde faz-se uso da seguinte expressão $u{\otimes }_{s}v=\frac{1}{2}\left(u\otimes v+v\otimes u\right)$, no cálculo das integrais, válida para quaisquer vetores u e v, sendo a expressão um produto tensorial resultando num tensor simétrico; $t_{n+1}^e$ são as forças ao longo do contorno externo (definidas na equação (25)) e sendo y o vetor das coordenadas de um ponto genérico do EVR.

Note que o segundo membro da equação (38) é nulo se forem adotadas flutuações periódicas ao longo do contorno do EVR, pois neste modelo tem-se b(y,t) = 0 sobre o domínio Ω_μ.

Por outro lado, considerando-se as equações (21) e (22), após a discretização do EVR em elementos, para uma iteração i de um incremento n of time (sendo Δt_n = t_n+1 – t_n), o tensor constitutivo tangente macroscópico ou homogeneizado C^ep é obtido a partir da equação (20), a qual pode ser escrita como: (ver mais detalhes em [¹⁵]):

(39) $$C_i^{ep}=\frac{\partial {\sigma }^i}{\partial {\epsilon }^i}=\frac{\frac{1}{V_{\mu }}{{\displaystyle \int }}_{{\text{Ω}}_{\mu }^h}\partial f_y\left({\epsilon }_{n+1}+{\nabla }^S{\tilde{u}}_{\mu }^i\right)dV}{\partial {\epsilon }^i}$$

onde ${\nabla }^S{\tilde{u}}_{\mu }^i$ é obtido com a solução do problema de equilíbrio do EVR definido em (27).

Por sua vez, a equação (39) pode ser dividida da seguinte maneira:

(40) $$C_i^{ep}=C_i^{ep(Taylor)}+{\tilde{C}}_i^{ep}$$

onde $C_i^{ep(Taylor)}$ é definido como o operador tangente do modelo de Taylor (obtido adotando-se ${\nabla }^S{\tilde{u}}_{\mu }^i=0$), sendo calculado pela média volumétrica do tensor constitutivo microscópico como segue:

(41) $$C_i^{ep\left(Taylor\right)}=\frac{1}{V_{\mu }}{\int }_{{\text{Ω}}_{\mu }}{\frac{\text{d}{\sigma }_{\mu }}{d{\epsilon }_{\mu }}|}_{i}dV=\frac{1}{V_{\mu }}{\int }_{{\text{Ω}}_{\mu }}{D}_{\mu }^{i}dV={\displaystyle {\sum }_{p=1}^{N_p}\frac{V_p}{V_{\mu }}{D}_{\mu }^{p\left(i\right)}}$$

onde D_μ é o tensor tangente microscópico e N_p o número de fases definidas no EVR.

A outra parte ${\tilde{C}}_i^{ep}$ da equação (40) representa a influência da flutuação dos deslocamentos no valor do tensor homogeneizado (ver mais detalhes em [¹⁵]). Considerando-se um discretização em N_e elementos, a mesma é definida por:

(42) $${\tilde{C}}_i^{ep}=\frac{\frac{1}{V_{\mu }}{\int }_{{\Omega }_{\mu }^h}\partial g_y({\nabla }^s{\tilde{u}}_{{\mu }_i})dV}{\partial {\epsilon }^i}=-\frac{1}{V_{\mu }}{G}_R^i{K}_R^{i^{-1}}{G}_R^{i^T}$$

ondeG_R é definido de acordo com o modelo em multi-escala; para flutuações periódicas é dado por: $G_R^i={\left[G_p+G_m{G}_I\right]}^i$, sendo G definida como: $G={\displaystyle \sum _{e=1}^{Ne}{D}_{\mu }^e{B}_e{V}_e}$.

Note que a parte sólida do EVR pode ser composta por diferentes fases, cada uma podendo ter diferentes propriedades elásticas e serem regidas por diferentes modelos constitutivos. Assim, no caso de se ter um EVR uniforme, ou seja, composto por apenas uma fase, o campo de flutuações do mesmo será nulo e, portanto, a resposta em multi-escala será igual àquela obtida através de uma análise não-linear convencional.

2.3 ALGORITMO

É importante notar que para calcular os vetores elásticos referentes ao problema de placa (como, por exemplo, os momentos elásticos ΔK_n definidos na equação (9), necessita-se de um módulo de Young E e um coeficiente de Poisson ν. Como a análise é feita em multi-escala, onde a microestrutura possui fases com diferentes propriedades elásticas, a fim de melhorar a taxa de convergência do processo iterativo do macro-contínuo (placa), esses valores (E eν da placa) serão adotados como a média volumétrica dos respectivos valores na micro-estrutura, ou seja:

(43) $$E={\displaystyle {\sum }_{p=1}^{N_p}\frac{V_p}{V_{\mu }}{E}_p}$$

(44) $$v={\displaystyle {\sum }_{p=1}^{N_p}\frac{V_p}{V_{\mu }}{v}_p}$$

Observe que as equações (43) e (44) não são usadas para resolver o EVR, apenas para calcular os vetores elásticos da placa. No algoritmo a seguir, i (sendo i ≥ 0) representa uma iteração do macro-contínuo enquanto i_EVR(sendo i_EVR ≥ 0)está relacionado com o processo iterativo do EVR. O processo incremental-iterativo necessário para alcançar o equilíbrio da placa em um incremento n é:

1. Cálculo do incremento de momentos elásticos ΔK_n (definida na equação (9)).

2. Para cada nó da placa:

   • 2.1

     Incremento de deformação $\Delta {\epsilon }_n^0$ é obtido com a lei de Hooke (equação (3)).

   • 2.2

     Um EVR é definido para cada nó da placa e resolvido seguindo o procedimento:

     • 2.2.1

       As macro-deformações $\text{Δ}{\epsilon }_{\text{n}}^{\text{i}}$ são impostas ao EVR .

     • 2.2.2

       Para i_EVR = 0, as flutuações de deslocamentos são nulas e, portanto, os deslocamentos no EVR (equação (21)) são: ${\left\{\text{Δ}{u}_{\mu }^t\right\}}_n^0={\left[\text{Δ}\epsilon \right]}_n^i{\left\{\begin{array}{c}{x}_1\\ x_2\end{array}\right\}}_y$.

     • 2.2.3

       Solução do problema de equilíbrio do EVR (processo iterativo i_EVR ≥ 0:

       • a)

         Valores nodais das forças ${\left\{\text{F}\right\}}_{\text{n}+1}^{{\text{i}}_{\text{EVR}}}$ (equação (29)). Para isso, os seguintes valores devem ser computados para cada elemento finito e:

         • –

           vetor do incremento de deformação de tentativa: ${\left\{\text{Δ}{\epsilon }_{\mu }^{\text{e}\left(\text{t}\right)}\right\}}_{\text{n}}^{{\text{i}}_{\text{EVR}}}={\text{B}}_{\text{e}}{\left\{{\text{Δu}}_{\mu }^{\text{e}\left(\text{t}\right)}\right\}}_{\text{n}}^{{\text{i}}_{\text{EVR}}}$,

         • –

           vetor do incremento de tensão de tentativa : ${\left\{\text{Δ}{\sigma }_{\mu }^{\text{e}\left(\text{t}\right)}\right\}}_{\text{n}}^{{\text{i}}_{\text{EVR}}}$, usando lei de Hooke,

         • –

           vetor de tensão de tentativa: ${\left\{{\sigma }_{\mu }^{e\left(t\right)}\right\}}_{n+1}^{i_{EVR}}={\sigma }_n^e+{\left\{\text{Δ}{\sigma }_{\mu }^{\text{e}\left(\text{t}\right)}\right\}}_{\text{n}}^{{\text{i}}_{\text{EVR}}}$.

         • –

           verifica o modelo constitutivo, obtendo a tensão ${\left\{{\sigma }_{\mu }^e\right\}}_{n+1}^{i_{EVR}}$ e o tensor constitutivo ${\left[D_{\mu }^e\right]}_n^{i_{EVR}}$ (todas as variáveis são atualizadas tendo como base os valores do incremento anterior (ver em [²⁵])

       • b)

         O vetor ${\left\{F_R\right\}}_{n+1}^{i_{EVR}}$ (equação (32)) é obtido de acordo com o modelo em multi-escala.

       • c)

         Verifica-se a convergência: $\frac{\sqrt{F_R{}^T{F}_R}}{\sqrt{F_P{}^T{F}_P}}\le tol$, onde F_p são as forças nodais no lado ${\text{Γ}}_i^{+}$ (ver figura (2)) e tol é a tolerância adotada para o processo iterativo. Se há convergência, segue no passo (2.2.4); se não há, o processo continua no passo d.

       • d)

         Obtém a matriz de rigidez do EVR (equação (30)) e a matriz ${\left[K_R\right]}_n^{i_{EVR}}$ (equação (32)).

       • e)

         Resolve o sistema dado pela equação (32) e obtém as correções de flutuações ${\left\{\delta {\tilde{u}}_R\right\}}_n^{i_{EVR}+1}$

       • f)

         Obtém o vetor total ${\left\{\delta \tilde{u}\right\}}_n^{i_{EVR}+1}$.

       • g)

         Calcula o novo incremento de deslocamento a ser imposto ao EVR: ${\left\{\text{Δ}{u}_{\mu }^t\right\}}_n^{i_{EVR}+1}={\left[\text{Δ}\epsilon \right]}_n^i\left\{\begin{array}{c}{x}_1\\ x_2\end{array}\right\}+{\left\{{\Delta \tilde{\text{u}}}_{\mu }\right\}}_n^{i_{EVR}+1}$, onde ${\left\{\text{Δ}{\tilde{u}}_{\mu }\right\}}_n^{i_{EVR}+1}={\left\{\text{Δ}{\tilde{u}}_{\mu }\right\}}_n^{i_{EVR}}+{\left\{\delta \tilde{u}\right\}}_n^{i_{EVR}+1}$ e retorna ao passo (a).

     • 2.2.4

       Obtém o tensor constitutivo homogeneizado ${\left[C^{ep}\right]}_n^i$ (equação (40))

     • 2.2.5

       Obtém a tensão homogeneizada ${\sigma }_{n+1}^{i+1}$ (equação (38)), onde o segundo termo é nulo para flutuações periódicas; calcula incremento de tensão $\text{Δ}{\sigma }_n^i=\text{Δ}{\sigma }_{n+1}^{i+1}-{\sigma }_{n-1}$ do macro.

     • 2.2.6

       O incremento de força de membrana $\Delta N_n^i=t\Delta {\sigma }_n^i$ é obtido.

3. Verifica-se o equilíbrio pela equação (13) do macro-contínuo.

   • 3.1

     Se o resíduo R_N não for nulo, a matriz tangente (equação (16)) é atualizada e as correções $\delta \text{Δ}{\epsilon }_n^{i+1}$ (equação (15)) calculadas. Atualiza-se o incremento de deformações $\delta \text{Δ}{\epsilon }_n^{i+1}$ (equação (14)) a ser aplicado na próxima iteração e retorna ao passo 2.2 para começar a nova iteração i+1.

   • 3.2

     Se o critério de convergência do macro é nulo, de acordo com a tolerância adotada, continua no passo 4.

4. Calculam-se os valores nodais das forças de membrana plásticas (equação (11)) e os valores nodais de deslocamentos, forças e esforços internos (equações (8) e (10)).

5. Retorna ao passo 1 e começa novo incremento n+1.

3. RESULTADOS E DISCUSSÃO

Com o intuito de verificar as potencialidades da formulação proposta para a análise chapas compostas por materiais heterogêneos, onde a microestrutura influencia o comportamento macromecânico, dois tipos de materiais ilustram as aplicações deste trabalho. Na primeira análise, são considerados materiais homogêneos, onde uma modelagem constitutiva convencional é empregada. Logo em seguida, inclusões com propriedades elásticas mais rígidas são incorporadas à matriz, onde a formulação apresentada neste trabalho é empregada. Na segunda análise, materiais que apresentam microestruturas porosas são utilizados, onde procura-se evidenciar a capacidade da formulação proposta em recuperar o comportamento macromecânico da estrutura influenciada pela modificação da distribuição de vazios e aumento da porosidade.

Nos exemplos numéricos descritos a seguir, considerou-se como macroestrutura, a chapa de alumínio tracionada lateralmente de dimensão (60cm x 30cm) e 1 cm de espessura, como mostrado na figura (3a). Para a discretização da placa, utilizou-se apenas um quarto da mesma, devido à simetria, sendo definidos 16 elementos e 36 nós no contorno com 24 células no domínio (ver figura (3b)). Como condição de contorno, foi definida uma tração p_n = 121,5N/mm aos nós 19 a 23; o lado correspondente aos nós 6 a 18 foi considerado livre enquanto prescreveu-se p_s = 0 e u_n = 0 para os lados correspondentes aos nós 1 a 5 e 24 a 36, sendo n e s as direções normal e tangencial, respectivamente. Ressalta-se que os valores dos deslocamentos não têm mudanças significativas quando se considera uma malha mais refinada para o macrocontínuo, ou seja, a convergência dos resultados foi obtida com a malha definida na figura (3b).

Figura 3
a) Geometria da chapa b) Discretização da chapa

Para essa placa foram feitos diferentes tipos de análises, a fim de demonstrar a diferença de resultados obtidos, de acordo com o tipo de análise. Considerando-se uma análise não-linear convencional, que equivale a adotar um EVR uniforme (com apenas uma fase), foram considerados dois diferentes modelos constitutivos para modelar o comportamento do material: Von Mises e Mohr-Coulomb, apesar de que por se tratar de uma placa de alumínio, o critério de Von Mises é mais indicado, pois o Mohr-Coulomb é mais apropriado para representar o comportamento de materiais frágeis. O objetivo de também se adotar o Mohr-Coulomb neste exemplo, é apenas mostrar que os resultados mudam, de acordo com o critério adotado, o que ressalta a importância de se adotar um critério apropriado para o material. Nessas análises foram adotadas as seguintes propriedades para o material: coeficiente de Poisson v = 0,2; módulo de Young E = 70GPa; considerando o modelo de Von Mises com strain hardening, adotou-se a tensão de escoamento σ_y = 243 MPa e coeficiente da curva de encruamento K = 0,032E. Por outro lado, para o Mohr Coulomb, além dos valores de tensão de escoamento e o coeficiente da curva de encruamento serem os mesmos, foi definido o ângulo de dilatação de 20° e ângulo de fricção de 10°.

Considerando-se a análise em multi-escala, a fim de aumentar a resistência e rigidez do material, foram definidas várias inclusões elásticas dentro da matriz (veja figura 4) acarretando numa heterogeneidade do material. As propriedades da matriz são iguais àquelas definidas para o material na análise não-linear convencional. Para governar o comportamento não-linear do material da matriz também foram considerados os critérios de Von Mises e Mohr-Coulomb. Na discretização da microestrutura (EVR) foram definidos 520 elementos finitos traingulares e 293 nós como se pode ver na figura (4). As propriedades elásticas adotadas para as inclusões foram as seguintes: E = 200 GPa e v = 0,20. Como condição de contorno no EVR, adotaram-se flutuações de deslocamentos periódicas, além de considerar uma tolerância tol = 1 x 10⁻⁵ ao longo dos processos iterativos do EVR e da placa. Como já dito anteriormente, cada EVR representa um ponto do macro-contínuo e, portanto, suas dimensões não são importantes, porém deve-se ter a informação sobre a proporcionalidade entre as fases constituintes, assim como a localização de cada constituinte.

Figura 4
Discretização do EVR

O deslocamento do nó 21 na direção X₂ (ver Figura (3b)) em função do fator de carga β para as diferentes análises pode ser visto na figura (5). Em primeiro lugar, considerando o material homogêneo (análises não-lineares convencionais), nota-se a diferença entre as respostas a partir do momento que ocorre o processo de plastificação do material, evidenciando a importância da adoção do modelo constitutivo próprio para cada tipo de material. O modelo de Von Mises evidencia um processo de plastificação mais pronunciado, prórpio de materiais dúcteis.

Figura 5
Deslocamento do ponto 21 ao longo do processo incremental, considerando-se diferentes modelos constitutivos.

Por outro lado, considerando a análise da chapa constituída por materiais heterogêneos (análise multi-escala), inicialmente percebe-se o ganho de rigidez e resistência no comportamento macromecânico da chapa com a introdução de inclusões mais rígidas na microestrutura do material. Observa-se que o nível de carregamento que leva ao início do processo de plastificação no macrocontínuo não foi alterado no caso da consideração do modelo de Von Mises. Ou seja, a matriz influencia de forma proeminente a plastificação do material. Já no caso da matriz modelada com o modelo de Mohr-Coulomb, a resposta macromecânica apresentou um início mais tardio do processo de plastificação. Por outro lado, o ganho de resistência e rigidez observado com a inclusão de material mais rígido na microestrutura também se deve a admissão de perfeita aderência entre inclusão e matriz. Contudo, no trabalho desenvolvido em [⁴], mostrou-se que com o decorrer do processo de carregamento, é adicionado o processo de descolamento da interface matriz/inclusão ao processo de plastificação da matriz, o que gera uma perda de resistência e rigidez da microestrutura do material. Contudo, o reflexo desse processo de danificação da microestrutura na macroestrutura será abordado em futuros trabalhos.

O segundo exemplo de aplicação refere-se a materiais com microestrutura dúctil porosa. Para tanto, são realizadas análises multi-escala na chapa definida na figura (3) considerando-se vazios dentro da matriz do EVR, onde diferentes frações volumétricas de vazios concentradas no centro de cada EVR ou vazios distribuídos no interior do EVR são adotadas. A matriz é considerada com comportamento governado pelo modelo de Von Mises, cujos parâmetros são os mesmos do exemplo anterior.

Inicialmente investiga-se a influência de se considerar vazios distribuídos (figura (6)) ou, mantendo a mesma proporção de fração volumétrica, no caso vf=37%, vazios concentrados no centro do EVR (figura (7)). Para a modelagem do EVR com cinco vazios distribuídos foram adotados 344 elementos finitos triangulares, e no caso do EVR com a concentração dos vazios em seu centro, 488 elementos finitos triangulares foram utilizados.

Figura 6
Discretização do EVR com cinco vazios – vf = 37%.

Figura 7
Discretização do EVR com um vazio concentrado, mantendo a mesma proporção de vazios do EVR da Figura 6 – vf = 37%.

A figura (8) ilustra o gráfico fator de carga x deslocamento do ponto 21 da chapa (figura (3)), onde foram utilizados os EVRs das figuras (6) e (7). Observa-se que com a mesma fração volumétrica de vazios na microestrutura, a resposta macromecânica é alterada quando se considera os vazios concentrados ou distribuídos no EVR. A consideração de vazios distribuídos na microestrutura leva a um comportamento mais flexível da placa.

Figura 8
Deslocamento do ponto 21 ao longo do processo incremental, considerando-se EVRs das figuras 6 e 7.

Em seguida, a fim de observar como a rigidez do macrocontínuo muda em função da fração volumétrica de vazios definida na microestrutura, EVRs com diferentes frações volumétricas de vazios são considerados, vf= 5, 10 e 30%, cujas discretizações estão mostradas nas figuras (9), (10) e (11), onde 688, 508 e 404 elementos finitos triangulares foram empregados, respectivamente.

Figura 9
Discretização do EVR – vf= 5%

Figura 10
Discretização do EVR – vf=10%

Figura 11
Discretização do EVR – vf=30%

A figura (12) ilustra a relação entre o fator de carga aplicada na chapa e o deslocamento do ponto 21. Observa-se que o aumento da fração volumétrica de vazios na microestrutura leva a uma diminuição do valor do deslocamento e da tensão onde ocorre o início do processo de plastificação no macrocontínuo. Obviamente que a presença de uma porosidade maior no maior configura-se numa perda de rigidez e de resistência do macrocontínuo.

Figura 12
Deslocamento do ponto 21 do macrocontínuo ao longo do processo incremental, considerando-se EVRs com vf = 5, 10, 30 e 37%.

Por fim, pode-se afirmar que as respostas apresentadas neste trabalho são coerentes com respostas esperadas e relatadas em outros trabalhos, vide, por exemplo, [¹], [²], [⁸], [⁹], [¹⁴] e [²⁶]. Portanto, os resultados apresentados evidenciam que o desenvolvimento da modelagem proposta é promissor.

4. CONCLUSÕES

No presente trabalho, uma formulação para a análise de chapas utilizando o Método dos Elementos de Contorno foi acoplada com o Método dos Elementos Finitos com o objetivo de realizar análise multi-escala. Uma formulação não-linear do MEC foi usado para modelar a chapa, aqui considerada como o macro-contínuo, enquanto que a microestrutura foi modelada pelo MEF escrito em termos de flutuação de deslocamentos. O Princípio da Macro-homogeneidade de Hill-Mandel, assim como a hipótese de média volumétrica dos tensores de tensão e de deformação, foi utilizado para fazer a transição do micro para o macro. Condições de flutuação de deslocamento periódicas foram impostas ao EVR, as quais apresentavam forte heterogeneidade devido à presença de inclusões rígidas com diferentes propriedades elásticas da matriz ou presença de vazios. Apesar da forte heterogeneidade do meio, a solução numérica apresentou estabilidade e taxa quadrática de convergência, atingindo o equilíbrio com um pequeno número de iterações.

O acoplamento MEC/MEF mostrou-se como uma boa alternativa às formulações multi-escala que utilizam do MEF para simular o comportamento mecânico das macro e microestruturas, como a maioria dos trabalhos sobre o tema. A principal vantagem da utilização do MEC reside no fato do mesmo ser melhor aplicado a problemas que apresentam singularidade no macro-contínuo devido à processos de localização de deformação e processos de fraturamento na microestrutura. Este tipo de problema será abordado em trabalhos futuros, onde a falha do macro-contínuo pode ocorrer por conta do descolamento da interface inclusão/matriz ou mesmo pela a propagação de fissuras na matriz em materiais frágeis, como o concreto.

As respostas qualitativas apresentadas neste trabalho mostram que a modelagem proposta é uma ferramenta própria para a identificação do comportamento macromecânico de chapas compostas por materiais heterogêneos. Pode-se afirmar ainda que a influência das heterogeneidades e vazios da microestrutura na macroestrutura foi muito bem representada pela formulação proposta. Senso assim acredita-se que as verificações realizadas neste trabalho, de modo a colher de forma satisfatória os principais fenômenos dissipativos envolvidos, são um passo inicial importante antes de aplicações e comparações com respostas experimentais em estruturas compostas por materiais heterogêneos.

5. AGRADECIMENTOS

Ao Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq, pelo suporte financeiro fornecido durante a realização do trabalho.

6. BIBLIOGRAFIA

• ^[1]
  AZIZI, R., “Micromechanical modeling of damage in periodic composites using strain gradient plasticity”, Engineering Fracture Mechanics, v. 92, pp. 101–113, 2012.
• ^[2]
  GIUSTI, S. M., BLANCO, P. J., SOUZA NETO, E. A., et al.,“An assessment of the Gurson yield criterion by a computational multi-scale approach”, Engineering Computation, v. 26, pp. 281–301, 2009.
• ^[3]
  CAVANCANTE, M. A. A., KHATAM, H., PINDERA, M., “Homogenization of elastic-plastic periodic materials by FVDAM and FEM approaches - an assessment”, Composites Part B-Engineering, v. 42, pp. 1713–1730, 2011.
• ^[4]
  PITUBA, J. J. C., FERNANDES, G. R., SOUZA NETO, E. A., “Modeling of cohesive fracture and plasticity processes in composite microstructures”, Journal of Engineering Mechanics-ASCE, publicado online, DOI: 10.1061/(ASCE)EM.1943-7889.0001123, 2016.
  » https://doi.org/10.1061/(ASCE)EM.1943-7889.0001123
• ^[5]
  PITUBA, J. J. C., SOUZA NETO, E. A., “Modeling of unilateral effect in brittle materials by a mesoscopic scale approach”, Computers and Concrete, v. 15, n. 5, pp. 735–758, 2015.
• ^[6]
  GHOSH, S., LEE, K., MOORTHY, S., “Two scale analysis of heterogeneous elastic–plastic materials with asymptotic homogenisation and Voronoi cell finite element model”, Computer Methods in Applied Mechanics and Engineering, v. 132, pp. 63–116, 1996.
• ^[7]
  TERADA, k., HORI, M., KYOYA, T., KIKUCHI, N., “Simulation of the multiscale convergence in computational homogenization approaches”, International Journal of Solids and Structures, v..37, pp. 2285–2311, 2000.
• ^[8]
  KOUZNETSOVA, V., GEERS, M. G. D., BREKELMANS, W. A. M., “Multi-scale constitutive modelling of heterogeneous materials with a gradient enhanced computational homogenization scheme”, International Journal for Numerical Methods in Engineering, v. 54, pp. 1235–1260, 2002.
• ^[9]
  KOUZNETSOVA, V., GEERS, M. G. D., BREKELMANS, W. A. M., “Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy”, Computer Methods in Applied Mechanics and Engineering, v. 193, pp. 5525–5550, 2004.
• ^[10]
  NEMAT-NASSER, S., HORI, M., Micromechanics: Overall Properties of Heterogeneous Materials, Elsevier Science, Amsterdam, 1999.
• ^[11]
  PITUBA, J. J. C., FERNANDES, G. R. “An anisotropic damage model for concrete”, Journal of Engineering Mechanics-ASCE, v. 137, pp. 610–624, 2011.
• ^[12]
  FERNANDES, G. R., PITUBA, J. J. C., SOUZA NETO, E. A. “FEM/BEM formulation for multi-scale analysis of stretched plates”,Engineering Analysis with Boundary Elements,. v. 54, pp.47–59, 2015.
• ^[13]
  FERNANDES, G. R., PITUBA, J. J. C., SOUZA NETO, E. A., “Multi-Scale Modelling For Bending Analysis of Heteregeneous Plates by Coupling BEM AND FEM”, Engineering Analysis with Boundary Elements, v. 51, pp.1–13, 2015.
• ^[14]
  PERIC, D., SOUZA NETO, E. A., FEIJÓO, R., et al., “On Micro-to-Macro Transitions for Multiscale Analysis of Heterogeneous Materials: Unified Variational Basis and Finite Element Implementation”, International Journal for Numerical Methods in Engineering, v. 87, pp. 149–170, 2011.
• ^[15]
  SOUZA NETO, E. A., FEIJÓO, R. A., Variational foundations of multi-scale constitutive models of solid: Small and large strain kinematical formulation National Laboratory for Scientific Computing (LNCC/MCT), Brazil, Internal Research & Development Report No. 16, 2006.
• ^[16]
  SOMER, D. D., SOUZA NETO, E. A., DETTMER, W. G., et al., “A sub-stepping scheme for multi-scale analysis of solids”, Computer Methods in Applied Mechanics and Engineering, v. 198, pp. 1006–1016, 2009.
• ^[17]
  WATANABE, I., TERADA, K., SOUZA NETO, E. A., et al., “Characterization of macroscopic tensile strength of polycrystalline metals with two-scale finite element analysis”, Journal of the Mechanics and Physics of Solids, v. 56, pp. 1105–1125, 2008.
• ^[18]
  FERNANDES, G. R., DENIPOTTI, G. J., KONDA, D. H., “A BEM formulation for analysing the coupled stretching-bending problem of plates reinforced by rectangular beams with columns defined in the domain”, Computational Mechanics, v.45, pp.523 - 539, 2010.
• ^[19]
  FERNANDES, G. R., KONDA, D. H., “A BEM formulation based on Reissner’s hypothesis for analysing the coupled stretching-bending problem of building floor structures”, Engineering Analysis with Boundary Elements, v.36, pp.1377 - 1388, 2012.
• ^[20]
  FERNANDES, G. R., VENTURINI, W. S., “Non-linear boundary element analysis of floor slabs reinforced with rectangular beams”, Engineering Analysis with Boundary Elements, v.31, pp.721 - 737, 2007.
• ^[21]
  FERNANDES, G. R., VENTURINI, W. S., “Building floor analysis by the Boundary element method”, Computational Mechanics, v.35, pp.277 - 291, 2005.
• ^[22]
  FERNANDES, G. R., “A BEM formulation for linear bending analysis of plates reinforced by beams considering different materials”, Engineering Analysis with Boundary Elements, v.33, pp.1132 - 1140, 2009.
• ^[23]
  FERNANDES, G. R., KONDA, D. H., “A BEM formulation based on Reissner’s theory to perform simple bending analysis of plates reinforced by rectangular beams”, Computational Mechanics, v.42, pp.671 - 683, 2008.
• ^[24]
  FERNANDES, G. R., VENTURINI, W. S. “Non-Linear Boundary Element Analysis of Plates Applied to Concrete Slabs”, Engineering Analysis with Boundary Elements, v.26, pp.169 - 181, 2002.
• ^[25]
  FERNANDES, G. R., SOUZA NETO, E. A., “Self-consistent linearization of non-linear BEM formulations with quadratic convergence”, Computational Mechanics, v.52, pp.1125–1139, 2013.
• ^[26]
  SANTOS, W. F., FERNANDES, G. R., PITUBA, J. J. C., “Análise da influência dos processos de plasticidade e fratura no comportamento mecânico de microestruturas de Compósitos de Matriz Metálica”, Revista Matéria, v. 21, pp. 577–598, 2016.

Datas de Publicação

• Publicação nesta coleção
  2017

Histórico

• Recebido
  20 Jul 2016
• Aceito
  05 Dez 2016