Simulation of resistance and deformability of rocks via discrete elements with cemented particles

Resumo

O Método dos Elementos Discretos (DEM) aparece como um dos métodos mais atraentes para entender problemas geotécnicos gerais e especificamente o comportamento de rochas devido às suas vantagens únicas em lidar com processos de fraturamento. Neste trabalho, foram feitas modelagens de ensaios de compressão simples e de tração indireta relativos a uma rocha magmática, conhecida como granito Lac du Bonnet, com base em uma abordagem discreta, utilizando o DEM e um modelo específico de contato e de cimentação entre grãos. Concluindo, foi realizada uma análise paramétrica da influência do módulo de Young em algumas propriedades macroscópicas do modelo numérico da rocha, visando gerar contribuições no entendimento do comportamento de resistência e de deformabilidade de rochas duras.

Palavras chaves: Método dos Elementos Discretos, bonded particle model, rocha, cimentação, resistência

Abstract

The Discrete Element Method (DEM) appears as one of the most attractive methods for understanding general geotechnical problems and specifically the behavior of rocks due to their unique advantages in dealing with fracturing processes. In this work, unconfined compression and indirect tensile tests were performed for a magmatic rock, known as Lac du Bonnet granite, based on a discrete approach using DEM and a specific contact and cementation model between the particles. Concluding, a parametric analysis of the influence of Young's modulus on some macroscopic properties of the rock’s numerical model was carried out, with the objective of generating contributions as to the understanding of the deformability and strength behavior of hard rocks.

Keywords: Discrete Element Methods, bonded particle model, rock, cementation, strength

1. Introdução

A grande dificuldade na modelagem de massas rochosas, seja qual for o método numérico, é que a rocha é um material geológico natural e por isso suas propriedades físicas e de engenharia não podem ser estabelecidas ou controladas através de um processo de artificial fabricação tal como para metais ou materiais plásticos. As massas rochosas são em grande parte descontínuas, não homogêneas, anisotrópicas e não elásticas.

Com o auxílio do Método dos Elementos Discretos (DEM), e através do modelo BPM (Bonded Particle Model), pode-se realizar a modelagem micro e macromecânica do comportamento tensão-deformação da rocha, apresentado uma abordagem na qual a rocha é tratada como um conjunto de partículas com resposta micro-macroscópicas sob carregamento.

Embora as rochas fraturadas não se pareçam com materiais granulares soltos, sua constituição inclui uma matriz cimentante envolvendo cristais, no caso de rochas magmáticas, ou partículas, como as rochas sedimentares detríticas. Portanto, os modelos de materiais granulares podem ser aplicados para investigar o comportamento mecânico quase microscópico de rochas, assumindo que na escala micro ou quase microscópica os materiais podem ser aproximados como conjuntos de partículas de tamanhos muito pequenos ligados entre si por diferentes modelos de contatos ou efeitos de cimentação. Assim, o comportamento mecânico global pode ser avaliado através das contribuições coletivas das partículas em processos de carga ou descarga que exibem movimento, deslocamento, destruição de ligações e re-ligação, deslizamento e rotação interparticular.

Neste trabalho, assume-se que a rocha se comporta como um material cimentado de grãos complexos, em que o cimento é deformável e pode se quebrar. Tal modelo conceitual pode, em princípio, explicar importantes aspectos do comportamento mecânico das rochas. Vários modelos numéricos têm sido propostos os quais imitam esse sistema com diferentes níveis de fidelidade. Para isto é utilizado o modelo de partículas cimentadas (BPM) o qual reproduz diretamente este sistema e assim, exibe um amplo conjunto de comportamentos emergentes que representam fielmente os da rocha real. O BPM, combinado com outros modelos de contato, fornece tanto uma ferramenta científica para investigar os micro-mecanismos, quanto para produzir complexos comportamentos macroscópicos.

2. Metodo dos Elementos Discretos

O DEM é uma técnica numérica que modela um sólido como uma coleção de partículas em forma de discos ou esferas, nas versões mais simples. O método foi desenvolvido originalmente por Cundall [1] para resolver problemas de mecânica das rochas. Cundall e Hart [2] mostraram que o DEM é melhor na modelagem de um material descontínuo do que outros instrumentos numéricos, como o Método dos Elementos Finitos (FEM). As partículas DEM podem se deslocar independente uma da outra e interatuam por meio de seus contatos. A finalidade é determinar o movimento de cada uma das partículas mediante a aplicação das leis do movimento de Newton. As forças que impulsionam o movimento são originadas a partir de forças externas, como a gravidade, e das interações que ocorrem entre pares de partículas em contato. Ao contrário dos métodos numéricos tradicionais baseados na teoria do meio contínuo, o DEM tem a vantagem de modelar o dano e a fratura frágil simultaneamente, o que acontece como uma consequência genérica da separação das partículas.

O DEM especifica propriedades das partículas e dos contatos na micro-escala (Figura 1), usando parâmetros como a rigidez nos contatos, a resistência dos contatos, o coeficiente de atrito das partículas, etc. No entanto, na macro-escala, o material é descrito através de constantes elásticas, resistência, etc. Embora não exista uma teoria completa para predizer as propriedades na macro-escala a partir da micro-escala e vice-versa, é estabelecida uma relação entre as propriedades micro e macro como pré-requisito em qualquer aplicação do DEM na engenharia.

Figura 1. Modelo mola-amortecedor do contato entre partículas (Figura adaptada de Ferreira [3])

O’Sullivan [4] afirma que a simulação progride como uma análise no tempo de um sistema de partículas, mesmo que o sistema de interesse responda de maneira quase estática. Um modelo DEM particulado consiste em um sistema ideal de partículas livres para se mover, conectadas por molas, amortecedores e/ou blocos deslizantes que simulam as interações nos contatos. Os contatos podem se formar entre os corpos e, à medida que o sistema se deforma, esses contatos podem quebrar e novos contatos podem se formar. Normalmente, uma pequena quantidade de sobreposição é permitida no contato entre os corpos, e esta sobreposição é análoga à deformação que ocorre nos contatos entre os corpos reais. Os "modelos constitutivos de contato simples", portanto, são usados para relacionar as forças de contato entre os corpos com a sobreposição de contato.

O ciclo de cálculo do DEM é um algoritmo de passo de tempo (time stepping) que consiste na aplicação repetida da lei de movimento em cada partícula, e da lei força- deslocamento em cada contato, além de uma constante atualização na posição de paredes e esferas (Figura 2). Contatos, os quais podem existir entre duas partículas ou entre uma partícula e uma parede, são formados e quebrados automaticamente durante a simulação. No início de cada passo de tempo, o conjunto de contatos é atualizado a partir das posições conhecidas de paredes e partículas. A lei de força-deslocamento é então aplicada em cada contato, de forma que sejam atualizadas as forças de contato baseadas no movimento relativo entre duas entidades e o modelo constitutivo de contato. Então, a lei de movimento é aplicada em cada partícula para atualizar a sua velocidade e posição, baseadas nas resultantes de força e momento que emergem das forças de contato e de qualquer força de corpo atuante na partícula. Também são atualizadas as posições das paredes, baseadas na velocidade prescrita nas mesmas. Como esse procedimento é feito para um número grande de partículas, o método é bastante dispendioso computacionalmente, limitando a extensão das simulações. Caso se queira trabalhar com número elevado de partículas, deve-se reduzir a duração virtual da simulação. Por outro lado, se a intenção é simular grandes períodos de tempo, deve-se reduzir o número de partículas. Apesar disso, o número de partículas em uma simulação ainda é bastante limitado.

Figura 2. Ciclo básico de cálculo do DEM (modificado de Ayquipa [5])

2.1. Formulação do Método dos Elementos Discretos

De acordo com Pinto [6], a formulação estabelecida para o desenvolvimento do modelo pelo DEM adota uma série de hipóteses que permitem simplificar o problema real, descartando os aspectos menos significativos e permitindo estabelecer um modelo físico e matemático do problema em estudo. Entre as hipóteses estabelecidas encontram-se as seguintes: as partículas são consideradas corpos rígidos; o contato acontece num ponto ou numa área muito pequena da interface entre cada partícula, no contato existe uma tolerância de penetração entre as partículas, limitada a uma fração do tamanho das partículas; a magnitude da penetração é relacionada com a força de contato através da Lei Força – Deslocamento; todas as partículas têm uma geometria conhecida (em 2D geralmente são utilizados discos e, em 3D, esferas); a geração de um meio empregando elementos discretos deve ser aleatória e os tamanhos dos mesmos devem ser tratados de forma similar (posição e diâmetro dos elementos distintos são aleatórios).

Na formulação podem-se incluir elementos rígidos ou paredes, os quais podem ter velocidade, deslocamentos ou forças prescritas. Estas paredes dão as condições de contorno às partículas, e podem compactar um arranjo de partículas.

O tipo de contorno mais utilizado é a parede rígida. Esse tipo de contorno é feito, em geral, de superfícies planas ou curvas, descritas analiticamente. Essas paredes não possuem inércia e as forças de contato determinadas na interação partícula-contorno são utilizadas somente para atualizar as coordenadas da partícula. Como as forças atuantes na parede rígida não influenciam seu movimento, seu deslocamento é feito por um controle explícito de suas velocidades translacional e rotacional. Quando as paredes são movimentadas, deformações e forças são impostas ao sistema de partículas. Em geral, forças de contato não são desenvolvidas entre paredes rígidas que se interceptam [4].

O movimento de translação e rotação de partículas rígidas (esféricas ou cilíndricas) é descrita por meio da Segunda Lei de Newton da dinâmica do corpo rígido. Para o i-ésimo elemento, tem-se:

${m}_{i}\ast {\ddot {u}}_{i}={F}_{i}$	(1)
${I}_{i}\ast {\overset {\cdot }{\omega }}_{i}={T}_{i}$	(2)

onde ${\textstyle {u}_{i}}$ é o deslocamento do centro da partícula em relação a um sistema de coordenadas (inercial), ${\textstyle {\omega }_{i}}$ é a velocidade angular, ${\textstyle {m}_{i}}$ a massa do elemento (partícula), ${\textstyle {I}_{i}}$ o momento de inércia, ${\textstyle {F}_{i}}$ a força resultante e ${\textstyle {T}_{i}}$ o torque resultante.

Os vetores força ${\textstyle {F}_{i}}$ e o torque ${\textstyle {T}_{i}}$ que atuam na partícula são obtidos através da soma de todas as interações da partícula com o meio.

${F}_{i}=\sum _{c=1}^{{N}_{c}}{F}_{i}^{c}+{F}_{i}^{ext}+{F}_{i}^{damp}$	(3)
${F}_{i}=\sum _{c=1}^{{N}_{c}}{F}_{i}^{c}+{F}_{i}^{ext}+{F}_{i}^{damp}$	(4)
${F}_{i}=\sum _{c=1}^{{N}_{c}}{F}_{i}^{c}+{F}_{i}^{ext}+{F}_{i}^{damp}$	(5)
${F}_{i}=\sum _{c=1}^{{N}_{c}}\left({{x}_{i}}^{\left(c\right)}\ast {F}_{i}^{c}+{q}_{i}^{c}\right)+$${T}_{i}^{ext}+{T}_{i}^{damp}$	(6)

