A Boundary Element Approach for an Interface Visco-Damage Model Exposed to Cyclic Shear Load

Cracks that appear between firmly connected parts of civil engineering constructions may substantially modify the mechanical properties of such structures, therefore, development of numerical approaches for solving such problems is highly demanded. In this paper, a model for interface damage evolution is presented. It may simulate total damage at the interfaces for displacement or force cyclic loading with a constant amplitude, which may occur for amplitudes substantially smaller than for a purely increasing load. The interface is supposed to be made of a very thin adhesive layer which may provide a stress-separation law as in a cohesive zone model (CZM) under uploading conditions and which may also mimic fatigue processes, including a rate-dependent damage evolution, similarly to cyclic CZMs introduced in [1–3].

The present solution approach introduces a model working with energy evolution, which involves the stored elastic energy, energy of the external load and also dissipation due to damage. To simulate the fatigue processes, the dissipation potential has a rate-dependent form with hysteresis associated to damage, see also [4–6]. The CZMs in the context of the energy evolution were implemented in [7, 8] and a simple viscous term was also tested in [9]. In the present implementation, the viscous term and also the damage triggering term of dissipation are formulated as state dependent in order to provide a better control of the hysteretic model so that the unloading and reloading paths were non-coincident. Additionally, the model is capable of capturing varying transferable force with increasing damage and to keep the interface undamaged for low stresses.

For the proposed model, also a numerical approach is introduced. The solution evolution is approximated by an appropriate time stepping algorithm which is capable of providing a variational structure to the solved problem by two recursively solved minimisations with respect to variables of deformation separated from damage variables as it was devised in [10, 11]. Similar methodology was also described and implemented in the previous author’s works [7, 8, 12]. The variational character of the solved problem requires algorithms of energy minimisation. The properties of the proposed model permit one to use quadratic programming (QP) algorithms based on conjugate gradient schemes of [13, 14], or, if the functionals are not quadratic, a sequential modification of QP algorithms can be utilised [15]. Finally, the numerical solution for spatial discretisation includes the symmetric Galerkin boundary element method (SGBEM) [16, 17] with its multi-domain form of [18] to provide the displacement and stress solutions for the solids. The advantage of this approach in the present model is obvious as no nonlinear phenomena are considered in the interior of solids.

The paper is structured as follows. Section 2 describes the interface damage model and the relations which govern its evolution. Section 3 explains some aspects of the numerical solution and its particular implementation. Finally, Section 4 presents numerical results which demonstrate the properties of the model and reveal its characteristic response to cyclic shear type loading.

The model considers a domain $\Omega$ consisting of two or more connected parts. Let us consider a domain containing two subdomains $\Omega^\eta, \eta=A$, or $B$, bounded by the respective Lipschitz boundaries $\Gamma^\eta$ as shown in Figure 1. A common part of the boundaries (an interface) is denoted $\Gamma_{\mathrm{C}}=\Gamma^A \cap \Gamma^B$. Each of the boundaries is additionally split into disjoint parts according to boundary conditions: $\Gamma_{\mathrm{D}}^\eta$ where displacements are prescribed, $\Gamma_{\mathrm{N}}^\eta$ where tractions are given, i.e. $\Gamma^\eta=\overline{\Gamma_{\mathrm{D}}^\eta} \cup \overline{\Gamma_{\mathrm{N}}^\eta} \cup \overline{\Gamma_{\mathrm{C}}}$. The normal and tangential vectors to the boundaries are denoted $\mathbf{n}^\eta$ and $\mathbf{s}^\eta$, respectively.

The model is governed by energy evolution. Its state is described by two variables: the displacement field $\mathbf{u}$ and an internal variable $\zeta$ characterising the actual damage level so that $\zeta$ =1 pertains to the initial undamaged state and $\zeta$ = 0 reflects total damage when a crack appears.

The energy of the system describing the behaviour of the model includes the stored energy and also dissipated energy. The stored energy functional is defined as

