Análise em multi-escala do problema bidimensional de placas submetidas ao cisalhamento e considerando descolamento de fases

Multi-scale analysis of the two-dimensional problem of plates submitted to shear loadings and considering phase debonding

RESUMO

O presente trabalho trata da simulação do comportamento mecânico de microestruturas dúcteis reforçadas por inclusões elásticas e a repercussão de seu comportamento na resposta macromecânica de estruturas submetidas a carregamentos predominantes de cisalhamento. Para tanto, uma modelagem baseada em conceitos de homogeneização computacional utilizando um acoplamento direto entre escalas e do MEC/MEF é utilizado. Também, o fenômeno de descolamento de fase é levado em conta com o uso de modelo de fratura e contato onde a abertura da fratura por escorregamento é considerada. Por outro lado, a ductilidade da matriz metálica é representada pelo modelo de Von Mises. As análises numéricas são realizadas em duas etapas: uma primeira, em nível de microestrutura e, em seguida, é realizada a análise de uma chapa através do acoplamento direto entre as escalas. Contata-se coerência dos resultados encontrados nas micro e macro escalas com a utilização da modelagem proposta.

Palavras-chave
Modelagem multi-escala; homogeneização; elementos de contorno; descolamento de fase

ABSTRACT

The present work deals with numerical analyses on mechanical behavior of ductile microstructures reinforced by elastic inclusions and the consequences upon micromechanical behavior of structures submitted to shear loadings. For this proposal, a proposed modeling based on computational homogenization technique using full coupling between scales and BEM/FEM is used. Also, the phase debonding is taken into account by using cohesive fracture and contact models where the opening sliding phenomenon is considered. On the other hand, the ductility of the metallic matrix is represented by Von Mises model. Initially, the numerical analyses are performed only at microstructural level and then a full coupled multi-scale analysis of a plate is considered, showing coherence between the micro and macro scale results.

Keywords
Multi-scale modeling; homogenization; boundary elements; phase debonding

1. INTRODUÇÃO

A modelagem do comportamento mecânico de materiais ainda é uma tarefa desafiadora principalmente quando lhe damos com materiais cujo comportamento mecânico é fortemente influenciado pela heterogeneidade da microestrutura. Atualmente o desenvolvimento de novos modelos constitutivos auxiliando o entendimento do comportamento de novos materiais é um tema de interesse. A obtenção desses novos materiais com a manipulação da geometria e composição dos constituintes da microestrutura visa a melhoria nas características mecânicas, acústicas e/ou térmicas do macrocontínuo, o que requer formulações teóricas sofisticadas assim como o aperfeiçoamento de técnicas numéricas. De acordo com AZIZI [²], uma classe de materiais que merece destaque é o compósito de matriz metálica reforçado por inclusões rígidas (CMM), onde propriedades como alta rigidez e resistência e baixa densidade levam a um aumento de aplicação desses materiais nos dias atuais. Por outro lado, essas mesmas melhorias nas propriedades apresentam a desvantagem de redução da ductilidade do material quando o descolamento de fase acontece em nível microestrutural.

Dentro desse contexto, com o passar dos anos, modelos constitutivos fenomenológicos foram propostos baseados nos conceitos da Mecânica do Contínuo e Termodinâmica dos Processos Irreversíveis aplicados a análises macroscópicas. Isso foi desenvolvido porque características similares do comportamento mecânico de materiais são observadas em nível macroscópico, como por exemplo: elasticidade, viscosidade, deformação plástica, ruptura frágil e ruptura dúctil. Tais características são reflexos do que ocorrem na microestrutura dos materiais [¹]. Entretanto, nas últimas décadas, formulações teóricas levando em conta fenômenos dissipativos que ocorrem na microescala e que repercutem no comportamento mecânico do macrocontínuo têm sido propostas para tratar de processos complexos em materiais heterogêneos. Surgiram, então, as chamadas teorias multi-escala utilizando modelos baseados no Método dos Elementos Finitos (MEF) para uma ou mais escalas do problema, [²-¹²] e mais recentemente métodos sem malha [¹³], teoria de volume finito [¹⁴] e método dos elementos finitos generalizados [¹⁵] vêm sendo empregados na análise de uma ou mais escalas. Já formulações baseadas no Método dos Elementos de Contorno (MEC) vêm sendo desenvolvidas, principalmente pelos autores do presente trabalho [¹⁶-²³]. Já em se tratando de modelos aplicados apenas à microestrutura utilizando o MEC, cita-se o trabalho de YANG e QIN [²⁴] que obteve propriedades mecânicas em compósitos reforçados por inclusões circulares rígidas e de RODRIGUES et al. [²⁵] que desenvolveram um modelo para a análise elástica de compósitos reforçados por fibras curtas enquanto que ARAÚJO e GRAY [²⁶] trabalharam com microestruturas reforçadas por nano tubos de carbono.

De modo a apresentar a contribuição do presente trabalho, incialmente pode-se citar que FERNANDES et al. [¹⁹] desenvolveram um modelo multi-escala para a análise de chapas com acoplamento MEC e MEF, onde o MEC foi usado para modelar a macroestrutura enquanto que o MEF foi utilizado para a microestrutura que possuía inclusões ou vazios definidos dentro da matriz para simular o comportamento de materiais heterogêneos e, no caso de considerar inclusões, uma perfeita aderência entre fases foi assumida. SOUZA JÚNIOR et al. [¹⁸] obtiveram uma formulação baseada no MEC para flexão sendo a microestrutura modelada via MEF e ainda considerando perfeita aderência na zona de interface dos compósitos analisados. A consideração do fenômeno do descolamento de fase que ocorre na interface inclusão/matriz em microestruturas passou a ser estudada por FERNANDES et al. [¹⁷] que desenvolveu uma abordagem multi-escala puramente cinemática com o objetivo de simular processos de microfissuração em nível da microescala usando o MEF, onde elementos finitos de contato e fratura coesiva desenvolvidos em PITUBA et al. [⁷] foram adicionados na formulação [¹⁹] para interfaces entre matriz e inclusão em materiais como concreto e compósitos de matriz metálica, entre outros. Dentro desse contexto, o trabalho de FERNANDES et al. [¹⁷] tratou de estudar o comportamento mecânico de materiais compósitos dúcteis considerando o descolamento de fase quando submetidos a estados de tensões normais de tração e compressão, onde observou-se que a relação entre as rigidezes dos constituintes da microestrutura afeta fortemente o comportamento de estruturas compostas por esses materiais. Contudo, por se tratar de materiais compostos essencialmente por materiais dúcteis, onde a matriz se faz mais presente, o estudo do comportamento desse compósito e de estruturas compostas por ele quando submetidas a estados de cisalhamento é algo desejável. Portanto, a principal contribuição do presente trabalho é mostrar como o descolamento de fase, que ocorre em nível da microestrutura de materiais dúcteis, como o CMM, afeta a resposta macromecânica de estruturas submetidas a estados de cisalhamento e, para isso, utiliza-se uma modelagem multi-escala.

É importante observar que o modelo baseado no MEF utilizado para modelar o comportamento mecânico da microestrutura considerando o descolamento de fase, pode ser acoplado facilmente a qualquer modelo numérico, ou seja, ou MEC ou MEF pode ser usado para modelar o macrocontínuo. Em particular, o MEC já provou ser uma ferramenta numérica adequada para tratar com diferentes problemas da Engenharia, como por exemplo, problemas de placas que são estudadas neste trabalho, sendo particularmente recomendado para o cálculo dos efeitos de cargas concentradas e em linha, assim como avaliar valores de gradientes altos. Além disso, a mesma ordem de erro é esperada quando são calculados deslocamentos e forças, pois as tensões não são obtidas por diferenciação de funções de aproximação como em outros métodos numéricos. Além disso, quando fenômenos dissipativos são considerados, apenas a aproximação de forças dissipativas é requerida sobre o domínio, isto é, as forças elásticas não são aproximadas sobre o domínio, o que é uma vantagem do MEC, onde usualmente uma malha mais refinada não é necessária para a obtenção de bons resultados. Pensando em trabalhos futuros quando for considerada a localização de deformação que ocorre em nível microestrutural podendo levar à formação de uma fratura no macrocontínuo, o MEC pode ser especialmente conveniente.