onde ${\textstyle {F}_{i}^{ext}}$ e ${\textstyle {T}_{i}^{ext}}$ são forças externas, ${\textstyle {F}_{i}^{c}}$ é a força na partícula devido aos contatos pela interação com as demais partículas e os obstáculos, ${\textstyle {F}_{i}^{damp}}$ e ${\textstyle {T}_{i}^{damp}\,}$ são força e torque devido ao amortecimento do sistema, ${\textstyle {{x}_{i}}^{\left(c\right)}}$ é o vetor que liga o centro da partícula do i-ésimo elemento com o ponto do contato ${\textstyle c}$, ${\textstyle N_{c}}$ é o número de contatos e ${\textstyle {q}_{i}^{c}}$ é o torque devido ao rolamento ou rotação das partículas.

2.2 Discretização do meio

Uma das fases mais importantes em uma simulação com DEM é a discretização do meio, pois tem uma forte influência nos resultados obtidos nas simulações. Para modelar materiais granulares sem coesão, uma partícula de DEM representa um grão de material, permitindo o uso das propriedades mecânicas. No caso da simulação de materiais sólidos com partículas circulares (ou esféricas) e contatos cimentados, como o caso do BPM, a partícula DEM não tem um paralelo físico, e é apenas uma forma de discretização do meio. Neste último caso, a caracterização do conjunto de partículas está fortemente relacionada com os resultados obtidos e deve ser considerada no processo de calibração ou estimação dos parâmetros do modelo.

O empacotamento das partículas utilizado na simulação deve permitir a representação das propriedades mecânicas do material, como isotropia ou coeficiente de Poisson. Existem vários métodos para a geração de empacotamentos e condições iniciais, sendo eles divididos em três classes: dinâmicos, construtivos e algoritmos de rearranjo coletivo. Todos os métodos anteriormente mencionados geram pacotes aleatórios. Em Morales et al. [7], é utilizado o método de avanço frontal, para a simulação destas estruturas policristalinas, com este algoritmo o tamanho das partículas pode seguir qualquer distribuição. Este algoritmo é o mais eficiente dentre os atualmente existentes, porque permite obter empacotamentos com uma alta fração de volume, o que é muito importante no caso destas estruturas.

2.3 Relação micro-macro

A relação entre os parâmetros locais e as propriedades mecânicas pode ser definida usando quantidades diferentes que descrevem a configuração local do sistema granular. Para obter o comportamento macroscópico, é necessário o uso dessas grandezas que permitam a análise mecânica do sistema, considerando o comportamento no nível micromecânico e a caracterização do sistema das partículas, conforme ilustrado na Figura 3.

Figura 3. Esquema da relação micro-macro. Figura adaptada de Labra [8]

A finalidade destes tipos de trabalhos com DEM é estabelecer leis de escala que relacionem os parâmetros fenomenológicos de um material na macro-escala (módulo elástico, resistência à compressão) aos parâmetros da micro-escala (parâmetros que caracterizam o contato entre as partículas), e para estabelecer a dependência da resposta macro-escala (material atingindo à ruptura sob tração) em consequência dos fatores da micro-escala.

Autores, como Walton [9], Bathurst e Rothenburg [10], Chang e Misra [11,12], mostraram que a relação tensão-deformação pode ser derivada com base em uma abordagem microestrutural contínua para empacotamento aleatório. Chang e Misra [12] demonstraram que o tensor de módulos de elasticidade que relaciona a tensão- deformação pode assumir várias formas, dependendo da simetria das propriedades mecânicas do material granular, e que está intimamente relacionado com a estrutura de empacotamento. O estudo analisa a relação entre a simetria de propriedades mecânicas e a estrutura de empacotamento granular aleatório, usando uma abordagem baseada em micromecânica.

A maior parte das análises teóricas encontradas na literatura mostra que os parâmetros relacionados ao módulo de Young e o coeficiente de Poisson são a rigidez dos contatos (${\textstyle k_{n}}$, ${\textstyle k_{s}}$), a relação entre a rigidez normal e tangencial ${\textstyle \left({\frac {{k}_{n}}{{k}_{s}}}\right)}$, o raio da partícula ${\textstyle R}$, a porosidade (${\textstyle n}$), e a distribuição do tamanho das partículas.

Trabalhos mais recentes foram realizados no grupo de Geotecnia da Universidade de Brasília. Neves [13] estudou o comportamento de materiais granulares com o auxílio do software PFC2D, mostrando que na compressibilidade do meio, o formato dos grãos anguloso sobressai com relação aos outros fatores envolvidos, como tensão confinante e coeficiente de atrito entre partículas. O autor também conclui que a porosidade inicial das amostras é de grande importância para a determinação do comportamento do material granular quando cisalhado seja o meio formado por partículas arredondadas ou angulosas.

Albuquerque [14] demostrou como o ângulo de atrito das partículas influência nas respostas de resistência do meio macroscópico para um estudo em estado plano de deformações, observando que o atrito interparticular leva a um aumento do módulo de Young e do atrito macroscópicos e a uma diminuição do coeficiente de Poisson global.

Tedesco [15] considerou aspectos muito importantes na realização de ensaio de enrocamentos como a quebra dos grãos, módulo de Young e atrito das partículas. Os resultados mostraram que o principal parâmetro, que controla o comportamento dos ensaios de compressão e cisalhamento, é o módulo de Young da partícula. Tal fato deve-se à dependência entre o cálculo das forças nos contatos e esse parâmetro. Concluiu que para maiores valores de módulo de Young, maiores as tensões necessárias para deformar as amostras.

2.4 Formulação do modelo BPM (Bonded Particle Model)

O modelo BPM simula o comportamento mecânico de uma coleção de partículas rígidas esféricas ou de tamanhos não uniformes que podem ser unidas em seus pontos de contato e podem tanto transladar como rodar, movendo-se independentemente uma da outra. As partículas rígidas interagem apenas nos contatos suaves, que possuem rigidez finita normal e de cisalhamento. O comportamento mecânico deste sistema é descrito pelo elemento de cada partícula e pela força e momento atuando em cada contato. As leis do movimento de Newton fornecem a relação fundamental entre o movimento das partículas e as forças e os momentos resultantes desse movimento. As ligações suportam carga e podem quebrar.

A suposição da rigidez das partículas é razoável quando os movimentos ao longo das interfaces representam a maior parte da deformação em um material. A deformação de um conjunto de partículas empacotadas ou de um sistema granular, tal como areia, é bem descrita por esta suposição, porque a deformação resulta principalmente do deslizamento e rotação das partículas como corpos rígidos, e a abertura e intertravamento nas interfaces - não da deformação individual de partículas. A adição de ligações entre as partículas do sistema corresponde à adição de cimento real entre os grãos de uma rocha sedimentar, como no caso do arenito, ou de cimento natural entre os grãos de uma rocha cristalina, como o granito. A deformação de um conjunto de partículas ligadas, ou uma rocha, deve ser semelhante, e ambos os sistemas devem apresentar processos de formação de danos semelhantes sob carga crescente à medida que as ligações são quebradas progressivamente e ambos os sistemas evoluem gradualmente em direção a um estado granular. Se os grãos individuais ou outras características microestruturais são representados como aglomerados de partículas ligadas, então o esmagamento e a não homogeneidade do material em uma escala maior do que o tamanho do grão também podem ser abrangidos por este modelo [16,17,18].

2.4.1 Comportamento e parâmetros do grão-cimento

O modelo BPM imita o comportamento mecânico de uma coleção de grãos unidos por cimento. A discussão a seguir considera cada grão como uma partícula e cada entidade de cimento como uma ligação paralela (parallel bond). A força total e o momento que atuam em cada contato cimentado são compostos por uma força, ${\textstyle {F}_{i}}$, resultante de uma sobreposição de partícula-partícula, indicada na Figura 4 como o comportamento do grão, e uma força e um momento, ${\textstyle {\overline {F}}_{l}\,}$ e ${\textstyle {\overline {M}}_{l}}$, carregados pela ligação paralela e denotados como o comportamento do cimento. Estas quantidades contribuem para a força e momento resultante que atuam sobre as duas partículas em contato. A lei força-deslocamento para este sistema será agora descrita, primeiro para o comportamento do grão e depois para o comportamento do cimento. Observa- se que se uma ligação paralela não estiver presente em um contato, então somente a parte do comportamento força-deslocamento baseada no grão ocorre.

Figura 4. Comportamento força-deslocamento do sistema cimento-grão, modificado de Potyondy e Cundall [17]

A parte do comportamento força-deslocamento baseada no grão em cada contato é descrita pelos seguintes parâmetros (Figura 4): a rigidez normal e de cisalhamento, ${\textstyle {K}_{n}}$ e ${\textstyle {K}_{s}}$; o coeficiente de atrito, ${\textstyle \mu }$. As duas partículas em contato são assumidas como esferas no programa STAR CCM+, usado neste trabalho. Sempre que duas partículas se sobrepõem, é formado um contato no centro da região de sobreposição ao longo da linha que une os centros de partículas ${\textstyle {{x}_{i}}^{\left(c\right)}}$, na Figura 4, e duas molas lineares são inseridas (com rigidez derivadas das rigidezes das partículas, assumindo que as duas molas das partículas atuam em série), juntamente com um bloco deslizante na direção de cisalhamento.

O vetor de força de contato, ${\textstyle {F}_{i}}$ (que representa a ação da partícula A na partícula B), pode ser decomposto em componentes normais e de cisalhamento em relação ao plano de contato como:

${F}_{i}={F}^{n}\ast {n}_{i}+{F}^{s}\ast {t}_{i}$

(7)

onde ${\textstyle {F}^{n}}$ e ${\textstyle {F}^{s}}$ indicam as componentes de força normal e de cisalhamento, respectivamente, e ${\textstyle {n}_{i}}$ e ${\textstyle {t}_{i}}$ são os vetores unitários que definem o plano de contato. A força normal é calculada por:

${F}^{n}={K}^{n}\ast {U}^{n}$

(8)

onde ${\textstyle {U}^{n}}$ é a sobreposição e ${\textstyle {K}^{n}}$ sendo a rigidez normal do contato dada por:

${K}^{n}={\frac {{{k}_{n}}^{\left(A\right)}\ast {{k}_{n}}^{\left(B\right)}}{{{k}_{n}}^{\left(A\right)}+{{k}_{n}}^{\left(B\right)}}}$

(9)

com ${\textstyle {{k}_{n}}^{\left(A\right)}}$ e ${\textstyle {{k}_{n}}^{\left(B\right)}}$ sendo a rigidez normal de cada partícula.

A força de cisalhamento é calculada de forma incremental. Quando o contato é formado, ${\textstyle {F}^{s}}$ é zero. Cada aumento subsequente de deslocamento relativo de cisalhamento, Δ ${\textstyle {U}^{s}}$, produz um incremento da força de cisalhamento elástica, Δ ${\textstyle {F}^{s}}$, que é adicionado a ${\textstyle {F}^{s}}$ (depois que ${\textstyle {F}^{s}}$ for girada para considerar o movimento do plano de contato). O incremento da força elástica de cisalhamento é dado por:

