AN ENERGY BASED FORMULATION OF A QUASI-STATIC INTERFACE DAMAGE MODEL WITH A MULTILINEAR COHESIVE LAW

. A new quasi-static and energy based formulation of an interface damage model which provides interface traction-relative displacement laws like in traditional trilinear (with bilinear softening) or generally multilinear cohesive zone models frequently used by engineers is presented. This cohesive type response of the interface may represent the behaviour of a thin adhesive layer. The level of interface adhesion or damage is deﬁned by several scalar variables suitably deﬁned to obtain the required traction-relative displacement laws. The weak solution of the problem is sought numerically by a semi-implicit time-stepping procedure which uses recursive double minimization in displacements and damage variables separately. The symmetric Galerkin boundary-element method is applied for the spatial discretization. Sequential quadratic programming is implemented to resolve each partial minimization in the recursive scheme applied to compute the time-space discretized solutions. Sample 2D numerical examples demonstrate applicability of the proposed model.


1.
Introduction. Computational simulations of crack initiation and propagation are currently under intensive development. It is known that the traditional linearelastic fracture mechanics (LEFM), which assumes a specific square root stress singularity at the crack tip, is in general unable to predict crack initiation at stress concentrations, weak stress-singularities or in regions of uniform stress fields, because it assumes a pre-existing crack. Nevertheless, there exist several fracture models which avoid stress singularities, like cohesive zone models (CZM) [3,5,6,24,26,36] or variationally based interface damage and plasticity models [9,15,23,29,32,43]. These models allow the engineers to simulate initiation as well as propagation of cracks. Following [32], a new quasi-static and energy based formulation was recently developed in [40], where a scalar damage variable for describing damage evolution and crack initiation and propagation was suitably defined so that it provides the same interface traction-relative displacement laws as in several well-known CZMs.
(SQP) [28], depending on the actual form of the arising functionals, as it was also done in [40,43,44]. 2. A cohesive zone model with fracture mode sensitivity.

