Simulation of subcritical crack propagation in quase-brittle materials applying a version of the discrete element method formed by beams



==Resumo==

Há numerosas aplicações de interesse tecnológico onde materiais de comportamento quase frágil, como é o caso de materiais cimentícios, rochas e materiais compostos formados pela mistura de cerâmicas com outras fases, são submetidos a cargas oscilantes e sofrem degradação das suas propriedades. Desenvolver ferramentas de análise que permitam entender e prever os mecanismos que governam esta degradação é um problema aberto na engenharia moderna. Neste contexto, no presente trabalho, se utiliza uma versão do método dos elementos discretos (DEM) formado por barras para explorar as possibilidades do mesmo na simulação dos mecanismos de degradação que ocorrem em materiais quase frágeis quando submetidos à fadiga. Simulações sobre corpos de prova de geometria simples são apresentadas e vários aspectos deste problema são discutidos, entre eles: é explorada a relação entre os parâmetros micro e macromecânicos do modelo empregado e como o efeito de escala é capturado e influi nesta relação. Os resultados preliminares apresentados deixam em evidência a potencialidade da metodologia proposta para compreender os micromecanismos de dano que ocorrem nos materiais quase frágeis e também para predizer sua evolução.

==Abstract==

There are numerous applications, of technological interest, where materials of quasi‐brittle behavior, as in the case of cement‐based composite materials, rocks and composites formed by the mixture of ceramics with other phases, are subjected to oscillating loads and suffer degradation of its properties. To develop analysis tools that allow to understand and to predict this degradation governing mechanisms is an opened problem in modern engineering. In this context, in the present work, it is used a version of the discrete element method formed by beams to explore its possibilities in simulating the degradation mechanisms that occur in quasi‐brittle materials when subjected to fatigue. Simulations over simple geometry test subjects are presented and several aspects of this problem are discussed and among the problems: it is explored the relationship between micro and macro mechanic parameters of the used model and how the scale effect is captured and influences this relationship. The presented preliminary results show the potentiality of the proposed methodology to understand the micro mechanisms of damage that occur in quasi‐brittle materials and also to predict its evolution.

==Palavras‐chave==

Fatiga ;                             Materiais quase frágeis ;                             Metodo dos elementos discretos

==Keywords==

Fatigue ;                             Quasi‐brittle materials ;                             Discrete elements methods

==1. Introdução==

Na engenharia moderna, em várias situações de interesse tecnológico, materiais com comportamento quase frágil são submetidos à ação de cargas oscilantes, entre eles: estruturas convencionais como estradas, pontes ou edifícios construídos com concreto e concreto reforçado, estruturas feitas em rocha como galerias e túneis, e também estruturas realizadas com cerâmicas ou materiais compostos utilizados em dispositivos de engenharia de alta exigência mecânica.

O fenômeno de fadiga tem sido largamente estudado especialmente em metais, existindo uma tradição de muitos anos na sua avaliação, uma grande quantidade de informação experimental permite hoje contar com metodologias para avaliar e predizer o tempo em que defeitos podem propagar devido à fadiga. Isto se deve a ser o fenômeno da fadiga um dos fatores mais importantes que produz degradação de componentes metálicos.

Quando o material estudado é quase frágil, as características que governam a degradação são diferentes e muito menos estudos relacionados com fadiga existem neste caso.

Métodos que permitam prever o comportamento de materiais, levando em conta os diferentes tipos de não linearidades envolvidos, estão disponíveis em sistemas de elementos finitos comerciais, ferramenta mais utilizada na modelagem de estruturas, mas no que se refere à determinação da vida em fadiga, os métodos mais comuns consistem em introduzir regras empíricas ou semiempíricas, entre eles se destacam 2 procedimentos típicos, aqueles que permitem prever a nucleação de um defeito (descrita pela relação entre tensão última e número de ciclos aplicados, a lei clássica de Basquin apud Moura Branco [[#bib0200|[1]]] ) e os que se baseiam na propagação subcrítica de um defeito já nucleado (metodologias fundamentadas na lei proposta por Paris [[#bib0205|[2]]] ). Existe uma farta bibliografia que descreve as características destas 2 abordagens, aplicando‐as às particularidades próprias de cada material. Como exemplo de referência bibliográfica, se pode citar Suresh [[#bib0210|[3]]] .      

Metodologias de cálculo baseadas na mecânica dos meios contínuos têm tido uma sensível evolução nos últimos anos, tornando‐se a ferramenta mais versátil para modelar problemas de engenharia. Como exemplo para ilustrar este enfoque, podemos citar os trabalhos de Saintier et al. [[#bib0215|[4]]]  e Infante et al. [[#bib0220|[5]]] .      

Por outro lado, para simular a possibilidade de nuclear defeitos dentro do contínuo, 2 procedimentos podem ser citados: a utilização de elementos finitos com funções de interpolação estendidas, que permitem introduzir singularidades dentro do elemento, os chamados elementos finitos estendidos. Como bibliografia básica sobre este tipo de elementos podemos citar o livro de Mohammadi [[#bib0225|[6]]]  e como exemplo de aplicação relacionada a simular a propagação de fissuras em estruturas de concreto podemos citar o trabalho de Unger et al. [[#bib0230|[7]]] . Outra forma de encarar o problema é utilizar a técnica das interfaces coesivas combinadas ao método dos elementos finitos. Esta técnica foi implementada por Xu e Nedeelman [[#bib0235|[8]]]  em 1993 e basicamente consiste em introduzir uma lei de interface que permita o «descolamento» entre elementos. Como exemplo de aplicações utilizando esta técnica na simulação de fadiga em materiais quase frágeis pode‐se mencionar os trabalhos de Yang et al. [[#bib0240|[9]]] , Yamaguchi et al. [[#bib0245|[10]]] , Siegmund [[#bib0250|[11]]]  e Xu e Yuan [[#bib0255|[12]]] , em todos os casos a geração de fissuras e sua propagação devido a cargas oscilantes pôde ser prevista através da utilização desta técnica.      

Há métodos alternativos para modelar o comportamento de sólidos quase frágeis, nos quais fenômenos como a interação de «nuvens» de microfissuras e posterior localização, assim como outros efeitos, podem ser modelados com relativa facilidade dentro da mecânica do descontínuo, assim denominada por Munjiza [[#bib0260|[13]]] . Várias alternativas do DEM são plausíveis de ser utilizadas no contexto de materiais frágeis, cabendo destacar a respeito os trabalhos de Schlangen e Van Mier [[#bib0265|[14]]]  e Rinaldi [[#bib0270|[15]]] , nestes casos o DEM é utilizado para simular o processo de dano frente a cargas crescentes, no que diz respeito à simulação do fenômeno de fadiga em materiais frágeis é possível mencionar o review realizado por Hansen et al. [[#bib0275|[16]]]  sobre aplicações dos chamados «Bundle Models», sendo que entre as aplicações apresentadas figura a simulação de fadiga em materiais frágeis. Também se destaca o trabalho de Rinaldi et al. [[#bib0280|[17]]]  que aplica leis estabelecidas de propagação de fadiga a nível microscópico estudando qual é o comportamento macroscópico resultante.      

Neste contexto, no presente trabalho se pretende explorar a propagação subcrítica de defeitos nucleados em um modelo de material quase frágil, sob cargas oscilantes, sem embutir leis específicas que induzam comportamentos de fadiga determinados. As características do método utilizado são apresentadas na seção 3 deste trabalho.

==2. Lei de Paris==

Como apresentado por Anderson [[#bib0285|[18]]] , a velocidade de crescimento de fissuras em forma subcrítica pode ser relacionada ao quociente e a diferença do nível de tensões trativas na região onde a trinca se desenvolve. De forma geral pode‐se escrever que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{da}{dN}=f\left(S_{m\acute{a}x}/S_{min},S_{m\acute{a}x}-S_{min},h\ldots \right)</math>
|}
| style="text-align: right;" | ( 1)
|}

onde ''S''  representa a tensão trativa na região da trinca, ''h''  representa a função que leva em conta o histórico de tensões preexistentes e ''N''  representa o número de ciclos. A lei proposta por Paris  [[#bib0205|[2]]]  é um caso particular de (1) que se expressa como segue:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{da}{dN}=C\Delta K^m</math>
|}
| style="text-align: right;" | ( 2)
|}

onde ''C''  e ''m''  são constantes que dependem do material e da relação entre tensões ''S''<sub>''m'' á''x''</sub> /''S''<sub>min</sub>  adotada. Graficando esta lei em escala bilogarítmica se pode representar por uma reta como se ilustra na região II  da [[#fig0005|figura 1]] , onde se verifica que ''m''  seria o coeficiente angular e ''log C''  sua translação no eixo das ordenadas. Na proposta de Paris, apresentada na expressão (2), ''ΔK''  representa o incremento do fator de intensidade de tensões que, segundo a mecânica linear da fratura, pode ser escrito como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Delta K=\Delta \sigma \sqrt{\pi \alpha }\beta </math>
|}
| style="text-align: right;" | ( 3)
|}

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr1.jpg|center|284px|Lei de crescimento subcrítico típica de uma fissura por fadiga.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Lei de crescimento subcrítico típica de uma fissura por fadiga.

</span>
|}

onde ''Δσ''  indica o incremento das tensões num ciclo de carga, ''a''  indica o tamanho da fissura e ''β''  é um coeficiente adimensional que representa a forma da fissura e condições de contorno atuantes.      

Na [[#fig0005|figura 1]]  se apresenta em escala bilogarítmica uma curva típica de comportamento de propagação subcrítica de uma fissura já nucleada, sendo que ela só cresce ao se exercer sobre a trinca um ''ΔK''  superior a um valor mínimo, parâmetro que depende do material e da relação ''S''<sub>''min''</sub>''/S''<sub>''max''</sub> , superando este valor mínimo a trinca propaga em forma subcrítica, incrementando sua velocidade até entrar em regime na região II , modelada pela Lei de Paris. Quando a fissura vai chegando à sua condição crítica acelera seu crescimento na região III , que está fora da região abrangida pelo modelo de Paris. Outros modelos mais sofisticados, que permitem simular o comportamento em todos os estágios do crescimento da fissura, estão disponíveis na bibliografia, consultar a este respeito Anderson [[#bib0285|[18]]] .      

==3. Descrição do modelo numérico utilizado==

===3.1. Fundamentos===

A versão do DEM aqui empregado consiste na representação de um meio contínuo utilizando um arranjo de massas concentradas em nós interligados por elementos uniaxiais que têm rigidez equivalente à do contínuo que se deseja representar. A estratégia de discretização usa um módulo cúbico básico, construído com 20 barras conectadas a 9 nós, 8 externos e um interno central.

Na [[#fig0010|figura 2]]  se apresenta a topologia do arranjo adotado. Cada nó tem 3 graus de liberdade, sendo os deslocamentos nas direções das coordenadas. No caso de um material elástico isotrópico, a área transversal ''A''<sub>''i''</sub>  dos elementos longitudinais na representação discreta é dada pela seguinte expressão:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>A_i=\phi L^2</math>
|}
| style="text-align: right;" | ( 4)
|}

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr2.jpg|center|375px|Discretização do MED: a) célula cúbica básica, b) corpo prismático.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Discretização do MED: a) célula cúbica básica, b) corpo prismático.

</span>
|}

onde ''L''  é o comprimento dos lados do módulo cúbico básico. A área dos elementos diagonais é definida por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>A_d=\frac{2}{\sqrt{3}}\delta \phi L^2</math>
|}
| style="text-align: right;" | ( 5)
|}

Os comprimentos dos elementos longitudinais ''L''  e diagonais ''L''<sub>''d''</sub>  estão relacionados pela equação:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>L=\left(\frac{2}{\sqrt{3}}\right)L_d</math>
|}
| style="text-align: right;" | ( 6)
|}

Finalmente, para materiais isotrópicos os coeficientes ''ϕ''  e ''δ''  resultam:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{(9+8\delta )}{(18+24\delta )}</math>
|}
| style="text-align: right;" | ( 7)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\delta =9\nu (4-8\nu )</math>
|}
| style="text-align: right;" | ( 8)
|}

