Acessibilidade / Reportar erro

Formulação multi-escala para a análise de flexão de placas considerando processos dissipativos na microestrutura e acoplamento MEC/MEF

Multi-scale formulation for analysis of the plate bending problem considering dissipative processes in the microstructure and coupling BEM and FEM

RESUMO

Neste trabalho apresentam-se análises de flexão de placas compostas por materiais heterogêneos através de uma abordagem multi-escala. O macro-contínuo, representado neste trabalho pela placa, é modelado por uma formulação não-linear do Método dos Elementos de Contorno (MEC), que leva em conta o operador tangente consistente (CTO). A micro-escala é representada pelo EVR (Elemento de Volume Representativo), sendo seu problema de equilíbrio definido em termos de flutuação dos deslocamentos e solucionado através do Método dos Elementos Finitos (MEF), onde a hipótese de média volumétrica das tensões e deformações é adotada para se fazer a passagem do micro-contínuo para o macro-contínuo. A cada ponto do macro-contínuo, onde se necessita conhecer as tensões e o tensor constitutivo deve estar associado um EVR, onde se podem definir inclusões e/ou vazios no interior de uma matriz a fim de representar a micro-estrutura de um material heterogêneo. Nos exemplos numéricos são considerados diferentes EVRs com inclusões elásticas dentro de uma matriz, onde os modelos de Von Mises ou Mohr Coulomb são adotados, a fim de governar o comportamento do seu material. Consideram-se diferentes frações volumétricas para as inclusões a fim de verificar a influência na resposta homogeneizada da microestrutura e, consequentemente, no comportamento mecânico do macro-contínuo. Para solucionar o problema de equilíbrio do EVR devem-se adotar condições de contorno em termos de flutuações dos deslocamentos, que nos exemplos analisados no presente trabalho serão consideradas como periódicas.

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

ABSTRACT

Analyses of the bending problem of plates composed of heterogeneous materials are performed considering a multi-scale modelling. The macro-continuum, represented in this paper by the plate, is modelled by a nonlinear formulation of the boundary element method (BEM) taking into account the consistent tangent operator (CTO) in the iterative procedure required to solve the plate equilibrium problem. The micro-scale is represented by the RVE (Representative Volume Element) being its equilibrium problem solved in terms of dis-placements fluctuations by a Finite Element Formulation (FEM), where the volume averaging hypothesis of strain and stress tensors is used to make the micro-to-macro transition. To each point of the macro-continuum where the stresses and constitutive tensor have to be computed, a RVE must be assigned, where inclusions and voids can be defined inside the matrix in order to represent the microstructure of a heterogeneous material. In the numerical examples are considered different RVEs with elastic inclusions, while for the matrix the Von Mises or Mohr-Coulomb criteria can be adopted to govern its material behavior. Different volume fractions have been adopted for the inclusions in order to verify its influence on the homogenized response of the microstructure as well as on the mechanical behavior of the macrostructure. To solve the RVE equilibrium problem, boundary conditions in terms of displacement fluctuations have to be imposed, which for the numerical examples presented in this paper will be adopted periodic.

Keywords:
Multi-scale modelling; homogenization; boundary elements; plate bending

1. INTRODUÇÃO

A maioria dos materiais é heterogênea ao nível dos grãos ou da sua microestrutura. Além disso, a micro-estrutura de alguns materiais são manipulados adicionando novos constituintes, a fim de melhorar suas propriedades. Como qualquer heterogeneidade do material, assim como o processo de fissuração que ocorre na micro-escala afetam diretamente a resposta da estrutura [1[1] PITUBA, J. J. C., FERNANDES, G. R., “An anisotropic damage model for concrete”, Journal of Engineering Mechanics-ASCE, v. 137, n.9, pp. 610-624, 2011.4[4] MATALLAH, M., LA BORDERIE, C., “Inelasticity–damage-based model for numerical modeling of concrete cracking”, Engineering Fracture Mechanics, v. 76, pp. 1087-1108, 2009.], sua modelagem em diferentes escalas é muito importante para melhor representar o comportamento de materiais complexos [5[5] GAL, E., KRYVORUK, R, “Fiber reinforced concrete properties – a multiscale approach”, Computers and Concrete, An International Journal, v. 8, n.5, pp. 525-539, 2011.9[9] PITUBA, J. J. C., SOUZA NETO, E. A. “Modeling of unilateral effect in brittle materials by a mesoscopic scale approach”, Computers and Concrete, An International Journal, v. 15, n.5, pp. 735-758, 2015.]. Em muitos casos, a análise não-linear convencional não consegue representar de forma precisa o comportamento de tais estruturas.

Para se fazer a modelagem multi-escala deve-se definir no domínio do macro-contínuo pontos de interesse que são chamados de EVR (Elemento de Volume Representativo), os quais representam a microestrutura, ao nível dos grãos, do macro-contínuo na vizinhança infinitesimal do ponto [10[10] 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[15] SAAVEDRA-FLORES, E. I., SOUZA NETO, E. A., PEARCE, C., “A large strain computational multi-scale model for the dissipative behavior of wood cell-wall”. Computational Materials Science, v. 50(3): pp. 1202-1211, 2011.]. Podem-se definir vazios ou inclusões (elásticas ou elasto-plásticas) dentro da matriz do EVR, sendo a formação e propagação de micro-fissuras ou o processo de plastificação monitorados individualmente em cada EVR, considerandose as diferentes propriedades elásticas e modelo constitutivo adotado para cada fase do EVR. Assim, considerando-se uma estrutura sujeita a um carregamento qualquer, através do modelo na macro-escala, obtém-se para cada EVR o campo de deformações, que deverá ser imposto de forma constante a todos os pontos do EVR, como “carregamento”. Então, utilizando-se o modelo adotado na micro-escala, estuda-se o comportamento do material nas diferentes fases do EVR devido a essa deformação imposta. Em seguida, através de princípios de homogeneidade e conceito de média volumétrica, passa-se da micro-escala para a macro-escala e atualizam-se as tensões e a relação constitutiva para aquele ponto do macro. Com a relação constitutiva atualizada para todos os pontos da macro-escala (neste trabalho tratamos de placas), dá-se novo incremento de carga obtendo-se, através do modelo na macro-escala, novos campos de deformações a serem aplicados nos EVRs. Portanto, a análise na micro-escala alimenta aquela na macro-escala e vice-versa. Neste trabalho, a microestrutura é modelada pelo Método dos Elementos Finitos (MEF), através da formulação apresentada em [12[12] 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.].