where we have taken into account the bounds of the damage variable and the time dependent boundary condition on a part of boundary $\Gamma_{\mathrm{D}}^\eta$ which restricts admissible displacement field. The fourth-order tensor $\mathrm{C}^\eta$ includes elastic constants for the domain $\Omega^\eta$ and e is the small strain tensor. The interface integral defines stored energy of the adhesive, depending on the initial undamaged interface stiffnesses $k_{\mathrm{n}}$ and $k_{\mathrm{s}}$, and on the actual state of damage controlled by the degradation function $\Phi(\zeta)$ as in [7, 8], e.g. an exponential law of softening can be obtained by the function $\Phi(\zeta)=\frac{\varphi(a \zeta+\varepsilon)-\varphi(\varepsilon)}{\varphi(a+\varepsilon)-\varphi(\varepsilon)}$ with $\varphi(\zeta)=\exp \left(-\gamma_r^{-1}(\zeta)\right)$ and $\gamma_r$ being the upper incomplete gamma function $\gamma_r(x)=\frac{1}{\Gamma(r)} \int_x^{+\infty} t^{r-1} \mathrm{e}^{-t} \mathrm{~d} t, r$ is the principal parameter of the model, usually set to three to be consistent with common cohesive zone models. Other parameters $\alpha$ and $\varepsilon$, satisfying the conditions $\gamma_r\left(\frac{r-1}{2}\right)\it a \it 1,0 \lt \varepsilon \lt 1-a$ help to control the damage initiation and termination, especially in the numerical solution, theoretically they can be set to the values 1 and 0, respectively.

The terms [[$\cdot$]]$_{\mathrm{n}}$ and [[$\cdot$]]$_{\mathrm{s}}$, respectively, are the normal and tangential displacement gaps between opposite interface points, i.e. $[[ \mathbf{u} ]]=\mathbf{u}^A-\mathbf{u}^B$ (superscripts $A$ and $B$ refer to the bodies adjacent to the interface from both its sides). The last term with [[$\cdot$]]$_{\mathrm{n}}$, denoting the negative part of the relative normal displacement, reflects the normal-compliance contact condition being understood as a reasonable approximation of Signorini contact conditions due to asperities of the contacting surfaces.

The potential energy of external forces (acting only along a part of the boundary denoted $\Gamma_N^\eta$) is given by the relation

where, $\mathbf{f}^\eta$ are the prescribed tractions on the part of boundary $\Gamma_N^\eta$.

The dissipated energy is introduced by the (pseudo)potential which reflects the rate dependence and unidirectionality of the delamination process

where, $\alpha_1$ is related to the instant of damage initiation and its dependence on the current state of damage enables the endurance limit for fatigue life to be modified as long as damage approaches the total rupture; and $\alpha_2$ is visco-damage parameter which introduces hysteresis into the damage evolution and again its dependence on the current state of damage modifies the range of viscous effect. How these parameters affect the stress-separation diagram is shown in Figure 2, where the power relations were used for both $\alpha_1$ and $\alpha_2$ : $a_i(\zeta)=a_{i 0}\left(\zeta+\varepsilon_i\right)^{r_i}$, with $\varepsilon_i$ = 0.001 (set especially for the cases where $p_i<0$). The same power rules will also be used below.

The relations which govern the elastic state evolution in the solid $\Omega$ and damage evolution in the interface $\Gamma_C$ can be written in a form of nonlinear variational inclusions with initial conditions.

where, $\partial$ denotes (partial) subdifferential of a convex non-smooth function which can be replaced by Gateaux differential if the functional is sufficiently smooth, for e.g. $\mathcal{F}$. The initial condition for damage is usually $\zeta_0=1$, pertaining to the undamaged state.

The system provides also an energy balance relation which can be obtained by testing Eqn. (4), respectively, by $\dot{\mathbf{u}}$, and $\zeta$ as follows:

where, $\mathcal{D}(\zeta ; \dot{\zeta})=\left\langle\partial_{\dot{\zeta}} \mathcal{R}(\zeta ; \dot{\zeta}), \dot{\zeta}\right\rangle=\int_{\Gamma_{\mathrm{c}}}-a_1(\zeta) \dot{\zeta}+a_2(\zeta) \dot{\zeta}^2 \mathrm{~d} \Gamma$ denotes the overall dissipation rate. The change of stored energy and the energy dissipated during a specified time interval then sum to work done by the external force and displacement loads.