onde ''ν''  é o coeficiente de Poisson. Existe uma equivalência completa entre o modelo discreto e um contínuo isotrópico para ''ν''  = 0.25. Para ''ν''  ≠ 0.25 pequenas diferenças aparecem nos termos de corte. Mas é oportuno observar que não existe um sólido localmente isotrópico, como tão pouco existe o denominado meio contínuo. A isotropia nos sólidos é uma propriedade que resulta da distribuição aleatória da orientação dos elementos constituintes. Equações de movimento podem ser determinadas a partir de condições de equilíbrio de todas as forças que atuam sobre as massas nodais, resultando num sistema da forma (9):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>M_{\ddot{x}}+C_{\dot{x}}+F\left(t\right)-P(t)=0</math>
|}
| style="text-align: right;" | ( 9)
|}

No qual <math display="inline">\dot{x}</math>  e <math display="inline">\ddot{x}</math>  representam os vetores de velocidade e aceleração, respectivamente, ''M''  e ''C''  as matrizes de massa e amortecimento, e os vetores ''F'' (''t'' ) e ''P'' (''t'' ) as cargas nodais internas e externas, respectivamente. No modelo aqui utilizado, se considera a hipótese simplificativa de considerar o amortecimento ''C''  proporcional à massa ''M'' :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>C=MD_f</math>
|}
| style="text-align: right;" | ( 10)
|}

sendo ''D''<sub>''f''</sub>  uma constante vinculada ao coeficiente de amortecimento crítico, ''ξ''<sub>''n''</sub> , através da equação (11):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>D_f=\xi _n2\pi f_n</math>
|}
| style="text-align: right;" | ( 11)
|}

onde ''f''<sub>''n''</sub>  representa a frequência natural de vibração do modo ''n'' , expresso em [Hz]. O modo ''n''  é, em geral, o modo fundamental de vibração da estrutura.      

Quando as matrizes ''M''  e ''C''  são diagonais o sistema de equações (9) não é acoplado e pode ser integrado no domínio do tempo utilizando um método explícito de integração (método das diferenças finitas centrais). Se as coordenadas nodais são atualizadas a cada passo de integração, grandes deslocamentos podem ser considerados sem dificuldade. Em materiais elásticos lineares, a estabilidade do método de integração exige que o intervalo de tempo ''Δt''  não exceda o valor dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Delta t\leq \frac{0.6L}{C_p}</math>
|}
| style="text-align: right;" | ( 12)
|}

onde ''C''<sub>''p''</sub>  é a velocidade de propagação de ondas longitudinais (ondas P) em barras com deslocamento axial:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>C_p=\sqrt{E/\rho }</math>
|}
| style="text-align: right;" | ( 13)
|}