A formulação do MEC a ser usada na análise da macro-escala é descrita em [16[16] FERNANDES, G. R., SOUZA NETO, E. A., “Self-consistent linearization of non-linear BEM formulations with quadratic convergence”, Computational Mechanics. DOI: 10.1007/s00466-013-0867-2, v. 52, pp. 1125-1139, 2013.
https://doi.org/10.1007/s00466-013-0867-...
], onde se deve discretizar o contorno da placa em elementos nos quais os deslocamentos e esforços são aproximados para se obter a solução do problema. 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 aproximados os momentos iniciais ou inelásticos. Os momentos num ponto da placa são calculados integrando-se numericamente as tensões ao longo da espessura, usando um sistema gaussiano. Assim, para a se obter a solução em multi-escala, deve-se definir um EVR em cada ponto de Gauss definido ao longo da espessura e relativo a um determinado nó do macro-contínuo. Em suma, este trabalho é um desenvolvimento dos trabalhos anteriores dos autores procurando-se validar a formulação proposta com novas aplicações e emprego do modelo constitutivo de Mohr-Coulomb com vistas às aplicações futuras em compósitos como o concreto. Diante disso, são apresentados dois exemplos numéricos. No primeiro, simulações de estruturas homogêneas empregando-se modelos constitutivos fenomenológicos são comparadas com a modelagem proposta quando da inserção de heterogeneidade na microestrutura e sua repercussão não macroccontínuo. Já no segundo exemplo, adota-se uma inclusão no centro do EVR com diferentes frações de volume, a fim de verificar tanto a influência na resposta homogeneizada da microestrutura, quanto no comportamento mecânico da placa.

2. MATERIAIS E MÉTODOS

2.1 Formulação do MEC para modelar o macro-contínuo (problema não-linear de flexão de placas)

2.1.1 Relações Básicas

Neste trabalho, o problema do macro-contínuo é definido pelo problema não-linear de flexão de placas, que será modelado pelo Método dos Elementos de Contorno (MEC), através da formulação descrita em [16[16] FERNANDES, G. R., SOUZA NETO, E. A., “Self-consistent linearization of non-linear BEM formulations with quadratic convergence”, Computational Mechanics. DOI: 10.1007/s00466-013-0867-2, v. 52, pp. 1125-1139, 2013.
https://doi.org/10.1007/s00466-013-0867-...
], baseada nas hipóteses de Kirchhoff. Para definir o problema considere-se uma placa de espessura t, contorno externo Γ e domínio Ω, onde as variáveis são definidas segundo um sistema de coordenadas cartesianas onde x1 e x2 são as direções no plano da placa e x3 a direção transversal a esse plano. Assume-se que a placa suporta apenas cargas distribuídas g distribuídas na 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 (x˙)=dx/dt. Assim, as variáveis relacionadas ao problema de flexão de placas são: V˙n (taxa do esforço cortante equivalente); M˙n (taxa do momento de flexão); w˙ (taxa da deflexão); w˙,n (taxa da derivada direcional), 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 de flexão de placas serão omitidas aqui, mas as mesmas podem ser encontradas em [19[19] 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.]-[24[24] 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.].

As taxas dos momentos de flexão e de torção m˙ij são obtidas integrando-se as taxas de tensões σ˙ij ao longo da espessura th da placa, como segue:

(1) m ˙ i j = t h / 2 t h / 2 z σ ˙ i j d z

No contexto não-linear, as tensões (σ˙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 ε˙ije e plásticas ε˙ijp como segue:

(2) ε ˙ i j = ε ˙ i j e + ε ˙ i j p

Observe que, aplicando-se a lei de Hooke, as tensões σ˙ij são relacionadas às deformações elásticas, sendo ε˙ije=ε˙ijε˙ijp. Além disso, adotando-se as hipóteses de Kirchhoff, as deformações totais podem ainda serem escritas em termos das curvaturas w˙,ij:ε˙ij=x3w˙,ij. Assim, as taxas de momentos elásticos de tentativa m˙ije, podem ser escritos em termos das curvaturas totais w˙,ij:

(3) m ˙ i j e = D v w ˙ , k k δ i j + ( 1 v ) w ˙ , i j

onde δij é o delta de Kronecker, D = Eth3/(12(1 - v2)) é a rigidez à flexão da placa e v o coeficiente de Poisson.

Portanto, a taxa de momentos plásticos ou inelásticos, se outro tipo de modelo constitutivo for adotado, podem ser definidos como:

(4) m ˙ i j p = m ˙ i j e m ˙ i j

2.1.2 Equações Integrais

A equação integral do deslocamento para um ponto interno é obtida a partir do Teorema de Reciprocidade de Betti, que no caso de flexão de placas é dado por:

(5) Ω w , j k * m ˙ j k d Ω = Ω m j k * w ˙ , j k d Ω Ω w , j k * m ˙ j k ( p ) d Ω

onde os termos com índice * são relacionados ao problema fundamental.

Integrando-se a equação (5) por partes duas vezes, se obtém a conhecida representação de deflexão para um ponto de colocação q:

(6) K ( q ) w ( q ) = Γ ( V ˙ n * w ˙ M n * w ˙ n ) d Γ j = 1 N C R c j * w c j + j = 1 N C R ˙ c j w c j * + Γ ( V ˙ n w * M ˙ n w ˙ n ) d Γ + Ω g g ˙ w * d Ω Ω w , j k * m ˙ j k p d Ω

onde w*,V*eM* são valores fundamentais de deflexão, cortante equivalente e momentos de contorno; Nc é o número de cantos e Ωg a região carregada da placa. Para os valores dos termos livres K(q), ver em [18[18] 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.].

Note que a solução não-linear é obtida a partir de um processo incremental de carga. Assim considerando-se Δt = tn+1 - tn um passo de tempo, relacionado a um incremento de carga, o problema consiste em achar a solução no passo de tempo tn+1, sendo a solução já conhecida em tn. 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 curvaturas ΔW,ℓi é obtida derivando-se a equação (6) duas vezes, ou seja:

(7) Δ w , i = Γ ( V n , i * Δ w M n , i * Δ w n ) d Γ j = 1 N C R c j , i * Δ w c j + j = 1 N C Δ R c j w c j , i * + Γ ( Δ V n w , i * Δ M n w , n i * ) d Γ + Ω g Δ g w , i * d Ω g x i Ω w , j k * Δ m j k p d Ω

2.1.3 Equações Algébricas

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 dos momentos inelásticos (ou plásticos). Ao longo do contorno têm-se quatro valores definidos: Δw, Δw,n, ΔMn, e ΔVn, sendo dois desses valores dados como condição de contorno. Portanto, duas equações devem ser escritas para cada nó de contorno. Neste trabalho, será considerada uma equação de ΔW escrita no próprio nó e outra equação de ΔW escrita para um ponto externo bem próximo ao contorno. Esse esquema já foi adotado em diversos trabalhos, como por exemplo, [20[20] 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[23] FERNANDES, GR., KONDA, DH. “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.], gerando bons resultados. Além disso, nos cantos são definidas duas variáveis, sendo uma delas dada como condição de contorno. Portanto, uma equação de ΔW também deve ser escrita em cada canto da placa. 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 [22[22] FERNANDES, G. R., VENTURINI, W. S. “Building floor analysis by the Boundary element method”, Computational Mechanics, v. 35, pp. 277 - 291, 2005.]).

