An Algorithm to compute damage from load in composites

Abstract

We present a new method to model fracture of concrete based on energy minimisation. The concrete is considered on the mesoscale as composite consisting of cement paste, aggregates and micro pores. In this first step, the alkali-silica reaction is taken into account through damage mechanics though the process is more complex involving thermo-hygro-chemo-mechanical reaction. We use a non-local damage model that ensures the well-posedness of the boundary value problem (BVP). In contrast to existing methods, the interactions between degrees of freedom evolve with the damage evolutions. Numerical results are compared to analytical and experimental results and show good agreement.

Keywords

Concrete / damage / prediction / modelling / energy minimisation / ASR

Introduction

Concrete can be considered as two or three phase composite depending whether or not the interface between the cement paste and the aggregates are considered or not. Computational modeling of damage in concrete materials at the mesoscopic scale is challenging. The first challenge originates from the complex microstructure. The second challenge is related to the complex damage process of concrete at the mesoscopic scale. Besides other, a specific problem of this damage process is the alkali-silica-reaction. In this phenomenon, the aggregates locally dissolve in contact with the cement pore solution to form an expansive gel which induces cracking in the inclusions and the paste. Although this process is a highly complex thermo-hygro- chemo-mechanical problem, we focused on the mechanical degradation due to the formation of the gel pockets. To solve this problem a novel numerical technique for attributing damage was developed.

Rate independent materials undergoing strain softening manifest bands of infinite strain when strong ellipticity of the linearized boundary value problem (BVP) is lost. Though a solution can be obtained when the region of strain softening is restricted to a set of Lebesgue measure zero, difficulties occur in numerical analysis since the damage localizes in a single element leading to mesh dependent results under mesh refinement. The difficulties that emerge with the material instabilities due to the ill-posedness of the BVP can be avoided by regularising the governing equation so that the BVP remains well posed. Typical regularisation techniques are higher order continuum models such as gradient-enhanced models and polar theories (the most well-known theory is probably the Cosserat continuum), non-local models, viscous models and cohesive zone models. All these regularisation techniques introduce a characteristic length into the discretisation that can absolutely be physically motivated, see Ref. [1,2].

Cohesive zone models are frequently applied in the context of partition of unity (PU) enriched methods such as GFEM [3-5], XFEM and certain mesh-free methods [6-8]. PU enriched methods require algorithms to describe and track the crack surface and are complex when many cracks need to be modeled.

Gradient-enhanced models are typically described by differential equations that contain higher order spatial derivatives. However, such frameworks pose problems at interfaces between materials, as they implicitly or explicitly introduce weighing functions which are geometrically based. Moreover, they require higher order spatial derivatives. Ref. [2] refers gradient-enhanced models as weakly non-local models.

The introduction of a viscosity can also restore the well-posedness of the BVP. It can be regarded as introducing higher order time derivatives, similar to the gradient models [2]. However, often a too large amount of viscosity is needed to keep the BVP well posed since viscous models restore the well-posedness of the BVP only in a certain time scale.

Strongly non-local models are models of the integral type, see Fig. 2. In non-local formulations, the stress at a given material point does not only depend on that local strain but also on the strain of the neighboring points. This effect is obtained by averaging certain internal variables through a smoothing function in a given neighborhood, Fig. 2. However, the introduction of a smoothing function imposes the scale of the material, and thus the lower limit on element size: beyond a certain point, refining the mesh cannot significantly improve the quality of the solution. We propose a third way of distributing damage, based on the computation of successive equilibriums.

The relaxation of the elements where the failure criterion is the highest causes the value of the criterion to become lower in these elements as well as in its neighborhood. However, the elements in the neighborhood will see their criteria lower less than in the relaxed elements, until such point that the maximum criterion value is found in more than an element. In turn, the damage is incremented in the set of elements with maximum criterion until such point that the set is increased. Such a process causes the damage to be spread even without an explicit weighing function. However it requires the damage state to converge monotonically to its final equilibrium value with small enough increments that this process of damage “flooding” can occur. Monotonic convergence is guaranteed since the damage can only increase monotonically.

This process minimises dissipated energy. Damage state of the material correspond to energy spent in the transformation. These energy costs allows shifting between successive material states. From a given state, physically, one can only move to a state which is of lower elastic energy (i.e. larger dissipated energy). When minimal-energy solutions are sought, from a given state, for a given energy increment used in transforming the material, only the states accessible through continuously decreasing energy are thermody namically possible [9]. This process is reproduced in the damage flooding method outlined above, as only the minimum domain where damage can be incremented is affected.