${\textstyle \Delta {F}^{s}=-{k}^{s}\ast \Delta {U}^{s}}$

(10)

onde ${\textstyle {k}^{s}}$ é a rigidez de cisalhamento de contato dada por:

${k}^{s}={\frac {{{k}_{s}}^{\left(A\right)}\ast {{k}_{s}}^{\left(B\right)}}{{{k}_{s}}^{\left(A\right)}+{{k}_{s}}^{\left(B\right)}}}$

(11)

com ${\textstyle {{k}_{s}}^{\left(A\right)}}$ e ${\textstyle {{k}_{s}}^{\left(B\right)}}$ sendo a rigidez de cisalhamento das partículas. O incremento do deslocamento relativo durante o passo de tempo ${\textstyle \Delta t}$ é dado por ${\textstyle \Delta {U}_{i}}$ = ${\textstyle {V}_{i}}$Δt, onde ${\textstyle {V}_{i}}$ é a velocidade de contato

${V}_{i}={\left({{\overset {\cdot }{x}}_{i}}^{\left(c\right)}\right)}_{B}-{\left({{\overset {\cdot }{x}}_{i}}^{\left(c\right)}\right)}_{A}$

(12)

${V}_{i}=\left({{\overset {\cdot }{x}}_{i}}^{\left(B\right)}+{e}_{ijk}\ast {{w}_{j}}^{\left(B\right)}\ast \left({{x}_{k}}^{\left(C\right)}-\right.\right.$$\left.\left.{{x}_{k}}^{\left(B\right)}\right)\right)-\left({{\overset {\cdot }{x}}_{i}}^{\left(A\right)}+\right.$$\left.{e}_{ijk}\ast {{w}_{j}}^{\left(\left(A\right)\right)}\ast \left({{x}_{k}}^{\left(C\right)}-\right.\right.$$\left.\left.{{x}_{k}}^{\left(A\right)}\right)\right)$

(13)

${\textstyle {\overset {\cdot }{x}}_{l}}$ e ${\textstyle {w}_{j}}$ são as velocidades de translação e de rotação da partícula, respectivamente, ${\textstyle {e}_{ijk}}$ é o símbolo de permutação e (c) é o ponto de contato entre duas partículas. O vetor de incremento relativo de cisalhamento é:

$\Delta {U}_{i}^{s}={V}_{i}^{s}\ast \Delta t=\left({V}_{i}-{V}_{i}^{n}\right)\ast \Delta t$

(14)

Se ${\textstyle {U}^{n}\leq 0}$ (existe uma abertura), então as forças normais e de cisalhamento são ajustadas para zero; caso contrário, o deslizamento é considerado calculando o coeficiente de atrito do contato

$\mu =min\left({\mu }^{\left(A\right)},{\mu }^{\left(B\right)}\right)$

(15)

com ${\textstyle {\mu }^{\left(A\right)}}$ e ${\textstyle {\mu }^{\left(B\right)}}$ sendo os coeficientes de atrito de partículas. A força de cisalhamento é limitada a:

${F}^{s}=\mu \ast {F}^{n}$

(16)

A parte do comportamento força-deslocamento baseada no cimento em cada contato cimentado é descrita pelos seguintes parâmetros que definem uma ligação paralela (parallel bond), ilustrados na Figura 4: rigidez normal e de cisalhamento por unidade de área ${\textstyle {\overline {k}}^{n}}$ e ${\textstyle {\overline {k}}^{s}}$, resistência normal e cisalhante do material cimentante, ${\textstyle {\overline {\sigma }}_{c}}$ e ${\textstyle {\overline {\tau }}_{c}}$ ; e o multiplicador do raio da ligação ${\textstyle {\overline {\lambda }}}$ , de tal forma que o raio da ligação paralela é:

${\overline {R}}={\overline {\lambda }}min\left({R}^{\left(A\right)},{R}^{\left(B\right)}\right)$

(17)

com ${\textstyle {R}^{\left(A\right)}}$ e ${\textstyle {R}^{\left(B\right)}}$ sendo os raios das partículas. Uma ligação paralela aproxima o comportamento mecânico de um cimento elástico frágil que une as duas partículas ligadas. As ligações paralelas estabelecem uma interação elástica entre estas partículas, que atua em paralelo com a porção do comportamento força-deslocamento baseada no grão; assim, a existência de uma ligação paralela não impede o deslizamento. As ligações paralelas podem transmitir a força e o momento entre as partículas, enquanto os grãos só podem transmitir força. Uma ligação em paralelo pode ser vista como um conjunto de molas elásticas uniformemente distribuídas sobre uma secção transversal retangular em 2D e uma secção transversal circular em 3D, situada sobre o plano de contato e centrada no ponto de contato. Essas molas se comportam como uma viga cujo comprimento ${\textstyle {\overline {L}}}$, na Figura 4, se aproxima ao zero para ter um comportamento muito parecido com o das juntas.

A força total e os momentos transmitidos pela ligação paralela são denotados por ${\textstyle {\overline {F}}_{i}}$ e ${\textstyle {\overline {M}}_{i}}$, respectivamente, que representam a ação da ligação sobre a partícula B. Os vetores de força e momento podem ser decompostos em componentes normais e de cisalhamento em relação ao plano de contato da seguinte forma:

${\overline {F}}_{i}={\overline {F}}^{n}\ast {n}_{i}+{\overline {F}}^{s}\ast {t}_{i}$	(18)
${\overline {M}}_{i}={\overline {M}}^{n}\ast {n}_{i}+{\overline {M}}^{s}\ast {t}_{i}$	(19)

onde ${\textstyle {\overline {F}}^{n}}$, ${\textstyle {\overline {F}}^{s}}$ e ${\textstyle {\overline {M}}^{n}}$, ${\textstyle {\overline {M}}^{s}}$ denotam as forças e momentos axial e de cisalhamento, respectivamente, e ${\textstyle {n}_{i}}$ e ${\textstyle {t}_{i}}$ são os vetores unitários que definem o plano de contato (para o modelo em 2D, o momento de torção, ${\textstyle {\overline {M}}^{n}}$ = 0 e o momento de flexão ${\textstyle {\overline {M}}^{s}}$ atua na direção fora do plano). Quando a ligação paralela é formada, ${\textstyle {\overline {F}}_{i}}$ e ${\textstyle {\overline {M}}_{i}}$ têm valores iniciais igual a zero. Cada subsequente deslocamento relativo e incremento de rotação (Δ ${\textstyle {U}^{n}}$,Δ ${\textstyle {U}^{s}}$, Δ ${\textstyle {\theta }^{n}}$, Δ ${\textstyle {\theta }^{s}}$ com Δ ${\textstyle {\theta }_{i}}$=( ${\textstyle {{w}_{i}}^{\left(B\right)}{{w}_{i}}^{\left(A\right)}}$)Δt) produz um incremento da força elástica e momento que são adicionados aos valores atuais (depois que os vetores de componente de cisalhamento foram rodados para ter em conta o movimento do plano de contato). Os incrementos de força e momento elásticos são dados por:

$\Delta {\overline {F}}^{n}={\overline {k}}^{n}\ast A\ast \Delta {U}^{n}$	(20)
$\Delta {\overline {F}}^{s}=-{\overline {k}}^{s}\ast A\ast \Delta {U}^{ns}$	(21)
$\Delta {\overline {M}}^{n}=-{\overline {k}}^{s}\ast J\ast \Delta {\theta }^{n}$	(22)
$\Delta {\overline {M}}^{s}={-{\overline {k}}}^{n}\ast I\ast \Delta {\theta }^{s}$	(23)

Onde $A$, $I$ e $J$ são a área, momento de inércia e momento de inércia polar da secção transversal da ligação paralela, respectivamente. Estas quantidades são dadas por

${\textstyle A=\pi \ast {\overline {R}}^{2}}$ para 3D	(24)
${\textstyle I={\frac {1}{4}}\ast \pi \ast {\overline {R}}^{4}}$ para 3D	(25)
${\textstyle J={\frac {1}{2}}\ast \pi \ast {\overline {R}}^{4}}$ para 3D	(26)

A tensão normal e de cisalhamento máximas que atuam no contorno da ligação paralela são calculadas a partir da teoria da viga:

${{\overline {\sigma }}_{c}}^{max}={\frac {-{\overline {F}}^{n}}{A}}+{\frac {\left|{\overline {M}}^{s}\right|\ast {\overline {R}}}{I}}$	(27)
${{\overline {\tau }}_{c}}^{max}={\frac {-{\overline {F}}^{s}}{A}}+{\frac {\left|{\overline {M}}^{n}\right|\ast {\overline {R}}}{J}}$	(28)

Se a tensão de tração máxima excede a resistência à tração (${\textstyle {{\overline {\sigma }}_{c}}^{max}\geq {\overline {\sigma }}_{c}}$) ou a tensão de cisalhamento máxima excede a resistência ao cisalhamento (${\textstyle {{\overline {\tau }}_{c}}^{max}\geq {\overline {\tau }}_{c}}$), então a ligação paralela rompe e é removida do modelo juntamente com a sua força, momento e rigidez.

A forma simplificada do BPM representa o cimento, utilizando ligações de contato em vez de ligações paralelas. Uma ligação de contato aproxima o comportamento físico de uma substância de cimento extremadamente pequena aderida entre às duas partículas ligadas. A ligação de contato comporta-se essencialmente como uma ligação paralela de raio zero. Assim, uma ligação de contato não tem um raio ou rigidez normal e de cisalhamento, tal como em uma ligação paralela, e não pode resistir a um momento de flexão ou se opor ao rolamento; porém pode resistir a uma força que atua no ponto de contato. Também não é permitido o deslizamento quando uma ligação de contato está presente. A ligação de contato é definida pelos dois parâmetros de resistência à tração e de cisalhamento, $R_{n}$ e $R_{t}$, expressados em unidades de força. Quando o componente correspondente da força de contato excede qualquer um destes valores, a ligação de contato rompe, e apenas é considerada a parte baseada no grão do comportamento força-deslocamento.

2.4.2. Caracterização dos micro-parâmetros

Em geral, o modelo BPM é caracterizado pela quantidade, forma, distribuição do tamanho, empacotamento do grão e micropropriedades do grão-cimento. Cada um desses itens influencia o comportamento do modelo. A densidade do grão, $\rho$, não afeta o comportamento quasi-estático, mas é incluída para completá-lo. Com o software STAR CCM+, os pontos centrais são os grãos esféricos compostos por partículas individuais. Os diâmetros das partículas satisfazem a uma distribuição uniforme de tamanhos de partículas delimitada por $D_{min}$ e $D_{max}$. A construção do empacotamento, ou a conectividade do conjunto de partículas ligadas, é controlada pela razão ${\textstyle \left({\frac {{D}_{max}}{{D}_{min}}}\right)}$. Assim como também a distribuição utilizada para a geração do empacotamento. Para uma razão fixa, variando $D_{min}$, o tamanho da partícula absoluta é alterado, sem afetar a construção do empacotamento.

Tal caracterização separa os efeitos da construção do empacotamento e do tamanho de partícula no comportamento do material e identifica claramente $D_{min}$ como a escala de comprimento de controle do material. O último item que caracteriza o modelo BPM são as micro-propriedades do cimento-grão, expressas pelos seguintes conjuntos de parâmetros:

$\left\{{E}_{c},\left({\frac {{k}_{n}}{{k}_{s}}}\right),\mu \right\}$ micropropriedades do grão	(29)
$\left\{{\overline {E}}_{c},{\overline {\lambda }},\left({\frac {{\overline {k}}^{n}}{{\overline {k}}^{s}}}\right),{\overline {\sigma }}_{c},{\overline {\tau }}_{c}\right\}$ micropropriedades do cimento	(30)

onde ${\textstyle {E}_{c}}$ e ${\textstyle {\overline {E}}_{c}}$ são os módulos de Young do grão e do cimento, respectivamente, ${\textstyle \left({\frac {{k}_{n}}{{k}_{s}}}\right)}$ e ${\textstyle \left({\frac {{\overline {k}}^{n}}{{\overline {k}}^{s}}}\right)}$ , são, nesta ordem, as razões entre as rigidezes normal e ao cisalhamento dos grãos e cimento; ${\textstyle {\overline {\lambda }}}$ é o multiplicador de raio utilizado para definir o raio de ligação paralela definido na equação (17); ${\textstyle \mu }$ é o coeficiente de atrito do grão, e ${\textstyle {\overline {\sigma }}_{c}}$ e ${\textstyle {\overline {\tau }}_{c}}$ são as resistências à tração e ao cisalhamento do cimento, respectivamente. Na análise em seguida, os módulos do grão e do cimento estão relacionados com as rigidezes normais correspondentes, de modo que a rigidez das partículas e da ligação em paralelo são atribuídas como:

${\textstyle {k}_{n}=2\ast t\ast {E}_{c}}$ t=1 2D	(31)
${\textstyle {k}_{n}=4-R\ast {E}_{c}}$ 3D	(32)
${k}^{s}={k}_{n}({\frac {{k}_{n}}{{k}_{s}}})$	(33)
${\overline {k}}^{n}={\overline {E}}_{c}({R}^{\left(A\right)}+{R}^{\left(B\right)})$	(34)
${\overline {k}}^{s}={\overline {k}}^{n}({\frac {{\overline {k}}^{n}}{{\overline {k}}^{s}}})$	(35)

onde ${\textstyle R}$ é o raio da partícula. A deformabilidade de um material elástico linear isotrópico é descrita por duas constantes elásticas. Essas quantidades são propriedades decorrentes do modelo BPM e não podem ser especificadas diretamente. É possível, porém, relacionar os módulos de grão e cimento, ${\textstyle {E}_{c}}$ e ${\textstyle {\overline {E}}_{c}}$, respectivamente, para as rigidezes normais de modo a formar o material de cada contato como uma viga elástica com as suas extremidades nos centros de partículas, como se mostra na Figura 5. A rigidez axial de uma viga, é ${\textstyle K=AE/L}$, onde ${\textstyle A}$, ${\textstyle E}$ e ${\textstyle L}$ são, respectivamente a área da seção transversal, módulo e comprimento. Para o comportamento baseada no grão, tem-se:

${\frac {{k}_{n}}{2}}={\frac {\left(Lt\right)\ast {E}_{c}}{L}}\rightarrow {E}_{c}={\frac {{k}_{n}}{2t}},\quad \,t=$$1\quad \,2D$	(36)
${\frac {{k}_{n}}{2}}={\frac {\left({L}^{2}\right)\ast {E}_{c}}{L}}\rightarrow {E}_{c}=$${\frac {{k}_{n}}{2L}}={\frac {{k}_{n}}{4R}},\quad \,\,3D$	(37)

Assumindo que ${\textstyle {k}_{n}}$ = ${\textstyle {{k}_{n}}^{\left(A\right)}}$= ${\textstyle {{k}_{n}}^{\left(B\right)}}$ para o comportamento baseado no cimento, tem-se:

${\overline {k}}^{n}\ast A={\frac {A\ast {\overline {E}}_{c}}{L}}={\frac {A\ast {\overline {E}}_{c}}{{R}^{\left(A\right)}+{R}^{\left(B\right)}}}\rightarrow {\overline {E}}_{c}=$${\overline {k}}^{n}\ast \left({R}^{\left(A\right)}+{R}^{\left(B\right)}\right)$

(38)

Figura 5. A relação entre o módulo de elasticidade e a rigidez normal do sistema grão-cimento, (modificado- Itasca [19])

O módulo do grão depende do tamanho da partícula, e para conseguir um módulo de grão constante, as rigidezes das partículas devem ser dimensionadas com o raio da partícula. Esta análise não produz uma relação entre o coeficiente de Poisson e as rigidezes da partícula no nível micro, entretanto, um coeficiente de Poisson macroscópico será observado, e seu valor será afetado, pelo empacotamento do grão e pelas relações, ${\textstyle \left({\frac {{k}_{n}}{{k}_{s}}}\right)}$ e ${\textstyle \left({\frac {{\overline {k}}^{n}}{{\overline {k}}^{s}}}\right)}$. Para uma forma e um empacotamento de grãos fixos, o aumento dessas proporções aumenta o coeficiente de Poisson [17].

2.5. Características do software utilizado

Para a execução das simulações numéricas, optou-se pelo software STAR CCM+ da companhia CD-adapco. O STAR CCM+ é um software para análises e simulações multi-físicas e soluciona problemas de Mecânica dos Fluidos Computacional (CFD), com o auxílio dos últimos avanços em tecnologia CAE. O STAR CCM+ inclui preparação do modelo CAD, um gerador de malha (meshing) e o pós-processador.

O STAR CCM+ usa uma abordagem de célula hospedeira para detectar colisões. Duas partículas podem colidir apenas se estiverem na mesma célula. STAR CCM+ usa o mais rápido entre os seguintes algoritmos de colisão para detectar colisões em qualquer célula particular: o algoritmo NTC (No Time Counter) descrito por Schmidt [20] ou o algoritmo proposto por O’rourke [21]. Também estão disponíveis algoritmos baseados em malhas. Quanto aos modelos de contato implementados no software, que totalizam 13, estão presentes as formulações linear e de Hertz-Mindlin, além de outros que contemplam partículas de concreto, rocha, comportamentos viscosos, inelásticos e modelos definidos pelo usuário, conseguindo, assim, misturar vários modelos clássicos. O amortecimento numérico utilizado é do tipo não viscoso, e atua diretamente sobre as forças das partículas.

No STAR CCM+, o passo de tempo máximo permitido para uma partícula DEM é limitado pela suposição de que a força que atua sobre uma partícula só é afetada pelos vizinhos imediatos da partícula durante um único passo de tempo.

3. Metodologia da simulação

Para determinar os parâmetros microscópicos que reproduzirão da melhor forma os aspectos macroscópicos da amostra virtual foi seguida uma metodologia para aplicar os fundamentos do DEM. Primeiramente, gerou-se a amostra de rocha por diferentes métodos de discretização com a finalidade de aperfeiçoar a geração do empacotamento da amostra. Posteriormente, realizaram-se simulações numéricas em DEM discutindo aspectos referentes à modelagem do ensaio de compressão simples, tais como a geometria do problema, as propriedades das partículas e paredes, condições de contorno, quantidade de partículas, taxa de carregamento, etc.

3.1. Modelo de contato e parâmetros da partícula

Dentre os 13 tipos diferentes de modelos de contato com que conta o STAR CCM+, optou-se pelo uso do modelo Bonded Particle representado na Figura 6. Neste modelo as rigidezes das molas são calculadas pela formulação apresentada na secção 2.4. A escolha deste modelo se embasou em dois critérios: o primeiro deles é o fato de representar a resistência do meio mediante as ligações dos contatos, o que propicia a modelagem de rochas de forma mais realista; o segundo é a presença de elementos dissipadores de energia no contato, sendo possível diminuir a magnitude do amortecimento numérico utilizado no sistema. Para esse modelo são necessários os seguintes parâmetros: densidade, módulo de Young, coeficiente de Poisson, coeficiente de atrito, resistência normal e tangencial das ligações e coeficientes de restituição normal e tangencial.

Figura 6. Modelo utilizado nas simulações numéricas. O cimento da ligação está representado como um cilindro em 3D

3.2. Modelo geométrico e discretização da malha

A geometria tridimensional do sistema começou a ser definida no software STAR CCM+ por meio das paredes (limites) do modelo numérico. Primeiramente, definiu-se uma geometria cilíndrica de 60 mm de altura e 30 mm diâmetro com a finalidade de representar o ensaio real com maior precisão, estabelecendo-se uma relação 2:1 (Figura 7a).

Posteriormente foi realizada uma discretização (malha) da geometria de contorno. Dentre os vários métodos de discretização foi escolhido o trimmed mesh (Figura 7b), com base em várias simulações feitas no software, a partir das quais observou-se que, para o DEM, este é o método mais preciso, com convergência mais rápida.

(a)	(b)

Figura 7. (a) Modelo geométrico do cilindro. (b) Modelo géometrico discretizado mediante o método trimmed mesh

3.3. Geração da amostra e definição do ensaio virtual

Após definir os limites do modelo, foi necessário gerar as partículas que formarão parte do meio. No processo de geração da amostra, foi utilizado um método de geração aleatória de partículas tal que não exista superposição entre partículas. Em seguida, os raios das partículas são aumentados até seu valor final e o sistema é liberado para atingir o equilíbrio estático com atrito zero. Os diâmetros das partículas satisfazem a uma distribuição uniforme de tamanhos delimitados por ${\textstyle {D}_{min}}$e ${\textstyle {D}_{max}}$, para assegurar um empacotamento inicial razoavelmente denso. Neste procedimento, primeiramente foi escolhida uma quantidade de 50000 partículas. Embora o processo de geração aleatória seja muito rápido, independentemente da quantidade de partículas, para cálculos posteriores como os da compressão, 50000 partículas é um número significativamente alto, o que torna as modelagens muito complicadas e demoradas. Por isso, decidiu-se diminuir a 10% a quantidade de partículas

A amostra final foi gerada com 5.000 partículas, utilizando diâmetros aleatórios com uma distribuição normal, como apresentado na Figura 8. Os parâmetros numéricos utilizados para essa última etapa de geração estão descritos na Tabela 1.

Figura 8. Amostra final e distribuição do tamanho de diâmetro das partículas

Ao final da etapa da geração, executam-se tantos passos de tempo quanto forem necessários para que o empacotamento se estabilize, gerando assim um bom arranjo de partículas, com a finalidade de diminuir o número de partículas flutuantes (as quais possuem menos de três contatos) de forma a obter uma rede mais densa de contatos.

Tabela 1. Parâmetros utilizados na geração das amostras

Micropropriedades das partículas
Parâmetro	Símbolo	Valor	Unidade
Altura da amostra	L	60	mm
Diâmetro da amostra	W	30	mm
Número de partículas	N	5000	-
Diâmetro mínimo da partícula	${D}_{min}$	1.91	mm
Relação dos diâmetros das partículas	$\left({\frac {{D}_{max}}{{D}_{min}}}\right)$	1.252	-
Velocidade das partículas	${V}_{p}$	-5 em Z	m/s
Coeficiente de fricção estático partícula- partícula	$\mu$	0.5	-