Por fim, o artigo está apresentado da seguinte forma: a formulação estendida para tratar de descolamento de fase dentro de uma abordagem multi-escala cinemática e o modelo adotado de fratura coesiva serão apresentados na seção 2, assim como a formulação do MEC não-linear desenvolvida para a análise de do problema de chapas e o algoritmo desenvolvido são brevemente descritos. Já na seção 3, alguns resultados numéricos relativos aos processos dissipativos em EVRs são discutidos de modo a evidenciar como o comportamento observado no nível da microestrutura interfere no comportamento do macrocontínuo utilizando-se de análise multi-escala. Finalmente, na seção 4 são apresentadas algumas conclusões baseadas nos resultados apresentados neste trabalho.

2. MATERIAIS E MÉTODOS

2.1 Formulação do MEC para solução do problema de equilíbrio do macrocontínuo

O macrocontínuo é representado por uma placa de contorno externo Γ e domínio Ω, onde se tem definido o problema bidimensional sendo as forças (ṗ_n e ṗ_s) e os deslocamentos (${\dot{u}}_n \mathrm{e} {\dot{u}}_s$) definidos segundo o sistema local de coordenadas (n,s), sendo n e s, respectivamente, as direções normal e tangencial ao contorno da placa. Essa formulação é discutida com maiores detalhes em [¹⁶] e [²⁷] e será aqui colocada de forma resumida. Admite-se que o material possa apresentar deformações residuais, que nesse trabalho serão assumidas como deformações plásticas ${\dot{\epsilon }}_{ij}^p$. Com isso, todas as variáveis são escritas em termos de taxas, sendo (ẋ) = dx/dt, e as deformações são divididas na parcela elástica ${\dot{\epsilon }}_{ij}^e$ e na plástica ${\dot{\epsilon }}_{ij}^p$ como segue:

$${\dot{\epsilon }}_{ij} = {\dot{\epsilon }}_{ij}^e + {\dot{\epsilon }}_{ij}^p \mathrm{i}, \mathrm{j}= 1, 2$$ (1)

onde as componentes de deformação são dadas no sistema de coordenadas cartesiano (x₁, x₂) definido no plano da placa.

Aplicando-se a lei de Hooke, as forças internas elásticas de membrana ${\dot{N}}_{ij}^e$ são relacionadas à deformação total ${\dot{\epsilon }}_{ij}$ e, portanto, podem ser escritas em termos das derivadas dos deslocamentos totais como segue:

$${\dot{N}}_{ij}^e = \mu \left[{\dot{\mathrm{u}}}_{\mathrm{i},j} + {\dot{u}}_{\mathrm{j},i}\right] + \frac{2\mu v}{1-2v\text{'}}{\dot{u}}_{k,k}{\delta }_{\mathrm{ij}} \mathrm{i}, \mathrm{j}, \mathrm{k}= 1, 2$$ (2)

onde δ_ij é o delta de Kronecker; ν' = ν/(1 + ν) 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.

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

$${\dot{N}}_{ij}^p = {\dot{N}}_{ij}^e - {\dot{N}}_{ij}^{\mathrm{⬚}}$$ (3)

onde as forças ${\dot{N}}_{ij}^{\mathrm{⬚}}$ são calculadas a partir das tensões δ_ij, que numa análise não-linear convencional são calculadas de acordo com o modelo constitutivo adotado, mas numa abordagem em multi-escala, essas tensões são calculadas após a solução do EVR, fazendo-se a homogeneização do campo de tensão no EVR.

A equação integral do deslocamento para um ponto interno é obtida a partir do Teorema de Reciprocidade de Betti (ver mais detalhes em [¹⁹]), a qual após se fazer a integração por partes, resulta na seguinte representação integral:

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

Onde, i representa a direção da carga fundamental, $u_n^{*}, u_s^{*}, p_n^{*} \mathrm{e} p_s^{*}$ são valores fundamentais das forças e deslocamentos no plano da placa; ḃ_n e ḃ_s são as forças que atuam em Ω_b; K(q) = 1 e K(q) = 1/2, respectivamente, para pontos internos e pontos sobre o contorno não coincidente com um canto.

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,ℓ pode ser obtida derivando-se a equação (4) (ver maiores detalhes em [²⁸]). Então, a partir dos valores de Δu_i,ℓ, as forças elásticas $\Delta N_{jk}^{e\left(MEC\right)}$ podem ser calculadas através da equação (2).

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. 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 [¹⁹]):

$$\Delta X = \Delta L + R_N\Delta N^P$$ (5)

onde o vetor ΔX contêm as incógnitas do contorno, 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.

Para se obter a solução não-linear, deve-se ainda escrever três equações das forças de membrana elásticas $\Delta N_{jk}^{e\left(MEC\right)}$ nos nós das células, cuja equação algébrica pode ser escrita como:

$$\Delta N^{e\left(MEC\right)} = \Delta K + S_N\Delta$$ (6)

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 nas equações (3), (5) e (6) 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 (6) 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 [¹⁶]):

$$\Delta N^{MEC} = C_N\Delta \epsilon - \Delta K_N - S_N\Delta N^P + \Delta N$$ (7)

onde os valores nodais ΔN^P são calculados localmente através das forças obtidas com o tensor constitutivo elástico (C_N) e o incremento de deformações nodais (Δε) e as forças que satisfazem o modelo constitutivo (ΔN), ou seja, são dados por:

$$\Delta N^P = C_N\Delta \epsilon - \Delta N$$ (8)

Para um incremento n, considera-se que a placa encontra-se em equilíbrio se a seguinte equação for satisfeita:

$$\Delta K_n = \Delta {N_n}^{MEC} = 0$$ (9)

onde os vetores ΔN_n e ΔN_n^MEC são definidos, respectivamente, nas equações (6) e (7). Portanto, a equação (9) de resíduos da placa pode ser definida como:

$$R_{\mathit{N}}\left(\Delta {\epsilon }_n\right)\mathit{ }\mathit{=}\mathit{ }\mathit{2}\Delta K_{\mathit{n}}\mathit{ }\mathit{-}\mathit{ }{C}_{\mathit{N}}\Delta {\epsilon }_{\mathit{n}}\mathit{ }\mathit{+}\mathit{ }{S}_{\mathit{N}}\mathit{ }(C_N\Delta {\epsilon }_n - \Delta N_n)\mathit{ }\mathit{-}\mathit{ }\Delta N_{\mathit{n}}\mathit{ }\mathit{=}\mathit{ }\mathit{0}$$ (10)

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 (10) não for satisfeita, a mesma 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 (10). Nesse processo iterativo, o incremento de tentativa $\mathrm{\Delta }{\epsilon }_n^{i+1}$ a ser aplicado na iteração (i+1) é obtido adicionado as correções $\mathrm{\delta \Delta }{\epsilon }_n^{i+1}$, ou seja, $\mathrm{\Delta }{\epsilon }_n^{i+1} = \Delta {\epsilon }_n^i + \delta \mathrm{\Delta }{\epsilon }_n^{i+1}$, sendo $\mathrm{\delta \Delta }{\epsilon }_n^{i+1}$ calculadas linearizando-se a equação (10), através da seguinte equação:

$\delta \mathrm{\Delta }{\epsilon }_n^{\mathrm{i}+1}=-{\left[\frac{\partial {\mathrm{R}}_{\mathrm{N}}\left(\mathrm{\Delta }{\epsilon }_n^i\right)}{\partial \mathrm{\Delta }{\epsilon }_{\mathit{n}}^{\mathit{i}}}\right]}^{-1}{\mathrm{R}}_{\mathrm{N}}\left(\mathrm{\Delta }{\epsilon }_{\mathit{n}}^{\mathit{i}}\right)$ (11)