In this article, we present a method based on strongly non-local methods. Once the non-local problem posed, it is solved by the application of a Polak-Ribière or modified Newton-Raphson method. In classical methods, the stress train relationship does not change. Even if it is non-local, the interactions between the different degrees of freedom in the problem are constant. In the proposed algorithm, the sets of interacting degrees of freedom are computed ‘dynamically’ as the state of the material evolves: it is therefore a non-local damage model solved using a modified newton method.

In the first part of the article, the behavior laws considered are described as well as the geometry and boundary conditions of the problems considered. In the second part, the proposed algorithm is described and a proof of convergence is proposed. In the third part, solutions to simple problems obtained with classical algorithms and with the proposed one are compared. Finally, some results obtained for the complex case for which the algorithm was designed are shown.

Governing Equations

The equilibrium equation in the absence of body forces is given by

(1)\nabla ·\sigma =\mathrm{0}\mathrm{o}\mathrm{r}\mathbf{f}=\mathrm{0},

where σ is the Cauchy stress tensor and f is the vector containing internal and external forces. In small strain elasticity, the stress-strain relationship is given by

(2)\mathit{\sigma }=\mathbf{E}:є,

with compatibility conditions

(3){є}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(u_{i,j}+u_{j,i})={\nabla }_s\otimes \mathbf{u},

where E is the fourth-order elasticity tensor, ϵ is the strain tensor and u denotes the displacement field. For a scalar damage model, Eq. (2) is modified taking into account the degradation of the material

(4)\sigma =\underset{E_{\mathrm{e}\mathrm{f}\mathrm{f}}}{\underbrace{(\mathrm{1}-d)\mathbf{E}:є}},

where d is a damage as a function of the current state and history of the material, monotonically increasing from 0 to 1. Substituting Eq. (4) into Eq. (1), we obtain

(5)\mathbf{f}=\nabla ({\mathbf{E}}_{\mathrm{e}\mathrm{f}\mathrm{f}}{\nabla }_{\mathit{s}}\mathbf{u}).

Those equations can also be expressed in weak form by using e.g. the method of weighted residual. Integrating over the material domain Ω and choosing test functions δu in the subspace of kinematically admissible functions  and trial functions u in the function space , we obtain:

(6)\mathrm{f}={\displaystyle \underset{\mathrm{\Omega }}{\int }(\nabla ({\mathbf{E}}_{\mathrm{e}\mathrm{f}\mathrm{f}}{\nabla }_s\mathbf{u})\delta \mathbf{u})\mathrm{d}\mathrm{\Omega }\forall \{\mathrm{u}\in u,\delta \mathbf{u}\in \}}.

Integrating by parts, and choosing the test function in the same subspace as u, we obtain

(7)\mathbf{f}={\displaystyle \underset{\mathrm{\Omega }}{\int }({\nabla }^T\mathbf{u}({\nabla }_s{\mathbf{E}}_{\mathrm{e}\mathrm{f}\mathrm{f}})\delta \mathbf{u}+{\nabla }_s\mathbf{u}{\mathbf{E}}_{\mathrm{e}\mathrm{f}\mathrm{f}}{\nabla }_s^T\delta \mathbf{u})}\mathrm{d}\mathrm{\Omega }\forall \{\mathbf{u},\delta \mathbf{u}\in {\mathscr{H}}^{\mathrm{1}}\}.

Finally, if {\nabla }_s{\mathbf{E}}_{\mathrm{e}\mathrm{f}\mathrm{f}} is negligible over Ω, the usual weak form of linear elasticity is derived.

(8)\mathbf{f}={\displaystyle \underset{\mathrm{\Omega }}{\int }(\nabla \mathbf{u}{\mathbf{E}}_{\mathrm{e}\mathrm{f}\mathrm{f}}{\nabla }^T\delta \mathbf{u})\mathrm{d}\mathrm{\Omega }\forall \{\mathbf{u},\delta \mathbf{u}\in {\mathscr{H}}^{\mathrm{1}}\}},

with {\mathscr{H}}^{\mathrm{1}} a first-order Sobolev space.