3.4 Preparação do ensaio virtual

A partir do empacotamento gerado apresentado na Figura 8. E para preparar a amostra para o ensaio virtual, optou-se por utilizar um método específico do software STAR CCM+ chamado de “overset meshing” (sobreposição de malhas). Neste método, primeiramente é gerado um contorno muito maior que o tamanho da amostra. Este contorno foi escolhido também como um cilindro, pois o overset mesh requer que os elementos ao serem discretizados sejam compatíveis com a malha de contorno (trimmed mesh). Com este método foi discretizado um cilindro externo de dimensões 150 mm de altura e 150 mm de diâmetro; dentro deste cilindro também é discretizado um disco de dimensões 30 mm de altura e 80 mm de diâmetro que é utilizado como “placa de compressão”, como exibido na Figura 9a. As partículas são criadas num arquivo diferente de modelagem como foi explicado na seção 3.3. Os diâmetros e as posições (${\textstyle x,y,z}$) das partículas deste corpo de prova são extraídos numa tabela para depois serem importados num arquivo onde se tem o overset meshing, como é mostrado na Figura 9b.

(a)

(b)

Figura 9. (a) Overset mesh dos elementos discretizados. (b) Modelo geométrico final discretizado mediante o método overset mesh

3.5 Taxa de carregamento

As respostas obtidas com uma simulação via DEM e um material real são sensíveis à taxa de carregamento. Na maioria dos casos, escolhe-se uma taxa de carregamento que seja suficientemente pequena para garantir a solução quase estática. A resposta quase estática é obtida quando o carregamento for suficientemente lento para que as partículas se adaptem à redistribuição das forças não lineares que acompanham cada caso.

Esta resposta pode ser obtida simulando o ensaio sob deformação controlada, onde as velocidades dos muros são redefinidas como zero após cada evento não linear até atingir um novo estado de equilíbrio. Desta forma, pode-se traçar uma curva tensão-deformação completa. Pode-se garantir que uma condição quase-estática foi obtida, realizando-se o ensaio com diferentes taxas de deformação com velocidades constantes e mostrando que a resposta é a mesma para todas as velocidades abaixo de um valor de velocidade crítico.

Usualmente os ensaios de compressão reais são realizados no regime quase-estático, onde as velocidades são suficientemente pequenas para que os parâmetros de inércia não afetem os resultados. Assim, a taxa de deformação tem que ser definida com um valor o mais baixo possível. No entanto, para obter uma resposta com um tempo computacional razoável usando o STAR CCM+, é necessário sacrificar este valor usando uma maior taxa de deformação. Quando se aumenta a taxa de deformação, as curvas da resposta tensão- deformação apresentam oscilações e efeitos dinâmicos que podem afetar o comportamento do ensaio.

4. Resultados e análises

A rocha padrão para efeitos de comparação é um granito, conhecido na literatura como Lac du Bonnet granite, o qual foi testado em laboratório sob compressão simples e tração indireta por Martin [22] e Martin e Chandler [23]. A Tabela 2 mostra um resumo estatístico das principais propriedades macroscópicas obtidas por estes autores com uma bateria de ensaios experimentais que inclui 81 amostras em alguns ensaios.

Tabela 2. Resultados experimentais de Martin [22]

Propriedade	Lac du Bonnet granite
E (GPa)	69 $\pm$ 5.8 ($m=81$)
$\nu$	$0.26\pm 0.04\,\,(m=81)$
${\textstyle {q}_{u}}$ (MPa)	$200\pm 22\,(m=81)$
${\textstyle {\sigma }_{ci}}$ (MPa)	${\textstyle 90+{\sigma }_{3}}$, ${\textstyle {\sigma }_{3}<30}$
${\textstyle {\sigma }_{cd}}$ (MPa)	150, ${\textstyle {\sigma }_{3}=}$$0$
${\textstyle {\sigma }_{t}}$ (MPa)	$9.3\pm 1.3(m=39)$
${\textstyle \varphi }$ (MPa)	59
c (MPa)	30

A qualidade dos dados publicados tem sido explorada em diversos trabalhos numéricos, como em Potyondy e Cundall [17], onde estes ensaios foram simulados usando os softwares PFC2D e PFC3D. A seguir são reproduzidos os resultados experimentais, usando modelos de contato e cimentação adequados e disponíveis no software STAR CCM+, sendo finalmente comparados com outros resultados obtidos da literatura.

4.1 Simulação do ensaio de compressão simples

A partir das condiciones finais apresentadas na secção 3.4, iniciou-se a bateria de ensaios de compressão simples na amostra virtual. A codificação desses ensaios foi feita com o uso de uma parede móvel superior utilizada para aplicar as cargas. O procedimento é semelhante a um ensaio convencional, em que o corpo-de-prova é carregado por meio da aplicação de velocidade constante na parede superior (Figura 10). O ensaio de deformação controlada aplica uma deformação axial numa sequência de estágios.

Figura 10. Amostra virtual e Amostra real

Decidiu-se utilizar o modelo de contato linear junto com o BPM. Os valores dos parâmetros utilizados são apresentados na Tabela 3.

Tabela 3. Parâmetros microscópicos utilizados com modelo linear/BPM

Micropropriedades das partículas
Parâmetro	Símbolo	Valor	Unidade
Diâmetro mínimo da partícula	${D}_{min}$	1.91	mm
Relação dos diâmetros das partículas	${\frac {{D}_{max}}{{D}_{min}}}$	1.252	-
Densidade	$\rho$	2630	kg/m³
Modulo de Young das partículas	${E}_{c}$	70	GPa
Coeficiente de atrito da partícula	$\mu$	0.5	-
Relação entre a rigidez normal e cisalhante no contato em BMP	$\left({\overline {k}}^{n},{\overline {k}}^{s}\right)$	2.5	-
Relação entre a rigidez normal e cisalhante na partícula	$\left({k}_{n},{k}_{s}\right)$	2.5	-
Coeficiente normal de restituição	${C}_{n\,rest}$	0.3	-
Coeficiente tangencial de restituição	${C}_{t\,rest}$	0.3	-
Micropropriedades da ligação
Multiplicador do raio da ligação paralela	${\overline {\lambda }}$	1	-
Resistencia normal da ligação paralela, media	${\overline {\sigma }}_{c,mean}$	175	MPa
Resistencia normal da ligação paralela, desvio	${\overline {\sigma }}_{c,std}$	40	MPa
Resistencia cisalhante da ligação paralela cisalhante, media	${\overline {\tau }}_{c,\,mean}$	175	MPa
Resistencia cisalhante da ligação paralela, desvio	${\overline {\tau }}_{c,std}$	40	MPa
Velocidade da compressão	${\nu }_{comp}$	0.2	m/s
Implicit time step	-	1*10-4	seg
Time scale do DEM solver	-	0.2	-

As razões entre a rigidez normal e cisalhante dos grãos e do cimento foram as mesmas para reduzir o número de parâmetros livres. O valor da razão, ${\textstyle {k}_{n}/{k}_{s}=2.5}$ foi escolhido após algumas tentativas de modo a reproduzir o coeficiente de Poisson macroscópico observado para a rocha nos ensaios experimentais. Em seguida, o módulo dos grãos (${\textstyle E_{c}=70}$ GPa) foi escolhido após algumas tentativas para reproduzir o módulo de Young macroscópico.

Os valores médios de resistência normal e cisalhante do cimento também foram definidos iguais entre si ( ${\textstyle {\overline {\sigma }}_{c,\,mean}}$, ${\textstyle {\overline {\tau }}_{c,mean}}$=175 MPa) e controlam o valor global da resistência à compressão não confinada da amostra. Os desvios padrões das resistências normal e cisalhante também foram assumidos iguais e o valor adotado de 40 MPa foi definido de modo a corresponder à tensão de iniciação da fissura (σ_ci ver discussão abaixo). O coeficiente de atrito das partículas parece afetar apenas a resposta pós-pico e não está claro para qual valor deve ser calibrado. Assim, ${\textstyle \mu }$ = 0,5 é usado como um valor supostamente razoável.

Com esses parâmetros em nível microscópico são obtidos os resultados para a simulação numérica dos ensaios de compressão simples. Um resultado representativo após a ruptura é ilustrado na Figura 11 e o aspecto geral da superfície de ruptura é fisicamente compatível com o que se espera experimentalmente. Acredita-se que a "explosão" de algumas partículas se deva à rápida liberação de energia devido ao passo de tempo e amortecimento adotados para viabilizar as análises numéricas.

(a)

(b)

Figura 11. Modelo DEM com contato linear/BPM: (a) Seção transversal no meio da amostra, ressaltando o rompimento da mesma; (b) Representação do rompimento da amostra em 3D

Uma segunda análise dos resultados dessa bateria de ensaios se faz por meio das curvas σ versus ε, como ilustrado na Figura 12. Deste ensaio são determinados o módulo de Young e o coeficiente de Poisson macroscópicos do meio.

Figura 12. Curva de Tensão vs Deformação

As constantes elásticas macroscópicas são calculadas utilizando os incrementos de tensão e deformação ocorridos entre o início do ensaio e o ponto em que foi obtida metade da tensão de pico. Para o material 3D, o módulo de Young é calculado por:

$E={\frac {\Delta {\sigma }_{A}}{\Delta {\sigma }_{A}}}$

(39)

O coeficiente de Poisson é calculado a partir da deformação lateral e da deformação axial como:

$\nu =-{\frac {\Delta {\epsilon }_{L}}{\Delta {\epsilon }_{A}}}={\frac {1}{2}}\ast (1-{\frac {\Delta {\epsilon }_{v}}{\Delta {\epsilon }_{A}}})$

(40)

em que ${\textstyle {\epsilon }_{L}}$ é a deformação lateral e ∆ ${\textstyle {\epsilon }_{V}}$ = ∆ ${\textstyle {\epsilon }_{x}}$ + ∆ ${\textstyle {\epsilon }_{y}}$ + ∆ ${\textstyle {\epsilon }_{z}}$ é a deformação volumétrica.

4.2 Simulação do ensaio brasileiro

De forma semelhante aos ensaios de compressão simples, executou-se os ensaios de tração indireta, também conhecido como ensaio de tração por compressão diametral e ensaio brasileiro. O modelo geométrico consta de um cilindro de diâmetro 60 mm e altura 30 mm (Figura 13) e a geração da amostra é semelhante à do ensaio de compressão simples.

Figura 13. Ensaio brasileiro. Perspectiva 3D do modelo DEM e vista frontal

As micro-propriedades utilizadas são as mesmas apresentadas na Tabela 3. O teste foi executado movendo a parede rígida superior a uma velocidade constante (0,00002 m/s) até atingir a ruptura. O aspecto geral após a ruptura é ilustrado na Figura 14. Observa-se que as fissuras causadas durante o ensaio são distribuídas ao longo do plano diametral de aplicação da carga, como determina a análise teórica deste problema. A amostra DEM foi dividida por trincas induzidas por tração, como mostrado na Figura 14(c), e em boa concordância com o padrão de observado no laboratório e mostrado na Figura 14(d).