Para se obter a solução não-linear, devem-se ainda escrever 3 equações dos momentos elásticos de tentativa ΔMjke nos nós das células. Essas equações são obtidas considerando-se a equação (3), onde as curvaturas 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 inserir as condições de contorno, pode ser escrito como (ver mais detalhes em [16[16] FERNANDES, G. R., SOUZA NETO, E. A., “Self-consistent linearization of non-linear BEM formulations with quadratic convergence”, Computational Mechanics. DOI: 10.1007/s00466-013-0867-2, v. 52, pp. 1125-1139, 2013.
https://doi.org/10.1007/s00466-013-0867-...
]):

(8) Δ X = Δ L + R M Δ m P

onde o vetor ΔX contêm as incógnitas do contorno, ΔL representa a parte elástica dessas incógnitas, RM representa as correções devido ao incremento de momentos inelásticos.

Escrevendo-se as equações algébricas de ΔMjke para todos os nós de células, obtém-se (ver mais detalhes em [16[16] FERNANDES, G. R., SOUZA NETO, E. A., “Self-consistent linearization of non-linear BEM formulations with quadratic convergence”, Computational Mechanics. DOI: 10.1007/s00466-013-0867-2, v. 52, pp. 1125-1139, 2013.
https://doi.org/10.1007/s00466-013-0867-...
]):

(9) Δ M e = Δ K + S M Δ m P

onde ΔK é a solução elástica dada em termos de incremento de momentos, SM expressa o efeito do incremento de momentos inelásticos ΔmP.

Observe que ΔmP, ΔMe e Δm definidos, respectivamente, nas equações (4), (3) e (1) são calculados localmente para um determinado ponto, isto é, são obtidos levando-se em conta apenas seu incremento de curvatura ΔW,ij e de tensões Δσij (obtidas após a solução do EVR). Por outro lado, os incrementos de momentos elásticos ΔMe definidos na equação (9) são calculados levando-se em conta os incrementos de momentos inelásticos ΔmP de todos os nós da placa. Neste contexto, pode-se também definir a equação algébrica do incremento de momentos ΔM (ver mais detalhes em [16[16] FERNANDES, G. R., SOUZA NETO, E. A., “Self-consistent linearization of non-linear BEM formulations with quadratic convergence”, Computational Mechanics. DOI: 10.1007/s00466-013-0867-2, v. 52, pp. 1125-1139, 2013.
https://doi.org/10.1007/s00466-013-0867-...
]):

(10) Δ M = C M Δ χ Δ K S M ( Δ m P ) + Δ m

onde os valores nodais ΔmP são dados por:

(11) Δ m P = C M Δ χ Δ m

sendo Δχ o incremento de curvaturas nodais da placa e CM a matriz elástica obtida a partir da equação (3).

2.1.4 Equações de Equilíbrio e Operador Tangente Consistente

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

(12) Δ K n Δ M n = 0

onde os vetores ΔKn e ΔMn são definidos, respectivamente, nas equações (9) e (10).

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

(13) R M ( Δ χ n ) = 2 Δ K n C M Δ χ n + S M ( C M Δ χ n Δ m n ) Δ m n = 0

Após aplicar à placa o incremento de curvaturas Δχn e obter a distribuição de tensões Δσij ao longo da espessura da placa, para todos os nós de células, se a equação (13) não for satisfeita, haverá um resíduo RM, isto é o incremento não será elástico. Neste 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 Δχni é conhecido. O próximo incremento de tentativa Δχni+1 na iteração (i+1) é obtido adicionado as correções δΔχni+1, ou seja:

(14) Δ χ n i + 1 = Δ χ n i + δ Δ χ n i + 1

onde as correções δΔχni+1 são calculadas linearizando-se a equação (13), sendo dadas por:

(15) δ Δ χ n i + 1 = [ R M ( Δ χ n i ) Δ χ n i ] 1 R M ( Δ χ n i )

onde RM(Δχni)Δχni é o operador tangente consistente obtido derivando-se a equação (13), resultando em:

(16) R M ( Δ χ n i ) Δ χ n i = S M ( C M n e p ( i ) C M ) + C M + C M n e p ( i )

Na equação (16), CMnep(i) é a matriz que contém as relações constitutivas [Cmep]ki (que relaciona momentos e curvaturas) referentes a todos nós de células. Para um nó k, [Cmep]ki é calculada integrando-se o tensor constitutivo [Cep]ki (que relaciona tensões e deformações) ao longo da espessura da placa th como segue:

(17) [ C m e p ] k i = m k i Δ χ k i = t h / 2 t h / 2 ( x 3 ) 2 [ C e p ] k i d x 3

sendo [Cep]ki=(σki)εki, que na análise em multi-escala é calculado após resolver o problema de equilíbrio do EVR; na análise não-linear convencional ele é calculado de acordo com o modelo constitutivo adotado.

Observe que os valores de momentos Δmij definidos na equação (1), assim como o tensor [Cmep]ki definido na equação (17), são calculados numericamente usando a fórmula e a quadratura de Gauss. Para isso, devem ser definidos pontos de Gauss ao longo da espessura da placa (veja [21[21] 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.], [22[22] FERNANDES, G. R., VENTURINI, W. S. “Building floor analysis by the Boundary element method”, Computational Mechanics, v. 35, pp. 277 - 291, 2005.], [23[23] FERNANDES, GR., KONDA, DH. “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.]). Assim, após o cálculo das curvaturas nodais na placa (equação (14)), o incremento de deformações Δεijg referente a um ponto de Gauss de um determinado nó de célula pode ser obtido. Então, esse incremento de deformações é imposto ao respectivo EVR e o incremento de tensão Δσijg obtido após resolver o problema de equilíbrio do EVR discutido adiante. Finalmente, a partir da distribuição de tensões ao longo da espessura, os momentos Δmijsão obtidos considerando-se a equação (1), que numericamente resultam em:

(18) Δ m i j t 2 4 g = 1 N g Δ σ i j ( g ) ξ g ω g i,j = 1 , 2

onde Ng é o número de pontos de Gauss ao longo da espessura, ξge ωg são, respectivamente, a coordenada adimensional e fator ponderador do ponto de Gauss.

O critério de convergência adotado para o problema de equilíbrio da placa é dado por RMTRM/ΔKMTΔKMtol, onde tol é a tolerância adotada no processo iterativo.

2.1.5 Definição do EVR e Passagem do Micro para o Macro-Contínuo

Seja o macro-contínuo representado na figura (1), que neste trabalho é definido pela placa, onde x é um ponto qualquer do macro-contínuo. Na análise em multi-escala, cada ponto x é representado pelo EVR (Elemento de Volume Representativo) (veja [10[10] 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.] - [12[12] 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.]), 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 dimensão, discretização e definição dos nós no EVR são completamente independentes das dimensões ou da 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 microestrutura

Note na figura (1) que o EVR pode ser composto por vazios (domínio Ωμv) e por partes sólidas (domínio ΩμS), sendo Ωμ = Ωμv ∪ Ωμ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.

Neste trabalho, o EVR também é considerado como meio contínuo e, portanto o conceito de tensão permanece válido a nível microscópico. Assim, a tensão microscópica pode ser escrita em termos de deformação como: σμ(y, t) = fyμ(y, t)), onde fy é 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, para um instante qualquer t, tem-se:

(19) ε ( x , t ) = 1 V μ Ω μ ε μ ( y , t ) d V
(20) σ ( x , t ) = 1 V μ Ω μ σ μ ( y , t ) d V

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:

(21) ε μ ( y , t ) = s u μ ( y , t )

onde ∇S é o operador gradiente simétrico.

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

(22) C e p ( x , t ) = σ ( x , t ) ε ( x , t ) = 1 V μ Ω μ S σ μ ( y , t ) d V ε ( x , t ) = 1 V μ Ω μ S f y ( ε μ ( y , t ) ) d V ε ( 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 (20) e (22). Note que a tensão macroscópica σ (equação (20)) e o tensor constitutivo do macro-contínuo Cep(equação (22)) são obtidos a partir de seus respectivos campos no EVR. Portanto, a definição de vazios ou de inclusões com diferentes materiais no interior do EVR afeta diretamente os valores dessas quantidades macroscópicas, mudando a resposta numérica da estrutura.

2.1.6 Campo de Deslocamento no EVR

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

(23) u μ ( y , t ) = ε ( x , t ) y + u ˜ μ ( y , t )

onde a parcela ε(x, t)y varia linearmente, sendo obtida multiplicando-se a deformação constante ε imposta pelo macro pelas coordenadas do ponto y do EVR; a parcela u˜μ é 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 u˜μ nulo. Analogamente, a deformação microscópica pode ser escrita como:

(24) ε μ ( y , t ) = ε ( x , t ) y + ε ˜ μ ( y , t )

onde ε é a deformação homogênea imposta pelo macro-contínuo e ε˜μ é o campo das flutuações de deformações, sendo definida como: ε˜μ(y,t)=Su˜μ.

A condição para se ter um campo de flutuações de deslocamentos cinematicamente admissível é que u˜μK˜μ*, sendo K˜μ* o espaço vetorial dos campos de flutuações de deslocamentos minimamente restringidos e cinematicamente admissíveis, definido como (ver mais detalhes em [12[12] 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.]):

(25) K ˜ μ * { v / Ω μ v S n d A = 0 }

onde faz-se uso da seguinte expressão usv=12(uv+vu), no cálculo das integrais, válida para quaisquer vetores u e v.

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

(26) ε ˙ μ ( y , t ) = S u ˙ μ = ε ˙ ( x , t ) + ε ˜ μ ( y , t ) u ˜ μ K ˜ μ

onde K˜μ é o conjunto dos campos de flutuações de deslocamentos microscópicos cinematicamente admissíveis, sendo K˜μK˜μ*.

2.2 Formulação do MEF para modelar a micro-escala

A microestrutura será modelada pelo Método dos Elementos Finitos, sendo a formulação desenvolvida detalhadamente nos trabalhos [10[10] 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.] - [12[12] 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.] e será apresentada aqui de forma resumida.

2.2.1 Equação de Equilíbrio do EVR

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 te = te (y, t) o campo de forças externas de superfície que atuam no contorno ∂Ωμ do EVR. O PTV estabelece que o EVR estará em equilíbrio se e somente se a seguinte equação variacional for satisfeita a cada instante t:

(27) Ω μ S σ μ ( y , t ) : S η d V Ω μ S b ( y , t ) . η d V + Ω μ v σ μ ( y , t ) : S η d V Ω μ v b ( y , t ) . η d V Ω μ t e ( y , t ) . η d A = 0 η V μ

onde 𝕍μ é o espaço de deslocamentos virtuais que satisfaz 𝕍μ = K˜μ 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 [12[12] 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.]), considerando-se ainda as equações (21) e (24) e escrevendo-se a tensão em termos das deformações, σμ = fyμ), onde fy é o tensor constitutivo, a equação (27) pode ser reescrita da seguinte forma em termos de flutuação de deslocamentos:

(28) Ω μ S f y ( ε ( x , t ) + S u ˜ μ ( y , t ) ) : S η d V = 0 η V μ

Finalmente, a formulação fica completa com a escolha do espaço 𝕍μ, 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 [12[12] 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.], 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. Porém, neste trabalho, apenas será adotada a condição de contorno onde se impõem flutuações periódicas ao longo do contorno. Em [17[17] 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.] todas as condições de contorno são discutidas de forma detalhada.

Assim, o problema de equilíbrio de EVR consiste em encontrar o campo de flutuações de deslocamentos u˜uVμ tal que, a cada instante t, a equação (28) seja satisfeita. 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 poΔtn = tn+1 - tn e discretização h, sendo u˜u(n+1)=u˜u(n)+Δu˜u(n):

(29) G h n + 1 = Ω μ h B T f y ( ε n + 1 + B u ˜ μ ( n + 1 ) ) d V = 0

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

Se o incremento n é não-linear, a equação (29) é resolvida aplicando-se o método de Newton-Raphson, que consiste em encontrar as correções de flutuações δu˜μi+1 na iteração i + 1, tal que:

(30) F i + K i δ u ˜ μ i + 1 = 0

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

(31) F i = Ω μ h B T f y ( ε n + 1 + B u ˜ μ i ) d V = e = 1 N e B e T σ μ e ( i ) V e
(32) K i = [ Ω μ h B T D μ i B d V ] = e = 1 N e B e T D μ e ( i ) B e V e

onde Ve é o volume de um elemento e qualquer e Dμe é o tensor constitutivo tangente do elemento e, definido como:

(33) D μ i = ( d f y d ε μ | ε μ = ε n + 1 + B u ˜ μ i + 1 )

Após o cálculo das correções δu˜μi+1 pela equação (30), o próximo campo de flutuação de deslocamentos de tentativa a ser considerado na iteração i + 1 referente à microestrutura, é dado por: u˜μi+1=u˜μi+δu˜μi+1.

2.2.2 Tensão Homogeneizada

A tensão homogeneizada é calculada pela equação (20), considerando-se que o EVR é composto de vazios e parte sólida (Ωμ = Ωμv ∪ ΩμS), resulta em:

(34) σ = σ ( x , t ) = 1 V μ Ω μ S σ μ ( y , t ) d V + 1 V μ Ω μ v σ ( y , t ) d V

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

(35) σ n + 1 = 1 V μ [ Ω μ h t n + 1 e s y d A Ω μ s ( h ) b n + 1 s y d V ]

onde faz-se uso da seguinte expressão usv=12(uv+vu), no cálculo das integrais, válida para quaisquer vetores u e v; tn+1e são as forças ao longo do contorno externo (definidas na equação (31)) e y é o vetor das coordenadas de um ponto genérico do EVR.

2.2.3 Condição de Contorno Imposta ao EVR

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 [17[17] 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.] e [25[25] 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, n.3, pp. 577 - 598, 2016.]. Note-se que cada condição de contorno leva a uma resposta numérica diferente, definindo-se assim, diferentes modelos em multi-escala.