onde $-\left(\partial {\mathrm{R}}_{\mathrm{N}} \left(\Delta {\mathrm{\epsilon }}_{\mathrm{n}}^{\mathrm{i}}\right)\right)/\left(\partial \Delta {\mathrm{\epsilon }}_{\mathrm{n}}^{\mathrm{i}}\right)={\mathrm{S}}_{\mathrm{N}} \left({\mathrm{C}}_{\mathrm{N}\left(\mathrm{n}\right)}^{\mathrm{ep}\left(\mathrm{i}\right)} - {\mathrm{C}}_{\mathrm{N}}\right) + {\mathrm{C}}_{\mathrm{N}} + {\mathrm{C}}_{\mathrm{N}\left(\mathrm{n}\right)}^{\mathrm{ep}\left(\mathrm{i}\right)}$ é o operador tangente consistente obtido derivando-se a equação (10); $({\mathrm{C}}_{\mathrm{N}\left(\mathrm{n}\right)}^{\mathrm{ep}\left(\mathrm{i}\right)}$ é a matriz que contém os tensores elasto-plásticos ${\left[{\mathrm{C}}_{\mathrm{N}}^{\mathrm{ep}}\right]}_{\mathrm{k}}^{\mathrm{i}}$ para todos os nós de células. Para um nó k, ${\left[{\mathrm{C}}_{\mathrm{N}}^{\mathrm{ep}}\right]}_{\mathrm{k}}^{\mathrm{i}}$ é dado por: ${\left[{\mathrm{C}}_{\mathrm{N}}^{\mathrm{ep}}\right]}_{\mathrm{k}}^{\mathrm{i}} = \frac{\partial {\mathrm{N}}_{\mathrm{k}}^{\mathrm{i}}}{\partial \Delta {\mathrm{\epsilon }}_{\mathrm{k}}^{\mathrm{i}}} = \mathrm{t}{\left[{\mathrm{C}}_{\mathrm{ep}}\right]}_{\mathrm{k}}^{\mathrm{i}}, \mathrm{sendo} {\left[{\mathrm{C}}_{\mathrm{ep}}\right]}_{\mathrm{k}}^{\mathrm{i}} = \frac{\partial \left({\mathrm{\sigma }}_{\mathrm{k}}^{\mathrm{i}}\right)}{\partial {\mathrm{\epsilon }}_{\mathrm{k}}^{\mathrm{i}}}$, isto, é, ${\left[{\mathrm{C}}_{\mathrm{ep}}\right]}_{\mathrm{k}}^{\mathrm{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 podem ser calculadas e a equação (10) verificada. Se ela não for satisfeita, as correções $\delta \Delta {\epsilon }_n^{i+1}$ (equação 11) têm que ser calculadas continuando o processo na iteração i+1.

2.2 Modelagem da micro-escala utilizando o MEF

A micro-escala representa a micro-estrutura do material e é definida pelo EVR (Elemento de Volume Representativo), cujo problema de equilíbrio é resolvido pelo Método dos Elementos Finitos. Essa formulação do MEF é discutida mais detalhadamente nos trabalhos [⁴,¹⁹-²⁰ e ²⁹], sendo neste trabalho abordados apenas seus conceitos principais. A formulação apresentada aqui é uma extensão daquela discutida em BORGES et al. [³⁰], onde além dos elementos finitos triangulares definem-se elementos retangulares nas interfaces entre matriz e agregados, os quais são regidos por uma lei de contato ou de fratura, a fim de simular o descolamento que ocorre entre matriz e agregado durante o processo de fissuração do material.

2.2.1 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 onde se deseja obter o vetor das tensões e o tensor constitutivo a fim de verificar a equação de equilíbrio da placa. Na análise em multi-escala, cada ponto x é representado pelo EVR (Elemento de Volume Representativo) (ver ⁴, ¹⁹-²⁰ e ²⁹]), onde define-se V_μ como seu volume, Ω_μ como seu domínio e ∂Ω_μ como seu contorno, sendo y um ponto qualquer do mesmo. O EVR pode ser composto por vazios (domínio ${\Omega }_{\mu }^{\nu }$) e por partes sólidas (domínio ${\Omega }_{\mu }^{\nu }$), sendo ${\Omega }_{\mu } = {\Omega }_{\mu }^{\nu } \cup {\Omega }_{\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.

Figura 1
Representação de uma chapa sujeita a um carregamento de cisalhamento e EVR composto por matriz e inclusões com possibilidade de descolamento de fase.

Note que a microestrutura do material é heterogênea quando se têm inclusões e/ou vazios definidos no domínio da matriz, ou seja, se o EVR for composto por uma única fase ele se torna homogêneo e a resposta numérica resulta igual a de uma análise não-linear convencional. O campo de deslocamentos u_µ do EVR é dividido da seguinte maneira:

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

onde a parcela εy varia linearmente, sendo obtida multiplicando-se a deformação constante ε imposta pelo macrocontínuo pelas coordenadas do ponto y do EVR; a parcela ũ_µ é denotada flutuação de deslocamentos e representa a variação de deformação no EVR, isto é, se a deformação no EVR é constante, que ocorre quando o EVR é homogêneo, tem-se ũ_µ nulo. Analogamente, a deformação microscópica pode ser escrita em termos de taxas como:

$$\dot{{\epsilon }_{\mu }}\left(y, t\right) = {\nabla }^S\dot{{\tilde{u}}_{\mu }} = \dot{\epsilon }\left(x, t\right) +\dot{\widetilde{{\epsilon }_{\mu }}}\left(y, t\right)$$ (13)

onde $\dot{\epsilon }\left(x, t\right)$ é a deformação homogênea imposta pelo macrocontínuo e $\dot{{\epsilon }_{\mu }}\left(y, t\right)$ é a taxa do campo das flutuações de deformações, sendo definida como: $\dot{{\epsilon }_{\mu }}\left(y, t\right) = {\nabla }^S\dot{{\tilde{u}}_{\mu }}$, onde ∇^s é o operador gradiente simétrico.

2.2.2 Modelagem do deslocamento entre matriz e inclusão

O modelo de fratura coesiva proposto em PITUBA et al. [⁷] será brevemente apresentado. O modelo em questão é utilizado para representar o processo de descolamento de fase que ocorre ao redor das inclusões do CMM. De acordo com PITUBA et al. [⁷], os materiais compósitos dúcteis com inclusões e/ou poros podem usualmente ser considerados como um meio elasto-plástico durante o processo de carregamento, porém em regime último de carga estes materiais podem apresentar microfissuração que governa o processo de descolamento na interface matriz/inclusão, além de propagação intensa de microfissuração na região de matriz levando a uma perda de rigidez do material que pode eventualmente entrar em colapso. No presente trabalho, o modelo de fratura coesiva proposto por PITUBA et al. [⁷] será utilizado para simular o processo de descolamento de fase que ocorre na microestrutura de compósitos, ou seja, a modelagem multi-escala proposta aqui não tem condições de ser utilizada para simular o problema de possível localização de deformação na microestrutura que pode gerar uma fratura na macroestrutura, sendo necessários alguns aperfeiçoamentos como descrito em TORO et al. [³¹].

Em PITUBA et al. [⁷], SANTOS e PITUBA [⁸] e SANTOS et al. [³²], algumas simulações numéricas foram realizadas para evidenciar o processo de colapso de microestruturas de materiais compósitos. Contudo, no presente trabalho os elementos finitos de fratura coesiva e contato são incorporados à formulação MEF/MEC multi-escala para a análise chapas [¹⁶] compostas por CMM.

Neste modelo constitutivo foi adotada uma lei coesiva em modo misto resultando em um deslocamento de abertura efetivo δ:

$$\delta = \sqrt{{{\beta }_{ef}}^2{\delta }_S^2 + {\delta }_n^2}$$ (14)

onde δ_s e δ_n são os deslocamentos de abertura normal e de escorregamento, respectivamente. Já o parâmetro β_ef assume valores que servem como ponderadores das contribuições de δ_s e δ_n. Por outro lado, assumindo que o potencial de energia livre φ depende de δ, i. e.,

$$\varphi = e{\sigma }_C{\delta }_C\left[1 - e^{\left[-\left(1+\frac{\delta }{{\delta }_C}\right)\right]}\right]$$ (15)

onde σ_c a resistência à tração na direção normal à fratura e δ_c é um deslocamento característico de abertura que indica uma abertura crítica. Já a lei coesiva é então dada por:

$${\mathit{t}}_{\mathbf{e}\mathbf{f}} = \frac{t_f}{\delta }\left({{\beta }_{ef}}^2{\mathit{\delta }}_{\mathbf{S}} + {\delta }_n\mathit{n}\right)$$ (16)

onde t_ef é o vetor de tensão coesiva ao longo da fratura, t_f é uma variável escalar representando a tensão efetiva, n representa o vetor normal à fratura e δ_s é o vetor de deslocamento de abertura por escorregamento localizado na superfície da fratura, onde neste modelo é adotado que a resistência ao escorregamento é independente da direção conferindo uma característica de isotropia à superfície da fratura.

Para os casos em que a fratura estiver em processo de abertura, a lei de carregamento é dada por:

$$t_f = \frac{\partial \varphi }{\partial \delta } = {\sigma }_C{e}^{-\frac{\delta }{{\delta }_C}} \mathrm{se} \delta = {\delta }_{max} \mathrm{e} \dot{\delta } \ge 0$$ (17)

Já para casos de descarregamento é adotado um comportamento elástico linear e a lei passa a ser dada por:

$$t_f = \frac{t_{fmax}}{{\delta }_{max}}\delta \mathrm{se} \delta < {\delta }_{max} \mathrm{ou} \dot{\delta } < 0$$ (18)

As variáveis t_fmax e δ_max são valores máximos para a tensão efetiva e deslocamento de abertura efetiva, respectivamente, cujas atualizações se dão durante o processo de carregamento. Assim, β_ef, σ_c e δ_c são parâmetros do modelo coesivo.

Por outro lado, antes do processo de nucleação da fratura, é assumida a existência de uma rigidez elevada entre as faces da fratura situadas entre os elementos finitos representativos da matriz e da inclusão, de modo que tal rigidez impeça a penetração das faces. Isso é conseguido com a inserção de um fator de penalidade (λ_p) na formulação do modelo sendo, portanto, mais um parâmetro do modelo coesivo e de contato. A lei passa a ser dada para esses casos por:

$$t_f ={\lambda }_p\delta \mathrm{se} {\lambda }_p\delta \le {\sigma }_C$$ (19)

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

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. Levando-se em conta o Princípio de macro-homogeneidade de Hill-Mandel (ver mais detalhes em [²⁹]), chega-se à equação de equilíbrio do EVR, que consiste em encontrar um campo de tensões que esteja autoequilibrado e que pode ser escrita da seguinte forma em termos de flutuação de deslocamentos:

$${\int }_{{\Omega }_{\mu }^S}{f}_y\left(\epsilon \left(x, t\right) + {\mathit{\nabla }}^s{\tilde{u}}_{\mu }\left(y, t\right)\right)dV = 0$$ (20)

Após a discretização do EVR em elementos, a equação de equilíbrio para um passo de tempo e discretização h, sendo ${\tilde{u}}_{\mu \left(n + 1\right)}^{⬚} = {\tilde{u}}_{\mu \left(n\right)}^{⬚} + \Delta {\tilde{u}}_{\mu \left(n\right)}^{⬚}$ é definida como:

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

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 \left(n + 1\right)}^{⬚} = {\tilde{u}}_{\mu \left(n\right)}^{⬚} + \Delta {\tilde{u}}_{\mu \left(n\right)}^{⬚}$ tal que a cada instante t_n+1 a equação (21) seja satisfeita. Se no incremento n+1 a equação (15) não é satisfeita, a mesma é 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ção i+1, tal que $\Delta {\tilde{u}}_{\mu }^{i + 1} = \Delta {\tilde{u}}_{\mu }^i + \delta {\widetilde{\Delta u}}_{\mu }^{i + 1}$, sendo as correções obtidas a partir da seguinte equação:

$$F^i+ K^i \delta {\widetilde{\Delta u}}_{\mu }^{i + 1} = 0$$ (22)

onde K é 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:

$$F^i= {\int }_{{\Omega }_{\mu }^h}{B}^T{f}_y\left({\epsilon }_{n+1} + B{\tilde{u}}_{\mu }^i\right)dV = \underset{e=1}{\overset{Ne}{\Sigma }} B_e^T {\sigma }_{\mu }^{e\left(i\right)}{V}_e$$ (23)

$$K^i= \left[{\int }_{{\Omega }_{\mu }^h}{B}^T{D}_{\mu }^{i}BdV\right]\underset{e=1}{\overset{Ne}{\Sigma }} B_e^T D_{\mu }^{e\left(i\right)}{B}_e{V}_e$$ (24)

onde V_e é o volume de um elemento e qualquer e é o tensor constitutivo tangente do elemento e, definido de acordo com o modelo constitutivo adotado para a fase a qual esse elemento pertence.

Finalmente, a formulação fica completa 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. Além disso, o sistema (22) pode ser reduzido de acordo com essas condições de contorno, 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 (22) pode ser escrita na seguinte forma reduzida:

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

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). Essas condições de contorno não serão tratadas aqui, para mais detalhes ver [²⁷].

Observe que 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 do macrocontínuo, a resposta de qualquer material pode ser modelada por essa condição.

2.2.4 Tensões e Tensor Constitutivo Homogeneizado

Após alcançar convergência do processo iterativo do EVR, as tensões e o tensor constitutivo referente ao macrocontínuo devem ser obtidos aplicando-se técnicas de homogeneização. Assim, 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 macrocontí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:

$$\sigma \left(x, t\right) = \frac{1}{V_{\mu }}{\int }_{{\Omega }_{\mu }}{\sigma }_{\mu }\left(y, t\right) dV$$ (26)

$$\sigma \left(x, t\right) = \frac{1}{V_{\mu }}{\int }_{{\Omega }_{\mu }}{\sigma }_{\mu }\left(y, t\right) dV$$ (27)

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

Discretizando o EVR em elementos finitos, a equação (27) pode ser escrita da seguinte forma num passo de tempo Δt_n = t_n+1 - t_n relacionado a um determinado incremento de carga (ver mais detalhes em [³⁰]):

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

onde faz-se uso da seguinte expressão u ⊗_Sν = 1/2 (u ⊗ ν + ν ⊗ u), 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 (23)) e sendo y o vetor das coordenadas de um ponto genérico do EVR.

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