The derivation of the discrete form is done by choosing a suitable space in which to approximate the displacement. Typically the Lagrange polynomial spaces {\ell }^i are used in finite elements because they simplify post-processing. A compact functional basis is associated to the degrees of freedom of the elements. If h_i is the i^th basis function, the discrete problem corresponding to Eq. (8) in element e is [10]:

(9){\displaystyle \underset{i,j}{\cup }\left({\displaystyle {\int }_e({\nabla }_s{h}_i){\mathbf{E}}_{\mathrm{e}\mathrm{f}\mathrm{f}}{({\nabla }_s{h}_j)}^T\mathrm{d}e}\right)\mathbf{u}=\mathbf{f}}

where \underset{i,j}{\cup } denotes the assembly of sub-matrices.

Non-local constitutive model

In non-local models, the constitutive model is expressed in terms of a non-local relationship instead of the local one:

(10)\sigma =\mathbf{E}(\widetilde{є}):є,

here \widetilde{є} is a non-local expression of є:

(11)\widetilde{є}={\displaystyle \frac{\mathrm{1}}{V}}{\displaystyle \underset{\omega }{\int }\alpha (\xi )є\mathrm{d}\omega },

with \alpha (\xi ) being a smoothing function and ω being the domain with \alpha (ϵ)\ne \mathrm{0}. The non-locality of the formulation guarantees that the damage pattern will not converge to band of single elements, but rather to a band of a fixed width which is governed by the decay property of the smoothing function [11]. The formulation of α in this case is exposed below in the outline of the algorithm. It takes the shape of a Heaviside function over a ’dynamically’ determined region.

Damage simulation algorithm

Energy minimisation

In the quasi-static case, dynamic effects such as inertia are neglected.

At equilibrium the elastic energy variation in the sample must be equal to the energy needed to transform the material.

(12)\mathrm{\Delta }{E}_{\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}}=\mathrm{\Delta }{E}_{\mathrm{i}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}.

This is equivalent to the irreversibility criterion not being met anywhere in the domain, and ΔE_elastic is maximum.

We assume loading steps during which the load is monotonously growing and such that initial load causes only elastic deformations. The material is tested elastically with the final load and locally reaches its irreversibility limit. An irreversible transformation of the material properties in this domain is applied. After the relaxation, global equilibrium has not been reached in general: for the specified boundary conditions, the material may need to undergo further irreversible transformation. However, the transformation is such that there exists a set of boundary conditions between the test conditions and the initial conditions such that the material is at equilibrium. Therefore, although the material is not at equilibrium during the iteration, it follows a state path composed of potential equilibria.

At any step of the process, ΔE_{irreversible transformation} is the sum of all energies from the successive irreversible transformations of the material. As the variation in elastic energy is automatically found from the load and the current material state, if an equilibrium is found, the error in the energy release from the irreversible transformations is at most equal to one step of relaxation. Thus it can be made to be arbitrarily small.

As damage only softens the material, the elastic energy stored can only diminish as the algorithm progresses. Thus the convergence to the solution is monotonic. The alternative is the minimisation of this energy using explicit energy functionals, which might not be experimentally measurable, or computationally expensive to calculate [12].

Algorithm formulation

The fracture criterion  is normally computed from a function of the element state s and compared to a critical value, c:

(13)\{\begin{array}{ll}\mathrm{1}& \mathrm{i}\mathrm{f}(\mathbf{s})>c\\ \mathrm{0}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}.

This function is smooth and monotonously increasing with s. It is commonly normalized R → [0, 1] such that

(14)\nu =\{\begin{array}{ll}\mathrm{1}-\left|{\displaystyle \frac{c}{(\mathbf{s})}}\right|& i\mathrm{f}(\mathbf{s})>c\\ \mathrm{0}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}.

The strain limit criterion is a simple example. If the strain є is larger than є_c, the material fails. The normalized criterion reads \mathrm{1}-{\displaystyle \frac{{є}_c}{є}}. The update of the elements can be done using a non-local damage model, in which case all elements in a ball centered on the element with the maximum value will be updated [13,14]. The algorithm is given by:

Algorithm

For each element, define a neighborhood in which to compare the values of the criterion. h is the characteristic distance defined in the constitutive law.

1. Compute the elastic deformation from the load.

2. For each element, compute the criterion using Eq. (14) as a real normalized value, for example, if the fracture criterion is a maximum strain є_c, compute:

(15)=\mathrm{1}-\left|{\displaystyle \frac{{є}_c}{є}}\right|.