O modelo com flutuações periódicas no contorno é adequado para descrever o comportamento de materiais que têm 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. Observe-se que nos EVRs definidos na figura (2) cada lado Γi+ corresponde a um lado igual e oposto Γi-, sendo ni+a a direção normal ao contorno Γi+ e ni- a direção normal ao contorno Γi-, com ni+ = - ni-. Assim, para cada ponto y+ pertencente ao contorno Γi+ existe um ponto correspondente y- do contorno Γi-.

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

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

(36) u ˜ u ( y + , t ) = u ˜ μ ( y , t ) { y + , y } Ω μ

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

(37) V μ = { u ˜ μ K ˜ μ * / u ˜ μ ( y + , t ) = u ˜ μ ( y , t ) { y + , y } Ω μ }

A fim de satisfazer o princípio de Hill-Mandel as forças te devem ser anti-periódicas no contorno, isto é: te(y+, t), = -te (y-, t) e as forças de volume b(y, t) devem ser nulas no domínio do EVR. Portanto, a tensão homogeneizada é dada pela equação (35), sendo o segundo termo nulo. Neste modelo o sistema definido em (30) pode ser escrito da seguinte forma:

(38) { F p F m F I } i + [ K p p K p m K p I K m p K m m K m I K I p K I m K I I ] i { δ u ˜ p δ u ˜ m δ u ˜ I } μ i + 1 = 0

onde os índices p, m e I definem, respectivamente, pontos y+, y- e internos.

Considerando que δu˜m=δu˜p, o sistema dado pela equação (38) pode ser reduzido e representado da seguinte forma:

(39) { F p + F m F I } i + [ K p p + K p m + K m p + K m m K p I + K m I K I p + K I m K I I ] i { δ u ˜ p δ u ˜ I } μ i + 1 = 0

isto é, neste 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.

A equação (39) pode ainda ser escrita na seguinte forma geral:

(40) { δ u ˜ R } i + 1 = [ K R i ] 1 { F R } i

onde os vetores {δu˜R}i+1 e {FR}i assim como a matriz [KR]i são definidos de acordo com o modelo em multi-escala (ou condição de contorno). Para flutuações periódicas no contorno, eles podem facilmente ser obtidos da equação (39).

2.2.4 Tensor Constitutivo Homogeneizado

Considerando-se a equação (24), após a discretização do EVR em elementos, para uma iteração i de um incremento n de tempo (sendo Δtn = tn + 1 - tn), o tensor constitutivo tangente macroscópico ou homogeneizado Cep é obtido a partir da equação (22), a qual pode ser escrita como: (ver mais detalhes em [12[12] 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.]):

(41) C i e p = σ i ε i = 1 V μ Ω μ h f y ( ε n + 1 + S u ˜ μ i ) d V ε i

onde Su˜μi é obtido com a solução do problema de equilíbrio do EVR definido na equação (29).

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

(42) C i e p = C i e p ( T a y l o r ) + C ˜ i e p

onde Ciep(Taylor) é definido como o operador tangente do modelo de Taylor (obtido adotando-se Su˜μ(n+1)=0), sendo calculado pela média volumétrica do tensor constitutivo microscópico como segue:

(43) C i e p ( T a y l o r ) = 1 V μ Ω μ h f y ( ε n + 1 ) d V ε i = 1 V μ Ω μ h f y ( ε n + 1 ) ε μ d V = 1 V μ Ω μ h D μ i d V = p = 1 N p V p V μ D μ p ( i )

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

Note-se 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.

A outra parte C˜iep da equação (42) representa a influência da flutuação dos deslocamentos no valor do tensor homogeneizado (ver mais detalhes em [10[10] 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.]). Considerando-se uma discretização em Ne elementos, a mesma é definida por:

(44) C ˜ i e p = 1 V μ Ω μ h g y ( S u ˜ μ i ) d V ε i = 1 V μ G R i K R i 1 G R i T

onde GR é definido de acordo com o modelo em multi-escala; para flutuações periódicas é dado por: GRi = [Gp + Gm G1]i, sendo G definida como: G=e=1NeDμeBeVe.

2.2.5 Algoritmo

É importante notar que para se calcular o vetor de momentos elásticos ΔKn (definido na equação (9)), um módulo de Young E e um coeficiente de Poisson v devem ser definidos para o macro-contínuo. 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 v da placa) serão adotados como a média volumétrica dos respectivos valores na microestrutura, ou seja:

(45) E = p = 1 N p V p V μ E p
(46) v = p = 1 N p V p V μ v p

onde Np é o número de fases definidas no EVR.

Observe-se que as equações (45) e (46) 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 iEVR (sendo i ≥ 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 ΔKn (definida na equação (9)).

  2. Para cada nó da placa:

    • 2.1

      Incremento de curvatura Δχn0 é obtido com a lei de Hooke (equação (3)).

    • 2.2

      Um EVR é definido para cada ponto de Gauss definido ao longo da espessura da placa. O EVR é resolvido seguindo o procedimento:

      • 2.2.1

        As macro-deformações são impostas ao EVR: Δεni=x3Δχni.

      • 2.2.2

        Para iPVE = 0, as flutuações de deslocamentos são nulas e, portanto, os deslocamentos no EVR (equação23) são: {Δuμt}n0=[Δε]ni{x1x2}y.

      • 2.2.3

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

        1. Valores nodais das forças {F}n+1iRVE (equação (31)). Para isso, os seguintes valores devem ser computados para cada elemento finito e:

          • vetor do incremento de deformação de tentativa: {Δεμe(t)}niRVE=Be{Δuμe(t)}niRVE,

          • vetor do incremento de tensão de tentativa {Δσμe(t)}niRVE, usando a lei de Hooke,

          • vetor de tensão de tentativa: {Δσμe(t)}n+1iRVE+1=σne+{Δσμe(t)}niRVE.

          • verifica-se o modelo constitutivo, obtendo-se a tensão {σμe}n+1iRVE+1 e o tensor constitutivo [Dμe]niRVE (todas as variáveis são atualizadas tendo como base os valores do incremento anterior (ver em [26[26] SOUZA NETO, E. A., PERIC, D., OWEN, D. R. J. Computational methods for plasticity: Theory and applications, Wiley, Chichester, 814 pp, 2008.])).

        2. O vetor {FR}n+1iRVE (equação (40)) é obtido de acordo com o modelo em multi-escala.

        3. Verifica-se a convergência:

          FRTFRFPTFPtol, onde Fp são as forças nodais no lado Γ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.

        4. Obtêm-se a matriz de rigidez do EVR (equação (32)) e a matriz [KR](n)iREV (equação (40));

        5. Resolve-se o sistema dado pela equação (40) e obtém as correções de flutuações {δu˜R}niREV+1;

        6. Obtém-se o vetor total {δu˜}niREV+1.

        7. Calcula-se o novo incremento de deslocamento a ser imposto ao EVR:

          {Δuμt}niREV+1=[Δε]ni{x1x2}y+{Δu˜u}niREV+1, onde {Δu˜μ}niREV+1={Δu˜μ}niREV+{δu˜}niREV+1 e retorna-se ao passo (a).

      • 2.2.4

        Obtém-se o tensor constitutivo homogeneizado [Cep]ni (equação (20)).

      • 2.2.5

        Obtém-se a tensão homogeneizada σn+1i+1 (equação (35)), onde o segundo termo é nulo para flutuações periódicas; calcula-se o incremento de tensãoΔσni=Δσn+1i+1σn1 do macro.

      • 2.2.6

        Obtêm-se o incremento de momento Δmni (equação (18)) e o tensor constitutivo [Cmep]ni (equação (17)).

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

    • 3.1

      Se o resíduo RM não for nulo, a matriz tangente (equação (16)) é atualizada e as ções δΔχni+1 (equação (15)) são calculadas. Atualiza-se o incremento de curvatura Δχni+1 (equação (14)) a ser aplicado na próxima iteração e retorna-se 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 dos momentos plásticos (equação (11)) e os valores nodais de deslocamentos, forças e esforços internos (equações (8) e (10)).

  5. Retorna-se ao passo 1 e aplica-se novo incremento n + 1.

3. RESULTADOS E DISCUSSÃO

Como exemplo inicial de aplicação da formulação proposta, adota-se uma placa quadrada cujos lados medem 1m (ver figura (3a)), sendo a mesma apoiada nas quatro bordas e submetida a uma carga uniformemente distribuída de 1 N/cm2. Além disso, adotou-se para a mesma 4 cm de espessura e a discretização definida na figura (3a), onde têm-se 100 nós no contorno, 25 pontos internos, 72 células e 48 elementos no contorno. Considerando-se essa placa são realizadas diferentes análises numéricas.

Figura 3
a) Discretização da placa: 100 nós no contorno, 72 células b) Discretização do EVR com inclusões elásticas.