$$C^{ep}\left(x, t\right) = \frac{\partial \sigma \left(x, t\right)}{\partial \epsilon \left(x, t\right)} = \frac{{\displaystyle \frac{1}{V_{\mu }}}{\int }_{{\Omega }_{\mu }}\partial {\sigma }_{\mu }\left(y, t\right)dV}{\partial \epsilon \left(x, t\right)} = \frac{{\displaystyle \frac{1}{V_{\mu }}{\int }_{{\Omega }_{\mu }}\partial f_y\left({\epsilon }_{n+1}+{\nabla }^s{\tilde{u}}_{\mu }^i\right)dV}}{\partial \epsilon \left(x, t\right)}$$ (29)

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

A equação (29) pode ser dividida da seguinte maneira:

$$C_i^{ep} = C_i^{ep\left(Taylor\right)} + {\overline{C}}_i^{ep}$$ (30)

onde $C_i^{ep\left(Taylor\right)}$ é 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:

$$C_i^{ep\left(Taylor\right)} = \frac{1}{V_{\mu }}{\int }_{{\Omega }_{\mu }} {\overline{)\frac{d{\sigma }_{\mu }}{d{\epsilon }_{\mu }}}}_i dV = \frac{1}{V_{\mu }}{\int }_{{\Omega }_{\mu }}{D}_{\mu }^{i}dV = \underset{p=1}{\overset{N_p}{\Sigma }}\frac{V_p}{V_{\mu }}{D}_{\mu }^{p\left(i\right)}$$ (31)

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 (30) 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:

$${\tilde{C}}_i^{ep} = \frac{{\displaystyle \frac{1}{V_{\mu }}}{\int }_{{\Omega }_{\mu }^h}\partial g_y\left({\nabla }^S{\tilde{u}}_{{\mu }_i}\right)dV}{\partial {\epsilon }^i} = -\frac{1}{V_{\mu }}{G}_R^i{K_R^i}^{-1}{G_R^i}^T$$ (32)

onde G_R é definido de acordo com o modelo em multi-escala (ver mais detalhes em [²⁰]) e é uma forma reduzida do tensor G definido como: $G = \underset{e=1}{\overset{Ne}{\Sigma }}{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 (6)), 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 macrocontí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:

$$E = \underset{p=1}{\overset{N_p}{\Sigma }}\frac{V_p}{V_{\mu }}{E}_p$$ (33)

$$\nu = \underset{p=1}{\overset{N_p}{\Sigma }}\frac{V_p}{V_{\mu }}{\nu }_p$$ (34)

Observe que as equações (33) e (34) 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 macrocontí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 (6)).

• 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 $\Delta {\epsilon }_n^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 são: ${\left\{\Delta u_{\mu }^t\right\}}_n^0{\left[\Delta \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) Os valores nodais das forças ${\left\{F_R\right\}}_{\mathbf{n}\mathbf{+}\mathbf{1}}^{{\mathbf{i}}_{\mathbf{EVR}}}$ (Eq. 22) são calculados somando-se a contribuição do elemento finito triangular ${\left\{F^{\mathrm{elem}}\right\}}_{\mathbf{n}\mathbf{+}\mathbf{1}}^{{\mathbf{i}}_{\mathbf{RVE}}}$ com aquela referente aos elementos de contato ${\left\{F^{\mathrm{ef}}\right\}}_{\mathbf{n}\mathbf{+}\mathbf{1}}^{{\mathbf{i}}_{\mathbf{RVE}}}$. Para isso, seguem-se os seguintes passos:

  • (a.1) Para cada elemento finito triangular elem, calculam-se:

    • vetor do incremento de deformação de tentativa: ${\left\{{\mathrm{\Delta \epsilon }}_{\mu }^{\mathrm{elem}\left(t\right)}\right\}}_{\mathbf{n}}^{{\mathbf{i}}_{\mathbf{RVE}}} = {\mathbf{B}}_{\mathbf{e}}{\mathbf{\left\{}{\mathbf{\Delta u}}_{\mathbf{\mu }}^{\mathbf{elem}\left(t\right)}\mathbf{\right\}}}_{\mathbf{n}}^{{\mathbf{i}}_{\mathbf{RVE}}}$,

    • vetor do incremento de tensão de tentativa ${\mathbf{\left\{}{\mathbf{\Delta u}}_{\mathbf{\mu }}^{\mathbf{elem}\left(t\right)}\mathbf{\right\}}}_{\mathbf{n}}^{{\mathbf{i}}_{\mathbf{RVE}}}$ usando lei de Hooke,

    • vetor de tensão de tentativa: ${\left\{{\mathbf{\sigma }}_{\mathbf{\mu }}^{\mathbf{elem}\left(t\right)}\right\}}_{\mathbf{n}}^{{\mathbf{i}}_{\mathbf{RVE}}}\mathbf{ }\mathbf{=}\mathbf{ }{\mathbf{\sigma }}_{(n-1)}^{\mathbf{elem}\mathbf{\left(}\mathbf{t}\mathbf{\right)}}\mathbf{ }\mathbf{+}\mathbf{ }{\{\Delta {\sigma }_{\mu }^{\mathrm{elem}\left(t\right)}\}}_{\mathbf{n}}^{{\mathbf{i}}_{\mathbf{RVE}}}$,

    • verifica o modelo constitutivo, obtendo a tensão ${\left\{{\mathbf{\sigma }}_{\mathbf{\mu }}^{\mathbf{elem}\left(t\right)}\right\}}_{\mathbf{n}}^{{\mathbf{i}}_{\mathbf{RVE}}}$ e o tensor constitutivo .

  • (a.2) Para cada elemento finito de contato ou fratura coesiva, ef, calcula-se:

    • O vetor de tentativa para o salto de deslocamento: ${\left\{{\mathrm{\delta }}_{\mathrm{\mu }}^{\mathrm{ef}\left(\mathrm{t}\right)}\right\}}_{\mathrm{n}}^{\left({\mathrm{i}}_{\mathrm{RVE}}\right)} = {\mathrm{N}}_{\mathrm{e}}{\left\{\Delta {\mathrm{u}}_{\mathrm{\mu }}^{\mathrm{ef}\left(\mathrm{t}\right)}\right\}}_{\mathrm{n}}^{\left({\mathrm{i}}_{\mathrm{RVE}}\right)}$, onde N_e é a função de forma relacionada ao ponto de Gauss.

    • Os modelos de contato e fratura são verificados, obtendo os valores nodais das forças ${\left\{{\mathrm{F}}_{\mathrm{ef}}^{\mathrm{int}}\right\}}_{\mathrm{n}}^{{\mathrm{i}}_{\mathrm{RVE}}}$ e a matriz do tensor constitutivo ${\left\{{\mathrm{K}}_{\mathrm{ef}}\right\}}_{\mathrm{n}}^{{\mathrm{i}}_{\mathrm{RVE}}}$ que são calculados em função se a fratura está em modo de abertura em carregamento, em descarregamento ou se permanece fechada.

  • b) O vetor ${\left\{{\mathrm{F}}_{\mathrm{R}}\right\}}_{\mathrm{n}+1}^{{\mathrm{i}}_{\mathrm{EVR}}}$ (equação (25)) é obtido de acordo com o modelo em multi-escala.

  • c) Verifica-se a convergência. 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 (18)), somando-se a parcela dos elementos triangulares elem e dos elementos de contato ou fratura ef; então, obtém-se a matriz ${\left\{{\mathrm{K}}_{\mathrm{R}}\right\}}_{\mathrm{n}}^{{\mathrm{i}}_{\mathrm{EVR}}}$ (equação (25)), de acordo com o modelo multi-escala.

  • e) Resolve o sistema dado pela equação (25) e obtém as correções de flutuações incógnitas ${\left\{\mathrm{\delta }\Delta {\tilde{\mathrm{u}}}_{\mathrm{R}}\right\}}_{\mathrm{n}}^{{\mathrm{i}}_{\mathrm{EVR}}+1}$.

  • f) Obtém o vetor total ${\left\{\mathrm{\delta }\Delta {\tilde{\mathrm{u}}}_{\mathrm{R}}\right\}}_{\mathrm{n}}^{{\mathrm{i}}_{\mathrm{EVR}}+1}$.

  • g) Calcula o novo incremento de deslocamento a ser imposto ao $\mathrm{EVR}:{\{\Delta u_{\mu }^t\}}_n^{i_{EVR}+1} = {\left[\Delta \epsilon \right]}_n^i \left\{\begin{array}{c}{x}_1\\ x_2\end{array}\right\} + {\left\{\Delta {\tilde{u}}_{\mu }\right\}}_n^{i_{EVR}+1}$, onde ${\left\{\Delta {\tilde{u}}_{\mu }\right\}}_n^{i_{EVR}+1} = {\left\{\Delta {\tilde{u}}_{\mu }\right\}}_n^{i_{EVR}} +{\left\{\delta \widetilde{\Delta 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 (30))

• 2.2.5 Obtém a tensão homogeneizada ${\sigma }_{n + 1}^{i + 1}$ (equação (28)); calcula incremento de tensão $\Delta {\sigma }_n^i = \Delta {\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 (10) do macrocontínuo.

• 3.1 Se o resíduo R_N não for nulo, a matriz tangente (equação (11)) é atualizada e as correções $\delta \Delta {\epsilon }_n^{i+1}$ calculadas. Atualiza-se o incremento de deformações $\delta \Delta {\epsilon }_n^{i+1}$ (equação (11)) 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 (6)) e os valores nodais de deslocamentos, forças e esforços internos (equações (5) e (7)).

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

2.4 Metodologia empregada nos exemplos numéricos

Com o intuito de investigar a influência do descolamento de fase no comportamento estrutural de placas submetidas a carregamentos tangenciais, análises numéricas são realizadas em nível microestrutural assim como em nível macroestrutural dentro de uma abordagem multi-escala. Como já dito, na formulação proposta, cada ponto do macrocontínuo é representado por um Elemento de Volume Representativo (EVR) cuja resposta homogeneizada define os tensores constitutivo e de tensão macroscópicos. O EVR representa o material em nível microestrutural que para simular os fenômenos dissipativos sobre o domínio do EVR, modelos constitutivos próprios devem ser adotados para cada fase constituinte do EVR. Logo, para tratar da ductilidade presente no comportamento mecânico da matriz de CMM, o modelo de Von Mises foi assumido enquanto que inclusões elásticas foram inseridas na matriz. Com o objetivo de mostrar a importância de considerar o descolamento de fases no processo de fraturamento do CMM, em algumas análises são assumidas inclusões perfeitamente aderidas para efeito de comparação de resultados. Observa-se que quando o descolamento de fase é considerado, um modelo de fratura coesiva para meios dúcteis desenvolvido em [⁷] deve ser assumido em conjunto com elementos finitos de contato e fratura dispostos ao redor das inclusões. Este modelo considera a separação tangencial e normal da zona coesiva.

Logo, para mostrar como o comportamento mecânico da microestrutura submetidas a estados de deformação predominantes de distorção influencia a solução do problema do macrocontínuo, inicialmente análises numéricas apenas em nível microestrutural são apresentadas. Em seguida, análises multi-escala, considerando micro e macro, são realizadas. Portanto, a seção 3 é dividida em duas subseções: análise numérica de EVRs e análise numérica multi-escala de estruturas, especificamente chapas submetidas a carregamentos de cisalhamento e compostas por compósitos de matriz metálica.

Todas as simulações foram realizadas com o emprego de da homogeneização computacional sob condições de estado plano de tensões, mais especificamente de cisalhamento puro, e regime de pequenas deformações. Além disso, os EVRs são quadrados com dimensões L x L e espessura L/10. A formulação multi-escala apresentada neste trabalho é objetiva em relação às dimensões do EVR, ou seja, qualquer valor pode ser adotado para L. A tensão média é obtida pela imposição de macrodeformação imposta ao EVR e subsequente solução do seu problema de valor de contorno inicial de acordo com a condição de contorno assumida em termos de flutuações de deslocamento que nos casos estudados aqui foi considerada a flutuação periódica. Além disso, uma tolerância de 10^-5 foi assumida para a verificação da convergência em ambos os procedimentos iterativos: o relacionado ao macrocontínuo e o requerido para atingir o equilíbrio do EVR.

3. RESULTADOS e discussão

3.1 Análise numérica de EVRs de compósitos de matriz metálica submetidos a estados de cisalhamento

Inicialmente, a resposta homogeneizada relacionada aos EVRs consideradas na análise multi-escala realizadas no item 3.2 serão discutidas a fim de mostrar a coerência do comportamento mecânico delas e a resposta da estrutura composta pelos EVRs analisados. Como explicitado anteriormente, procura-se avaliar o emprego da modelagem proposta para os casos em que EVRs e, num emprego amplo, estruturas, ambos estão submetidos a estados predominantes de cisalhamento. Ainda para discutir alguns pontos sobre as propriedades mecânicas dos materiais, principalmente sobre a zona de interface, os parâmetros dos modelos constitutivos são adotados considerando os trabalhos de AZIZI [²], BORGES et al. [³⁰], BORGES e PITUBA [³³] e KIM e AL-RUB [³⁴]. A Figura 2 mostra, de forma genérica, a modelagem em elementos finitos usada no nível microestrutural assim como os modelos constitutivos adotados para cada fase (matriz, inclusões e zona de interface, se o descolamento de fase é considerado).

Figura 2
Visualização ilustrativa amplificada de uma microestrutura: modelagem em elementos finitos e modelos constitutivos correspondentes.

Para o primeiro exemplo numérico, um EVR (denominado de EVR1) com fração volumétrica (v_f) de 37% e inclusões aleatoriamente disposta na matriz é considerado (Figura 3). Diferentes dimensões são adotadas para as inclusões com o objetivo de definir uma forte heterogeneidade para o material. O EVR1 foi discretizado com 540 elementos finitos triangulares, sendo adotado para a matriz o modelo de Von Mises com endurecimento isótropo contendo as seguintes propriedades: v = 0,2;, E = 20 GPa e σy = 0,7MPa. Para definir a lei de endurecimento na matriz foram adotados os seguintes valores para (ε^p, σ_y), onde ε^p é a deformação plástica efetiva: (0; 0,7MPa); (0,6; 5 MPa), ou seja, um coeficiente da curva de encruamento K=7,167 MPa. Já para as inclusões elásticas foram considerados v = 0,35; E = 40 GPa. Para a análise onde o descolamento de fase é considerado, 95 elementos finitos coesivos e contato foram definidos nas interfaces entre matriz e inclusões, sendo seus parâmetros dados por: δc = 6 ×10^-5 m, σc = 10 MPa, β_ef = 0,707 e λp = 3 × 10⁶.

Figura 3
.EVR para compósito de matriz metálica com 5 inclusões (EVR1).

Como o trabalho trata do descolamento de fase em microestruturas dúcteis, o EVR1 foi analisado para o caso de estados de cisalhamento, onde a deformação macroscópica imposta (ε = [ε_x; ε_y; ε_xy]) relacionada a estados de cisalhamento, ε = [0,0000001; 0,0000001; 0,001], foi aplicada em 20 incrementos sendo a resposta homogeneizada mostrada na figura 4 considerando inclusões perfeitamente aderidas e descolamento de fase.

Figura 4
Respostas numéricas homogeneizadas para o EVR1 com 5 inclusões para deformação macroscópica imposta ε = [0,0000001; 0,0000001; 0,001].

Como pode ser observado na figura 4, a consideração do descolamento de fase não muda os resultados numéricos em ambos os casos porque a resistência à tração (σc = 10 MPa) é alta se comparada ao nível máximo de tensão que ocorre em pontos pertencentes ao domínio do EVR1 (0,8 MPa) levando a uma não nucleação de fratura nas interfaces da microestrutura. Observe que a resistência à tração é um critério para abertura de fratura diretamente relacionada com a abertura por escorregamento entre as faces da fratura pelo parâmetro β_ef.

Para o segundo exemplo numérico, um EVR (denominado EVR2) apresentando v_f = 10% com apenas uma inclusão posicionada no centro da microestrutura é considerado (Figura 5). Para a discretização foram adotados 580 elementos finitos triangulares, assim como 16 elementos finitos coesivos e de contato se o descolamento de fase é considerado. O modelo de Von Mises com endurecimento isótropo (K = 0,032E) foi adotado para a matriz, assim como as seguintes propriedades: v = 0,2; E = 70GPa e σ_y = 243MPa. Para a região de interface foram assumidos os mesmos valores de parâmetros do modelo de fratura coesiva e contato do exemplo anterior. Para a inclusão elástica foram adotados: v = 0,2 e E = 200 GPa.

Figura 5
EVR para microestrutura de compósito com matriz dúctil com uma inclusão rígida centrada (EVR2)

Com os mesmos argumentos do exemplo anterior, um estado de deformação macroscópica é imposto ao EVR2 para tratar de estados predominantes de cisalhamento definido por ε = [0,0000001; 0,0000001; 0,01]. A figura 6 ilustra a resposta homogeneizada pela aplicação do estado de deformação em 20 incrementos considerando inclusões perfeitamente aderidas à matriz e descolamento de fase. Como neste segundo exemplo numérico os módulos de elasticidade adotados para a matriz e inclusões são muito maiores que os valores do exemplo anterior, o nível de tensão no EVR2 atinge o critério de nucleação de fratura na zona de interface. Portanto, o descolamento de fase é um importante influenciador no comportamento mecânico do EVR2 mesmo em estágios iniciais de carregamento, como visto na figura 6.

Figura 6
Resposta numérica homogeneizada para EVR2 com uma inclusão para a deformação macroscópica imposta )

De acordo com os resultados numéricos apresentados nesta seção, pode-se concluir que os parâmetros adotados para o modelo de fratura coesiva e contato, assim como para a relação entre as propriedades elásticas da matriz e inclusões, têm um papel fundamental na simulação do comportamento mecânico da microestrutura do compósito. Mais detalhes sobre a sensibilidade paramétrica do modelo de fratura coesiva e contato podem ser encontrados em SANTOS et al. [³²].

3.2 Análise numérica de estruturas submetidas a estados de cisalhamento e compostas por compósito de matriz metálica

Nesta seção são apresentadas análises multi-escala de chapas submetidas a carregamentos de cisalhamento, sendo utilizados os EVRs definidos na seção anterior, assim como os parâmetros dos modelos constitutivos para a matriz, inclusões e zona de interface (se considerado o processo de descolamento de fase). De posse dos resultados das microestruturas (EVR1 e EVR2) analisadas no item anterior, procura-se verificar no presente item a repercussão na resposta da macroestrutura. Como já dito anteriormente, cada EVR representa um ponto do macrocontí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.

Considerou-se como macroestrutura, a chapa de alumínio de dimensão (60cm x 30cm) e 1 cm de espessura, como mostrado na figura (7a), onde é aplicada em suas extremidades uma carga tangencial, a fim de simular um caso de cisalhamento puro. 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 7b). Como condição de contorno, foi definida uma tração ps = 20 N/cm aos nós 19 a 23; o lado correspondente aos nós 6 a 18 e aquele referente aos nós 24 a 36 foram considerados livres enquanto prescreveu-se u_n = u_s = 0 para o lado correspondente aos nós 1 a 5, sendo n e s as direções normal e tangencial ao contorno, 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 (7b).

