A computational model of interaction between material and interface cracks

A quasi-static model for numerical solution of initiation and propagation of cracks along interfaces or inside materials is developed. The two types of cracks are modelled by the material damage theory with two independent damage parameters introduced. For cracks at the interface, in fact represented by contact of construction components, cohesive or adhesive contact is considered, for which several computational relationships based on energetic formulation exist. Accordingly, the appropriate modelling of bulk damage also includes energy consideration. In terms of cracks, it leads to so called diffuse cracks. The computational approach is referred to as phase field models. These will cause damage in a very narrow band representing the actual crack. The computational analysis provides stress-strain quantities and the damage variables to simulate both interface and material cracks. The proposed mathematical approach has a variational form based on an energetic formulation looking for a kind of weak solution. The solution is approximated by a time stepping procedure, a finite element code, and it utilizes quadratic programming algorithms.


Introduction
Cracks may appear in any part of a civil engineering construction. Such cracks may significantly change mechanical or other properties of the structure. Thus, developing new efficient algorithms for analysing and predicting is highly demanding. The cracks may be of various nature and basically can appear either inside the materials or along material interfaces. If a material crack forms first, it may hit an interface and later propagate along it, and vice versa, a primary interface crack may kink into the material. The present proposed approach enables to study both cases.
The present approach models cracks in terms of damage mechanics where the state of the structure is described by an internal variable [1] which reflects the damage state, e.g. density of microcracks or microvoids, and its limit value means a new crack formation.
As other nonlinear phenomena also damage and fracture can be treated variationally. Previous works of author [2,3] provided several such possibilities of solving problems with interface cracks in a quasi-static manner. To include material cracks into them, it was desirable to keep the same energetic philosophy of models' generalisation [4]. Therefore, a damage phase-field model [5][6][7] was considered to simulate such cracks. These cracks are also called diffuse as they are related to a scalar material damage variable which should by at least a continuous function so that the crack seams to be blurred, nevertheless the phasefield models can control the width of such a haze.
Numerical solution of these evolution models then, besides standard time stepping and finite element discretisation, includes optimisation algorithms, too. The simplest case then considers convex even quadratic functionals for whose minimisation there exist quadratic programming algorithms [8], which can be straightforwardly adapted for the general convex case [9]. Therefore, we first describe the model in terms of energies and governing relations for their evolution, and then we show some details for the model discretisation. Finally, we document the behaviour of the proposed approach by two simple examples.

Description of the model
We consider a domain Ω which consists of two or more bonded parts. An example of two such parts Ω and Ω is shown in Fig. 1, where their respective boundaries are denoted Γ and Γ . The common part of the boundaries -an interface -is Γ C . Each of the boundaries is additionally split into disjoint parts according to boundary conditions: Γ D where displacements are prescribed by known function , Γ N where tractions are given by . The behaviour of the model is governed by energy state evolution. This state is described by displacement field , a gap of interface displacement ⟦ ⟧, and two internal parameters and varying between 1 and 0: the former determines the damage state from the phase-field damage model in the bulk domains so that = 1 corresponds to undamaged state and = 0 pertains to a crack. Similarly, the latter parameter characterises the state of the interface, considered as a thin adhesive layer, with = 1 corresponding to the undamaged state and = 0 pertaining to an interface crack. The energy includes the stored energy which may be considered as follows: for an admissible displacement field satisfying the displacement boundary conditions on Γ D and admissible damage parameters and , lying in the interval [0,1], otherwise the value of would be infinite. The introduced parameters include stiffness matrix of the material of the domains, stiffness of the adhesive layer and the compressive stiffness g , to replace contact condition by a penalisation term. Further, c is fracture energy for domain cracks, c I is fracture energy of the interface cracks, is a characteristic length parameter for the phasefield model (determining the width of the blurred crack) and the parameter 0 introduces residual stiffness to avoid degeneration of the totally damaged material in the numerical for accumulating energy of the crack in Griffith-like models is replaced by a generalised Abrosio-Tortorelli [10] functional (the term containing c ) of the phase-field damage model diffuse crack. The interface term includes elastic energy of the damageable interface in a simplified way of the approaches from [2,3], in order to be of the same form as for the domain cracks. The energy of the crack formation process is usually dissipated from the system. While in the present model we included it into the stored energy, it is still necessary to guarantee a unidirectional character of the damage process at least by the assumption that there is no additional dissipated energy expressed as a potential (̇,) = 0, provided that ̇≤ 0, ̇≤ 0, otherwise it would be set to infinity.
Finally, if there are some external forces, their energy should be added to the system, too. In Fig. 1 we considered forces applied at a part of boundary Γ N , their energy includes the functional