2.1.
Description of the interface model. First, we describe a model which includes a possibility of having multilevel damage for cohesive crack evolution and fracture mode sensitivity. Only two 2D bodies will be considered for the sake of simplicity of explanation and notation, though the theory can be generalized to consider several 2D or 3D bodies. The bodies are supposed to occupy domains Ω η ⊂R 2 (η=A, B) with bounded Lipschitz boundaries ∂Ω η =Γ η , Figure 1. Let n η Figure 1 The used notation for two bonded domains.
denote the unit outward normal vector defined at smooth parts of Γ η , let t η denote the unit tangential vector such that it defines anti-clockwise orientation of Γ η .
The cohesive contact zone Γ C is defined as the common part of Γ A and Γ B , i.e. Γ C =Γ A ∩Γ B . The Dirichlet and Neumann boundary conditions are defined on the outer boundary parts, respectively, prescribing displacements as u η =g η D (t) on Γ η D and tractions as p η =f η N (t) on Γ η N at a time instant t, cf. (2b),(2c) below. It is naturally assumed that Γ η D and Γ η N are disjoint parts of Γ η , disjoint also with Γ C , and Γ η =Γ D η ∪Γ N η ∪Γ C . Additionally, the Dirichlet part is considered to lie far from the contact boundary, i.e. Γ D η ∩Γ C =∅.
The difference across Γ C of the functions defined in Ω A and Ω B will be denoted by [[·]]. In particular, the jump (difference) of displacements on the contact boundary Γ C means [[u]]=u A | Γ C −u B | Γ C . Further, we will use the jump of the normal and tangential displacements u n = u · n B = − u · n A = −u B | Γ C · n B − u A | Γ C · n A , u t = u − u n n B (1) We also use the convention that the 'dot' stands for the partial time derivative ∂ ∂t .
where C η and D η are the positively definite symmetric 4th-order tensors of elastic and viscous moduli, respectively, in the domain Ω η . The multilevel damage CZM is represented by the traction-relative displacement relations (2d) and the complementarity conditions for interface damage parameters ζ i (2e). The interface is considered as an infinitesimally thin layer of an adhesive, whose elastic moduli are introduced by the matrix K. Then, K is the positive-definite 2nd-order tensor of interface stiffness, K : 2×2 , depending on the actual state of damage expressed by the value of the multilevel damage parameter ζ=(ζ 1 , ζ 2 , . . . ζ m d ) with given number of damage levels m d : ζ=1=(1 1 , 1 2 , . . . 1 m d ) corresponds to an undamaged state, and ζ=0=(0 1 , 0 2 , . . . 0 m d ) represents total rupture. Each component of the matrix function K is defined by a function, increasing in each ζ i , with K(1) being the initial (undamaged) stiffness of the adhesive. Also, it is supposed that (K(ζ)v) ·v is a convex function of ζ for any v. It should be noted that in Γ C the problem is usually expressed in local coordinate system defined by normal and tangential components as in (2d). K is expressed in this coordinate system, too, and it is frequently simplified to a diagonal matrix. The diagonal terms of such K are convex functions.
The parameter k G is the compression stiffness of the normal-compliance contact allowing interpenetration of the bodies (we introduced the notation x − = min(0, x) for any real x), but for large values of k G being rather realistic replacement of the Signorini contact which can be efficiently implemented by QP.
The complementarity conditions in (2e) govern the evolution of the multilevel damage parameter ζ. The Lagrange parameter λ was introduced to guarantee the lower bound of ζ.
The interface fracture energy G c , required to break a unit length (or area in 3D) of the interface, is considered to depend on the current state of displacement jump [[u]] through the displacement fracture-mode-mixity angle ψ defined as In the present multilevel damage case it is split into m d parts so that For the fracture energy, any phenomenological law describing its dependence on the mode-mixity angle can be used, e.g. that of Hutchinson and Suo [14], for each damage level (2h) The parameters G iIc and G iIIc express the partial fracture energies at the i-th damage level in the pure Mode I and in the pure Mode II, respectively.
The normal and tangential components of the traction defined in (2c), are denoted by σ n and σ t in (2d), respectively, as interface stresses. The evolution BVP (2) holds for the instants t∈[0, T ] within a fixed time range T . We will further consider an initial BVP by prescribing the initial conditions at time t=0 for the displacement and the damage parameter: It should be noted that initiating the damage at its maximal allowed value does not require to define the condition ζ≤1 in (2e) as it follows by the condition . ζ≤0. It is noteworthy that the evolution BVP (2), when written in a weak formulation, takes the form of a nonlinear evolution governed by a stored energy functional E and a nonsmooth potential of dissipative forces R. In the Biot-equation form, it reads as a system of the differential inclusions for u=(u A ,u B ) and g D =(g A D ,g B D ), where the symbol ∂ refers to partial subdifferentials relying on the convexity of R(u; ·, . ζ), R(u; . u, ·), E (g D ; u, ·) and E (g D ; ·, ζ). It involves the functionals E , R in the specific form: and a linear functional F (t) acting on u as Note that in (4b) we denoted the split of the functional R according to the degree of homogeneity of its terms with respect to the rates of the arguments (the subscript index). The dissipation rate then reflects the different homogeneity degrees as occurs in (7) below. As follows from (3b) in view of (4c), the energy released due to progressing of interface damage (∂ ζ E ) equals the interface fracture energy (∂ . ζ R 1 ). This is a distinguished feature of the present formulation in comparison with others previously proposed potential based CZMs.
If we want to write the energy balance for the system (3), we have to use a shift of u to have homogeneous boundary conditions suitable for testing (3). Let u D define an appropriately smooth extension of the boundary condition g D . Since we assume that the interface Γ C is far from the Dirichlet boundary Γ D , we can restrict to [[u D ]]=0 on Γ C so that such a shift does not affect (4c). Then, ] on Γ C , and satisfies δ .
Now, testing (5a) by  ζ ) is allowed. The sum of such tests provides an energy inequality which, using some technical details derived in [31,39], can be formulated in the present case as the following equality: Let us also note that some integrals are written rather formally as e.g.
. ζ does not have to be sufficiently smooth. Integrating it over the time interval [0; T ] with the initial conditions (2i) renders the following energy balance, by substituting u S according to its definition, which expresses the energy conservation in terms of the original displacement u and damage ζ: the sum of the energy stored in the system at the time T and the energy dissipated during the time span [0; T ] due to interface crack growth and due to the bulk viscosity equals to the energy stored at the system at the time 0 plus the work made by the external boundary loadings g η D and f η N during the same time interval.