Figura 7
a) Geometria da chapa submetida a um carregamento de cisalhamento b) Visão ampliada da discretização de um quarto da chapa em células internas e elementos de contorno

O deslocamento do nó 21da chapa, na direção X₁, (ver Figura 7b) em função do fator de carga β para as diferentes análises pode ser visto na figura 8 considerando o material definido pelo EVR 1 e na figura 9 adotando o material definido pelo EVR 2. A Tabela 1 resume os dados das propriedades mecânicas dos EVRs utilizados nas análises.

Figura 8
Deslocamento no ponto 21 da chapa ao longo do processo incremental considerando-se o EVR 1, com 5 inclusões

Figura 9
Deslocamento no ponto 21 da chapa ao longo do processo incremental considerando-se o EVR 2, com 1 inclusão

Tabela 1
Propriedades mecânicas dos materiais dos EVRs 1 e 2

EVR 1 – Disposição aleatório de inclusões
Matriz				Inclusão
ν	E (GPa)	σy (MPa)	K (MPa)	ν	E (GPa)	Vf (%)
0,2	20	0,7	7,167	0,35	40	37
EVR 2 – Disposição centrada de inclusão unitária
Matriz				Inclusão
ν	E (GPa)	σy (MPa)	K (MPa)	ν	E (GPa)	Vf (%)
0,2	70	243	2240	0,2	200	10