3. If the elements within the neighborhood which have exceeded their criterion have a total surface, respectively volume, superior to that of the disc, respectively the ball, defined by h: Find all elements such that:

• the element is at a distance less than h from the locus of the maximum calculated criterion at this step in this neighborhood

• the element exceeds its criterion

4. Increase the damage in all found elements

5. Iterate form the new material state.

Treating elements in neighborhoods is done for the sake of numerical efficacy. It is reasonable to assume that when damage occurs throughout the material small, local damage increments will not affect the solution, especially if the material considered is geometrically complex. Indeed, composites under load have very complex stress patterns. The stress patterns are locally influenced by the inclusions and pores. This means that the influence of the weakening of an element is negligible beyond a distance in the order of the median diameter of the inclusions.

Condition for convergence

Solving for damage can be done using different numerical methods. The Newton family of solvers use a fixed point strategy: the material state is updated according to a stress-strain relationship at each step of the iteration until no update of the material state is required. In arc-length solvers, the boundary condition is updated together with the material state to find the equilibrium.

The method proposed here includes the solver in the expression of the algorithm: it functions like a Newton method where the descent direction is imposed. Imposing the direction ensures that the algorithm converges monotonically. Newton methods converge for small enough increments in the load, but here, the descent radius is a user specified parameter, in the form of the defined damage increment. The stability condition of the algorithm is derived below. Because the descent direction is chosen by the algorithm, it is necessary to verify that this does not affect the convergence conditions of Newton methods.

The algorithm is always stable if the second derivative of the compliance is positive (see appendix). This condition is equivalent to requiring that the compliance be strictly convex, which is the usual requirement for the convergence of a Newton method. Therefore, choosing the descent direction in the way proposed in this paper does not adversely affect the convergence of the Newton solver.

Convergence of the error

Due to the monotonic nature of the convergence, the error is easily derived. The damage is incremented in a set of elements such that at convergence, none of them is beyond its yield surface. At no point during the iteration was damage incremented in an element which did not reach its criterion. Therefore, an upper bound for the error in terms of the damage state of the elements is exactly one damage increment per element damaged during the iteration.

(16)

The area damaged {}_{\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{d}} depends on the geometric setup, but also on the load step, and can therefore be made arbitrarily small. Further, the damage increment δd is chosen freely and can also be made arbitrarily small. Of course, choosing a small damage increment yields a penalty in terms of computational efficiency.

Comparison with a classic non-local approach

To compare the results of the algorithm proposed in this paper to classical non-local methods coupled to a Newton solver, we have implemented a test setup with a square sample with two notches under direct tension. The damage was calculated using both methods. The same constitutive law was used in both cases, elastic perfectly plastic, with the same characteristic length set to 10. The setup is illustrated on Fig. 4.

The non-local smoothing function used in the classical non-local methods was a Gaussian:

(17)\overline{\sigma }={\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\pi }}{\int }}{\displaystyle \underset{\mathrm{0}}{\overset{\mathrm{2}}{\int }}\ell {\sigma }_{\mathrm{V}\mathrm{o}\mathrm{n}\mathrm{M}\mathrm{i}\mathrm{s}\mathrm{e}\mathrm{s}}(r)e^{-\frac{r^{\mathrm{2}}}{{\ell }_{\mathrm{m}\mathrm{a}\mathrm{t}}}}\mathrm{d}r\mathrm{d}\theta }.}

The solution was obtained using a modified (fixed point) Newton method. The resulting stress-strain curve is plotted on Fig. 5. A comparison of the damage patterns can be seen on Fig. 6.

The stress-strain responses obtained from the two methods are identical. This is despite the implicit smoothing functions being different. Further, at each step of the iteration of the proposed algorithm, only a subset of elements have their mechanical properties recalculated whereas at Newton step of the classical method, all elements have their properties recalculated.

The damage patterns are nearly identical, but the pattern obtained with the classical approach is slightly less sharp. This can be explained by the smoothing function used, which spreads damage more uniformly than the proposed algorithm. The smoothing function is arbitrarily defined, and imposes a shape for the spread of the damage. To maintain isotropy, this shape is necessarily circular. On the contrary, the proposed algorithm follows the flow of the damage, and the shape is defined by the fracture criterion. This does not induce anisotropic behavior, because it follows the contour curves defined by the deformation of the material at each step.