Numerical solution to the model described above requires both time discretisation and spatial discretisation. In what follows, we separately provide some properties of both. Let us stress that due to SGBEM being used to calculate elastic response of the solids, the numerical procedures are expected to be expressed in terms of boundary data, unlike the stored energy in Eqn. (1).

The time discretisation which ensures convergence of the numerical solution includes decoupling of the system (4) by a fractional-step method: separation of variables then may be done to guarantee separate convexity of the energy functional $\mathcal{E}$ and additive splitting of the dissipation potential $\mathcal{R}$.

The discretisation scheme is introduced by a semi-implicit algorithm with a fixed time step $\tau$ such that the solution is obtained at the instants $t^k=k \tau$ for $k=1, \ldots, T / \tau$ and denoted $\mathbf{u}^k$ for displacements and $\zeta^k$ for the damage variable. In order to obtain such an algorithm from Eqn. (4), the rate of damage is approximated by the finite difference $\dot{\zeta} \approx \frac{\zeta^k-\zeta^{k-1}}{\tau}$. The differentiation with respect to the rate is accordingly replaced by the differentiation with respect to $\zeta^k$, so that from Eqn. (4) we obtain

$\tau \partial_{\zeta^k} \mathcal{R}\left(\zeta^{k-1} ; \frac{\zeta^k-\zeta^{k-1}}{\tau}\right)+\partial_{\zeta^k} \mathcal{E}\left(t^k ; \mathbf{u}^k, \zeta^k\right) \ni 0$ (6b)

which are solved separately to calculate $\mathbf{u}^k$ and $\zeta^k$, respectively. Let us notice that for characterising the state we used $\zeta^{k-1}$ in $\mathcal{R}$, which is quite reason able while we consider nonzero hysteretic term $\alpha_2$.

The separation of variables provides a variational structure to the solved problem, containing two minimisations in each time step: the first minimisation with respect to displacements of

provides $\mathbf{u}^k$, and the second minimisation of

$\mathcal{H}_\zeta^k(\zeta)=\mathcal{E}\left(t^k ; \mathbf{u}^k, \zeta\right)+\tau \mathcal{R}\left(\zeta^{k-1} ; \frac{\zeta-\zeta^{k-1}}{\tau}\right)$ (7b)

provides $\zeta^k$.

The minimisation system (7) is solved recursively for $k=1, \ldots, T / \tau$, starting from $k=1$ with initial conditions $\mathbf{u}^0=\mathbf{u}_0, \zeta^0=\zeta_0$ from Eqn. (4). In fact, the functional $\mathcal{E}$ is quadratic with respect to $\mathbf{u}$ so that various efficient numerical QP algorithms can be used in the solution, see e.g. [14]. On the other hand, with respect to $\zeta$ it is quadratic if the function $\Phi(\zeta)$ is quadratic. Nevertheless, our present choice is not covered so that the QP algorithm should be applied sequentially as in [7].

The SGBEM algorithm is used to calculate the elastic strain energy in Eqn. (7a), in fact from Eqn. (1), and it provides the complete boundary-value solution from the given boundary data.

Therefore, it is convenient to change the domain integral to boundary integrals. First, we introduce the linear mapping $\mathcal{S}:(\mathbf{g}, \mathbf{f}, \mathbf{w}) \mapsto \mid \mathbf{u}$ with $\mathbf{u}$ being the weak solution to the transmission boundary value problem:

$\mathbf{u}^\eta=\mathbf{g}^\eta \text { on } \Gamma_{\mathrm{D}}^\eta \text {, }$, (8b)

$\mathbf{p}^\eta=\mathbf{f}^\eta \quad$ with $\quad \mathbf{p}^\eta=\left(\mathrm{C}^\eta \mathrm{e}\left(\mathbf{u}^\eta\right)\right) \cdot \mathbf{n}^\eta \quad$ on $\Gamma_{\mathrm{N}}^\eta$, (8c)