De uma maneira geral, como pode-se observar nas figuras. 8 e 9, os deslocamentos finais obtidos para a estrutura com o emprego dos dois EVRs foram similares, por volta de 0,03cm. Porém a chapa onde considerou-se o EVR 2 suportou uma carga cerca de 10 vezes maior que aquela suportada quando adotou-se o EVR 1, o que é coerente, pois a inclusão do EVR 2 tem um módulo de elasticidade cinco vezes maior que aquele do EVR 1 e sua matriz foi definida com um módulo de elasticidade 3,5 vezes maior e uma tensão de escoamento 350 vezes maior, isso tudo levando a níveis de tensões elevadas desencadeando todo o processo de descolamento de fase mostrando coerência das respostas micromecânicas (figura 6) e macromecânicas (figura 9). Outro aspecto importante a ser observado nos gráficos da figura 8 é que quando se adotou o EVR 1, a chapa não apresentou diferença significativa no seu comportamento mecânico quando foram considerados elementos de fratura e contato nas interfaces entre matriz e inclusão. Isso é coerente com o comportamento apresentado na análise da microestrutura (figura 4), onde a consideração do descolamento de fase na microestrutura não apresenta ser um fator importante no fenômeno dissipativo devido à relação entre as propriedades dos materiais na microestrutura (inclusões e matriz) o que leva a um nível de tensão na microestrutura insuficiente para iniciar um processo de nucleação na zona de interface.

Por fim, a diferença apresentada na figura 9 foi bem significativa nas análises onde se adotou o EVR 2, pois a consideração de elementos de fratura e contato reduziu consideravelmente a resistência e rigidez da chapa. Verifica-se também que as chapas analisadas apresentam um comportamento mecânico com endurecimento gradual [³⁵], difícil de definir uma tensão de plastificação homogeneizada representativa do material da chapa submetida a estado de cisalhamento, devido à presença da matriz dúctil tanto com a consideração do fenômeno de descolamento de fase quanto com a consideração de perfeita aderência entre as fases.

4. CONCLUSÕES

Este trabalho tratou da aplicação de um modelo multi-escala desenvolvido em FERNANDES et al. [¹⁷] e estendido no presente trabalho para a incorporação de efeitos de cisalhamento no processo de descolamento de fase que ocorre na microestrutura e sua repercussão na resposta macroscópica homogeneizada. Para a análise crítica da formulação foram analisadas placas submetidas a carregamentos de cisalhamento e compostas por materiais heterogêneos considerando descolamento de fase, onde o Método dos Elementos de Contorno e o Método dos Elementos Finitos foram adotados para modelar o macro e microcontínuos, respectivamente. Neste modelo, que é objetivo em relação às dimensões do EVR, a resposta macroscópica é definida por campos homogeneizados de tensão e do tensor constitutivo sobre o domínio do EVR (Elemento de Volume Representativo) que representa o material em nível microestrutural. A microestrutura do EVR é definida por uma matriz onde inclusões podem ser definidas de maneira a melhorar as características mecânicas do material ou vazios que podem ser considerados para modelar o comportamento de materiais porosos, apesar que este segundo caso não foi tratado no presente trabalho. Para governar o comportamento do material da região de matriz foi adotado o modelo elastoplástico de von Mises enquanto que as inclusões foram assumidas como tendo comportamento elástico. Afim de capturar o processo de descolamento de fase que ocorre, por exemplo, em compósitos de matriz metálica (CMM) submetidos a estados de cisalhamento, elementos finitos coesivos e de contato (PITUBA et al. [⁷]) foram definidos na zona de interface ao redor das inclusões.

As respostas qualitativas obtidas se mostraram bastante satisfatórias e coerentes evidenciando a boa descrição do comportamento microestrutural dos EVRs compostos por matriz metálica reforçados com inclusões elásticas e sua repercussão na resposta macroscópica de estruturas compostas por esses materiais. Os exemplos numéricos mostraram que se o nível de tensão atinge o critério de nucleação de fratura na zona de interface, o fenômeno de descolamento de fase assume um papel importante no comportamento macroestrutural. Portanto, quando o fenômeno de descolamento de fase é levado em conta em materiais que apresentam ductilidade sob cisalhamento, a melhoria das propriedades mecânicas do material não depende apenas da fração volumétrica de inclusões e quantidade de inclusões inseridas na matriz, mas também da relação entre as propriedades mecânicas da matriz. Por sua vez, SANTOS et al. [³²] e FERNANDES et al. [¹⁷] identificaram que a quantidade de inclusões intensifica a influência do fenômeno de descolamento de fase tanto em microestruturas submetidas a estados predominantes de tração quanto na sua repercussão na resposta de macrocontínuos tracionados.

Observe que nos exemplos numéricos discutidos neste trabalho não se considera o caso de materiais onde o problema de localização de deformação na microestrutura pode levar à nucleação na macroestrutura, como por exemplo, em materiais quase-frágeis. Este tipo de problema será tratado em trabalhos futuros onde um critério deve ser adotado adequadamente para definir quando a nucleação de fratura ocorre no macrocontínuo. Além disso, para problemas de localização, a presente formulação multi-escala pode ser estendida adequadamente com o objetivo de assegurar a equivalência de energia entre a microestrutura e macrocontínuo de modo a efetuar corretamente a transição entre escalas para esse caso, vide TORO et al. [²²] e SANCHEZ et al. [³⁶]. Observa-se que os CMMs podem apresentar localização de deformação se o processo de faturamento propagar também na matriz, porém não de forma tão intensa como em materiais quase-frágeis que leva a um processo de nucleação e propagação de fratura no macrocontínuo. Contudo, este caso não é considerado nos exemplos numéricos onde o comportamento da matriz é governado por um modelo elastoplástico sem considerar o processo de dano/faturamento. Ainda mesmo se as análises numéricas fossem realizadas com a presente formulação utilizando controle de deslocamentos ou métodos de comprimento de arco, a possibilidade de surgir um comportamento global com softening poderia levar a uma inconsistência energética [³¹] e consequentemente, respostas numéricas inválidas. Todos esses aspectos são alvos de desenvolvimento da presente formulação para a futura aplicação em materiais como o concreto, por exemplo.