Figura 14. Ensaio brasileiro: (a) Iniciação das trincas de tração 3D. (b) Início da ruptura. (c) Padrão de ruptura observado na simulação via DEM. (d) Padrão de ruptura observado no teste de laboratório

A curva força-tempo obtida na simulação do teste brasileiro é mostrada na Figura 15. Observa-se que a força aumenta quase linearmente até um valor de pico. Esta curva permite determinar a resistência à tração do material modelado com o método de elementos discretos.

Figura 15. Curva força versus tempo no ensaio brasileiro

A resistência à tração é calculada mediante a seguinte fórmula:

${\sigma }_{t}={\frac {2\ast {F}_{f}}{\pi \ast D\ast L}}$

(41)

onde $F_{f}$ é a força de pico que atua sobre a placa diametral, e ${\textstyle D}$ e ${\textstyle L}$ são o diâmetro e a altura, respectivamente, do corpo de prova.

De posse dos valores de resistência à tração e de resistência à compressão é possível determinar os parâmetros macroscópicos de coesão e ângulo de atrito do material, com base no critério de ruptura de Mohr-Coulomb. O critério de Mohr-Coulomb é um modelo linear empírico que tem sido adotado para diferentes rochas e solos. Este critério pode ser expresso em termos de tensões principais como:

${\sigma }_{1}=2\ast c\ast {\sqrt {{k}_{p}}}+{k}_{p}\ast {\sigma }_{3}$

(42)

onde

${k}_{p}={\frac {1+\mathrm {sin} \,\varphi }{1-\mathrm {sin} \,\varphi }}$

(43)

Com as tensões de ruptura obtidas através do ensaio de compressão simples onde ${\textstyle {\sigma }_{1}=q_{u}}$ e ${\textstyle {\sigma }_{3}}$ =0 e para o ensaio brasileiro para o qual ${\textstyle {\sigma }_{1}}$ = 3* ${\textstyle {\sigma }_{3}}$ e ${\textstyle {\sigma }_{3}}$ = − ${\textstyle {\sigma }_{t}}$, tal como ilustrado na Figura 15, é possível determinar os parâmetros de coesão (${\textstyle c}$) e ângulo de atrito (${\textstyle \varphi }$) do material, da seguinte forma:

${k}_{p}={\frac {{q}_{u}-3\ast {\sigma }_{3}}{{\sigma }_{t}}}$

(44)

onde obtém-se o valor da coesão como;

$c={\frac {{q}_{u}}{2\ast {\sqrt {{k}_{p}}}}}$

(45)

e do ângulo de atrito:

$\varphi =\mathrm {{sin}^{-1}} \,\left({\frac {{k}_{p}-1}{{k}_{p}+1}}\right)$

(46)

A Figura 16 mostra a envoltória de resistência obtida considerando o critério de ruptura de Mohr-Coulomb. Estes valores serão apresentados na Tabela 4.

Figura 16 Envoltória de resistência de Mohr-Coulomb

4.3. Apresentação e discussão dos resultados quantitativos

Nesta seção discute-se quantitativamente os resultados das análises numéricas tanto do ensaio de compressão simples quanto de tração indireta, comparando-se com os valores obtidos experimentalmente e numericamente por outros autores. A Tabela 4 repete os dados experimentais de Martin [22] e resultados numéricos obtidos por Potyondy e Cundall [17] com o programa PFC e acrescenta na última coluna os dados obtidos neste trabalho com o programa STAR CCM+.

Em cada coluna da Tabela 4 os valores dos parâmetros macroscópicos são mostrados em termos de média mais ou menos o desvio padrão. Foram usadas $m=81$ amostras em alguns resultados experimentais, $m=10$ amostras nas simulações numéricas com o PFC e $m=6$ amostras nas simulações numéricas com o STAR CMM+. Ressalta-se que a aleatoriedade dos resultados numéricos advém do próprio modelo BPM que adota um input em forma de distribuição normal (média e desvio padrão) para os parâmetros microscópicos da cimentação nas direções normal e tangencial com mostrado na Tabela 3. Entretanto, o pacote de partículas usado em cada simulação é sempre o mesmo, mudando apenas a distribuição espacial de resistência da cimentação.

Tabela 4 Macropropriedades dos modelos e do Lac du Bonnet granite

Propriedade	Lac du Bonnet granite	PFC2D model, Dmed=0.72 mm	PFC3D model, Dmed=1.53 mm	STAR CCM+3D model, Dmed=2.12 mm
E (GPa)	69 ${\textstyle \pm }$ 5.8 (m=81)	70.9 ${\textstyle \pm 0.9\,(m=10)}$	69.2 ${\textstyle \pm }$ 0.8 ${\textstyle (m=10)}$	69.07 ${\textstyle \pm }$ 1.38 ${\textstyle (m=6)}$
${\textstyle \nu }$	0.26 ${\textstyle \pm }$ 0.04 (m=81)	0.237 ${\textstyle \pm }$ 0.011 ${\textstyle (m=10)}$	0.256 ${\textstyle \pm }$ 0.014 ${\textstyle (m=10)}$	0.255 ${\textstyle \pm }$ 0.0145 ${\textstyle (m=6)}$
${\textstyle {q}_{u}}$ (MPa)	200 ${\textstyle \pm }$ 22 (m=81)	199.1 ${\textstyle \pm }$ 13 ${\textstyle (m=10)}$	198.8 ${\textstyle \pm }$ 7.2 ${\textstyle (m=10)}$	189.42 ${\textstyle \pm }$ 3.6 ${\textstyle (m=6)}$
${\textstyle {\sigma }_{ci}}$ (MPa)	90 ${\textstyle +}$${\sigma }_{3}$, ${\textstyle {\sigma }_{3}}$<0	71.8 ${\textstyle \pm }$ 21.8 ${\textstyle (m=10)}$	86.6 ${\textstyle \pm }$ 11 ${\textstyle (m=10)}$	93.08 ${\textstyle \pm }$ 3.92 ${\textstyle (m=6)}$
${\textstyle {\sigma }_{cd}}$ (MPa)	150, ${\textstyle {\sigma }_{3}}$=0	NA ${\textstyle (m=10)}$	190.3 ${\textstyle \pm }$ 7.5 ${\textstyle (m=10)}$	130.28 ${\textstyle \pm }$ 5.87 ${\textstyle (m=6)}$
${\textstyle {\sigma }_{t}}$ (MPa)	9.3 ${\textstyle \pm }$ 1.3 (m=39)	44.7 ${\textstyle \pm }$ 3.3 ${\textstyle (m=10)}$	27.8 ${\textstyle \pm }$ 3.8 ${\textstyle (m=10)}$	38.49 ${\textstyle \pm }$ 8.87 ${\textstyle (m=6)}$
φ (deg)	59	29.5 ${\textstyle \pm }$ 4.8 ${\textstyle (m=10)}$	32.1 ${\textstyle \pm }$ 2.4 ${\textstyle (m=10)}$	17.78 ${\textstyle \pm }$ 12.2 ${\textstyle (m=6)}$
c (MPa)	30	58.5 ${\textstyle \pm }$ 8.5 ${\textstyle (m=10)}$	55.1 ${\textstyle \pm }$ 4.2 ${\textstyle (m=10)}$	70.64 ${\textstyle \pm }$ 20.35 $(m=6)$

Comparando os resultados numéricos com os experimentais, percebe-se a excelente concordância para os parâmetros obtidos do ensaio de compressão simples, principalmente para o módulo de Young ($E$), o coeficiente de Poisson ($\nu$) e a resistência à compressão simples ($q_{u}$). Isto não é uma surpresa justamente porque os valores dos parâmetros microscópicos de módulo de Young ($E_{c}$) e razão de rigidez ($k_{n}/k_{s}$) dos grãos e de resistência da cimentação, apresentados na Tabela 3, foram calibrados justamente para tentar reproduzir os parâmetros macroscópicos observados no granito de Lac du Bonnet. A comparação dos demais parâmetros (${\textstyle {\sigma }_{ci}}$ e ${\textstyle {\sigma }_{cd}}$), obtidos do ensaio de compressão simples requer um pouco mais de detalhamento teórico.

A relação tensão-deformação de uma rocha sujeita a testes de compressão pode ser classificada em cinco regiões, de acordo com vários estudos [24-28]. Estas regiões são apresentadas na Figura 17 e descritas a seguir:

1. Região de fechamento das fissuras. O estágio inicial da curva tensão- deformação exibe deformabilidade não linear, devido ao fechamento de microfissuras pré-existentes na amostra, e depende da densidade inicial da fissura e da sua geometria.

2. Região de deformação elástica. Uma vez que as fissuras pré-existentes estão fechadas, as deformações axial, lateral e volumétrica aumentam linearmente com a tensão axial. Assim, a deformação elástica da amostra de rocha intacta domina esta parte das curvas tensão-deformação.

3. Início do fissuramento e região de crescimento de fissuras estáveis. Esta região é caracterizada pelo início de novas fissuras e pela propagação controlada das fissuras pré-existentes. Alguns autores [22,23] relataram que esta região começa para um nível de tensão igual a 30-50% do pico de tensão, e referiram-se a essa tensão, como tensão de início da fissura ${\textstyle {\sigma }_{ci}}$. Além disso, no ponto final desta região pode-se identificar a transição da deformação volumétrica de compressão para dilatação.

4. Fissura de dano e região de crescimento instável. Esta região foi primeiramente definida para identificar onde a estrutura do espécime da rocha muda, devido à dilatância volumétrica que é gerada. A densidade de fissuras aumenta rapidamente e pode atingir cerca de sete vezes o valor das fissuras pré-existentes [29]. Este aumento imediato resulta na ocorrência da inversão da deformação volumétrica em consequência de a deformação lateral crescer tão rapidamente que excede a deformação axial. Tal mudança ocorre frequentemente em um nível de tensão axial da ordem de 70-85% da tensão de pico, e é referida como a tensão de dano ${\textstyle {\sigma }_{cd}}$.

5. Região pós-pico. O início desta região é marcado pelo pico de tensão ${\textstyle {\sigma }_{f}}$ que é usado para estabelecer a envoltória da curva de resistência à ruptura. Nesta fase, inicia-se a superfície de fratura macro e ocorre ruptura da amostra.

Figura 17 Relação tensão - deformação obtida experimentalmente a partir do ensaio de compressão simples no granito Lac du Bonnet [23]

Os três níveis de tensão característicos, ou seja, a tensão de iniciação da fissura ${\textstyle {\sigma }_{ci}}$, a tensão de dano ${\textstyle {\sigma }_{cd}}$ e o pico de tensão ${\textstyle {\sigma }_{f}}$, são muito importantes na relação tensão- deformação, mas é muito difícil detectá-los em laboratório, exceto pelo pico de resistência. O resultado experimental de Martin & Chandler [23] para uma das amostras ensaiadas, representativa da média, e mostrado na Figura 17. Na curva com abscissa para a direita e ordenada para cima tem-se a tradicional relação entre a tensão axial e a deformação axial. Nesta porção identifica-se claramente uma resistência de pico em torno de 210 MPa e a deformação axial na ruptura no pico em torno de 0,36%. No geral a curva é praticamente linear com uma leve curva (para cima) no início e uma curvatura (para baixo) antes do pico.

