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 [1] ) e os que se baseiam na propagação subcrítica de um defeito já nucleado (metodologias fundamentadas na lei proposta por Paris [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 [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. [4] e Infante et al. [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 [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. [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 [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. [9] , Yamaguchi et al. [10] , Siegmund [11] e Xu e Yuan [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 [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 [14] e Rinaldi [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. [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. [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 [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:

${\frac {da}{dN}}=f\left(S_{m{\acute {a}}x}/S_{min},S_{m{\acute {a}}x}-S_{min},h\ldots \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 [2] é um caso particular de (1) que se expressa como segue:

${\frac {da}{dN}}=C\Delta K^{m}$

( 2)

onde C e m são constantes que dependem do material e da relação entre tensões S_{m áx} /S_min adotada. Graficando esta lei em escala bilogarítmica se pode representar por uma reta como se ilustra na região II da 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:

$\Delta K=\Delta \sigma {\sqrt {\pi \alpha }}\beta$

( 3)

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

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 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_min/S_max , 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 [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 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_i dos elementos longitudinais na representação discreta é dada pela seguinte expressão:

$A_{i}=\phi L^{2}$

( 4)

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

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

$A_{d}={\frac {2}{\sqrt {3}}}\delta \phi L^{2}$

( 5)

Os comprimentos dos elementos longitudinais L e diagonais L_d estão relacionados pela equação:

$L=\left({\frac {2}{\sqrt {3}}}\right)L_{d}$

( 6)

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

$\phi ={\frac {(9+8\delta )}{(18+24\delta )}}$

( 7)

$\delta =9\nu (4-8\nu )$

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

$M_{\ddot {x}}+C_{\dot {x}}+F\left(t\right)-P(t)=0$

( 9)

No qual ${\textstyle {\dot {x}}}$ e ${\textstyle {\ddot {x}}}$ 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 :

$C=MD_{f}$

( 10)

sendo D_f uma constante vinculada ao coeficiente de amortecimento crítico, ξ_n , através da equação (11):

$D_{f}=\xi _{n}2\pi f_{n}$

( 11)

onde f_n 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:

$\Delta t\leq {\frac {0.6L}{C_{p}}}$

( 12)

onde C_p é a velocidade de propagação de ondas longitudinais (ondas P) em barras com deslocamento axial:

$C_{p}={\sqrt {E/\rho }}$

( 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. [19] . Mais detalhes sobre a formulação do DEM podem ser encontrados em Kosteski et al. [20] .

3.2. Modelo constitutivo

Riera e Rocha [21] e Rocha et al. [22] adaptaram a lei constitutiva para materiais quase frágeis proposta por Hilleborg [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 figura 3 , que permite a consideração dos efeitos irreversíveis de nucleação e propagação de fissuras.

Figura 3. Lei constitutiva bilinear adotada no modelo utilizado.

A área embaixo da curva força versus deformação (área do triângulo OAB na 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 figura 3 são:

• Força, F : é a força agindo sobre o elemento, função da deformação longitudinal ɛ .
• Rigidez do elemento, EA_i : dependendo se o elemento é normal ou diagonal diferente tipo de rigidez é considerada (EA_l ou EA_d ), as áreas A_l e A_d são definidas nas expressões (4) e (5).
• Comprimento do módulo cúbico, célula básica do modelo (L ).
• Tenacidade G_f : a energia dissipada na fratura por unidade de área.
• Área equivalente de fratura das barras ${\textstyle A_{i}^{f}}$ , 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 é:

$\Gamma =G_{f}\Lambda =G_{f}L^{2}$

( 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 figura 2 (a).

Desta forma é possível escrever:

$\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}$

( 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 ${\textstyle \left[\left(2/{\sqrt {3}}\right)\delta \right]}$ 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_A permite estabelecer a equivalência entre Γ dado pela expressão (14) e Γ_DEM 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:

$G_{f}L^{2}=G_{f}({\frac {22}{3}}c_{A})L^{2}$

( 16)

E determinar que c_A  = 3/22. Finalmente, as áreas equivalentes de fratura para as barras diagonais e longitudinais ficam definidas como segue:

${\begin{array}{c}A_{l}^{f}=c_{A}L^{2}={\frac {3}{22}}L^{2}\\A_{d}^{f}=c_{A}{\sqrt {3}}L_{d}^{2}={\frac {4}{22}}L^{2}\end{array}}$

( 17)

Este resultado está fundamentado em assumir certa relação entre os modos de fratura I , II e III , definida quando consideramos o fator ${\textstyle \left[\left(2/{\sqrt {3}}\right)\delta \right]}$ 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. [24] . Se introduz a aleatoriedade no modelo de 2 formas:

• Considerando a tenacidade do material G_f 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. [25] e Riera et al. [26] .

A versão do modelo de barras utilizado no presente trabalho foi proposta por Riera [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 [28] , [29] ), a geração e posterior propagação de um sismo (Dalguer et al. [30] , [31] ), o estudo do efeito de escala em concreto (Rios e Riera [32] ) e em buchas de rocha (Miguel et al. [33] ; Iturrioz et al. [34] ), no cálculo dos parâmetros de fratura em problemas estáticos e dinâmicos (Kosteski [35] ), no estudo da resistência dos materiais frágeis sob altas taxas de deformação (Riera et al. [36] ), e na simulação de ensaios de emissão acústica (Iturrioz et al. [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 figura 4 se apresenta a configuração geométrica do modelo e na 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 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.

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

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 ɛ_r /ɛ_p  = 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_f , um campo aleatório caracterizado por um valor médio e coeficiente de variação especificado na 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. [37] . Para todos os modelos estudados foi utilizado um amortecimento proporcional à massa, caracterizado por D_f  = 500. Um estudo mais aprofundado sobre a influência do amortecimento pode ser encontrado em Kosteski [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. [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 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.

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

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:

$\sigma _{m{\acute {e}}dia}={\frac {100n{\acute {o}}s\times 165N}{0,75\times 0,0075}}\cong 2,93MPa$

( 18)

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

Na 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.

Figura 6. Curva força aplicada versus tempo no modelo implementado.

4.1. Resultados

Na 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. [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.

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

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. [38] . Para poder perceber com clareza até onde acontece a propagação subcrítica da fissura, se apresenta na 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 figura 7 . Da análise da 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.

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

No trecho situado entre (t/tmáx = 0,80) e (t/tmáx = 0,82), na 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 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 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 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 [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 [39] , adequada à geometria do corpo modelado:

${\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}}$

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

$log\left({\frac {da}{dN}}\right)=log\left(C\Delta K^{m}\right)=m\cdot log\left(\Delta K\right)+log\left(C\right)$

( 20)

4.3. Determinação de da/dN utilizando a distribuição espacial dos elementos que atingem a sua deformação crítica ɛ_p

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 ɛ_p , versus as coordenadas do baricentro da barra respectiva. O mesmo procedimento foi aplicado utilizando a deformação de ruptura, ɛ_r .

Na figura 9 , o gráfico gerado, onde se visualiza a coordenada x do baricentro de cada barra no tempo em que se atinge ɛ_p (pontos vermelhos) e ɛ_r (pontos azuis), observa‐se que a porção inferior da fileira de pontos vermelhos representa as primeiras barras a atingirem ɛ_p . 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 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 figura 6 .

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).

Se poderia argumentar que os instantes de ruptura, ɛ_r , 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 ɛ_p , e a propagação da trinca. Conhecendo a velocidade de propagação do trinca, na 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 figura 1 .

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

Na 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:

$log\left({\frac {da}{dN}}\right)\cong 0.22log\left(\Delta K\right)-4.37$

( 21)

Da onde se obtêm:

${\begin{array}{c}m=0.22\\log\left(C\right)=-4.37\\C=4.26\times 10^{-5}\end{array}}$

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

${\frac {da}{dN}}=4.26\times 10^{-5}\Delta K^{0.22}$

( 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 [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, ɛ_p , reduzida para 80% do valor base para evitar danos fora da área de interesse induzida pelas condições de contorno aplicadas. Na figura 11 , verifica‐se um dos modelos gerados, onde a faixa de elementos enfraquecidos é indicada com um marco em azul.

Figura 11. Região de elementos com ɛp reduzido a 80% do valor base.

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 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_f  = 1.2, utilizando nas barras a lei constitutiva bilinear descrita anteriormente, com valores de ɛ_p  = 8.28 × 10⁻⁰⁴ e ɛ_r  = 5.0 × 10⁻⁰² , o que implica num valor de k_r  = 60. A força média aplicada a cada nó carregado, conforme 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 tabela 2 .

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 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 [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 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.

Figura 12. Diferentes configurações para os 3 modelos analizados.

Percebe‐se, na 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 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.

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

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

Figura 14. Curvas a versus N para os 3 modelos simulados.

Na figura 15 se apresentam os resultados obtidos por Siegmund [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, δ₀ , que é proporcional à discretização dos modelos. Os números associados às curvas correspondem à razão h_s /δ₀ (altura normalizada). Quanto maior o número associado a uma curva do gráfico, maior o corpo que esta curva representa.

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 hs /δ0 , (Siegmund [11] ).

Ao comparar os gráficos exibidos na figura 14 e na 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 figura 13 , onde a propagação subcrítica da macrofissura acontece.

Na 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 figura 16 b.

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.

Dentro da região microfissurada (4), indicada na 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 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 figura 17 o comprimento b de cada modelo é expresso em função do comprimento l₀ , correspondente a 0,375 m

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

Na 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 figura 18 está ilustrado o gráfico gerado por Siegmund [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, δ₀ , que é proporcional a discretização dos modelos.

Figura 18. Resultados de Siegmund quanto à distribuição de dano (Siegmund [11] ).

Na 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 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 figura 18 , têm concordância qualitativa com as curvas observadas na 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 [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 [11] rompe muito antes do que acontece nos modelos de dimensões 12 e 37,5 H da figura 18 .

A partir da comparação entre os gráficos ilustrados na figura 14 e na figura 15 e os ilustrados na figura 17 e na 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 [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 [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

1. [1] C. Moura Branco; Fadiga de Estruturas Soldadas, Fundação Calouste Gulbenkian; Lisboa, Portugal (1984)
2. [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
3. [3] S. Suresh; Fatigue of Materials; Cambridge University Press, UK (2000)
4. [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
5. [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
6. [6] S. Mohammadi; Extended Finite Element for Fracture Analysis of Structures; Blackwell Publishing, Iran (2008)
7. [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
8. [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
9. [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
10. [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
11. [11] T. Siegmund; Cyclic Crack Growth and Length Scales. School of Mechanical Engineering; Purdue University, Indiana, USA (2007)
12. [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
13. [13] Munjiza, A., (2009) Guest editorial From: Engineering Computations. International Journal for Computer‐Aided Engineering and Software, 26/6.
14. [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
15. [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
16. [16] A. Hansen; Failure processes in elastic fiber bundles; Reviews of Modern Physics (2010) http://dx.doi.org/10.1103/RevModPhys.82.499
17. [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
18. [18] T.L. Anderson; Fracture Mechanics–Fundamentals and Applications, Department of Mechanical Engineering; Texas A&M University, USA (1994)
19. [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
20. [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
21. [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
22. [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)
23. [23] Hillerborg, A., A Model for Fracture Analysis. Cod. LUTVDG/TVBM 300‐51‐8, (1971).
24. [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
25. [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)
26. [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.
27. [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
28. [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
29. [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
30. [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
31. [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
32. [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
33. [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
34. [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
35. [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)
36. [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
37. [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
38. [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‐012‐9684‐4
39. [39] E.E. Gdoutos; Fracture Mechanics; an Introduction, Springer, Xanthi, Grécia (2005)