Neste exemplo objetiva-se mostrar que a formulação proposta é capaz de capturar os diferentes comportamentos da estrutura com a manipulação da microestrutura heterogênea. Para tanto, inicialmente, os resultados numéricos são obtidos adotando-se uma análise não linear convencional, onde um material homogêneo é utilizado a fim de comparar com aqueles referentes a uma análise não-linear multi-escala. Nesse segundo caso, a microestrutura do material é modelada inserindo inclusões mais rígidas na matriz, o que aumenta a rigidez da estrutura.

Para a análise não-linear convencional adotou-se que a estrutura é composta por um material que obedece o modelo de Von Mises com haderning isotrópico com as seguintes propriedades: módulo de elasticidade E = 2000 kN/cm2, coeficiente de Poisson v = 0,2; tensão de escoamento de 70 N/cm2 e módulo de encruamento K= 430 N/cm2. Já na análise multi-escala, o material tem sua microestrutura representada por um EVR onde cinco inclusões elásticas são inseridas na matriz metálica, cujas propriedades são: v = 0,35 e E = 4000 kN/cm2, enquanto que para a matriz do EVR é definida com as mesmas propriedades do material homogêneo da análise não-linear convencional. Condições periódicas foram utilizadas na modelagem do EVR. Na discretização do EVR adotaram-se 520 elementos triangulares e 293 nós (veja figura (3b)) e como condição de contorno no mesmo foi considerada flutuações periódicas. Vale ressaltar que na formulação aqui desenvolvida, cada EVR representa um ponto do macro e, portanto, suas dimensões não são importantes, porém é fundamental o conhecimento da proporcionalidade entre as fases constituintes e a localização dos constituintes. Ou seja, os resultados numéricos não mudam em função das dimensões adotadas para os lados do EVR.

Esse tipo de geometria de EVR pode representar materiais como o concreto, onde inclusões elásticas mais rígidas (agregados) são adicionadas à matriz (argamassa), de comportamento mais flexível. Também se pode pensar na modelagem de materiais compósitos com matriz metálica reforçada (CMM), [25[25] 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, n.3, pp. 577 - 598, 2016.]. Contudo, focando em materiais como o concreto, para modelar o comportamento da matriz, além do modelo de Von Mises, é adotado também o modelo de Mohr-Coulomb, a fim de mostrar que os resultados mudam de acordo com o modelo constitutivo adotado, o que ressalta a importância de adotar um modelo constitutivo que melhor represente o comportamento do material. No entanto, para modelar o comportamento da argamassa, que é um material frágil, o modelo de Mohr-Coulomb seria mais indicado por considerar uma limitação de resistência para estados predominantes de tração. Para o modelo de Mohr Coulomb foram adotados os mesmos parâmetros já citados para o modelo de Von Mises, além de ser definido o ângulo de dilatação de 20° e ângulo de fricção de 10°, o que leva a uma regra não-associativa, [9[9] PITUBA, J. J. C., SOUZA NETO, E. A. “Modeling of unilateral effect in brittle materials by a mesoscopic scale approach”, Computers and Concrete, An International Journal, v. 15, n.5, pp. 735-758, 2015.].

A figura (4) apresenta o deslocamento do ponto central da placa ao longo do processo incremental, considerando-se os diferentes tipos de análises. Observa-se que os resultados das análises multi-escala estão bem próximos e os resultados mostram um comportamento mais rígido da estrutura, devido à consideração de inclusões elásticas na microestrutura do material. Os resultados da análise não-linear convencional foram bem diferentes a partir do ponto que houve plastificação, o que era esperado devido a menor rigidez do material homogêneo. Focando no caso das análises multi-escala, a estrutura evidenciou processos de plastificação devido ao emprego dos modelos constitutivos mencionados, entretanto, com a consideração de perfeita aderência entre agregado e matriz, sendo os agregados mais rígidos que a matriz, foi evidenciada uma rigidez maior da estrutura mesmo com o processo de plastificação da matriz, o que mostra que a inserção de inclusões rígidas na microestrutura do material leva a um ganho de resistência e rigidez. Também percebe-se que a inserção de inclusões rígidas na matriz limita o comportamento dúctil do material quando comparado às análises convencionais, principalmente no caso do emprego do modelo de Von Mises (abreviado como VM na figura (4)). Ainda discutindo as análises convencionais, com o modelo de Mohr Coulomb (abreviado como MC na figura (4)) a estrutura resiste a maiores carregamentos, porém apresenta menores deslocamentos caracterizando materiais frágeis enquanto o modelo de Von Misses apresenta um grande aumento no valor do deslocamento com pequeno incremento de carga, caracterizando o comportamento de materiais mais dúcteis.

Figura 4
Deslocamentos no ponto central da placa considerando-se diferentes tipos de análises.

Analisando agora a placa definida na figura (5) para representar o macro-contínuo, o EVR considerado é composto por uma matriz metálica, de comportamento elasto-plástico, onde se define uma inclusão elástica, cuja área será modificada a fim de demonstrar como o aumento de sua área afeta diretamente a resistência e rigidez do material e, consequentemente, o comportamento da estrutura. Esse tipo de EVR pode representar materiais como o CMM, onde inclusões são adicionadas ao material a fim de melhorar suas propriedades elásticas ligadas à rigidez e à resistência.

Figura 5
a) Condições de Contorno b) Discretização da placa.

A figura (5a) ilustra a placa analisada, onde os lados menores são considerados livres enquanto que os outros dois lados são engastados. Para a definição da geometria da placa adotou-se: espessura th=4,0 cm, a=200,0 cm e b=100,0 cm. Como carregamento, considerou-se uma carga g=3,5 kN/cm2 uniformemente distribuída sobre a área central definida na figura, a fim de representar uma carga concentrada de 787,5 kN aplicada no centro da placa. Na discretização da placa foram adotados 36 elementos no contorno e 144 células no domínio como definido na figura (5b).

Para modelar o comportamento elasto-plástico da matriz será adotado o critério de Von Mises, sendo adotado: encruamento K=2000 kN/cm2, coeficiente de Poisson v = 0,3; Módulo de Elasticidade E = 20000 kN/cm2, tensão de escoamento σy = 40,0 kN/cm2. Por outro lado, as seguintes propriedades elásticas serão adotadas para as inclusões: v = 0,3 e E = 160000 kN/cm2. Serão considerados três diferentes tamanhos para a inclusão: o primeiro EVR possui uma fração de volume de inclusão de 5%, onde foram adotados 720 elementos finitos triangulares e 393 nós na discretização (veja figura (6)); para o segundo EVR adotou-se fração de volume de inclusão de 10%, com 584 elementos finitos triangulares e 325 nós (figura (7)); finalmente, para o terceiro EVR foi utilizada uma fração volumétrica de inclusão de 30%, com 572 elementos finitos triangulares e 319 nós (figura (8)). Condições periódicas foram utilizadas na modelagem do EVR.

Figura 6
Discretização do EVR com fração de volume de 5% para a inclusão
Figura 7
Discretização do EVR com fração de volume de 10% para a inclusão
Figura 08
Discretização do EVR com fração de volume de 30% para a inclusão

No que segue são apresentadas as respostas dos EVRs em termos de tensão homogeneizada na direção x quando submetidos a uma deformação macroscópica, tendo como objetivo a simulação de estados predominantes de tração e de compressão, evidenciando a diferença de comportamento mecânico da microestrutura heterogênea. Para o caso de estados predominantes de tração foi imposta a deformação macroscópica ε = [εx; εy; γxy] = [0,005; −0,000695; 0,0]. Já para o caso de estados predominantes de compressão, foi imposta a deformação macroscópica ε = [εx; εy; γxy] = [-0,05; 0,000695; 0,0]. Nos dois casos, as deformações macroscópicas foram aplicadas em 10 incrementos utilizando uma tolerância de 10-8 para verificação do processo incremental-iterativo. Condições periódicas foram utilizadas na modelagem do EVR.

No caso de estados predominantes de compressão, a figura (9) ilustra as respostas obtidas. Observa-se que o aumento da fração volumétrica de inclusão (reforço da matriz metálica) confere um aumento da rigidez e resistência da microestrutura do material. Isso é refletido no comportamento mecânico da macroestrutura ilustrada na figura (10). Observa-se ainda que o aumento da fração volumétrica de inclusão também evidencia sua influência no comportamento elástico do EVR, onde a plasticidade da matriz fica menos evidente. Em particular, a resposta do EVR com 5% de inclusão mostra um comportamento mecânico elástico até 42 kN/cm2 (420 MPa) quando há o início dos processos de plastificação da matriz. Já o aumento da fração volumétrica da inclusão mostra valores menores de tensão para início da plastificação, isso porque a presença da inclusão rígida faz com que o nível de tensão no EVR seja aumentado.

Figura 9
Comparação dos resultados de Tensão homogeneizada de compressão na direção x (σx) do EVR
Figura 10
Deslocamento no ponto central da placa ao longo do processo incremental.

versus deformação macroscópica de encurtamento imposta na direção x (εx) para os casos de fração volumétrica de inclusão de 5%, 10% e 30%. Os valores foram considerados positivos para melhor visualização.

No caso de estados predominantes de tração, as respostas obtidas foram exatamente as mesmas obtidas para estados predominantes em compressão. Isso é devido ao uso do modelo de Von Mises que considera que os materiais têm o mesmo comportamento mecânico em tração e compressão.

A figura (10) apresenta as respostas obtidas com o emprego das análises multi-escala da placa, sendo analisado o deslocamento no ponto central da placa com a evolução do fator de carga. Observa-se que o ganho de rigidez demonstrado nas análises da microestrutura (figura (9)) é refletido no comportamento mecânico do macro-contínuo. É interessante constatar como a manipulação da microestrutura, inserindo inclusões com maior rigidez na matriz elasto-plástica, leva a um comportamento desejável da estrutura, que em nosso caso, procurou-se aumentar a rigidez da estrutura com a manipulação do material homogêneo, conferindo-lhe maior rigidez e, em contrapartida, heterogeneidade. Portanto, com a formulação multi-escala é possível analisar estruturas compostas por materiais complexos em sua microestrutura recorrendo a modelos constitutivos simples em formulação e com poucos parâmetros, sem a necessidade de recorrer a modelos constitutivos fenomenológicos com alto grau de complexidade. Em suma, os resultados apresentados neste trabalho estão em concordância com respostas esperadas, as quais também estão relatadas em [7[7] NGUYEN, V.P., LLOBERAS VALLS, O., STROEVEN, M., et al, “Homogenization-based multiscale crack modelling: from micro-diffusive damage to macro-cracks”, Computer Methods in Applied Mechanics and Engineering, v. 200, n.9–12, pp. 1220–1236, 2011.], [10[10] 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.], [11[11] 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 Computations, v. 26, n.3, pp. 281-301, 2009.] e [14[14] 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.], tanto em análises desacopladas da microestrutura apenas, como também em análises envolvendo a micro e macroestrutura de forma totalmente acoplada.

4. CONCLUSÕES

Neste trabalho foi apresentada a formulação multi-escala para análise de flexão de placas, onde o MEC foi usado para modelar a macro-escala e o MEF adotado para modelar a micro-escala. Neste modelo cada nó de célula da placa é representado por um EVR, cujo problema de equilíbrio deve ser resolvido para se obter a tensão e o tensor constitutivo do nó da placa através de um processo de homogeneização.

Foram apresentadas análises numéricas, onde os resultados obtidos com a análise não-linear convencional foram comparados com aqueles obtidos com as respostas do emprego da formulação multi-escala. Adotaram-se diferentes modelos constitutivos para representar o comportamento do material, tendo definidas inclusões elásticas no EVR a fim de enrijecer a estrutura. O modelo de Von Mises é ideal para representar o comportamento de materiais dúcteis, os quais apresentam valores consideráveis de deformação plástica antes de atingir o estado de colapso, sendo o modelo de Mohr Coulomb mais indicado para materiais frágeis, que rompem sem apresentar deformações plásticas muito consideráveis. Contudo, o emprego desses modelos, associados à geometria e à fração volumétrica dos constituintes da microestrutura, levam à capacidade de obtenção de respostas numéricas de estruturas compostas por materiais heterogêneos ao se utilizar a formulação multi-escala.