$\left.\begin{array}{l}[[ \mathbf{u} ]]=\mathbf{w} \\ \mathbf{p}^A+\mathbf{p}^B=0\end{array}\right\} \quad$ on $\Gamma_{\mathrm{C}}$, (8d)

where, $\mathbf{w}$ is supposed to be prescribed displacement gap at $\Gamma_C$ and $\mathbf{p}$ denotes boundary (and interface) tractions.

The functional $\mathcal{H}_u^k$ in Eqn. (7a) is converted, using also Eqn. (1), to $\mathcal{H}_{\mathrm{w}}^k$ as

where, $\mathbf{u}^\eta\left(\mathbf{g}\left(t^k\right), \mathbf{f}\left(t^k\right), \mathbf{w}\right)$ and $\mathbf{p}^\eta\left(\mathbf{g}\left(t^k\right), \mathbf{f}\left(t^k\right), \mathbf{w}\right)$ are, respectively, displacement and tractions obtained by Eqn. (8) through the (generalised) Poincaré-Steklov operator as described in [19]. They are obtained by the SGBEM taking the known displacement gap $\mathbf{w}$ from the minimisation procedure in the same way as in [18, 19]. Once all the boundary data (displacements and tractions) are obtained from the solution of the SGBEM code, the pertinent energy functional in Eqn. (9) can be calculated.

The implementation of SGBEM uses the weighted system of BIEs written in the following block form, cf. [18, 20]:

with weighting functions $\psi$ and $\phi$ defined by the element shape functions used for approximation of displacement $\mathbf{u}$ and traction $\mathbf{p}$. The introduced matrices are defined as follows:

$\mathbf{K}^\eta=\left(\begin{array}{cccc}-\mathbf{U}_{\mathrm{DD}}^\eta & \mathbf{T}_{\mathrm{DN}}^\eta & -\mathbf{U}_{\mathrm{DC}}^\eta & \mathbf{T}_{\mathrm{DC}}^\eta \\ \mathbf{T}_{\mathrm{ND}}^{\eta^*} & -\mathbf{S}_{\mathrm{NN}}^\eta & \mathbf{T}_{\mathrm{NC}}^{\eta^*} & -\mathbf{S}_{\mathrm{NC}}^\eta \\ -\mathbf{U}_{\mathrm{CD}}^\eta & \mathbf{T}_{\mathrm{CN}}^\eta & -\mathbf{U}_{\mathrm{CC}}^\eta & \omega^\eta \frac{1}{2} \mathbf{I}_{\mathrm{CC}}^\eta+\mathbf{T}_{\mathrm{CC}}^\eta \\ \mathbf{T}_{\mathrm{CD}}^{\eta^*} & -\mathbf{S}_{\mathrm{CN}}^\eta & \omega^\eta \frac{1}{2} \mathbf{I}_{\mathrm{CC}}^\eta+\mathbf{T}_{\mathrm{CC}}^{\eta^*} & -\mathbf{S}_{\mathrm{CC}}^\eta\end{array}\right), \quad \omega^\eta= \begin{cases}-1 & \text { if } \eta=\mathrm{A}, \\ 1 & \text { if } \eta=\mathrm{B},\end{cases}$

The above equations use the (weakly singular) Kelvin fundamental solution $U_{i j}^n(x, y)$, and the associated derivatives obtained by the differential traction operator used also in Eqn. (8c)-the strongly singular function $T_{i j}^n(x, y)$ and the hypersingular function $S_{i j}^n(x, y)$. The compact form of BIEs in Eqn. (10) uses the notation

where, $\omega$ stands for $\phi$ or $\psi$, while $v$ stands for $u$ or $p$, further, $q$ and $r$ stand for $\mathrm{D}, \mathrm{N}$, and C , and eventually $Z^\eta$ stands for $U^\eta, T^\eta, T^{\eta^*}$ or $S^\eta$, and where the inner integral can be regular, weakly singular, Cauchy principal value or Hadamard finite part integral. In the previous relations, $\mathbf{I}^\eta$ denotes the identity operator with the subscripts and superscripts specifying the part of the boundary where it is restricted to.

We use the model to find solutions for various amplitudes and various maxima of the same load type applied to the domain shown in Figure 3 (see also [6–8]). The model formulation is implemented in a MATLAB computer code.