2.2.
An engineering insight to the multilinear cohesive model. The most simple case of the introduced damage model is obtained by representing the interface by only one 'spring'. This simplified model can be useful for analysing the condition which triggers the damage process in the 'spring' under consideration. The energy functionals (with no bulk integrals, no viscosity, no compression) from (4) reduce where The function K is a convex function with K(0)=0 and K(1)=k, k being the stiffness of the undamaged spring. Notice that K(ζ) is reduced to a scalar function in this simple case. The activation criterion for initiation of the i-th level of the interface damage can be derived from the condition (3b) or from (2e). The pertinent complementarity conditions read where we suppose u≥0, representing elongation of the spring. Also the constraint ζ≥0 should be taken into account, while the upper bound ζ≤1 is satisfied naturally if the initial condition satisfies it and .
ζ≤0. The situation can be simplified by considering K separated to m d functions all Φ i being convex and increasing. Then, e.g. (9) is also split into m d uncoupled complementarity problems It means that reaching the damage activation point for ζ i and continuing increasing the load, ζ i starts to decrease in time and the inequality Eq. (11) 1 holds as equality which provides the relation for ζ i . The spring stress σ is the derivative of E in (8) with respect to u, we thus obtain the stress-displacement curve as: for inequality in (11) 1 prior to the damage initiation at i-th level, or inequality in (11) 1 after reaching lower bound for ζ i , respectively. Here, also (−1) denotes the inverse of the pertinent function. As we require Φ i to be a convex increasing function, the inverse exists.
As we intend to obtain a trilinear cohesive model, or generally a multilinear CZM, we start with a bilinear one and then we separate the multilinear CZM in a series of bilinear ones. As it was introduced in [32,45] where the parameter β determines the slope of the softening part of CZM, very large β can be used to approximate an elastic-brittle interface model [20,35]. The stress-displacement relation obtained by differentiation, as recalled in (12), is then governed by the following relation: where the lower bound u 0 is the damage activation point: for all smaller values the relation is linear σ=ku as ζ=1, and where the upper bound u c is associated to the total damage with ζ=0: for all greater values the spring stress σ vanishes. It also should be noted that the maximal allowed stress is expected to appear for u 0 and due to (14) equals σ max = 2Gc kβ 1+β . Also the three parameters required by the model can be easily obtained from the stress-displacement relation as The stress-displacement curve is plotted in Figure 2(a). The graph also reflects the fact that upon unloading the process is irreversible. If now the stress-displacement relation is given by a broken line as shown in Figure 2(b) for positive u, it can be seen as a sum of the functions of the shape from Figure 2(a) (for positive u), so that the stiffness K is given by the relation (10) with all Φ i given by (13): The model then requires knowledge of all G ic , k i and β i . Two basic ways how such separation can provide these parameters are shown in Figure 3: either the damage in the variable ζ i is triggered just at the moment when ζ i−1 reaches zero (u 0i =u ci−1 ), as in Figure 3(a) named Mode 'subsequent', or all ζ i are triggered at once, but they have different u ci , as in Figure 3(b), named Mode 'at once'. The relations are clear from geometrical description in the pictures and from the relation (15). In the case of Mode 'subsequent' we obtain and for Mode 'at once' we have where, if necessary, we consider u m d =u c , u m d +1 =u c +1 and σ m d =σ m d +1 =0. The Mode 'subsequent' was previously discussed in [13,38] whereas the Mode 'at once' in [8], see also [12] for a different partition of fracture energy in a particular case of bilinear softening. Noteworthy the proposed model of combined damage is not restricted only to the above described multilinear CZMs, but can be used for any other model, for example, for those studied in [40].
In a real interface, both normal and tangential displacements and stresses appear. The damage may depend on the direction of the load, forming, after total rupture, a crack in Mode I, Mode II (generally in 3D, also mode III), or in a mixed mode. To include such a possibility into the present simplified model, the relations in (8) are modified so that K is a diagonal matrix expressed in the local coordinate system defined by the normal and tangential vectors providing the functionals where G c (u) reflects the fracture-mode dependence as e. g. in (2h). The stresses then similarly to (12) can be written in the form with ζ i =1, or ζ i being a solution of the system 1 2 K n ζi (ζ)u 2 n + 1 2 K t ζi (ζ)u 2 t =G ic (u), or ζ i =0, using the same three respective cases as described below Eq. (14), only the second one is a bit complicated if K uses different functions Φ for normal and tangential components. Making such a simplification for K (i.e. assuming identical functions Φ for both components), and denoting u n =|u| cos ψ, u t =|u| sin ψ, we may also write G ic (u) in the form G iIc g(ψ), see (2f). In compression (u n <0) we may consider no contribution to damage due to the term with k G in (4a). Then, for example, the normal stress component in (20) can be written as if we use K n and K t in the form of (10) with the same Φ i . As we have seen in relations (12) and (14), such a relation leads to a graph in the form of Figure 3, if we consider constant ψ in (21). Thus, for a fixed fracture-mode-mixity angle ψ we still may obtain a multilinear stress-displacement relation. An example of the total stress-displacement relations is shown in Figure 4.