Esta última não linearidade fica mais perceptível quando se plota a tensão axial versus a deformação lateral (curva com a abscissa para a esquerda). Porém é na curva de deformação volumétrica versus deformação axial que se distingue mais claramente as cinco regiões descritas anteriormente. A zona I de fechamento de trincas pré-existentes ocorreu até cerca de 50 MPa; a zona II de linearidade elástica até cerca de 90 MPa, quando se inicia a formação de trincas por cisalhamento das ligações; a zona III de crescimento estável das trincas vai até cerca de 165 MPa; e a zona IV de dano instável até o pico, é marcada claramente pela dilatância na curva de variação volumétrica.

Na Figura 18 mostra-se o resultado da simulação numérica, fazendo um contraponto com a Figura 17. Acrescentou-se ainda uma nova curva em vermelho, com ordenada à direita, que mostra a evolução do número de ligações quebradas versus o avanço das deformações axiais. O módulo de Young global, o pico de resistência e o início das fissuras mostram uma boa concordância com os resultados experimentais, dentro da faixa variação dos resultados numéricos e reais. Pode-se observar que a ruptura macroscópica da rocha real é causada pela iniciação e propagação de microfissuras que podem ser correlacionadas com a ruptura progressiva das ligações à medida que se aumentam as forças externas. Quando um contato rompe, as partículas que o formavam ficam livres e deslizam umas em relação às outras, dependendo apenas da força gerada pelo atrito entre elas. Desta forma, a facilidade de deslizamento destas partículas sem ligação no contato faz com que a deformação dos contatos próximos seja maior, podendo romper progressivamente, até gerar uma zona de ruptura generalizada e um mecanismo de separação da amostra em partes distintas, uma que permanece praticamente como um corpo rígido e outra que desliza e se fragmenta progressivamente até a falha final.

Figura 18 Diagramas tensão-deformação obtidos numericamente do ensaio de compressão simples DEM, mostrando também o número de ligações quebradas

O processo de ruptura progressiva da amostra simulada via DEM também pode ser classificado em cinco fases distintas, como foi descrito por Shen [30]:

Fase I: Praticamente nenhuma carga axial é absorvida enquanto o pistão desce inicialmente. Esta fase inicial no modelo numérico corresponde à acomodação da placa superior até mobilizar um número razoável de partículas para aplicar uma carga que possa ser considerada uniforme. Também podem ocorrer pequenas acomodações de partículas no interior do corpo de prova não se percebe a curva côncava típica do fechamento inicial de trincas em amostras reais, uma vez que não foram incluídas trincas (ou ligações rompidas) pré-existentes na amostra numérica.

Fase II: Não são gerados danos na amostra, uma vez que o número de ligações quebradas permanece nulo nesta fase. Nesta fase, a capacidade de deformação elástica das cimentações desempenha um papel dominante, porque as relações entre a tensão e as deformações lateral, volumétrica e axial são todas lineares. A zona de linearidade elástica foi observada até cerca de 90 MPa, como no ensaio real.

Fase III: O início das quebras das ligações cimentícias indica o começo da terceira região e a tensão correspondente é referida como tensão de iniciação da fissura, ${\textstyle {\sigma }_{ci}}$, nas simulações via DEM. Como resultado da geração de fissuras, a relação tensão versus deformação lateral começa a apresentar não linearidade. Além disso, a deformação volumétrica atinge o seu valor máximo no final desta fase, devido à taxa crescente de deformação lateral que excede a deformação axial.

Fase IV: O nível de tensão axial no momento da reversão da deformação volumétrica marca o início desta fase, e é chamado de tensão de dano por fissuras, ${\textstyle {\sigma }_{cd}}$. Este ponto está em torno de 130 MPa no modelo numérico. Além disso, pode ser observado pela inclinação da curva do número de ligações quebradas-deformação axial, que as fissuras se propagam rapidamente, o que está em boa concordância com a principal característica da região de fissura instável em rochas reais. Nota-se claramente a reprodução da dilatância também na simulação numérica.

Fase V: É definida pelo o início do comportamento pós-pico. Após o pico, as respostas de deformação lateral e de deformação volumétrica exibem uma dilatância contínua, até atingir uma resistência residual determinada por atrito. A quebra de ligação cessa. Não foi identificado um estado crítico com estabilização das deformações, pois a amostra rompe fisicamente com um comportamento frágil.

Quanto ao ensaio de tração indireta por compressão diametral, o valor médio calculado de resistência (${\textstyle {\sigma }_{t}=38.5}$ MPa) é 4.1 vezes superior ao valor médio experimental (9,3 MPa). Isto ocorre porque os parâmetros micoscópicos foram calibrados tendo como principal preocupação reproduzir os valores macroscópicos do ensaio de compressão simples. Para simplificar o número de combinações durante a calibração, assumiu-se a mesma razão de rigidez normal e cisalhante para os grãos e para o contato, bem como a mesma resistência normal e cisalhante para a cimentação.

A relação entre a resistência à compressão e a resistência à tração do granito de Lac du Bonnet era da ordem de 21,5 vezes (200/9.3), enquanto que a razão calculada foi de apenas 4,9 vezes (189/38.5). Este valor se deve principalmente ao fator de se ter utilizado o mesmo valor (175 MPa) de resistência microscópica em tração e em cisalhamento para a cimentação. Desta forma estima-se que a resistência à tração da cimentação deveria ser reduzida tentativamente por um valor da ordem de 5 vezes para algo em torno de 35 MPa.

Os valores de resistência à tração obtidos com o programa STAR-CCM+ são compatíveis com os valores numéricos determinados por Potyondy e Cundall [17], que também não logrou bons resultados pelo mesmo motivo. As diferenças entre as análises numéricas devem- se à adoção de diferentes diâmetros médios para os grãos. É conhecida na literatura a limitação do modelo BPM para reproduzir valores compatíveis de resistência a compressão e a tração, quando estes diferem da ordem de 20 vezes. Isto se deve ao fato de que as ligações não transmitem momentos após a quebra [31]. Para sanar este problema, outros modelos foram propostos recentemente, tal como FJM (Flat Joint Model) ou a utilização de clusters de partículas com diferentes valores de resistência de cimentação.

O elevado valor numérico da resistência à tração tem implicações diretas no cálculo dos parâmetros globais de resistência ao cisalhamento, coesão e ângulo de atrito, calculados de acordo com as equações 5.7 e 5.8, com base apenas em ensaios de tração e de compressão. Altos valores de resistência à tração levam a círculos de Mohr com maiores diâmetros, uma vez que a tensão desvio é quatro vezes o valor de ${\textstyle {\sigma }_{t}}$ no ensaio de tração indireta. Isto implica em maiores valores de coesão e menores valores de ângulo de atrito quando conjugado com os valores do ensaio de compressão simples. Com base apenas nos dados experimentais de ${\textstyle q_{u}=200}$ MPa e ${\textstyle {\sigma }_{t}=9.3}$ MPa, pode-se estimar uma coesão $c=23$ MPa e um ângulo de atrito de ${\textstyle \varphi =64}$°. Estes valores são compatíveis com os valores obtidos de ensaios triaxiais e mostrados na Tabela 4, mas o atrito observado parece um tanto elevado.

4.4 Análise paramétrica da rigidez

Os resultados macroscópicos da simulação numérica dependem diretamente dos parâmetros físicos e geométricos dos grãos que constituem as partículas. Em especial os módulos de Young dos grãos (${\textstyle E_{c}}$) são usados diretamente para calcular os parâmetros de rigidez dos contatos. O tamanho dos grãos também tem uma grande influência nos valores globais de resistência e rigidez da amostra, uma vez que o raio da partícula afeta diretamente a área da seção transversal da cimentação, através do parâmetro multiplicador do raio da ligação paralela (${\textstyle {\overline {\lambda }}}$), como mostrado na equação (17). Desta forma nesta seção mostra-se uma análise paramétrica do efeito de ${\textstyle E_{c}}$.

Nas Figuras 19(a) e (b), mostra-se como o efeito de uma micro-propriedade (neste caso, o módulo de Young das partículas) atua sobre uma macro-propriedade (neste caso, o valor máximo da resistência a compressão simples (${\textstyle q_{u}}$) do material, e sua deformação no ponto de ruptura).

(a)

(b)

Figura 19. (a) Curva de resistência a compressão simples versus Módulo de Young das partículas. (b) Curva de deformação máxima versus Módulo de Young das partículas

Observa-se que a resistência última (${\textstyle q_{u}}$) aumenta com o aumento do módulo de Young das partículas. Para reproduzir o valor experimental (${\textstyle q_{u}=200\pm 22}$ MPa), o módulo da partícula da rocha analisada deve estar na faixa de 50 GPa (ver Figura 18(a)) até cerca de 100 GPa (estimado). O valor adotado de ${\textstyle E_{c}=70}$ GPa é bem ajustado.

Já na Figura 18(b), percebe-se que deformação axial necessária para atingir a resistência de pico diminui rapidamente à medida que a rigidez das partículas aumenta. Isto se explica pelo modelo de contato adotado, uma vez que as rigidezes normal e cisalhante da amostra aumentam diretamente com os módulos elásticos equivalentes das partículas em contato. Neste caso logrou-se um bom ajuste em forma de potência, mostrado na figura, com um alto coeficiente de determinação (${\textstyle R^{2}=0.98}$).

5. Conclusões

Como apresentado neste trabalho, uma das principais vantagens da abordagem com elementos discretos é que o comportamento constitutivo empírico complexo pode ser substituído por uma lógica mais física e fundamental de contato direto e cimentação entre partículas. Neste contexto, quando a resistência à tração ou ao cisalhamento da zona cimentante é excedida, as ligações entre as partículas rompem. Uma vez que as ligações são quebradas, o campo de deslocamento numa formulação de elemento discreto é controlado pela geometria de partículas e parâmetros do modelo de contato suave (interpenetrações). Este conceito permite simular a inicialização e propagação de danos localizados numa amostra de rocha, reproduzindo adequadamente o fenômeno de ruptura. O processo de dano evolui progressivamente e naturalmente a cada ruptura de contato, o que leva a uma redistribuição das tensões até que haja uma densidade crítica de quebra de ligações que permita que uma superfície de ruptura se propague através de todo o corpo.

Entretanto, os modelos discretos via DEM para materiais granulares apresentam três dificuldades intrínsecas: (1) geração do estado inicial de empacotamento das partículas e distribuição de porosidade; (2) determinação da deformabilidade das próprias partículas; e (3) a necessidade de calibração dos parâmetros para modelos de partículas contra experimentos. O efeito da deformação das partículas, até mesmo simulando a quebra das mesmas, pode ser introduzido nas abordagens DEM, mas o custo computacional aumenta consideravelmente, já que o número de partículas para problemas práticos quase sempre é muito grande. O ponto mais crítico é a necessidade de calibração de parâmetros de partículas, tais como distribuição de tamanho, rigidez de contato ou parâmetros de ligação (atrito, coesão, força de ligação, etc.). Uma vez que esses parâmetros não são propriedades geotécnicas convencionais e não podem ser facilmente medidos na escala necessária em laboratório, a maneira usual de superar essa dificuldade é simular o comportamento macroscópico medido das rochas, como as curvas de tensão-deformação obtidas a partir de testes laboratoriais usando amostras de rocha de tamanho padrão. Os parâmetros em nível microscópico são, então, ajustados em aproximações sucessivas até reproduzir as curvas experimentais. Além dos parâmetros físicos também surgem vários parâmetros de controle não físicos, como o amortecimento numérico, de modo a viabilizar o tempo computacional da simulação. A depender dos modelos de contato e de cimentação adotados podem surgir várias combinações e o processo de calibração pode demandar mais tempo do que uma aplicação em escala real.