The material properties of both layers in our example are the same: $E$ = 70 GPa, $v$ = 0.35. The upper layer is loaded either by a time-dependent displacement load with the maximum $g_{\text {max }}=v t_0$ or by a time-dependent force load with the maximum $f_{\text {max }}=p t_0$, the respective mean values $m_g, m_f$ and amplitudes $a_g, a_f$ may vary, introducing a parameter $\rho$ such that $a_g=\rho g_{\text {max }}$, $a_f=\rho f_{\text {max }}$. It is also shown in Figure 3, where as an example, the displacement loading function uses $\rho=0.25$ and the force loading function uses $\rho=0.5$. The numerical values of the other loading parameters used in tests are $v=0.1753 \mathrm{~mm} \mathrm{~s}^{-1}, p=100 \mathrm{MPa} \mathrm{s}^{-1}$. For $t_0$, various values are used to obtain various upper limits of the load.

The characteristics of the interface needed for functionals in Eqns. (1) and (3) (and explanatory text below them) are $\alpha_{10}=21.5 \mathrm{~J} \mathrm{~m}^{-2}, \alpha_{20}=100 \mathrm{~J} \mathrm{~s} \mathrm{~m}^{-2}$, the exponents $r_1, r_2$ have various values in the numerical tests to show their influence on the solution, the initial interface stiffnesses are $k_{\mathrm{n}}=2 k_{\mathrm{s}}=3.88 \mathrm{TPa} \mathrm{m}^{-1}$, and $k_{\mathrm{g}}=14 \mathrm{TPa} \mathrm{m}^{-1}$. Finally, for the degradation function $\Phi(\zeta)$, we used $r=3, \alpha=0.99, \varepsilon=0.005$.

The used boundary element discretisation is uniform such that the mesh size (the element length) is $h=2.5 \mathrm{~mm}$, the time-stepping is performed in steps of $\tau=12.5 \mathrm{~ms}$.

The first series of results presents a global response of the layered structure to both types of loading. As a measure, we defined the stiffness $K$ which relates the average force and average displacement applied or enforced at the same locus of load application at the right face of the upper layer. We first fix the parameters $\rho, r_1, r_2$ and let vary the maximum of the load introduced by $t_0$. The graphs in Figure 4 show how the stiffness decreases due to damage evolution. The larger $t_0$ is, the smaller number of cycles $N$ is needed to form an interface crack. At the end of the loading, the stiffness falls down to zero as the layers separate totally. The stiffness decreases faster for force loading which is in close connection to the faster damage (and crack) propagation in this case as we will see below.

The influence of the amplitude of the load, with the maximum of the load kept fixed, is shown in Figure 5. Here, the parameters $r_1$ and $r_2$ are the same as before and maximum of the load corresponds to $t_0=0.4 \mathrm{~s}$. The parameter $\rho$, which controls the amplitude, varies. In any type of loading, it is clearly seen that intermediate values of $\rho$ provide almost the same curve, which is in agreement with fatigue behaviour usually displayed by Haigh diagram, see [21]. For $\rho$ close to unity the damage propagates not only around the maximum, but also close to the minimum (negative) of the load which is relevant in the present shear type of loading.

Finally, for the displacement type of loading, we also intend to show how the parameters $r_1$ and $r_2$ modify the stiffness changes in Figure 6. The increasing parameter $r_1$ causes the endurance limit to decrease as damage progresses so that the number of cycles needed to make corresponding damage also decreases. The increasing parameter $r_2$ causes the decrease of hysteretic effect and thus decreases the number of required cycles to cause a crack propagation. Here, we also used negative values which may cause an opposite result. However, it is important to see that it is necessary to adjust also the parameter $\varepsilon_2$ to a positive value, for the viscosity not to be infinite as $\zeta$ approaches zero.