3.
Numerical examples. The numerical procedures for the solution of the problems with cohesive contact model obeying the multilinear law described above are tested by solving various sample problems. Although the trilinear model is mostly used in these examples, also a more complicated case of tetralinear model is employed. Though, the model was described with viscosity in the bulk domains, the rate independent calculation without viscosity can be used as a reasonable approximation of the proposed model when the vanishing-viscosity concept [30,33] is taken into account. In this way the interface model is tested first neglecting viscosity, and only afterward some viscosity is included into the model to see its influence.
The numerical procedures used for solving the aforementioned problem includes time and spatial discretizations, as usual. In fact, these discretizations are equivalent to those described and used in [40], so only their basic features are sketched here. The time discretization uses a semi-implicit fractional-step algorithm in order to provide a suitable variational structure to the solved discretized problem [22,23,32,43] by separating the displacement and damage variables and recursive application of minimization with respect to the separated quantities. The procedures are developed so that the solved problem is formulated in terms of the boundary and interface data, with the spatial discretization carried out by SGBEM [41]. The whole model expressed in terms of the displacement variables allows for implementation of QP algorithms, based on the conjugate gradient schemes applied to problems with bound constraints [10]. As it can be guessed from Section 2.1, the functionals are not quadratic with respect to ζ. The optimization with respect to the damage variables thus uses QP sequentially (SQP).
3.1. Verification of the model by a two-squares problem. This simple example serves for verification of the stress-displacement relation at a selected point of the interface when solved by the computer code using SGBEM and SQP.
The geometry and load arrangement correspond to the problem of uniaxial tension, which tests the satisfaction of the stress-displacement relation of the used CZM. The plane strain problem layout is shown in Figure 5 The interface fracture and stiffness properties are characterized by the fracture energy in Mode I G c =7.5J m −2 , by the initial normal stiffness k n =20GPa m −1 and by the normal compression stiffness parameter k G =20TPa m −1 . As the load is tensile, the fracture-mode-mixity angle as well as the tangential stiffness of the adhesive are irrelevant. For characterizing the interface traction-relative displacement relation in the cohesive contact zone, the trilinear model defined by (16) with m d =2 is used. The parameters are set according to the two forms introduced in Figure 3, defined also by the relations (17) and (18). The characteristics of the traction-relative displacement law are shown in Figure 5(b), which in the Mode 'subsequent' lead to and in the Mode 'at once' to Notice that in any case G c =G 1c +G 2c and k n =k n1 +k n2 .
In this example, the BE mesh is uniform so that each face of a square subdomain of length a 1 is divided into 16 elements. The time step is τ =1ms.
The vertical load g (prescribed displacements) applied in time steps is the following function of time (24) up to the total time 0.4s, see also Figure 5(c). The bottom face of the structure is fixed. The applied load causes tension in the structure and in particular across the interface, and the stiffness of the interface is set to E/a 1 in order to have approximately the stiffness of the square bodies. Thus, there is no sudden change of displacement and the whole interface deforms in the same way.
The computed stress-and damage-displacement relations are shown in Figure 6. Due to variation of the loading, when the load is decreased, the damage of the interface stops evolving (the horizontal segments in ζ distribution), and the stressdisplacement relation obeys a linear relationship with a reduced stiffness until the previously achieved maximum is reached again. As loading is decreased three times, the previously described behaviour appears also three times, though the last one appears after the total damage of the interface.
The numerical data for the normal traction can be compared with those shown in Figure 5(c). They coincide very well, following the expected behaviour of the interface damage process. We can also see that the stress and displacement variations are independent of the way of defining the multilevel damage parameters.
3.2. Double cantilever beam. A typical example for analyzing a Mode I crack is the double cantilever beam (DCB). The dimensions used here are adopted from [13] and are shown in Figure 7(a). The material parameters of the beam are E=176.6GPa and ν=0.34 with no viscosity.
The initial length of the crack, propagating in the midplane, is ini . The bilinear CZM is considered. The considered interface traction-relative displacement law is shown in Figure 7(b). The split of the fracture energy provides G 1c =7.745kJm −2 and G 2c =1.34kJm −2 in the Mode 'subsequent'. Such a split corresponds to crack initiation at the crack tip and subsequent crack propagation with a large process zone, as explained e.g. for fibre composites by fibre bridging, see [13].
The mesh of boundary elements is slightly refined in the part of the interface Γ C where the crack is supposed to propagate, so in this zone the element length is 1mm and far from this region the elements are 2mm long. The time step is τ =20ms  g=3.30mm σ n ζ 1 ζ 2 Figure 8 Deformations of DCB, the damage evolution and normal stress distribution in the partially cracked interface at selected time instants corresponding to prescribed displacement g, ini is the initial crack length, cf. Figure 7.
and propagates in Mode I, the plots of normal stresses in the enclosed graph show, after the crack initiation phase, an approximately self-similar distribution near the crack tip, the zone where ζ 1 develops. Outside of this region the stresses are more or less constant. The influence of ζ 2 is hardly observable here (i.e. positive tensile stresses to the left of the propagating wave of normal stresses), but it will be demonstrated in Figure 9.
Although the normal stress distribution in the fracture process zone is observable in snapshots in Figure 8, the plot of the stress-relative displacement curve at a chosen point of arising crack documents even better the expected behaviour of the stress distribution, see Figure 9. The right graph focuses on the variation of Normal stress-relative displacement graphs at the interface point The damage range is kept the same in the right and left part, only the range for the normal stress σ n is changed in the right picture.
the normal stress after the crack is initiated. The curves also document how the expected stress-relative displacement relation of the proposed bilinear CZM formulation is followed in a chosen point. They also show, that the damage parameter ζ 2 begins to vary after ζ 1 has reached zero and that it also depends linearly on the displacement jump. Finally, Figure 10 presents the total applied force at the part of boundary of the lower beam, where the displacement load is applied.
F 1 F 2 Figure 10 The applied forces for DCB calculated at the place where the vertical displacement is imposed.
In addition to the expected vertical force also a smaller horizontal force appears due to fixing both ends of the beams in the horizontal direction. Nevertheless, it is the vertical force component which naturally reflects the structural behaviour for a propagating crack in the opening mode. Its maximum is closely related to the initiation of the crack propagation as can be reasoned by comparing the results in Figures 8 and 10. 3.3. A mixed mode crack with a tetralinear cohesive zone model. This example shows a more complicated CZM which was obtained by a combination of a trilinear and a trapezoidal CZMs. A mixed mode bending of a beam, shown in Figure 11(a), is studied. Figure 11 The mixed mode beam, cf. [40]; (a) the problem layout: =120mm, C =92mm, ini =8mm, h=20mm, w=2mm, s=2mm, (b) used stressdisplacement law in the cohesive zone: u 1 =2u 0 , u 2 =3u 0 , u c =4u 0 , σ 2 = 1 8 σ 0 and σ 0 =7.5 MPa, where u 0 =0.01 mm in normal component; in the tangential component either the same value (no dependence on mode-mixity) or u 0 =0.04 mm (mixed-mode dependent: G IIc =4G Ic ).
Initially, there is a crack of the length ini in front of the cohesive zone with the following characteristics: the fracture energy in Mode I G Ic =159.4J m −2 , and the maximum cohesive normal stress in pure tension σ n,max =7.5MPa, which uniquely defines the normal stiffness of the undamaged adhesive k n . The stress-relative displacement relation is supposed to obey the graph of Figure 11(b) in the cohesive zone, which can be defined in terms of parameters k i , G ic and β i calculated from (17) in the Mode 'subsequent' with m d =3. The initial stiffness in the tangential direction k t is the same as in the normal direction.
The crack propagates in a mixed mode so that two choices for the parameters in G ic from (2h) obtained by splitting G c by (2g) are considered: first, G IIc =G Ic (no sensitivity to fracture mode mixity) and, second G IIc =4G Ic (fracture mode sensitive). Both can be identified by an appropriate choice of u 0 defined in the caption to Figure 11. The elastic properties of the blocks are E=70GPa and ν=0.35, which naturally define the elastic stiffness tensor C. The viscosity moduli in D are considered to be D=τ r C with either τ r =10ms or τ r =0 (no viscosity).
The reference discretization is uniform, the BE mesh size (element length) is The vertical load g (again in terms of displacements) is applied incrementally in time steps so that the maximum applied displacement is g max =0.2mm and it is reached at a constant velocity 1mm s −1 .
The right face of the structure is fixed and the form of the interface crack is influenced by supporting the lower left corner of the bottom beam.
First, the total force response of the structure is observed. Figure 12 shows the evolution of reaction obtained in the simple support constraint with respect to the prescribed displacement g.  Figure 12 The total reaction forces for the mixed-mode beam calculated at the simple support constraint: (a) observing the influences of the fracture mode mixity (G IIc =4G Ic or G IIc =G Ic ) and viscosity (solid lines for no viscosity τ r =0, dashed lines for τ r =10ms), (b) changes of the solution for various discretizations. In Figure 12(a), the curves correspond to two aforementioned options of mixedmode sensitivity distinguished by various ratios of Mode II and Mode I fracture energies, either (G IIc =4G Ic or G IIc =G Ic ), and to two choices of the time-relaxation parameter τ r . The influence of the sliding mode in the crack propagation is clearly observed by later appearing of a softening effect in the force evolution for the case with higher fracture energy in Mode II. Also the natural tendency of viscosity to eliminate sudden changes in overall structural response can be observed.
The other graphs, shown in Figure 12(b), document how the changes in discretization affect the variation of the vertical reaction in the case τ r =0 and G IIc =4G Ic . Unlike all other plots in this section, here, the reference time step is ten times greater τ =2.5ms, the mesh reference number is N =4. The other curves are pertinent to τ =5ms (N =2) and τ =1.25ms (N =8). Accordingly, also the BE meshes are adapted. A convergence to an a-priori unknown force evolution is evident.
Comparing both graphs in Figure 12, we can see the difference obtained by refinement of the time step by the factor 10, if we focus on the curve associated to G IIc =4G Ic and τ r =0 in the part (a) and N =4 in the part (b) which use the same spatial meshes but different time steps.
The overall response of the structure can be also documented by the evolution of the energy. In (7), the balance of the energy for the evolution problem was formulated. However, in the numerical solution we obtain rather an imbalance, similarly to [40,43], which can be seen in Figure 13. For the same meshes and time-steps as used previously and for both viscid and inviscid cases the evolution of the left-hand and right-hand sides of (7) after discretization is shown. Unlike  Figure 13 The energy evolution and fulfillment of the energy balance (7) for various discretizations calculated for the mixed-mode beam, G IIc =4G Ic : (a) no viscosity τ r =0, (b) viscosity with τ r =10ms.
that equation, the relation is considered only for a time-step from time t till t+∆t: the left-hand side term containing the stored energy E at the instant t+∆t and the dissipated energy D within the time step, the right-hand side containing the stored energy E at the instant t and the work by the external forces W within the same time step. Additionally, as the time steps are different for various meshes, a power function is plotted dividing the energy by the size of the time step, therefore the axis shows ∆E/∆t. A convergence of the right-hand side to the left-hand side in the energy balance relation is again evident for both viscid and inviscid cases. The magnified deformations (the used factor is 50) of the beam are shown in Figures 14 and 15. The difference between them is that Figure 14 compares deformations with various fracture-mode-mixity sensitivities and Figure 15 compares different bulk viscosities. Anyhow, in the evolution of damage parameters, we can see the same behaviour of the tetralinear CZM, where each subsequent ζ i , i=1, 2, 3 starts to develop only after the previous one is brought to zero according to the proposed Mode 'subsequent'.
The differences taking into account fracture-mode sensitivity or viscosity are clearly observed in damage and stress distributions. In both figures, the maximal difference is developed during the rapid propagation of the crack (corresponding also to the softening period observed in Figure 12). As the crack propagation is slowed down in the final part of the loading process because of the prevailing compression, the viscous solution catches the non-viscous one up. The difference in the fracture-mode sensitivity figure is seen in the whole loading process.
Finally, we tried to plot stresses at selected points of the interface as functions of produced displacement jumps during the loading process. Figures 16 and 17 provide an insight to these relations, respectively for the cases G IIc =4G Ic or G IIc =G Ic , including also surfaces for expected stress-relative displacement relations similar to those plotted in Figure 4.
In both figures, the changing fracture mode can be guessed, because the ratio between the tangential and normal displacement jumps evolve during loading and crack propagation. The points were chosen so that this ratio is different: the point at x 1 =60mm is located below the applied displacement load so that in the initial part of loading there is a compression there, which is not shown in the graphs, the point at x 1 =89mm is close to the initial crack tip so that starting from the very beginning of loading the normal stress is tensile and it is combined with the g=0.05mm σ n,4 σ t, 4 σ n,1 σ t,1 g=0.10mm σ n,4 σ t,4 σ n,1 σ t,1 x 1 =0mm x 1 =92mm g=0.15mm σ n,4 σ t,4 σ n,1 σ t,1 g=0.20mm σ n,4 σ t,4 σ n,1 σ t,1 Figure 14 Deformations of the beam, the damage evolution and stress distribution in the cracked interface at selected time instants corresponding to the prescribed displacement g: comparison of the cases G IIc =4G Ic or G IIc =G Ic (referenced respectively by indices 4 and 1) with no viscosity.
tangential stress. The third point at x 1 =72mm is between them and it is initially exposed only to tangential stress as it can be seen that the pertinent lines lie in the plane of [[u]] n =0. g=0.20mm τ r =0ms τ r =10ms ζ=1 ζ=0 ζ i,0 ζ i,10 Figure 15 Deformations of the beam, the damage evolution and stress distribution in the cracked interface at selected time instants corresponding to the prescribed displacement g: comparison of the cases τ r =10ms and τ r =0 (referenced respectively by indices 10 and 0) with G IIc =4G Ic .
for a specific cohesive zone model (CZM) with a trilinear traction-relative displacement law, which is widely used in engineering calculations, e.g., in concrete fracture modeling.
Eventually, a general formulation covering any multilinear traction-relative displacement formulation in CZMs was developed here. Such a general formulation requires more scalar damage variables, whose number m d is by one less than the number of straight branches of the traction-relative displacement curve.  Figure 17 Stress-relative displacement graphs at selected points of the interface: τ r =0, G IIc =G Ic .
Regarding further generalizations, the methodology introduced is not restricted only to piecewise-linear traction-relative displacement curves, many types of wellknown CZMs can be combined, e.g., the exponential Ortiz-Pandolfi CZM [24] and the Raous model [29], one only needs to decide when to initiate each particular model in the overall damage process.
Numerical tests for the multilinear CZMs, in particular trilinear and tetralinear, showed very satisfactory and promising results. Hence, this new formulation is expected to be used in the solution of more complicated engineering problems of fracture of concrete and composite materials in the future.