Discretisation
Numerical solution requires both a time stepping algorithm and a spatial discretisation. The latter is implemented by a standard finite element approach using triangular elements [11] and will not be discussed. The time discretisation by a semi-implicit fractional-step method relies on separate quadraticity of the proposed functional Eq. (1) in relation to each state variable which then guarantees a variational character of the solved problem.
We may choose a fixed time step for a time range [0, ] such that the solution is obtained at the instants = for = 1,2, … / and denoted for displacements and , for the damage variables. In order to obtain such an algorithm from Eq.

(4)
The separation of variables thus provided two minimisations which should be performed at each time step: the first minimisation with respect to displacements of ( ) = ( ; , −1 , −1 ) + ( ) (5) renders as its minimiser (the constraint for u is hidden in definition of ), and the second minimisation with constraints ≤ −1 and ≤ −1 of ( , ) = ( ; , , ) provides as a constrained minimiser. These two minimisations are solved recursively for = 1,2, … / . In fact, both functionals are quadratic so that various quadratic programming algorithms can be used in the numerical solution, see e.g. [8].

Examples
The proposed model will be tested in two problems. The first one is to demonstrate a comparison with a standardly used approach for phase-field damage from [5] and it is used as a benchmark but does not include an interface. Therefore, the other example demonstrates mutual interaction between an interface and a domain crack. The first example's geometry is shown in Fig. 2. It shows an asymmetric beam with an initial crack and three circular cavities. The parameters are taken from [5], it means elastic properties are = 20.8 GPa, ν = 0.3 and the fracture energy c = 1 kJm -2 . In relation to the mesh size, which is ℎ = 0.1 mm, the length parameter = 0.25 mm. There is an initial crack of the length = 1 mm. The displacement loading ( ) is linearly increasing by the velocity = 1 m ms -1 , where the time steps are chosen 10 ms for < 0.3 s and 0.1 ms for > 0.3 s. The crack propagation is shown in Fig. 3. Its form is the same as in the mentioned reference, which we intended to see. Of course, if the cavities were filled by a material, the crack could continue along the interface. Fig. 2. Geometry of the example according to [5] and the used mesh.
The other example is shown in Fig. 4. We consider three mutually bonded blocks, where only one of the interfaces, plotted thicker in the picture, is damageable and thus allowing cracking. The material parameters are set the same as in the previous case, only c = 0.1 kJm -2 , while that of the vertical interface is c I = 1 Jm -2 . The cracking along other interfaces is forbidden by setting the fracture energy very high. The initial stiffness of the damageable interface is 2 TPam -1 and its degradation function is quadratic, to be able to use only quadratic programming algorithms. Discretisation characteristics are: mesh size (min.) ℎ = 1 mm, time step 0.1 ms.

Fig. 3.
Crack propagation for the first example. Deformation is magnified two times and the snapshots correspond to the time instants = 0.3s, 0.3059s, 0.3078s, where the first instant corresponds to the moment of the damage triggering.
The pictures in Fig. 5 reveal the character of the fracture process. First, there appears an interface crack, which of course stops propagating when it reaches the bottom block. This crack does not immediately continue to this block as the energy release rate has not reached the fracture energy for the domain crack.
This continuing crack appears later and is spread in the vertical direction until it approaches the fixed bottom side of the bottom block when it kinks to the more or less horizontal direction. .

Conclusion
A computational model for problems, where appear both interface and domain cracks, has been introduced. The model is based on a damage theory and includes internal parameters controlling either the interface cracks as a corrupted thin adhesive layer or domain cracks as a diffuse crack of a phase-field damage model.
From the computational point of view, the proposed problem is solved by a semi-implicit time-stepping method, guaranteeing a variational structure for a numerical approximation of the model and then quadratic programming algorithms with spatial discretisation by finite elements.
Of course, there are some parameters in the model which should be suitably tuned to follow experimental observations. Nevertheless, it might be expected that the proposed approach turns to be useful in practical engineering calculations.
The author acknowledges support from the Ministry of Education of the Slovak Republic by the grant VEGA 1/0199/20 and VEGA 1/0374/19.