Após várias tentativas de geração do empacotamento inicial, conclui-se que a forma mais adequada para geração dos corpos de prova e condições de contorno dos ensaios pretendidos era acoplar as técnicas de “overset meshing” e “trimmed mesh”. No primeiro (overset meshing) gera-se um cilindro maior que controla a discretização de fato, o segundo método (trimmed mesh) discretiza o contorno cilíndrico que define a amostra do corpo de prova. Em seguida procedeu-se o empacotamento das partículas. O número final e o tamanho das partículas tiveram que ser ajustados para 5000 partículas com diâmetro médio de 2.12 mm de modo a viabilizar computacionalmente as análises. Vale a observação que, recentemente, novos métodos de empacotamento de partículas vem sendo desenvolvidos [32-34], os quais possuem o potencial de, no futuro próximo, tornar mais simples, preciso e prático o processo de geração de malhas de partículas para simulações numéricas utilizando o método dos elementos discretos.

O próximo passo importante foi a definição do modelo de contato entre as partículas. Tentativas iniciais com o modelo de Hertz-Mindlin não lograram uma boa reprodução dos resultados experimentais. Por fim decidiu-se pelo modelo de contato linear acoplado ao modelo de cimentação BPM. A calibração final se deu para um número limitado de parâmetros com base nos ensaios de compressão simples.

O modelo adotado conseguiu reproduzir adequadamente todos os estágios de evolução de dano em amostras de rochas sujeitas a ruptura por compressão simples. Neste aspecto, auxilia muito a visualização da evolução do número de partículas com ligações quebradas ao longo do ensaio. Identifica-se claramente a zona de progressão instável de danos por meio da aceleração das deformações laterais e tendência de dilatância com a evolução das quebras até atingir a resistência de pico. O comportamento pós-pico é típico de uma rocha frágil.

Em seguida procedeu-se a simulação do ensaio brasileiro de tração indireta usando os mesmos parâmetros calibrados do ensaio de compressão simples. Neste caso, consegue-se também uma boa representação qualitativa do mecanismo de inicialização e propagação de fissuras. Entretanto o valor numérico obtido para a resistência à tração da amostra foi cerca de 4 vezes superior em relação ao valor experimental. A incongruência no valor calculado de resistência a tração se refletiu numa sobrestimação no parâmetro de coesão da resistência ao cisalhamento da rocha e uma subestimação no valor de ângulo de atrito. Estas discrepâncias devem-se ao fator de se ter adotado o mesmo valor para os parâmetros microscópicos de resistência da cimentação sob tração e cisalhamento (175 MPa). Em retrospectiva, percebe-se que o parâmetro de resistência à tração do modelo de cimentação deveria ser minorado por um fator de 5 vezes de modo a manter a mesma relação entre resistência a compressão e tração macroscópica observada nos ensaios experimentais.

Pode-se concluir que o modelo de partículas cimentadas apresentado é uma ferramenta de pesquisa para auxiliar na compreensão de processos de ruptura frágeis. Foi demonstrado que proporcionou uma reprodução adequada do comportamento macroscópico das rochas de forma quantitativa e qualitativa. Além disso, o estudo oferece um procedimento de calibração direto que é extensível ao uso em outras simulações microestruturais baseadas em descontinuidades. Ele também permite estabelecer expressões de regressão para estimar as propriedades macroscópicas do material em função propriedades das partículas.

A modelagem de partículas ligadas ou cimentadas está em uma fase inicial de desenvolvimento. Atualmente se tem um modelo genérico para rocha intacta. Futuros modelos de partículas cimentadas devem ser construídos de acordo com tipos de rocha particulares, como granito, rocha calcária ou arenito.

Finalmente, fez-se uma análise paramétrica da influência do módulo de Young em algumas propriedades macroscópicas adotadas em modelos fenomenológicos do tipo elástico perfeitamente plástico. Desta análise conclui-se que a resistência a compressão da amostra aumenta e a deformação na ruptura diminui com o aumento da rigidez dos grãos.

Agradecimientos

Ao Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), à Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) e à empresa Furnas Centrais Elétricas S.A., pelo apoio financeiro.

Referências

[1] Cundall P.A. A computer model for simulating progressive, large-scale movements in blocky rock systems. In Symp. Int. Soc. Rock Mech., 2:132-150, 1971.

[2] Cundall P.A., Hart R.D. Numerical modeling of discontinua. Engineering Computations, 9:101-113, 1992.

[3] Ferreira S.T. Estudo e implementação do Métod dos Elementos Discretos em 3D. Rio de Janeiro: Instituto Militar de Engenaria. 2009.

[4] O' Sullivan C. Particulate discrete element modelling. London and New York: Spon Press, 2011.

[5] Ayquipa C.J. Simulação 3D pelo método dos elementos discretos de refluxo de material de sustentação de fraturas em poços de petróleo. Dissertação de Mestrado pelo Programa de Pós-Graduação em Engenharia Civil da PUC-Rio, 2008.

[6] Pinto C.N. Uso de elementos discretos na modelagem numérica da perfuração de poços de petróleo por brocas PDC. Dissertação de Mestrado, Pontifícia Universidade Católica do Rio de Janeiro, 2011.

[7] Morales I.P. Contributions to the generalization of advancing front particle packing algorithms. Numerical Methods in Engineering, 107:993-1008, 2016.

[8] Labras C.A. Advances in the development of the discrete element method for excavation processes. Tese de Doutorado pela Universitat Politécnica de Catalunya, 2012.

[9] Walton O.R. Numerical simulation of inelastic, frictional particle-particle interactions, in Particulate Two-Phase Flow. M.C. Roco, Ed., Butterworth–Heinemann, Stoneham, MA, p. 884-911, 1993.

[10] Bathurst R., Rothenburg L. Micromechanical aspects of isotropic granular assemblies with linear contact interactions. Journal of Applied Mechanics, ASME, 55(1):17-23, 1988.

[11] Chang C., Misra A. Theoretical and experimental study of regular packing of granules. Journal of Engineering Mechanics, ASCE, 115(4):704-720, 1989.

[12] Chang C., Misra A. Packing structure and mechanical properties of granulates. Journal of Engineering Mechanics, 116(5):1077-1093, 1990.

[13] Neves C.E.V. Comportamento de materiais granulares usando o método dos elementos discretos. Dissertação de Mestrado, Publicação G.DM-173/09, Departamento de Engenharia Civil, Universidade de Brasília, Brasília, DF, 166 p, 2009.

[14] Albuquerque B.C.P. Influência do ângulo de atrito da partícula nos parâmetros macroscópicos de resistência do meio no estado plano de deformações. Distrito Federal. Dissertação de Mestrado, Publicação G.DM-276/16, Departamento de Engenharia Civil e Ambiental, Universidade de Brasília, Brasília, DF, 87p, 2016.

[15] Tedesco B.M.M.S. Modelagem numérica de ensaios de enrocamento. Distrito Federal. Dissertação de Mestrado, Publicação G.DM-277/16, Departamento de Engenharia Civil e Ambiental, Universidade de Brasília, Brasília, DF, 87p, 2016.

[16] Robertson D., Bolton M. DEM simulations of crushable grains and soils. In: Kishino Y, editor. Powders and grains. Lisse. The Netherlands: Balkema. pp. 6–23. 2001.

[17] Potyondy D., Cundall P. A bonded-particle model for rock. International Journal of Rock Mechanics and Mining Sciences, 41(8):1329, 2004.

[18] Autio J., Wanne T., Potyondy, D. Particle mechanical simulation of the effect of schistosity on strength and deformation of hard rock. In: NARMS-TAC 2002. Toronto: University of Toronto Press. pp. 275–82, 2002.

[19] ITASCA. PFC2D 2.0 Particle Flow Code in Two Dimensions, Theory and Background. 2008.

[20] Schmidt D.P., Rutland C.J. A new droplet collision algorithm. Journal Comput Phys, 164:62–80, 2000.

[21] O’Rourke P.J. Collective Drop Effects on Vaporizing Liquid Sprays. Department of Mechanical and Aerospace Engineering, Princeton University, 1981.

[22] Martin C.D. The strength of massive Lac du Bonnet granite around underground openings. Tese de Doutorado, University of Alberta, 1993.

[23] Martin C.D., Chandler N. A. The Progressive Fracture of Lac du Bonnet Granite. International Journal of Rock Mechanics and Mining Sciences, 31(6):643–659, 1994.

[24] Tapponnier P., Brace W.F. Development of stress-induced microcracks in westerly granite. Int J Rock Mech Min Sci, 13:395–430, 1976.

[25] Goodman R.E. Introduction to Rock Mechanics. Wiley, New York, 1989.

[26] Lajtai E.Z., Carter B.J., Scott D.E.J. Mapping the state of fracture around cavities. Engineering Geology, 31, N. 3–4, 277–289, 1991.

[27] Gatelier N., Pellet F., Loret B. Mechanical damage of an anisotropic porous rock in cyclic triaxial tests. International Journal of Rock Mechanics and Mining Sciences, 39(3):335–354, 2002.

[28] Taheri A., Yfantidis N., Olivares C.S., Connelly B.J., Bastian T.J. An experimental study on degradation of mechanical properties of sandstone under different cyclic loadings. Geotechnical Testing Journal, 39(4):673-687, 2016.

[29] Hallbauer D., Wagner H., Cook N.G.W. Some observations concerning the microscopic and mechanical behavior of quartzite specimen in stiff, triaxial compression tests. Int J Rock Mech Min Sci, 10:713–26, 1972.

[30] Shen Z., Jiang M., Wan R. Numerical study of inter-particle bond failure by 3D discrete element method. International Journal of Rock Mechanics & Mining Sciences, 40(4): 175–189, 2016.

[31] Wu S., Xu X. A study of three intrinsic problems of the classic discrete element method using flat-joint model. Rock Mech Rock Eng, 49(5):1813-1830, 2015.

[32] Morales I.P., Valera R.R., Morfa C.A.R., Farias M.M. Dense packing of general-shaped particles using a minimization technique. Computational Particle Mechanics, 4(2):165–179, 2017.

[33] Morfa C.A.R., Cortés L.A., Farias M.M., Morales I.P.P., Valera R.R., Oñate E. Systemic characterization and evaluation of particle packings as initial sets for discrete element simulations. Computational Particle Mechanics, 5(3):319-334, 2018.

[34] Morfa C.A.R., Morales I.P.P., Farias M.M., Oñate E., Valera R.R., Casañas H.D.G. General advancing front packing algorithm for the discrete element method. Computational Particle Mechanics, 5(1):13-33, 2018.