The number of iterations to convergence as a function of the loading steps are plotted on Fig. 7. The increment is a displacement of 5 × 10^-3. The number of steps required for the algorithm described in this paper depends on the damage increments required for convergence. Therefore, if few elements need to be damaged, the algorithm compares favorably to a Newton method. Further, the update is only done on a subset of the elements, which can be advantageous if the matrix assembly is costly. Unlike classical methods, here the number of iterations is proportional to the amount of damage required to reach equilibrium. If a loading step causes a lot of damage, the number of iterations is large. The spike is caused by the formation of a large damaged zone.

The shape of the damage pattern can be affected by changing the fracture criterion, while keeping the stress-strain relationship constant. This effect is described in the next section.

Choice of fracture criterion

Setup

Fracture criteria relate energy dissipation to stress or strain configurations. They relate to microstructural material properties, and prescribe the preferred mode of fracture of materials. The same numerical experiments were run with different mesh refinements and fracture criteria. The sample was discretised with our in-house Delaunay mesher at four different mesh densities [13]. The goal of these experiments was to observe the effect of the choice of fracture criterion on the damage pattern, for materials with otherwise identical constitutive laws (elastic-perfectly plastic). Because the criteria differ, the critical equivalent stress values are not the same in each case.

The geometrical setup is per Fig. 8. This setup, with additional initial cracks is usually used as a benchmark for partition of unity (PU) crack propagation. It is interesting because of the inbuilt symmetry and mixed-mode failure it induces. Also, this experiment has a very high phase contrast between the phases (pores and matrix).

The material is linear elastic, with parameters E = 147.25 and ν = 0.3.

The stiffness tensor E is:

(18)\mathbf{E}={\displaystyle \frac{E}{\mathrm{1}-{\nu }^{\mathrm{2}}}}\left(\begin{array}{ccc}\mathrm{1}& {\displaystyle \frac{\mathrm{1}}{\nu }}& \mathrm{0}\\ {\displaystyle \frac{\mathrm{1}}{\nu }}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}-\nu }{\mathrm{2}}}\end{array}\right).

To demonstrate the convergence properties of the algorithm, the following modified von Mises criterion was used (σ_i is the i^th principal stress):

(19)(\mathbf{s})=\mathrm{1}-{\displaystyle \frac{{\sigma }_c\sqrt{\mathrm{2}}}{\sqrt{{({\sigma }_{\mathrm{0}}-{\sigma }_{\mathrm{1}})}^{\mathrm{2}}+{\sigma }_{\mathrm{0}}^{\mathrm{2}}+{\sigma }_{\mathrm{1}}^{\mathrm{2}}}}}.

A material which can only fail along shear planes was also tested:

(20)(\mathbf{s})=\mathrm{1}-{\displaystyle \frac{{\sigma }_c\sqrt{\mathrm{2}}}{\sqrt{{({\sigma }_{\mathrm{0}}-{\sigma }_{\mathrm{1}})}^{\mathrm{2}}}}}.

And finally a Mohr-Coulomb criterion:

(21)(\mathbf{s})=\mathrm{1}-{\displaystyle \frac{{\sigma }_c}{\max \left|{\sigma }_i\right|}},

with the critical value σ_c = 23 in all cases.

A displacement is imposed on the upper and lower boundaries, such that complete failure of the specimen does not occur. Here the imposed displacements are 1 and -1 and the computation is done in one step.

Results

Figure 9 shows the results from a set of simulations. Each row shows a different yield criteria (von Mises, Mohr-Coulomb, pure shear failure), within the rows the mesh refinement increases from left to right. The damage pattern is strongly affected by the choice of failure criterion, but not by the successive mesh refinements. This figure illustrates the ability of the algorithm to spread damage as non-local smoothed approaches without the need to specify explicitly a distribution function for the damage. The shape of the area damaged is determined from the choice of fracture criterion, and not from an a priori guess.

The algorithm avoids crack coalescence on a single band of elements, and does not exhibit pathological mesh dependency. Most importantly, the total energy dissipated during deformation converges with mesh refinement. The time to convergence was experimentally found to be proportional to (n^{\mathrm{1.72}}) where n is the number of unknowns. The time taken to find which element to damage is small compared to the time taken to solve the linear system of equations, thus the performance in terms of computing time of the algorithm depends on the choice of solver. When very coarse meshes are used, as the damage is constant per element, the connectivity of the damaged zone can be wrong. This is due to insufficient fineness of the mesh.