Next, we show the speed of the crack length increase when the crack propagates. Table 1 summarises the number of time steps $n_a$ needed to obtain cracks of specified lengths $a$ in the sense that $\zeta$ at the pertinent nodal point of the interface falls below 0.1 for particular parameters $r_1, r_2$ and $\rho$. With given $\tau$, the number of cycles $\Delta N$ which is needed to extend the crack from the length $a_1$ to $a_2$ is $\Delta N=\left(n_{a_2}-n_{a_1}\right) \tau /\left(4 a_f\right)$. In fact, we had solved the problem until we reached a crack along the whole interface.

The data from Table 1 can be used to calculate the parameters of the crack-growth power law known for fatigue processes, see [22]. Linearised regression is used to obtain $\Delta a / \Delta N=2.71 \times 10^8 \cdot g_{\text {max }}^{6.05}\left(R^2=0.9903\right)$ for the displacement loading. Here, the load maximum $g_{\text {max }}$ is applied instead of stress intensity factor used in the aforementioned Paris law. Similarly, the power relation for the force load reads $\Delta a / \Delta N=7.25 \times 10^{-8} \cdot f_{\max }^{5.17}\left(R^2=0.9946\right)$, where the maximum of the applied force is used as the independent variable.

The damage evolution in Figure 7 for one particular $t_0=0.3 \mathrm{~s}$ and the same parameters can be used to check some values from Table 1. Let us look, for example, at the force load. The crack reaches the length of $a=10 \mathrm{~mm}$ (which corresponds to $x_1=170 \mathrm{~mm}$) after $n_a=1610$ time steps. For time span $t_0$ we have 24 time steps, therefore $n_a$ in terms of load cycles makes $N=33.29$. Similarly, the length of $a=17.5 \mathrm{~mm}\left(x_1=162.5 \mathrm{~mm}\right)$ is reached at the step $n_a=1780$, which corresponds to $N=36.83$. Both these values can be guessed from the left picture in Figure 7.

Althought the evolution of deformation in Figure 8 and Figure 9 is not very interesting, at least, we can compare the two loading types, here. We used the results for the same $t_0=0.3 \mathrm{~s}$ as before. The figures show the maximal displacements at selected instants after $N$ cycles, starting from $N=1$, which are comparable. For the force load, the crack spreads along the whole interface after $N=41$ cycles, while in case of the displacement load, the structure withstands $N=118$ cycles. These facts can be compared with the damage propagation in Figure 7, too.

Finally, a local behaviour of the model can be represented by the traction-separation relation and evolution of damage at a point of the interface. Both of these relations are shown in Figure 10 and Figure 11. A point in the middle part of the interface has been chosen: $x_1=127.5 \mathrm{~mm}$. The graphs also show how the interface stiffness (the slope of the stress-separation relation) is decreasing during the loading process. The number of cycles to provoke the total damage ($\zeta=0$) is naturally different for displacement and force loads. In the case of the smallest load maximum, the cycles cannot be easily counted in these graphs due to their large number. Nevertheless, for the other maxima it is possible: the graph of $\zeta$ evolution provides at the instant $t_0=0.3 \mathrm{~s}$ the total time of rupturing at the selected point $t=24.85 \mathrm{~s}$ and $t=69.3 \mathrm{~s}$ for the force load and for the displacement load, respectively. These two values expressed in number of cycles $N$ are respectively $N=41.42$ and $N=115.5$ and can again be checked by the graphs in Figure 7.

A model for solving interface damage has been presented. Its parameters have been set to such values that the model behaviour is sensitive to cyclic loading. The presented results confirm that the model is capable of propagating the damage in the sense of fatigue life. There are several parameters in the model which influence some aspects of the solution, and which were also studied in the example.

From the computational point of view, to obtain the solution of an evolution problem, semi-implicit time discretisation was implemented leading to a numerical approximation model with a variational structure. The approximated model then utilises the method of sequential quadratic programming with a space discretisation by SGBEM. The choice of the latter mentioned numerical technique enables the formulation of the minimisation problems only in terms of interface variables such as displacement gap and interface damage, if all nonlinear processes in the model are accumulated only at the interface.

Of course, there are many parameters in the model, some of them were discussed in the paper some have been kept fixed. All of them should be appropriately adjusted to follow satisfactorily experimental observations. Nevertheless, it might be expected that the proposed approach turns out to be useful in practical engineering calculations.