Os resultados obtidos são coerentes, isto é, obteve-se resposta mais rígida com a análise em multi-escala se comparada à análise não-linear convencional e uma curva mais flexível quando se adotou o modelo de Von Mises para modelar a matriz da microestrutura. Além disso, é importante dizer que o programa mostrou ser estável, sendo a convergência do processo incremental-iterativo sempre alcançada com poucas iterações. Daí a importância da formulação tangente consistente apresentada neste trabalho, principalmente numa análise multi-escala, onde o custo computacional é um limitador de sua utilização para grandes problemas.

Por outro lado, os exemplos numéricos evidenciaram a capacidade da modelagem multi-escala em capturar a influência da modificação das propriedades da microestrutura na rigidez do macro-contínuo, evidenciando que a validação qualitativa da formulação proposta foi alcançada. Acredita-se que esse é um passo inicial antes de aplicações e comparações com respostas experimentais em estruturas compostas por materiais heterogêneos.

Por fim, os resultados apresentados até então são promissores. A formulação apresentada aqui está em aperfeiçoamento para se considerar o processo de fraturamento na microestrutura com o intuito de modelar o descolamento da interface inclusão/matriz em CMMs (ver os trabalhos [6[6] PITUBA, J. J. C., FERNANDES, G. R., SOUZA NETO, E. A., “Modeling of cohesive fracture and plasticity processes in composite microstructures”, Journal of Engineering Mechaniccs-ASCE, v. 142, n.10, DOI: 10.1061/(ASCE)EM.1943-7889.0001123, 2016.
https://doi.org/10.1061/(ASCE)EM.1943-78...
] e [25[25] 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, n.3, pp. 577 - 598, 2016.]) e sua influência no macro-continuo. Já no caso de materiais frágeis como o concreto, o surgimento e localização do processo de danificação na microestrutura e sua posterior transição para o macro-contínuo em forma de fraturamento, é uma linha de pesquisa ainda em aberto, onde alguns trabalhos foram publicados, [27[27] RODRIGUES, E. A., MANZOLI, O. L., BITENCOURT JR, L. A. G., et al., “2D mesoscale model for concrete based on the use of interface element with a high aspect ratio”, International Journal of Solids and Structures, v. 94-95, pp. 112-124, 2016.] e [28[28] TORO, S., SANCHEZ, P. J. S., BLANCO, P. J., et al., “Multiscale formulation for material failure accounting for cohesive cracks at the macro and micro scales”, International Journal of Plasticity, v. 76, pp. 75-110, 2016.].

5. AGRADECIMENTO

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

6. BIBLIOGRAFIA

  • [1]
    PITUBA, J. J. C., FERNANDES, G. R., “An anisotropic damage model for concrete”, Journal of Engineering Mechanics-ASCE, v. 137, n.9, pp. 610-624, 2011.
  • [2]
    PITUBA, J. J. C., “Anisotropic damage model on the effects of damage process due to shearing stress in concrete”, Acta Scientiarum: Technology, v. 35, n.2, pp. 227-236, 2013.
  • [3]
    PITUBA, J. J. C., “A damage model formulation: unilateral effect and RC structures analysis”, Computers and Concrete, An International Journal, v. 15, n.5, pp. 709-733, 2015.
  • [4]
    MATALLAH, M., LA BORDERIE, C., “Inelasticity–damage-based model for numerical modeling of concrete cracking”, Engineering Fracture Mechanics, v. 76, pp. 1087-1108, 2009.
  • [5]
    GAL, E., KRYVORUK, R, “Fiber reinforced concrete properties – a multiscale approach”, Computers and Concrete, An International Journal, v. 8, n.5, pp. 525-539, 2011.
  • [6]
    PITUBA, J. J. C., FERNANDES, G. R., SOUZA NETO, E. A., “Modeling of cohesive fracture and plasticity processes in composite microstructures”, Journal of Engineering Mechaniccs-ASCE, v. 142, n.10, DOI: 10.1061/(ASCE)EM.1943-7889.0001123, 2016.
    » https://doi.org/10.1061/(ASCE)EM.1943-7889.0001123
  • [7]
    NGUYEN, V.P., LLOBERAS VALLS, O., STROEVEN, M., et al, “Homogenization-based multiscale crack modelling: from micro-diffusive damage to macro-cracks”, Computer Methods in Applied Mechanics and Engineering, v. 200, n.9–12, pp. 1220–1236, 2011.
  • [8]
    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.
  • [9]
    PITUBA, J. J. C., SOUZA NETO, E. A. “Modeling of unilateral effect in brittle materials by a mesoscopic scale approach”, Computers and Concrete, An International Journal, v. 15, n.5, pp. 735-758, 2015.
  • [10]
    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.
  • [11]
    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 Computations, v. 26, n.3, pp. 281-301, 2009.
  • [12]
    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.
  • [13]
    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.
  • [14]
    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.
  • [15]
    SAAVEDRA-FLORES, E. I., SOUZA NETO, E. A., PEARCE, C., “A large strain computational multi-scale model for the dissipative behavior of wood cell-wall”. Computational Materials Science, v. 50(3): pp. 1202-1211, 2011.
  • [16]
    FERNANDES, G. R., SOUZA NETO, E. A., “Self-consistent linearization of non-linear BEM formulations with quadratic convergence”, Computational Mechanics DOI: 10.1007/s00466-013-0867-2, v. 52, pp. 1125-1139, 2013.
    » https://doi.org/10.1007/s00466-013-0867-2
  • [17]
    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.
  • [18]
    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.
  • [19]
    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.
  • [20]
    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.
  • [21]
    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.
  • [22]
    FERNANDES, G. R., VENTURINI, W. S. “Building floor analysis by the Boundary element method”, Computational Mechanics, v. 35, pp. 277 - 291, 2005.
  • [23]
    FERNANDES, GR., KONDA, DH. “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.
  • [24]
    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.
  • [25]
    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, n.3, pp. 577 - 598, 2016.
  • [26]
    SOUZA NETO, E. A., PERIC, D., OWEN, D. R. J. Computational methods for plasticity: Theory and applications, Wiley, Chichester, 814 pp, 2008.
  • [27]
    RODRIGUES, E. A., MANZOLI, O. L., BITENCOURT JR, L. A. G., et al, “2D mesoscale model for concrete based on the use of interface element with a high aspect ratio”, International Journal of Solids and Structures, v. 94-95, pp. 112-124, 2016.
  • [28]
    TORO, S., SANCHEZ, P. J. S., BLANCO, P. J., et al, “Multiscale formulation for material failure accounting for cohesive cracks at the macro and micro scales”, International Journal of Plasticity, v. 76, pp. 75-110, 2016.

Datas de Publicação

  • Publicação nesta coleção
    2017

Histórico

  • Recebido
    09 Jul 2016
  • Aceito
    02 Nov 2016
Laboratório de Hidrogênio, Coppe - Universidade Federal do Rio de Janeiro, em cooperação com a Associação Brasileira do Hidrogênio, ABH2 Av. Moniz Aragão, 207, 21941-594, Rio de Janeiro, RJ, Brasil, Tel: +55 (21) 3938-8791 - Rio de Janeiro - RJ - Brazil
E-mail: revmateria@gmail.com