Damaging the elements leads to a decrease in the elastic energy stored in the mesh compared to the undamaged elastic test solution. The loss can be measured in absolute terms or as a fraction relative to the original value. The elastic energy at convergence for the same meshes was studied. Both the relative and absolute energy losses converge as the number of elements is increased. There is a distinction between absolute and relative energy, because the absolute energy at the beginning of the iteration depends on the discretisation of the problem: the pores are not satisfactorily represented when the mesh is coarse.

Practical application: Micro-mechanical model of ASR

ASR is an important durability issue in concrete structures [15]. It is characterized by the reaction of the aggregates with the alkalies present in the cement pore solution. The product of the reaction is a hydrophilic gel which swells, causing the aggregates to crack. Cracking then extends to the paste. Ben Haha carried out experiments where the macroscopic expansion and the mechanical properties were measured and reported against measures from image analysis of the damage state of the aggregates [16]. The observations of Ben Haha were used as a basis to formulate a numerical model of asr. The mesostructure of concrete was modeled explicitly, as well as the gel pockets formed by the reaction, see [17] for a complete description of the setup. These simulations are meant to reproduce the real experimental conditions as closely as possible. Due to the large size of the problem, linear elements were used. Figure 11 compares the simulated and observed damage patterns.

Computed damage, swelling and loss of mechanical properties were com pared to the experimental results and were found to be in very good agreement. As the mechanical properties of the gel are unknown, its stiffness was fit as a fraction of that of C-S-H, the main hydrate in cement paste, whose chemical composition is close to the reported composition of asr gel [18]. This is the parameter α in Fig. 12.

Full-scale simulations using the real particle size distribution of aggregates and realistic placement of asr gel pockets were run. The elastic properties of the material were measured, as well as their critical stresses. The model and the experiments were directly compared. Damage modeling required an algorithm which provides detailed damage patterns. This was possible using the present algorithm, because it uses all the available precision from the mesh to produce damage paths.

In Fig. 12 damage is measured as the relative volume of all cracks in the aggregates to the total aggregate volume, in both the experiments and the simulation. It is termed “apparent reaction” to distinguish it from the “reaction” which is the relative volume of aggregate effectively transformed into gel by asr. Dunant and colleagues show that there is a strong correlation between the two in [17].

The damage algorithm produced realistic cracking patterns both in the aggregates and the paste, despite the strength and stiffness differences between those two phases. The experimentally observed relation between strains with respect to the damage state, which is due to the complex interactions between growing gel pockets partly healing the cracks, and the propagation and merging of cracks, was reproduced by the simulation.

The location of the damage, as well as the density of branching, etc., induced a loss of mechanical properties in he simulation comparable to the experiment. The initial softening of the samples is due to aggregate degradation. Further softening is caused by the cracks propagating in the cementitious matrix. This process, predicted by the model, was corroborated by microscopic observations.

Direct simulation of such damage patterns is possible using the proposed algorithm, because it converges to finely defined damage paths, and is applicable to composite problems.

Discussion

The algorithm outlined in this paper was developed to solve composite problems with explicit geometry. It captures the shape of damage domains without the need to specify a smoothing function. The material law is expressed as the fracture criterion as a function of the damage state. Normalization of the criterion allows the computation of mixed-material damage evolution.

The ASR example serves as a large-scale verification of the validity of the approach, by demonstrating the feasibility of such complex setups. Also, such simulations, of great interest to the community working on asr required this approach approach, due to the particular properties of the problem (changing microstructure, high crack density, complex geometrical layout).

The algorithm described hinges on the description of the material behavior as irreversibility occurs. Localization of the damage, in the sense that micro-defects coalesce into a macro-crack, is implicitly described by the damage state update: a cohesive law can be introduced in the material behavior there. This algorithm can thus also be used to make the transition between continuum and discontinuous techniques such as the XFEM: the introduction of a localized discontinuity in an element can be considered as a further damage state.

The damage increments can be adjusted to each problem to obtain the best performance. As a direct simulation, this algorithm cannot be guaranteed to be absolutely stable. If the increments are too large, the softening of the material at the local level will not be adequately captured. However, if the increments are too small, convergence will be very long to obtain. We experienced no unstable behavior in our experiments, however.

Crucially, this algorithm converges with mesh refinement. This is necessary when definite damage paths are sought, for example in cases where damage is both present throughout the material and highly localized. The ASR example showed how such properties are required for the validation of a microstructure-level degradation mechanism.