Na expressão (13), ''E''  representa o módulo de Young e ''ρ''  a densidade do material. No caso de materiais não lineares, a estabilidade do processo de integração pode ser verificada controlando a energia total no processo. A convergência de soluções utilizando o DEM em elasticidade linear e em problemas de instabilidade elástica foi verificada, entre outros, por Dalguer et al.  [[#bib0290|[19]]] . Mais detalhes sobre a formulação do DEM podem ser encontrados em Kosteski et al. [[#bib0295|[20]]] .      

===3.2. Modelo constitutivo===

Riera e Rocha [[#bib0300|[21]]]  e Rocha et al. [[#bib0305|[22]]]  adaptaram a lei constitutiva para materiais quase frágeis proposta por Hilleborg [[#bib0310|[23]]] , estendendo esta versão do DEM ao estudo de problemas de fratura em materiais quase frágeis por meio da relação constitutiva bilinear (ECR) apresentada na [[#fig0015|figura 3]] , que permite a consideração dos efeitos irreversíveis de nucleação e propagação de fissuras.

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr3.jpg|center|284px|Lei constitutiva bilinear adotada no modelo utilizado.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Lei constitutiva bilinear adotada no modelo utilizado.

</span>
|}

A área embaixo da curva força versus deformação (área do triângulo OAB na [[#fig0015|figura 3]] ) representa a densidade de energia necessária para fraturar a área de influência do elemento. Assim, para um determinado ponto P sob a curva, a área do triângulo OPC representa a densidade de energia elástica reversível armazenada no elemento, e a área do triângulo OAP define a densidade de energia dissipada pelo dano. Quando submetido à compressão se admite que o comportamento do material é linear elástico, ou seja, a falha por compressão é induzida por tração indireta.      

Os parâmetros do modelo constitutivo ilustrado na [[#fig0015|figura 3]]  são:
* Força, ''F'' : é a força agindo sobre o elemento, função da deformação longitudinal ''ɛ'' .                          
* Rigidez do elemento, ''EA''<sub>''i''</sub> : dependendo se o elemento é normal ou diagonal diferente tipo de rigidez é considerada (''EA''<sub>''l''</sub>  ou ''EA''<sub>''d''</sub> ), as áreas ''A''<sub>''l''</sub>  e ''A''<sub>''d''</sub>  são definidas nas expressões (4) e (5).                          
* Comprimento do módulo cúbico, célula básica do modelo (''L'' ).                          
* Tenacidade ''G''<sub>''f''</sub> : a energia dissipada na fratura por unidade de área.                          
* Área equivalente de fratura das barras <math display="inline">A_i^f</math> , onde o subíndice ''i''  será ''l''  para barras normais e ''d''  para diagonais: este parâmetro força a condição de que a energia dissipada na fratura do contínuo é equivalente à dissipada na fratura do modelo discreto. A energia dissipada pela fratura de uma porção de material de tamanho ''L''  × ''L''  × ''L''  devido à fratura paralela a uma de suas faces é:                          

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Gamma =G_f\Lambda =G_fL^2</math>
|}
| style="text-align: right;" | ( 14)
|}

A energia dissipada pela fratura num módulo típico do modelo discreto leva em conta a contribuição de 5 barras longitudinais (4 coincidentes com as arestas do módulo e uma interna) e 4 elementos diagonais, ver [[#fig0010|figura 2]] (a).      

Desta forma é possível escrever:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Gamma _{DEM}=G_f\left[4\cdot \left(0.25\right)\cdot c_A+c_A+4c_A\left(\frac{2}{\sqrt{3}}\delta \right)\right]L^2</math>
|}
| style="text-align: right;" | ( 15)
|}

O primeiro termo entre parênteses da expressão (15) leva em conta a contribuição dos elementos longitudinais localizados nas bordas do módulo cúbico. O segundo termo representa a barra longitudinal interna e o terceiro termo representa a contribuição dos 4 elementos diagonais. Neste terceiro termo, o fator <math display="inline">\left[\left(2/\sqrt{3}\right)\delta \right]</math>  corresponde à relação entre as rigidezes das barras diagonais e longitudinais, é possível chegar a esta relação fazendo o quociente entre as expressões (5) e (4), já vistas. O coeficiente ''c''<sub>''A''</sub>  permite estabelecer a equivalência entre ''Γ''  dado pela expressão (14) e ''Γ''<sub>''DEM''</sub>  pela expressão (15), considerando que ''ν''  = 0.25 e utilizando a expressão (8) para definir ''δ''  = 1.125, é possível simplificar a equação (15) como segue:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>G_fL^2=G_f(\frac{22}{3}c_A)L^2</math>
|}
| style="text-align: right;" | ( 16)
|}

E determinar que ''c''<sub>''A''</sub>  = 3/22. Finalmente, as áreas equivalentes de fratura para as barras diagonais e longitudinais ficam definidas como segue:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}
A_l^f=c_AL^2=\frac{3}{22}L^2\\
A_d^f=c_A\sqrt{3}L_d^2=\frac{4}{22}L^2
\end{array}</math>
|}
| style="text-align: right;" | ( 17)
|}

Este resultado está fundamentado em assumir certa relação entre os modos de fratura I , II  e III , definida quando consideramos o fator <math display="inline">\left[\left(2/\sqrt{3}\right)\delta \right]</math>  no terceiro termo da expressão (15). Outras relações entre os modos de fratura poderiam ser introduzidas no modelo com extrema facilidade levando em conta esta característica particular da microestrutura de cada material. Mais detalhes acerca desta propriedade do material se pode obter em Ayatollahi et al. [[#bib0315|[24]]] . Se introduz a aleatoriedade no modelo de 2 formas:
* Considerando a tenacidade do material ''G''<sub>''f''</sub>  como um campo aleatório com uma distribuição de probabilidades de Weibull tipo III  com média, desvio e comprimento de correlação determinados como parâmetros de entrada.                          
* Introduzindo um campo de perturbações na malha dos elementos, governada por uma distribuição de probabilidades normal, o que permite capturar melhor o comportamento de corpos de prova de materiais quase frágeis quando a ruptura é governada pelas tensões de compressão. Sobre este particular ver em Iturrioz et al. [[#bib0320|[25]]]  e Riera et al. [[#bib0325|[26]]] .                  

A versão do modelo de barras utilizado no presente trabalho foi proposta por Riera [[#bib0330|[27]]]  para determinar a resposta dinâmica de placas e cascas sob carga de impacto. O DEM tem sido usado com sucesso para resolver problemas de dinâmica estrutural, tais como cascas submetidas a carregamento impulsivo (Riera e Iturrioz [[#bib0335|[28]]] , [[#bib0340|[29]]] ), a geração e posterior propagação de um sismo (Dalguer et al. [[#bib0345|[30]]] , [[#bib0350|[31]]] ), o estudo do efeito de escala em concreto (Rios e Riera [[#bib0355|[32]]] ) e em buchas de rocha (Miguel et al. [[#bib0360|[33]]] ; Iturrioz et al. [[#bib0365|[34]]] ), no cálculo dos parâmetros de fratura em problemas estáticos e dinâmicos (Kosteski [[#bib0370|[35]]] ), no estudo da resistência dos materiais frágeis sob altas taxas de deformação (Riera et al. [[#bib0375|[36]]] ), e na simulação de ensaios de emissão acústica (Iturrioz et al. [[#bib0320|[25]]] ). Nas referências citadas é possível encontrar mais detalhes sobre a fundamentação teórica da formulação utilizada no DEM.      

==4. Simulação da propagação subcrítica de fissuras==

Para estudar a propagação subcrítica de fissuras em um corpo submetido a cargas oscilantes, foi desenvolvido um modelo utilizando o DEM, que consiste numa placa retangular com uma trinca reta que emana de uma de suas bordas laterais, com tensões aplicadas nas 2 faces em direção normal à fissura. Na [[#fig0020|figura 4]]  se apresenta a configuração geométrica do modelo e na [[#tbl0005|tabela 1]]  os parâmetros que definem o mesmo. Considerou‐se, para este estudo, estado plano de deformações, implementado no DEM utilizando somente um elemento na espessura e restringindo todos os deslocamentos na direção perpendicular ao plano xz, no qual a figura é exibida. Na [[#fig0020|figura 4]]  são indicadas as outras condições de contorno aplicadas na placa. A introdução de fissura inicial tem por objetivo localizar a propagação subcrítica do defeito na simulação realizada. O modelo de placa aqui utilizado foi submetido à tração, abrindo a trinca em modo I , com forças uniformemente distribuídas nas superfícies inferior e superior, de forma simétrica.

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr4.jpg|center|339px|Geometria do modelo estudado, onde b=0,75m (100 módulos) e h=0,3075m (41 ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Geometria do modelo estudado, onde ''b''  = 0,75 ''m''  (100 módulos) e ''h''  = 0,3075 ''m''  (41 módulos).                  

</span>
|}

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 1.

Parâmetros do modelo de elementos discretos

|-

! L (ver figura 2a)
| 0,0075m
|-

! Número de módulos na direção x
| 100
|-

! Número de módulos na direção y
| 1
|-

! Número de módulos na direção z
| 41
|-

! Coeficiente de Poisson [ν]
| 0,25
|-

! Massa específica [ρ]
| 2.400kg/m3
|-

! Módulo de Young [E]
| 35GPa
|-

! Deformação crítica, ɛp
| 2.18×10-4
|-

! Energia específica de fratura [Gf]
| 1.155N/m
|-

! Coeficiente de variação de Gf
| 5%
|}

Com os parâmetros considerados, a relação entre a deformação crítica e a deformação de ruptura será de ''ɛ''<sub>''r''</sub> /''ɛ''<sub>''p''</sub>  = 60.      

É importante salientar que, no presente trabalho, não se tem por objetivo caracterizar um material específico, mas sim explorar a capacidade da ferramenta em simular a propagação subcrítica de um defeito.

Nesta análise foi considerado, para definir a tenacidade do material, ''G''<sub>''f''</sub> , um campo aleatório caracterizado por um valor médio e coeficiente de variação especificado na [[#tbl0005|tabela 1]] . Se considera neste caso que o comprimento de correlação do campo é da ordem da discretização adotada, implementações para desvincular o nível de discretização do comprimento de correlação do campo aleatório estão explicadas em Puglia et al. [[#bib0380|[37]]] . Para todos os modelos estudados foi utilizado um amortecimento proporcional à massa, caracterizado por ''D''<sub>''f''</sub>  = 500. Um estudo mais aprofundado sobre a influência do amortecimento pode ser encontrado em Kosteski [[#bib0370|[35]]] .      

Para formar a macrofissura se debilitaram as barras correspondentes a um módulo do modelo, formas mais sofisticadas de criar trincas nos modelos do DEM podem ser vistas em Kosteski et al. [[#bib0385|[38]]] .      

Para aplicar o carregamento, foram selecionados os nós centrais de todas as células cúbicas localizadas nas faces superior e inferior da placa. Na [[#fig0025|figura 5]] , a forma de aplicação das forças é detalhada, assim como as restrições ao deslocamento aplicadas à superfície traseira (à direita) do modelo.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr5.jpg|center|373px|Esquema de aplicação de força e condições de contorno.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

Esquema de aplicação de força e condições de contorno.

</span>
|}

A variação temporal da carga aumenta gradualmente, primeiro atingindo uma magnitude fixa e depois oscilando em torno desse valor. Em cada nó carregado, o vetor de força atinge o valor médio de 165 N. Após entrada em regime, a amplitude do carregamento oscilatório assume valor constante de 46 N, o que gera força máxima de 188 N e mínima de 142 N por nó. O carregamento oscila conforme função senoidal. A tensão nominal agindo no modelo é calculada dividindo o somatório das forças aplicadas numa superfície da placa pela área correspondente, desconsiderando a trinca. Na Equação (18) é calculada a tensão nominal média:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\sigma _{m\acute{e}dia}=\frac{100n\acute{o}s\times 165N}{0,75\times 0,0075}\cong 2,93MPa</math>
|}
| style="text-align: right;" | ( 18)
|}

As restrições ao movimento aplicadas à superfície traseira da placa, indicadas na [[#fig0020|figura 4]] , correspondem à restrição na direção x aplicada aos nós centrais das células cúbicas constituintes daquela superfície.      

Na [[#fig0030|figura 6]]  se ilustra a função de força transferida a um elemento normal próximo a ponta do defeito inicial. Até a entrada em regime, em aproximadamente 0,7(t/tmáx), ocorrem 880 ciclos de carga com amplitude ascendente. Após entrada em regime, 1.020 ciclos de carga com amplitude constante são aplicados ao modelo até o tempo máximo simulado (t/tmáx) = 1.

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr6.jpg|center|375px|Curva força aplicada versus tempo no modelo implementado.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Curva força aplicada versus tempo no modelo implementado.

</span>
|}

===4.1. Resultados===

Na [[#fig0035|figura 7]]  são apresentadas 6 configurações obtidas durante todo o processo da propagação. Os elementos em cinza claro correspondem a barras sem danos, barras danificadas são indicadas em cinza escuro e barras rompidas não são apresentadas. O modelo constitutivo das barras não leva em conta plasticidade, ao serem descarregadas, as barras retornam ao seu comprimento original, de modo que não há deformação residual. Resultados obtidos com o DEM onde a lei constitutiva uniaxial é elastoplástica foram publicados em Kosteski et al. [[#bib0385|[38]]] . A região em cinza escuro precede o desenvolvimento da macrofissura, representando o «''cluster»''  de microfissuras característico na propagação de trincas em materiais quase frágeis.

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr7.jpg|center|357px|Configurações que ilustram o avanço subcrítico da fissura, onde t*= (t/tmáx) ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

Configurações que ilustram o avanço subcrítico da fissura, onde t<sup>*</sup>  =  (t/tmáx) × 100.                  

</span>
|}

Cabe salientar também que a direção de propagação das fissuras não é influenciada de forma significativa pela orientação da malha de elementos discretos. A este respeito ver Kosteski et al. [[#bib0385|[38]]] . Para poder perceber com clareza até onde acontece a propagação subcrítica da fissura, se apresenta na [[#fig0040|figura 8]]  o balanço energético durante o processo, onde variações abruptas nas formas de energia são associadas à instabilidade na propagação da fissura. Nesta mesma figura, se indicam os tempos nos quais acontecem as configurações apresentadas na [[#fig0035|figura 7]] . Da análise da [[#fig0040|figura 8]]  se observa que, até o momento em que a energia elástica é crescente, o processo de propagação da fissura pode ser considerado subcrítico (até aproximadamente t/tmáx = 0,93), após este tempo a dissipação da energia se acelera junto com a energia cinética. Isto indica que desde (t/tmáx = 0,70), quando se atinge a situação de regime, até (t/tmáx = 0,96), pode‐se considerar que a fissura propagou em forma subcrítica. Notar que, para facilitar a visualização, tanto a energia dissipada como elástica estão multiplicadas por fatores que reduzem seu valor em 50 e 25 vezes, respectivamente.

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr8.jpg|center|373px|Variação de energia cinética, elástica e de dano no tempo.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

Variação de energia cinética, elástica e de dano no tempo.

</span>
|}

No trecho situado entre (t/tmáx = 0,80) e (t/tmáx = 0,82), na [[#fig0040|figura 8]] , se observa que a energia cinética decai e aumenta a energia elástica (potencial), indicando acúmulo energético. Pode‐se perceber que a velocidade de crescimento da trinca varia pouco neste intervalo, esta última tendência pode ser observada também na [[#fig0035|figura 7]] , nos quadros 1 e 2.      

Entre (t/tmáx = 0,82) e (t/tmáx = 0,88), a taxa de crescimento da energia elástica diminui e aumenta a energia cinética e a dissipada, após (t/tmáx = 0,88) ocorre o contrário. Entre (t/tmáx = 0,92) e (t/tmáx = 0,94), novamente há um decréscimo de energia elástica acumulada acompanhado por uma aceleração do dano e aumento de energia cinética. Este comportamento «espasmódico» das energias também resulta em pequenas oscilações na forma da macrofissura como se evidencia nas configurações apresentadas na [[#fig0035|figura 7]] . Após (t/tmáx = 0,94), o crescimento subcrítico começa a se instabilizar. A energia elástica flutua e a energia cinética cresce abruptamente junto com a dissipada na fratura.      

Quando a energia dissipada deixa de crescer é sinal claro de que o corpo de prova rompeu. Entre o início do carregamento cíclico em regime, em (t/tmáx = 0,70), e a instabilização, em aproximadamente (t/tmáx = 0,96), ocorrem 880 ciclos de carga em regime.

Ao examinar a variação de energias no corpo, apresentada na [[#fig0040|figura 8]] , é possível concluir que a trinca propaga, de forma estável (subcrítica), desde (t/tmáx = 0,70), até à eminência da fratura, aos (t/tmáx = 0,96). Com o crescimento da fissura ocorrendo em forma subcrítica, busca‐se caracterizar a vida em fadiga do material segundo o modelo clássico de Paris [[#bib0205|[2]]] .      

===4.2. Caracterização da vida em fadiga===

Como citado anteriormente, o modelo clássico de Paris relaciona com uma função potencial a velocidade de propagação subcrítica de uma fissura versus o incremento do fator de intensidade de tensões ao qual o corpo fissurado está submetido. Na metodologia utilizada, o incremento do fator de intensidade de tensões, ''ΔK'' , é calculado através da equação extraída de Gdoutos  [[#bib0195|[39]]] , adequada à geometria do corpo modelado:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}
K_I=\sigma \sqrt{\pi a}\left[1.12-0.23\left(\frac{a}{b}\right)+10.55\left(\frac{a}{b}\right)^2\right.\\
\left.-21.72\left(\frac{a}{b}\right)^3+30.39\left(\frac{a}{b}\right)^4\right]
\end{array}</math>
|}
| style="text-align: right;" | ( 19)
|}

onde σ é a tensão que atua na direção e sentido de abrir a trinca, se a trinca não estivesse lá (equação 18), ''a''  é o comprimento da trinca, ''b''  é o comprimento do corpo e o termo entre colchetes equivale ao coeficiente de correção, ''β'' , apresentado na expressão (3). Conhecendo a mínima e máxima tensão aplicada, se obtêm ''ΔK''  a partir da Equação (19). Esta varia no tempo, aumentando na medida em que cresce a trinca.      

Obtendo a relação (''da'' /''dN'' ) versus (''ΔK'' ), a partir dos valores obtidos com o modelo de elementos discretos, e aplicando logaritmo em ambos os lados da expressão, se pode ajustar na região onde a curva bilogarítmica apresenta uma tendência linear o modelo de Paris, como se mostra a seguir:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>log\left(\frac{da}{dN}\right)=log\left(C\Delta K^m\right)=m\cdot log\left(\Delta K\right)+log\left(C\right)</math>
|}
| style="text-align: right;" | ( 20)
|}

====4.3. Determinação de ''da/dN''  utilizando a distribuição espacial dos elementos que atingem a sua deformação crítica ''ɛ''<sub>''p''</sub>====

Neste método de obtenção de ''da'' /''dN''  trabalha‐se com a hipótese de que a velocidade de propagação da fissura equivale à velocidade com que a área afetada pelo dano na ponta da fissura avança. Esta velocidade é medida computando o instante em que cada barra atinge a deformação crítica ''ɛ''<sub>''p''</sub> , versus as coordenadas do baricentro da barra respectiva. O mesmo procedimento foi aplicado utilizando a deformação de ruptura, ''ɛ''<sub>''r''</sub> .      

Na [[#fig0045|figura 9]] , o gráfico gerado, onde se visualiza a coordenada x do baricentro de cada barra no tempo em que se atinge ''ɛ''<sub>''p''</sub>  (pontos vermelhos) e ''ɛ''<sub>''r''</sub>  (pontos azuis), observa‐se que a porção inferior da fileira de pontos vermelhos representa as primeiras barras a atingirem ''ɛ''<sub>''p''</sub> . Este conjunto de pontos forma uma curva que é associada à velocidade de propagação do dano através da peça e, logo, à velocidade de propagação da trinca. Desse modo, 10 pontos são selecionados para gerar a curva ''a''  versus ''N'' , os pontos são indicados em verde na mesma  [[#fig0045|figura 9]] . O número de ciclos ''N''  é contado a partir do tempo normalizado de 70%, quando o carregamento entrou em regime, como já foi ilustrado na  [[#fig0030|figura 6]] .

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr9.jpg|center|301px|Obtenção de a versus N a partir da monitoração da deformação nas barras (número ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

Obtenção de ''a''  versus ''N''  a partir da monitoração da deformação nas barras (número de ciclos ''N''  contados a partir de quando o carregamento entra em regime (0,7[t/tmáx]), ver Figura 6).                  

</span>
|}

Se poderia argumentar que os instantes de ruptura, ''ɛ''<sub>''r''</sub> , representariam melhor a propagação da fissura, no entanto, é preciso considerar como ocorre o desenvolvimento do dano no modelo. O crescimento do defeito não se dá pela ruptura organizada de barras sucessivas, mas pelo avanço de um «''cluster»''  de barras afetadas pelo dano, dentro do qual elementos podem romper desordenadamente. Assim, justifica‐se a hipótese de que existe correlação entre a velocidade de propagação do início do dano, representada pela relação x versus tempo normalizado onde as barras atingem a deformação crítica ''ɛ''<sub>''p''</sub> , e a propagação da trinca. Conhecendo a velocidade de propagação do trinca, na [[#fig0050|figura 10]]  se apresenta a curva log ''(da/dN)''  versus log(''ΔK'' ) obtida. A linha de tendência modela o intervalo equivalente à zona II  de propagação, já apresentada na [[#fig0005|figura 1]] .

<span id='fig0050'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr10.jpg|center|262px|Curva log(da/dN) versus log(ΔK) utilizando os resultados simulados com DEM.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 10.

Curva log(''da'' /''dN'' ) versus log(''ΔK'' ) utilizando os resultados simulados com DEM.                  

</span>
|}

Na [[#fig0050|figura 10]] , a velocidade ''da'' /''dN''  está expressa em micrômetros por ciclo, o que translada a curva para cima do eixo x, salienta‐se, entretanto, que os seguintes cálculos são realizados com as unidades no SI (distância em metros e tensão em Pascal).      

A linha de tendência aplicada fornece a equação da reta:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>log\left(\frac{da}{dN}\right)\cong 0.22log\left(\Delta K\right)-4.37</math>
|}
| style="text-align: right;" | ( 21)
|}

Da onde se obtêm:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}
m=0.22\\
log\left(C\right)=-4.37\\
C=4.26\times 10^{-5}
\end{array}</math>
|}
| style="text-align: right;" | 
|}

O que permite caracterizar o material conforme a Lei de Paris:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{da}{dN}=4.26\times 10^{-5}\Delta K^{0.22}</math>
|}
| style="text-align: right;" | ( 22)
|}

Cabe salientar aqui que não se teve neste estudo interesse em modelar um material em particular e sim realizar um estudo preliminar sobre as possibilidades de o método escolhido simular o crescimento subcrítico de uma fissura devido à ação de cargas oscilantes.

==5. Estudo do efeito de escala==

Apresentam‐se seguidamente os resultados obtidos aplicando o modelo de elementos discretos apresentado ao estudo do efeito de escala na propagação subcrítica de macrofissuras num material quase frágil sujeito à fadiga. Siegmund [[#bib0250|[11]]] , ao estudar o mesmo efeito, aplicou o método de zona coesiva (MZC) e verificou diferenças na forma de propagação da fissura e na distribuição de dano para modelos de tamanhos diferentes. A geometria do modelo utilizado pelo autor é semelhante ao modelo utilizado neste trabalho, o que permite realizar comparação entre os resultados obtidos através do DEM e do MZC.      

Com o objetivo de estudar o efeito da escala na forma como um corpo se danifica e rompe, foram gerados modelos com 3 tamanhos diferentes, mantendo constantes a tensão e a intensidade de tensões solicitantes.

Para garantir a propagação da trinca em modo I , elementos contidos numa faixa central longitudinal dos modelos tiveram sua deformação crítica, ''ɛ''<sub>''p''</sub> , reduzida para 80% do valor base para evitar danos fora da área de interesse induzida pelas condições de contorno aplicadas. Na [[#fig0055|figura 11]] , verifica‐se um dos modelos gerados, onde a faixa de elementos enfraquecidos é indicada com um marco em azul.

<span id='fig0055'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr11.jpg|center|315px|Região de elementos com ɛp reduzido a 80% do valor base.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 11.

Região de elementos com ''ɛ''<sub>''p''</sub>  reduzido a 80% do valor base.                  

</span>
|}

Para o estudo do efeito de escala não se aplicou restrição ao deslocamento na direção x. Fora isso, as condições de contorno aplicadas aos modelos utilizados neste estudo são idênticas ao descrito na [[#fig0020|figura 4]] .      

Como no modelo apresentado anteriormente, a tenacidade do material é definida por um campo aleatório caracterizado por um valor médio e coeficiente de variação de 5%. Trabalha‐se com fator de falha ''R''<sub>''f''</sub>  = 1.2, utilizando nas barras a lei constitutiva bilinear descrita anteriormente, com valores de ''ɛ''<sub>''p''</sub>  = 8.28 × 10<sup>−04</sup>  e ''ɛ''<sub>''r''</sub>  = 5.0 × 10<sup>−02</sup> , o que implica num valor de ''k''<sub>''r''</sub>  = 60. A força média aplicada a cada nó carregado, conforme [[#fig0025|figura 5]] , é de 515 N e a amplitude de oscilação de 48 N. O efeito de escala é verificado através de 3 modelos com tamanhos diferentes, cujas dimensões são descritas na [[#tbl0010|tabela 2]] .

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabela 2.

Dimensões dos modelos estudados

|-

| 
| Maior
| Interm.
| Menor
|-

! Módulos na direção x
| 100 (0,75m)
| 75 (0,5625m)
| 50 (0,375m)
|-

! Módulos na direção y
| 1 (0,0075m)
| 1 (0,0075m)
| 1 (0,0075m)
|-

! Módulos na direção z
| 41 (0,3075m)
| 31 (0,2325m)
| 21 (0,1575m)
|}

Massa específica, coeficiente de Poisson, módulo de Young e o coeficiente de variação da tenacidade permanecem como descrito na [[#tbl0005|tabela 1]] . Note que, como se trabalha em estado plano de deformações, só se considera um módulo na direção da espessura e se restringem os deslocamentos dos nós na direção ''y'' .      

Siegmund [[#bib0250|[11]]]  avalia o efeito da escala ao comparar, para modelos de diferentes tamanhos, a velocidade de propagação das fissuras e o nível de dissipação energética devida ao dano ao longo do comprimento dos modelos. Neste trabalho, os mesmos aspectos são avaliados através das metodologias apresentadas, fundamentadas no DEM, e os resultados obtidos são comparados aos resultados do autor citado.      

Seguindo o procedimento descrito anteriormente, em cada um dos 3 modelos gerados a trinca é verificada visualmente e a seguir sua criticidade é avaliada através do balanço energético durante todo o processo. Os 3 modelos são submetidos à mesma tensão e intensidade de tensões. Trabalha‐se novamente com o tempo em forma normalizada.

Na [[#fig0060|figura 12]]  os 3 modelos gerados, em instantes diferentes da propagação. Para os 3 modelos, o número de ciclos é contado a partir do instante em que a carga oscilante entra em regime.

<span id='fig0060'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr12.jpg|center|320px|Diferentes configurações para os 3 modelos analizados.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 12.

Diferentes configurações para os 3 modelos analizados.

</span>
|}

Percebe‐se, na [[#fig0060|figura 12]] , que os modelos de 3 tamanhos têm as trincas propagando em modo I  ao longo de toda sua extensão.      

Novamente, para verificar com maior clareza até onde a propagação da fissura ocorre de forma subcrítica, é avaliada a variação das formas de energia (cinética, elástica e dissipada) em cada caso. Na [[#fig0065|figura 13]]  as variações energéticas em cada um dos 3 modelos são apresentadas. Salienta‐se que, para facilitar a visualização no gráfico, tanto a energia dissipada como a elástica estão multiplicadas por fatores que reduzem seu valor em 50 e 25 vezes, respectivamente. As regiões indicadas por elipses correspondem aos intervalos de tempo onde se avalia o nível de dano, posteriormente.

<span id='fig0065'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr13.jpg|center|284px|Balanço de energia para os corpos de 3 tamanhos.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 13.

Balanço de energia para os corpos de 3 tamanhos.

</span>
|}

O critério de avaliação utilizado assume que a trinca passa a propagar de forma crítica a partir do momento em que se verifica queda da energia elástica acumulada no corpo. Dessa forma, observa‐se que ocorre propagação subcrítica até (t/tmáx ≈ 0,94) para o corpo de tamanho menor (t/tmáx ≈ 0,91) para o corpo de tamanho intermediário e (t/tmáx ≈ 0,94) para o corpo de maior tamanho.

Na [[#fig0070|figura 14]]  as curvas ''a''  versus ''N''  são apresentadas, com o comprimento da trinca, ''a'' , normalizado em relação ao comprimento do corpo, ''b'' . A inclinação da curva correspondente ao corpo menor é ilustrada pela linha tracejada verde e a do corpo intermediário pela vermelha. Nessa figura se observa que, para corpo de maior tamanho, maior número de ciclos ocorre durante a propagação do defeito, até à ruptura.

<span id='fig0070'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr14.jpg|center|323px|Curvas a versus N para os 3 modelos simulados.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 14.

Curvas ''a''  versus ''N''  para os 3 modelos simulados.                  

</span>
|}

Na [[#fig0075|figura 15]]  se apresentam os resultados obtidos por Siegmund [[#bib0250|[11]]] . O autor simula corpos de diversos tamanhos, como se verifica na quantidade de curvas geradas. O tamanho do corpo e a variação do comprimento da trinca, ''Δa'' , são normalizados em função de um comprimento característico, ''δ''<sub>0</sub> , que é proporcional à discretização dos modelos. Os números associados às curvas correspondem à razão ''h''<sub>''s''</sub> /''δ''<sub>0</sub>  (altura normalizada). Quanto maior o número associado a uma curva do gráfico, maior o corpo que esta curva representa.

<span id='fig0075'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr15.jpg|center|258px|Resultados de Siegmund em termos da velocidade de propagação da macrofissura, o ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 15.

Resultados de Siegmund em termos da velocidade de propagação da macrofissura, o número do lado da curva indica a dimensão característica do modelo medida em ''h''<sub>''s''</sub> /''δ''<sub>0</sub> , (Siegmund [[#bib0250|[11]]] ).                  

</span>
|}

Ao comparar os gráficos exibidos na [[#fig0070|figura 14]]  e na [[#fig0075|figura 15]]  se pode perceber a tendência comum de trincas em corpos maiores necessitarem maior quantidade de ciclos até à ruptura. Verifica‐se também que a propagação de trincas em corpos menores pode ter um período de nucleação mais longo do que em corpos maiores. Através de mecanismos diferentes, ambos os métodos (DEM e MZC) captam tendências semelhantes quanto ao efeito de escala na forma como se propaga a macrofissura. Para avaliar o efeito de escala no nível de dissipação energética, devida ao dano, ao longo do comprimento dos modelos (direção x), foram gerados gráficos onde se verifica, num determinado instante, a coordenada x dos pontos centrais de cada barra e a energia dissipada pelas barras que se encontram com seu baricentro nesta posição. Para cada modelo o estado do dano foi avaliado em um instante dentro dos intervalos indicados por elipses na [[#fig0065|figura 13]] , onde a propagação subcrítica da macrofissura acontece.      

Na [[#fig0080|figura 16]] a se apresenta o gráfico gerado para o modelo maior (L = 0,75m), onde estão indicadas as regiões equivalentes ao defeito inicial (1), ao macrodefeito (2), a região microfissurada à margem da trinca (3) e a região microfissurada que precede a trinca (4). A localização de cada região no modelo é verificada na [[#fig0080|figura 16]] b.

<span id='fig0080'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr16.jpg|center|320px|a) Nível de dissipação energética ao longo da direção x do modelo maior, b) ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 16.

a) Nível de dissipação energética ao longo da direção x do modelo maior, b) regiões indicadas do modelo.

</span>
|}

Dentro da região microfissurada (4), indicada na [[#fig0080|figura 16]] , a ruptura das barras é parcial e seu nível de dissipação energética é menor quanto mais distantes estiverem da macrotrinca. Os pontos dentro da região 4 descrevem uma curva que é capturada e as curvas descritas pela região 4 de cada modelo são sobrepostas e visualizadas na [[#fig0085|figura 17]] , onde as coordenadas x e a quantidade de energia dissipada em cada modelo são vistos de forma normalizada, para que se possa comparar os níveis de dissipação energética em coordenadas correspondentes dos 3 modelos. Na [[#fig0085|figura 17]]  o comprimento ''b''  de cada modelo é expresso em função do comprimento ''l''<sub>0</sub> , correspondente a 0,375 m

<span id='fig0085'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr17.jpg|center|360px|Comparação entre níveis de dissipação energética.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 17.

Comparação entre níveis de dissipação energética.

</span>
|}

Na [[#fig0085|figura 17]]  se verifica que, em uma coordenada x equivalente nos 3 modelos, a proporção de energia dissipada é maior para modelos menores. De forma coerente, também se verifica que, para corpos menores, o comprimento da região microfissurada na ponta da trinca representa uma extensão maior do corpo em questão.      

Na [[#fig0090|figura 18]]  está ilustrado o gráfico gerado por Siegmund [[#bib0250|[11]]]  na análise do efeito de escala sobre o nível de dissipação energética ao longo do comprimento dos modelos. O autor expressa o tamanho de seus modelos em função de uma altura ''H'' . As coordenadas x são normalizadas em função de um comprimento característico, ''δ''<sub>0</sub> , que é proporcional a discretização dos modelos.

<span id='fig0090'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_735392864-1-s2.0-S0213131515000280-gr18.jpg|center|240px|Resultados de Siegmund quanto à distribuição de dano (Siegmund [11]).]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 18.

Resultados de Siegmund quanto à distribuição de dano (Siegmund [[#bib0250|[11]]] ).                  

</span>
|}

Na [[#fig0090|figura 18]]  se verifica que para os modelos de dimensões ''12''  e ''37,5 H''  as curvas de dissipação energética ao longo de x têm formato semelhante ao verificado na  [[#fig0080|figura 16]] , sendo possível, em ambas as figuras, identificar a região de dissipação energética correspondente à região da trinca.      

As curvas correspondentes aos modelos de dimensões ''12''  e ''37,5 H'' , na  [[#fig0090|figura 18]] , têm concordância qualitativa com as curvas observadas na [[#fig0085|figura 17]] , onde se verifica que para coordenadas x correspondentes, corpos de menores dimensões tem nível de dissipação energética mais alto.      

Deve‐se observar que as diferenças entre as escalas dos modelos simulados por Siegmund [[#bib0250|[11]]]  são muito maiores que as diferenças apresentadas neste trabalho. O autor citado consegue verificar que, para modelos de dimensões reduzidas, a ruptura ocorre de forma abrupta, sem haver propagação subcrítica de um macro defeito. Dessa forma, o modelo de menores dimensões simulado por Siegmund [[#bib0250|[11]]]  rompe muito antes do que acontece nos modelos de dimensões ''12''  e ''37,5 H''  da  [[#fig0090|figura 18]] .      

A partir da comparação entre os gráficos ilustrados na [[#fig0070|figura 14]]  e na [[#fig0075|figura 15]]  e os ilustrados na [[#fig0085|figura 17]]  e na [[#fig0090|figura 18]] , se pode verificar concordância qualitativa entre os resultados obtidos através das metodologias apresentadas neste trabalho, fundamentadas no DEM, e os resultados obtidos por Siegmund [[#bib0250|[11]]] , que utiliza o MZC.      

==6. Conclusões==

No presente trabalho se aplicou o DEM para simular a propagação de uma fissura dentro de um corpo de geometria simples constituído de material quase frágil. No transcurso do trabalho foi possível obter as seguintes conclusões:
* Método numérico apresentado se mostrou uma ferramenta adequada para simular a propagação de fissuras em regime subcrítico.
* Os testes realizados mostraram que a propagação subcrítica segue o comportamento previsto pela Lei de Paris.
* Estudo preliminar de escala realizado mostrou que os resultados obtidos são coerentes com os apresentados por Siegmund [[#bib0250|[11]]] , que realizou um estudo similar utilizando elementos finitos junto ao método das interfaces coesivas para simular a propagação subcrítica de uma macrofissura.  

==Agradecimentos==

Se agradece à Fundação CAPES (www.capes.gov.br/ ) pelo auxílio financeiro recebido durante a realização desta pesquisa.

==Referências==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0200'></span>
[[#bib0200|[1]]] C. Moura Branco; Fadiga de Estruturas Soldadas, Fundação Calouste Gulbenkian; Lisboa, Portugal (1984)</li>
<li><span id='bib0205'></span>
[[#bib0205|[2]]] P.C. Paris, M.P. Gomez, W.E. Anderson; A rational analytic theory of fatigue; The Trend in Engineering, 13 (1) (1961), pp. 9–14</li>
<li><span id='bib0210'></span>
[[#bib0210|[3]]] S. Suresh; Fatigue of Materials; Cambridge University Press, UK (2000)</li>
<li><span id='bib0215'></span>
[[#bib0215|[4]]] N. Saintier, T. Palin-luc, J. Bénabes, F. Cocheteux; Non‐local energy based fatigue life calculation method under multiaxial variable amplitude loadings. Elsevier; International Journal of Fatigue, 54 (2013), pp. 68–83</li>
<li><span id='bib0220'></span>
[[#bib0220|[5]]] V. Infante, J.M. Silva; Case studies of computational simulations of fatigue crack propagation using finite elements analysis tools. Elsevier; Engineering Failure Analysis, 18 (2011), pp. 616–624</li>
<li><span id='bib0225'></span>
[[#bib0225|[6]]] S. Mohammadi; Extended Finite Element for Fracture Analysis of Structures; Blackwell Publishing, Iran (2008)</li>
<li><span id='bib0230'></span>
[[#bib0230|[7]]] J.F. Unger, S. Eckardt, C. Könke; Modelling of cohesive crack growth in concrete structures with the extended finite element method. Elsevier; Comput Methods Appl Mech Engrg., 196 (2007), pp. 4087–4100</li>
<li><span id='bib0235'></span>
[[#bib0235|[8]]] X.P. Xu, A. Needleman; Void nucleation by inclusions debonding in a crystal matrix; Model Simul Mater Sci Engng., 1 (1993), pp. 111–132</li>
<li><span id='bib0240'></span>
[[#bib0240|[9]]] Z.J. Yang, J.F. Chen; Finite element modelling of multiple cohesive discrete crack propagation in reinforced concrete beams. Elsevier; Engineering Fracture Mechanics, 72 (2005), pp. 2280–2297</li>
<li><span id='bib0245'></span>
[[#bib0245|[10]]] T. Yamaguchi, T. Okabe, T. Kosaka; Fatigue simulation for Ti/GFRP laminates using cohesive elements. VSP; Advanced Composite Materials, 19 (2010), pp. 107–122</li>
<li><span id='bib0250'></span>
[[#bib0250|[11]]] T. Siegmund; Cyclic Crack Growth and Length Scales. School of Mechanical Engineering; Purdue University, Indiana, USA (2007)</li>
<li><span id='bib0255'></span>
[[#bib0255|[12]]] Y. Xu, H. Yuan; Computational analysis of mixed‐mode fatigue crack growth in quasi‐brittle materials using extended finite element methods. Elsevier; Engineering Fracture Mechanics, 76 (2009), pp. 165–181</li>
<li><span id='bib0260'></span>
[[#bib0260|[13]]] Munjiza, A., (2009) Guest editorial From: Engineering Computations. International Journal for Computer‐Aided Engineering and Software, 26/6.</li>
<li><span id='bib0265'></span>
[[#bib0265|[14]]] E. Schlangen, J.G.M. van Mier; Crack propagation in sandstone: Combined experimental and numerical approach; Rock Mech and Rock Eng, 28 (2) (1995), pp. 93–110</li>
<li><span id='bib0270'></span>
[[#bib0270|[15]]] A. Rinaldi; Advances In Statistical Damage Mechanics: New Modelling Strategies; G. Voyiadjis (Ed.), Damage Mechanics and Micromechanics of Localized Fracture Phenomena in Inelastic Solids, Springer (2011) CISM Course Series</li>
<li><span id='bib0275'></span>
[[#bib0275|[16]]] A. Hansen; Failure processes in elastic fiber bundles; Reviews of Modern Physics (2010) [http://dx.doi.org/10.1103/RevModPhys.82.499 http://dx.doi.org/10.1103/RevModPhys.82.499]</li>
<li><span id='bib0280'></span>
[[#bib0280|[17]]] A. Rinaldi, P. Peralta, D. Krajcinovic, Y.-C. Lai; Prediction of scatter in fatigue properties using discrete damage mechanics; Elsevier, International Journal of Fatigue, 28 (2006), pp. 1069–1080</li>
<li><span id='bib0285'></span>
[[#bib0285|[18]]] T.L. Anderson; Fracture Mechanics–Fundamentals and Applications, Department of Mechanical Engineering; Texas A&M University, USA (1994)</li>
<li><span id='bib0290'></span>
[[#bib0290|[19]]] L.A. Dalguer, K. Irikura, J.D. Riera, H.C. Chiu; The importance of the dynamic source effects on strong ground motion during the 1999 Chi‐Chi Taiwan, earthquake: Brief interpretation of the damage distribution on buildings; Bull Seismol Soc Am., 91 (2001), pp. 1112–1127</li>
<li><span id='bib0295'></span>
[[#bib0295|[20]]] L.E. Kosteski, I. Iturrioz, R.G. Batista, A.P. Cisilino; The truss‐like discrete element method in fracture and damage mechanics; Engineering Computations, 6 (2011), pp. 765–787</li>
<li><span id='bib0300'></span>
[[#bib0300|[21]]] J.D. Riera, M.M. Rocha; A note on velocity of crack propagation in tensile fracture; Revista Brasileira de Ciências Mecânicas, XII/3 (1991), pp. 217–240</li>
<li><span id='bib0305'></span>
[[#bib0305|[22]]] M.M. Rocha, J.D. Riera, N.J. Krutzik; Extension of a model that aptly describes fracture of plain concrete to the impact analysis of reinforced concrete, Int.; Conf. and Structural Mechanics in Reactor Technology, SMiRT 11, Trans. Vol. J., Tokyo, Japan (1991)</li>
<li><span id='bib0310'></span>
[[#bib0310|[23]]] Hillerborg, A., A Model for Fracture Analysis. Cod. LUTVDG/TVBM 300‐51‐8, (1971).</li>
<li><span id='bib0315'></span>
[[#bib0315|[24]]] M.R. Ayatollahi, M.R.M. Aliha, M.M. Hassani; Mixed mode brittle fracture in PMMA ‐ An experimental study using SCB specimens; Materials Science and Engineering, A‐Structural Materials Properties Microstructure and Processing, 417 (1‐2) (2006), pp. 348–356</li>
<li><span id='bib0320'></span>
[[#bib0320|[25]]] I. Iturrioz, J.D. Riera, L.F.F. Miguel; Introduction of imperfections in the cubic mesh of the Discrete Element Method, 22th International Conference on Structural Mechanics in Reactor Technology (SMiRT 22); San Francisco, CA USA (2013)</li>
<li><span id='bib0325'></span>
[[#bib0325|[26]]] Riera, J. D., Miguel, L. F. F., Iturrioz, I, (2013), «Study of Imperfections in the Cubic Mesh of the Truss‐Like Discrete Element Method (DEM)», International Journal of Damage Mechanics, submitted to publication.</li>
<li><span id='bib0330'></span>
[[#bib0330|[27]]] J.D. Riera; Local effects in impact problems on concrete structures, Proceedings Conference on Structural Analysis and Design of Nuclear Power Plants. Oct.; Porto Alegre, RS, Brasil (1984) Vol. 3, CDU 264.04:621.311.2:621.039</li>
<li><span id='bib0335'></span>
[[#bib0335|[28]]] J.D. Riera, I. Iturrioz; Discrete element dynamic‐response of elastoplastic shells subjected to impulsive loading; Communications in Numerical Methods in Engineering, 11, Wiley & Sons, U. K (1995), pp. 417–426</li>
<li><span id='bib0340'></span>
[[#bib0340|[29]]] J.D. Riera, I. Iturrioz; Discrete element model for evaluating impact and impulsive response of reinforced concrete plates and shells subjected to impulsive loading; Nuclear Engineering and Design, 179 (1998), pp. 135–144</li>
<li><span id='bib0345'></span>
[[#bib0345|[30]]] L.A. Dalguer, K. Irikura, J.D. Riera, H.C. Chiu; The importance of the dynamic source effects on strong ground motion during the 1999 Chi‐Chi, Taiwan, earthquake: Brief interpretation of the damage distribution on buildings; Bull Seismol Soc Am., 91 (2001), pp. 1112–1127</li>
<li><span id='bib0350'></span>
[[#bib0350|[31]]] L.A. Dalguer, K. Irikura, J.D. Riera; Simulation of tensile crack generation by three‐dimensional dynamic shear rupture propagation during an earthquake; J Geophys Res., 108 (B3) (2003), p. 2144</li>
<li><span id='bib0355'></span>
[[#bib0355|[32]]] R.D. Rios, J.D. Riera; Size effects in the analysis of reinforced concrete structures, Engineering Structures; Elsevier, Engineering Structures, 26 (2004), pp. 1115–1125</li>
<li><span id='bib0360'></span>
[[#bib0360|[33]]] L.F.F. Miguel, J.D. Riera, I. Iturrioz; Influence of size on the constitutive equations of concrete or rock dowels; Int J Numer Anal Meth Geomech, 32 (15) (2008), pp. 1857–2188 [http://dx.doi.org/10.1002/nag.699 http://dx.doi.org/10.1002/nag.699]</li>
<li><span id='bib0365'></span>
[[#bib0365|[34]]] I. Iturrioz, L.F.F. Miguel, J.D. Riera; Dynamic fracture analysis of concrete or rock plates by means of the Discrete Element Method; LAJSS, 6 (2009), pp. 229–245</li>
<li><span id='bib0370'></span>
[[#bib0370|[35]]] L.E. Kosteski; Aplicação do método dos elementos discretos formado por barras no estudo do colapso de estruturas; Tese de Doutorado, Universidade Federal do Rio Grande do Sul (2012)</li>
<li><span id='bib0375'></span>
[[#bib0375|[36]]] J.D. Riera, L.F.F. Miguel, I. Iturrioz; Strength of brittle materials under high strain rates in DEM simulations; Computer Modeling in Engineering & Sciences, 82 (2011), pp. 113–136</li>
<li><span id='bib0380'></span>
[[#bib0380|[37]]] B.V. Puglia, I. Iturrioz, J.D. Riera, L. Kosteski; Random field generation of the material properties in the truss‐like discrete element method, Mec. Comp; Cilamce‐Mecom 2010, XXIX (2010), pp. 6793–6807</li>
<li><span id='bib0385'></span>
[[#bib0385|[38]]] L.E. Kosteski, R. Barrios, I. Iturrioz; Crack propagation in elastic solids using the truss‐like discrete element method; Int J Fract (2012) [http://dx.doi.org/10.1007/s10704%E2%80%90012%E2%80%909684%E2%80%904 http://dx.doi.org/10.1007/s10704‐012‐9684‐4]</li>
<li><span id='bib0195'></span>
[[#bib0195|[39]]] E.E. Gdoutos; Fracture Mechanics; an Introduction, Springer, Xanthi, Grécia (2005)</li>
</ol>