Por fim, as análises apresentadas no presente trabalho são válidas de acordo com a presente formulação concluindo-se que a utilização de modelos constitutivos simples na microestrutura aliada à formulação multi-escala leva a análises mais realísticas da macroestrutura.

AGRADECIMENTOS

Ao Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq e Fundação ao Amparo à Pesquisa do Estado de Goiás, FAPEG, pelo suporte financeiro fornecido durante a realização do trabalho.

BIBLIOGRAFIA

• ¹
  PITUBA, J. J. C., FERNANDES, G. R. “An anisotropic damage model for concrete”, Journal of Engineering Mechanics-ASCE, v. 137, pp. 610-624, 2011.
• ²
  AZIZI, R., “Micromechanical modeling of damage in periodic composites using strain gradient plasticity”, Engineering Fracture Mechanics, v. 92, pp. 101-113, 2012.
• ³
  OLIVER, J., CAICEDO, M., ROUBIN, E.; et al Continuum approach to computational multiscale modeling of propagating fracture. Computer Methods in Applied Mechanics and Engineering, v. 294, pp. 384-427, 2015.
• ⁴
  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.
• ⁵
  NGUYEN, V.P., LLOBERAS VALLS, O., STROEVEN, M. et al On the existence of representative volumes for softening quasi-brittle materials – a failure zone averaging scheme. Computer Methods in Applied Mechanics and Engineering, 199: pp.3026–3036, 2010.
• ⁶
  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.
• ⁷
  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, DOI: 10.1061/(ASCE)EM.1943-7889.0001123, 2016.
  » https://doi.org/10.1061/(ASCE)EM.1943-7889.0001123
• ⁸
  SANTOS, W. F., PITUBA, J. J. C.. “Yield surfaces of material composed of porous and heterogeneous microstructures considering phase debonding”. Latin American Journal of Solids and Structures, v. 14, pp. 1387-1415, 2017.
• ⁹
  CAVALCANTE, M. A. A., MARQUES, S. P. C.. Microstructure effects in wavy-multilayers with viscoelastic phases. European Journal of Mechanics A, Solids, v. 64, pp. 178-185, 2017.
• ¹⁰
  Li L. X., Chen Y. L., Lu Z. C. Generalization of multi-scale finite element method to plane elasticity problems. Applied Mathematical Modelling, v. 39, n. 2, 642-653, 2014.
• ¹¹
  de Souza Neto, E. A., Blanco, P. J., Sánchez, P. J. et al An RVE-based multiscale theory of solids with micro-scale inertia and body force effects. Mechanics of Materials. v. 80, pp. 136-144, 2015.
• ¹²
  Edmans, B. D., Alfano, G., Bahai, H. Nonlinear multi-scale homogenization with different structural models at different scales. International Journal for Numerical Methods in Engineering. v. 94, pp. 355-373, 2013.
• ¹³
  Rodrigues, D. E. S., Belinha, J., Pires, F. M. A., et al Homogenization technique for heterogeneous composite materials using meshless methods, Engineering Analysis with Boundary Elements. v. 92, pp.73-89, 2018.
• ¹⁴
  Cavalcante, M. A. A., Pindera, M-J, Generalized FVDAM theory for elastic-plastic periodic materials. International Journal of Plasticity, v. 77, pp. 90-117, 2016.
• ¹⁵
  Patil, R. U., Mishra, B. K., Singh, I. V., A new multiscale XFEM for the elastic properties evaluation of heterogeneous materials. International Journal of Mechanical Sciences, v. 122, pp.277-287, 2017.
• ¹⁶
  FURTADO, A. S., PITUBA, J. J. C., FERNANDES, G. R.. “Análise da influência de microestruturas heterogêneas na resposta macromecânica do problema bidimensional de placas”. Materia-Rio de Janeiro, v. 22, pp. e11830, 2017.
• ¹⁷
  FERNANDES, G. R.; FURTADO, A. S.; PITUBA, J. J. C.; et al “Multiscale Analysis of Structures Composed of Metal Matrix Composites Considering Phase Debonding”. Journal of Multiscale Modeling, v. 8, n. 3&4, p. 1740004, 2017.
• ¹⁸
  SOUZA JÚNIOR, N. S., FERNANDES, G. R., PITUBA, J. J. C.. Formulação multi-escala para a análise de flexão de placas considerando processos dissipativos na microestrutura e acoplamento MEC/MEF, Matéria (UFRJ), v. 22, e11821, 2017.
• ¹⁹
  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.
• ²⁰
  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.
• ²¹
  BENEDETTI, I., ALIABADI, M. H. A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials. Computer Methods in Applied Mechanics and Engineering, v. 265, pp. 36-62, 2013.
• ²²
  FERNANDES, G. R., Ohland, G. A., VIEIRA, J. F.. A Boundary Element Formulation to perform elastic analysis of heterogeneous microstructures. ENGINEERING ANALYSIS WITH BOUNDARY ELEMENTS, v. 87, pp. 47-65, 2018.
• ²³
  Sfantos G.K., Aliabadi M.H., Multi-scale boundary element modelling of material degradation and fracture, Computer Methods in Applied Mechanics and Engineering. n. 196, pp. 1310–1329, 2007.
• ²⁴
  YANG, Q.-S.; QIN, Q.-H. Micro-Mechanical Analysis of Composite Materials by BEM. Engineering Analysis with Boundary Elements, v. 28, n. 8, pp. 919-926, 2004.
• ²⁵
  RODRIGUES, R.Q.; GALVIS, A.F.; SOLLERO, P.; et al Fast BEM Multi-Domain Approach for the Elastostatic Analysis of Short Fibre Composites. European Journal of Computational Mechanics, v. 26, n.5, pp. 525-540, 2017.
• ²⁶
  ARAUJO, F.C.; GRAY, L.J. Evaluation of Effective Material Parameters of CNT-Reinforced Composites via 3D BEM. Computer Modeling in Engineering & Sciences, v. 24, n. 2-3, pp. 103-121, 2008.
• ²⁷
  DORO, V. C., FERNANDES, G. R.. Análise não linear de chapas através de uma formulação do método dos elementos de contorno com convergência quadrática. Matéria (UFRJ), v. 21, pp. 27-48, 2016.
• ²⁸
  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.
• ²⁹
  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, n. 16, 2006.
• ³⁰
  BORGES, D. C., QUARESMA, W. M. G., FERNANDES, G. R., et al. Evaluation of a proposed model for concrete at mesoscopic scale. REVISTA IBRACON DE ESTRUTURAS E MATERIAIS, v. 10, pp. 1087-1112, 2017.
• ³¹
  TORO, S., SANCHEZ, P.J. BLANCO, P. J., SOUZA NETO, E. A., HUESPE, A.E; FEIJÓO, R. A. Multiscale formulation for material failure accounting for cohesive cracks at the macro and micro scales. Journal of Plasticity, v. 76, pp. 75-110, 2016.
• ³²
  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.
• ³³
  BORGES, D. C., PITUBA, J. J. C. “Analysis of quasi-brittle materials at mesoscopic level using homogenization model”. Advances in Concrete Construction, v. 5, pp. 221-240, 2017.
• ³⁴
  KIM, S. M., AL-RUB, R. K. A. Meso-scale computational modeling of the plastic-damage response of cementitious composites. Cement and Concrete Research, vol. 41, pp.339-358, 2011
• ³⁵
  HOSFORD, W. F., CADDELL, R. M., Metal Forming: Mechanics and Metallurgy Cambridge University Press, 4^th Edition, 1993.
• ³⁶
  SÁNCHEZ, P. J., BLANCO, P. J., HUESPE, A. E., et al Failure-oriented multi-scale variational formulation: micro-structures with nucleation and evolution of softening bands. Computer Methods in Applied Mechanics and Engineering, n. 257, pp. 221-247, 2013.

Datas de Publicação

• Publicação nesta coleção
  16 Set 2020
• Data do Fascículo
  2020

Histórico

• Recebido
  10 Jul 2018
• Aceito
  11 Nov 2019