The boundary conditions are assumed to be varying monotonously during a loading step. Thus, this algorithm is ill-suited to capture snap-back: at each step of the simulation, an equilibrium will be reached for the set boundary conditions if it exists. An extra loop in the algorithm would be required to keep the load at equilibrium. However, the introduction of this loop is made possible by the monotonic convergence property of the algorithm.

The implementation of the algorithm is well-suited for a space-time imple mentation: the damage can be expressed as a continuous function over the time slice and can be increased in a manner similar to the quasi-static case. This opens the possibility of computing damage in visco-elastic materials.

Conclusion

A direct algorithm for the computation of non-local damage in bundle of fiber materials has been presented. This algorithm can be applied to cases where load depends on the varying geometry of the sample, and can be used in composite cases. The convergence of the algorithm in terms of energy and domain definition was studied. The convergence was shown to be monotonic. Applications to simple and complex cases were presented, highlighting the possibilities of the algorithm to study durability issues, notably in cementitious composites.

Future work involves the extension of the implementation to XFEM crack initiation and propagation as well as the implementation of more complex damage models, using the same element-ranking approach to compute the locus of the crack initiation. This algorithm should also be implemented in space-time and tested in visco-elastic cases, using damage functions to form materials which are functionally graded in time.

The framework used in this paper, Amie, is also implemented in 3D and is available as Open Source. Interested readers should contact the authors.

References

[1] Bazant Z P, Belytschko T. Wave propagation in a strain-softening bar: Exact solution. Journal of Engineering Mechanics, 1985, 111(3): 381-389

[2] Bazant Z, Jirasek M. Non-local integral formulations of plasticity and damage: survey of process. Journal of Engineering Mechanics, 2002, 128(1): 1119-1149

[3] Babuška I, Melenk I. Partition of unity method. International Journal for Numerical Methods in Engineering, 1997, 40(4): 727-758

[4] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601-620

[5] Strouboulis T, Babuška I, Copps K. The design and analysis of the generalized finite element method.Computer Methods in Applied Mechanics and Engineering, 2000, 181(1-3): 43-69

[6] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous crack initiation, nucleation and propagation in statics and dynamics. Computational Mechanics, 2007, 40(3): 473-495

[7] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. International Journal of Solids and Structures, 2008, 75(5): 943-960

[8] Bordas S, Duot M, Le P. A simple a posteriori error estimator for the extended finite element method. Communications in Numerical Methods in Engineering, 2007 (In press)

[9] Stratonovitch R L. Derivation of irreversibility of thermodynamic processes from microscopic reversibility. Theoretical and Mathematical Physics, 1978, 36(1): 607-616

[10] Simo J C, Hughes T J R. Computational Inelasticity,vol. 7. Springer, 1998.

[11] Bažant Z P, Jirásek M. Nonlocal integral formulations of plasticity and damage: Survey of progress. Journal of Engineering Mechanics, 2002, 128(11): 1119

[12] Francfort G A, Marigo J J. Revisiting brittle fracture as an energy minimisation problem. Journal of the Mechanics and Physics of Solids, 1998, 46(8): 1319-1342

[13] Dunant C, Vinh P N, Belgasmia M, Bordas S, Guidoum A. Architecture tradeoffs of integrating a mesh generator to partition of unity enriched object-oriented finite element software. Revue Européenne de Mécanique Numérique, 2007, 16(2): 237-258

[14] Bordas S, Nguyen V P, Dunant C, Nguyen-Dang H, Guidoum A. An extended finite element library. International Journal for Numerical Methods in Engineering, 2007, 71(6): 703-732

[15] Ponce J M, Batic O R. Different manifestations of the alkali-silica reaction in concrete according to the reaction kinetics of the reactive aggregate. Cement and Concrete Research, 2006, 36(6): 1148-1156

[16] Ben Haha M, Gallucci E, Guidoum A, Scrivener K L. Relation of expansion due to alkali silica reaction to the degree of reaction measured by sem image analysis. Cement and Concrete Research, 2007, 37(8): 1206-1214

[17] Dunant C F, Scrivener K L. Micro-mechanical modelling of alkali silica-reaction-induced degradation using the amie framework. Cement and Concrete Research, 2010, 4(4): 517-525

[18] Kawamura M, Iwahori K. Asr gel composition and expansive pressure in mortars under restraint. Cement and Concrete Composites, 2004, 26(1): 47-56

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg