A BEM formulation applied in the mechanical material modelling of viscoelastic cracked structures

The present study aims at performing a mechanical analysis of 2D viscoelastic cracked structural materials using the Boundary Element Method (BEM). The mesh dimensionality reduction provided by the BEM and its accuracy in representing high gradient fields make this numerical method robust to solve fracture mechanics problems. Viscoelastic models address phenomena that provide changes on the mechanical material properties along time. Well-established viscoelastic models such as Maxwell, Kelvin–Voigt and Boltzmann are used in this study. The numerical viscoelastic scheme, which is based on algebraic BEM equations, utilizes the Euler method for time derivative evaluation. Therefore, the unknown variables at the structural boundary and its variations along time are determined through an ordinary linear system of equations. Moreover, time-dependent boundary conditions may be considered, which represent loading phases. The dual BEM formulation is adopted for modelling the mechanical structural behaviour of cracks bodies. Three examples are considered to illustrate the robustness of the adopted formulation. The results achieved by the BEM are in good agreement with reported data and numerical stability is observed.


Introduction
The application of structures composed by viscoelastic materials has increased in the last years. The structural design of mechanical components by composites and polymers, in mechanical engineering, and concrete, in civil engineering, has enlarged the application of such materials in engineering fields. Viscoelastic materials present elastic and viscous mechanical properties which cause relaxation, creep and hysteresis in structures (Zhu et al. 2011(Zhu et al. , 2013. In the context of mechanical modelling of viscoelastic materials, the analytical approaches available are restricted to a limited class of problems, in which particular boundary conditions and structural geometries are assumed (Oliveira and Leonel 2013) As a result, the robust mechanical modelling, which accounts for complex boundary conditions and structural geometries, requires the application of numerical techniques.
The Boundary Element Method (BEM) is a numerical technique widely applied in the literature for modelling fracture problems (Zhu et al. 2011(Zhu et al. , 2013Oliveira and Leonel 2013). Due to the mesh dimensionality reduction provided by the BEM, the remeshing procedures during the crack growth become a less complex task. Moreover, its accuracy in representing stress concentration make the BEM a recommendable numerical technique for solving fracture mechanics problems, especially into the viscoelastic domain. Classical BEM formulations solve viscoelastic problems using Laplace transform method coupled to the principle of elastic-viscoelastic correspondence (Liu and Antes 1997;Rizzo and Shippy 1971). This approach requires a convolutional relation between strain and stress tensors, which results into a space transformation (Lee and Westmann 1995).
A time-dependent BEM formulation was presented by Syngellakis andWu (2004, 2008) for the mechanical analysis of viscoelastic structures subjected to quasi-static and dynamic effects. In these studies, the principle of correspondence was applied to the determination of a fundamental solution, which was utilized for the fracture modelling of polymers. An integral formulation based on the Kelvin's fundamental solution was presented by Ashrafi and Farid (2009). Material properties were represented by time-dependent functions to analyze the viscoelastic mechanical behaviour of polymers.
The viscoelastic reciprocity principle was utilized by Cezario et al. (2011) to the development of a time-dependent BEM formulation. Such a formulation is defined considering the Stieltjes integral and material functions provided by the Boltzmann's model. The symmetric Galerkin BEM was applied by Perez- Gavilan and Aliabadi (2001) into dynamic analyses of viscoelastic structures in frequency domain. The mechanical analyses of stationary rolling contact problems considering linear-viscoelastic rollers were presented by Abascal (2004, 2006). This type of problem is efficiently solved by BEM once high stress and strain gradients are present at the contact surface, i.e., at the body's boundary. The applied formulation is based on the BEM and the correspondence principle neglecting all the inertial effects.
The mechanical analysis of viscoelastic materials by the BEM may be performed through alternative formulations, as proposed by Mesquita and Coda (2001, 2002, 2003. Such formulation utilizes integral equations written exclusively at the bodies' boundaries, which are obtained from the weighted residual technique. A time-dependent differential system of equations is achieved considering Kelvin-Voigt and Boltzmann models. Then, a proper time marching procedure is applied to solve the time-dependent problem. This formulation is simpler, as well as accurate. Therefore, this BEM formulation was adopted in the present study. Experimental and theoretical researches involving fracture analysis of viscoelastic materials have been observed in the literature. However, numerical approaches in this scientific field are relatively limited (Syngellakis 2002). In the context of BEM models applied into the viscoelastic fracture mechanics problems, it is worth mentioning (Sladek et al. 1984). In this work, a numerical approach based on the Laplace transform is presented to determine the crack opening displacements into a penny shaped crack. The mechanical material behaviour is represented by the Kelvin-Voigt model. (Sun and Hsiao 1988) applied a BEM approach to evaluate displacements and stresses fields surrounding a crack filled with failed. The mechanical analysis of viscoelastic anisotropic bodies containing holes, inclusions and initial defects were performed by Chen and Hwu (2011). In this work, the correspondence principle was used considering the sub-region BEM approach. The displacement discontinuity method (DDM) was utilized by Wang and Birgisson (2007), Birgisson et al. (2002Birgisson et al. ( , 2004 for mechanical analyses of asphalts. A DDM time-dependent integral approach is presented and applied into quasistatic cases. The functional of total potential energy was utilized by Lee and Kim (1995) to determine an expression for the strain energy released rate. Such functional was obtained using the direct time domain BEM.
As previously presented, a small amount of researches are available in the literature concerning the viscoelastic modelling of fracture mechanics problems by numerical approaches based on the BEM. Therefore, the present study aims to contribute into this scientific domain. A numerical model is presented, which is composed by the coupling of Maxwell, Kelvin-Voigt and Boltzmann viscoelastic models to the BEM algebraic equations. The main contribution of this study is the mechanical analysis of nonhomogeneous viscoelastic bodies in fracture conditions. The timedependent problem is solved using the explicit time marching process introduced in Mesquita and Coda (2001, 2002, 2003. Then, the one-step Euler method is adopted to approximate the required time derivatives. The mechanical behaviour of viscoelastic cracked materials is represented by the dual BEM formulation (Portela et al. 1992;Leonel and Venturini 2010;. This BEM approach is the most popular to model the mechanical behaviour of materials containing cracks. In the dual BEM, the singular integral representation is written along the boundary elements positioned at one crack surface, whereas the hyper-singular integral representation is applied along the opposite crack surface. The singular integral representation is applied along the entire external boundaries.
Two applications of single structural viscoelastic material are presented to illustrate the accuracy of the implemented BEM approach. In these applications, experimental and analytical results available in the literature are compared against the responses achieved by the BEM. In addition, structures composed of nonhomogeneous viscoelastic materials are also modelled, which is the main contribution of the present research. The nonhomogeneous structures are modelled using the sub-region technique. This BEM technique enforces compatibility of displacements and equilibrium of forces along the interfaces of all multiple bodies that compose the nonhomogeneous structure . One application concerning nonhomogeneous viscoelastic cracked structure is presented. This type of structure has not being properly addressed in the literature. Therefore, the application presented in this study serves as benchmark for future numerical researches addressing this subject.
It is worth mentioning that the responses achieved by the implemented BEM approach presents good agreement to references adopted in literature and numerical stability as well.

Review over the viscoelastic models
The mechanical behaviour of viscoelastic materials is idealized by simple elements utilized into the electrical circuit's theory (Tschoegl 1989). Such elements are dashpots and springs, which represent viscous and elastic material rheological characteristics, respectively. Then, the convenient composition of such elements lead to the representation of different viscoelastic models.
The Maxwell's viscoelastic model is represented by a series scheme composed by one dashpot and one spring, as illustrated in Fig. 1a.
The following group of equations governs this viscoelastic model: where r and e represent normal stress and normal strain, respectively. E is the Young's modulus and g the viscous coefficient. The superscripts visc and elast indicate viscous and elastic components, respectively. The dot over the variables represent the condition of time variation and t t a function of time. The Maxwell's viscoelastic model is governed by a differential equation obtained from Eq. (1). The compatibility condition presented in Eq. (1) has to be differentiated with respect to the time. Then, equilibrium condition and constitutive relation are applied to obtain the following representation: The total strain is obtained by integrating Eq.
(2) with respect to the time. Then: Equation (3) has to be integrated by parts, which leads to the following: in which D t t À s ð Þ¼ 1 E þ t t Às g is denominated creep function for viscoelastic Maxwell's model. It is important emphasizing that D(t t -s) is a linear function with respect to the time. Then, according to this model, the material flows indefinitely if a time-constant stress state is observed. Analogously, the relaxation function for Maxwell's viscoelastic model is defined as follows: The composition illustrated in Fig. 1b represents the Kelvin-Voigt's viscoelastic model. A parallel association of one dashpot and one spring characterizes such a viscoelastic model. The governing equations for Kelvin-Voigt model are the following: Such a viscoelastic model is governed by a differential equation, which is written by introducing compatibility condition and constitutive relation into the equilibrium condition, both of them presented in Eq. (6). These algebraic procedures lead to the following: The solution of the last differential equation is obtained as follows, for a given time-history stress state: The creep function for Kelvin-Voigt viscoelastic model is obtained by integrating parts of Eq. (8). Then: The creep function presented in Eq. (9) tends asymptotically to the elastic solution when t t tends to infinity, i.e., e(?) = r 0 /E. When this condition is observed, the total stress is supported by the spring element, which is directly evaluated by Eq. (7). The relaxation function for Kelvin-Voigt viscoelastic model is not available.
The viscoelastic Boltzmann's model is represented by a composition of dashpot and springs in series and parallel scheme as illustrated in Fig. 1c. The following equations govern this viscoelastic model: in which the superscript v/e indicates the values evaluated on the parallel sequence scheme of spring and dashpot.
The Boltzmann's viscoelastic model is mathematically expressed by a differential equation, which is obtained from the differentiation of the compatibility condition presented in Eq. (10) with respect to the time. Then, the constitutive relation and equilibrium conditions are applied to obtain the following representation: The strain on materials governed by such viscoelastic model is determined as follows: And the stress state by the following equation: cates the relaxation function for Boltzmann's viscoelastic model. This relaxation function shows that, during the external load application, both elastic and viscous strains are developed.

Integral equations of BEM
The integral equations required by the BEM in elastostatics are obtained from the differential representation of equilibrium, which is written in terms of displacements as follows: in which l represents the material shear modulus, t the material Poisson's ratio, u l indicate the displacements components and b l the body forces. The singular integral BEM equation is obtained by applying the Betti's theorem or weighted residual techniques. This equation is presented below disregarding the body forces.
in which u k and P k indicate displacements and tractions at the body's boundary, respectively. The free term c lk is equal to the Kroenecker operator multiplied by 0.5, for smooth contours. u lk * and P lk * represent the fundamental kernels for displacements and tractions, respectively (Portela et al. 1992).
It should be mentioned that Eq. (15) is capable to model the mechanical behaviour of 2D elastic bodies. However, the single application of this equation leads to the singularities over the final system of algebraic equations, when cracked bodies are considered. Thus, in such cases, different BEM approaches have to be adopted.
Among the BEM formulations available in the literature for this purpose, i.e., the mechanical analysis of cracked structures, the dual BEM formulation is the most popular (Oliveira and Leonel 2013a, b;Syngellakis and Wu 2008). The singular integral representation is applied into the discretization of one crack boundary in such BEM approach. In addition, the oppositely crack boundary is discretized by the hyper-singular integral representation. The entire external boundary is discretized by the singular integral representation.
The hyper-singular integral representation mentioned in the last paragraph is obtained from Eq. (15). Equation (15) is written in terms of displacements. Thus, by differentiating it, one obtains an integral representation written in terms of strains. Afterwards, Hooke's law is applied to obtain an integral representation based on stresses. Then, the equilibrium of Cauchy is applied to obtain an integral representation written in terms of tractions, which is the well-known hyper-singular integral representation. This integral equation is the following: where the terms S kij * and D kij * indicate the hyper-singular kernels, which are obtained from P lk * and u lk * , respectively (Portela et al. 1992;Hong and Chen 1988). f k represents the cosines of the normal exterior direction to the structural boundary.

Algebraic BEM equations for multi-domain analysis
To simulate the mechanical behaviour of solids composed by multi-domains, i.e., nonhomogeneous materials, the sub-region BEM technique has to be applied. In such approach, the body is divided into a finite amount of homogeneous sub-regions interconnected by interfaces.
BEM analyses involving singular and hyper-singular integral representations are performed using Eqs. (15) and (16). When multi-domains are considered, these equations have to be applied at each sub-domain individually. Then, the classical BEM system of algebraic equations is obtained for each sub-region, i, of the entire solid as follows: in which matrix H contains the integration kernels P lk * and S kij * , whereas matrix G contains the integration kernels u lk * and D kij * . Vectors U and P contain the displacement and traction values on the body boundary, respectively.
The global system of algebraic equations presented in Eq. (17) cannot be solved directly just by imposing the boundary conditions of the problem, because along the interfaces neither tractions nor displacements values are known. Therefore, it is required to enforce compatibility of displacements and equilibrium of forces along all interfaces. These conditions are written as follows: The compatibilities conditions, Eq. (18), coupled to the boundary conditions have to be imposed on the global system of equations. By performing a convenient change on the columns of H and G matrices, all known variables are placed at the right hand side of this algebraic system, whereas unknown variables, x, are placed at its left hand side . This system is presented as follows: Once f f g is the vector of known boundary values, the system is solved and the unknowns variables determined.

Integral formulations for viscoelastic materials based on the BEM
A brief review on the rheological models was presented in ''Review over the viscoelastic models''. In the present section, the differential equations that govern the viscoelastic models previously presented are applied in the determination of the BEM viscoelastic integral and algebraic representations. The formulation for the Maxwell viscoelastic model is presented in details. The formulations for the Kelvin-Voigt and Boltzmann viscoelastic models are shortly presented to avoid repetitive algebraic procedures.

Integral and algebraic representations for the maxwell viscoelastic model
To present the integral and algebraic representations for the Maxwell viscoelastic model, Eq. (2) has to be rewritten in terms of stresses as follows: Int J Adv Struct Eng (2017) 9:1-12 5 The last equation can be modified by assuming proportionality conditions between c and g. Such assumption, g = cE, was presented in (Mesquita and Coda 2001, 2002, 2003, which leads to: Equation (21) is also valid for 3D case. Therefore, it can be rewritten as follows: where C ijkl is the constitutive elastic tensor. The stress state in each material point is represented by Eq. (22) for Maxwell viscoelastic model. Then, such equation is applied on the equilibrium equation, r ij, j ?b i = 0, to obtain the BEM representation. For this purpose, the equilibrium equation has to be weighted by the fundamental solution of displacements, u * , which leads to: To support the algebraic developments that follow, the auxiliary variable below is defined: Bearing in mind that u ki,j * C ijmn = r kmn * , Eq. (24) assumes the following form: Equation (25) has to be integrated by parts, leading to: BEM formulations require a fundamental problem. In the context of the present study, the fundamental Kelvin problem is considered, i.e., r kmn,n * = -D(s, f)d km and r kl,l =b k D indicates the Dirac function. f and s are the field and source points, respectively. Thus, considering the Kelvin problem, the last equation assumes the following form: Due to a mathematical property of Dirac function, one writes the second integral of Eq. (27) as follows: As a result, Eq. (27) is rewritten in the following form: The result presented in Eq. (29) enables to write Eq. (23) as follows: To simplify the mechanical analysis, the body forces and its variations along time are assumed as nil. Therefore, Eq. (30) assumes the following form: The last equation is valuable for source points positioned at the body domain. Thus, convenient limits must to be carried out, as usual in BEM, to write Eq. (31) for source points positioned exclusively at the body boundary. Such limits lead to the equation presented below: The algebraic representation of Eq. (32) is performed by approximating displacements and traction by shape functions. High order polynomial functions may be adopted. Therefore, considering the approximation of these variables at the boundary, the algebraic representation for the Maxwell viscoelastic model is the following: The last equation is a differential equation in time domain, which is solved by forward finite differences technique. Therefore, a linear approximation is considered for evaluating the first time derivative terms. This procedure requires the division of the total analysis time into finite time steps, Dt. Bearing in mind that s is the actual time one has: The algebraic representation presented in Eq. (33) has to be rewritten considering the definitions introduced in Eq. (34). Thus, the final algebraic representation for the Maxwell viscoelastic model is obtained as follows: The last equation has to be solved for finite time increments to evaluate the unknown variables at the boundary. The stress state at the internal points is obtained after the determination of all unknown variables at the boundary. For this purpose, Eq. (31) has to be derived with respect to the source points, which results in the following: As presented in ''Review over the viscoelastic models'', the Maxwell viscoelastic model assumes that r. ij = r. ij e = C ijkl u. k,l . Therefore, Eq. (36) is rewritten as follows: where r Ã kil is the fundamental solution for stress (Portela et al. 1992;Hong and Chen 1988). As usual in BEM formulations, the last integral equation is evaluated by approximating displacements and tractions at the boundary by shape functions. This approximation enables obtaining the algebraic representation of Eq. (37), which is the following: To solve the differential equation above, finite difference technique is applied. As a result, the following algebraic representation is obtained: It is worth emphasizing that viscous and elastic stresses are equal in this viscoelastic model. Thus, Eq. (39) is sufficient for determining the stress state in materials governed by the Maxwell viscoelastic model.

Integral and algebraic representations for the Kelvin-Voigt viscoelastic model
To develop the algebraic and integral equations that represent the mechanical behaviour of materials governed by the Kelvin-Voigt viscoelastic model, mathematical procedures similar to the presented in the last sub-section have to be applied. To avoid the introduction of repetitive matter to that presented in the last sub-section, such procedures are omitted. Thus, the unknown variables at the boundary are evaluated by the following algebraic equation: The total stress state components at the internal points are achieved by the following algebraic representation: As presented in ''Review over the viscoelastic models'', the Kelvin-Voigt viscoelastic model considers the contribution of elastic and viscous portions on the total stress state. Therefore, the elastic stresses are evaluated as follows: On the other hand, the viscous stresses are determined by the following expression: Integral and algebraic representations for the Boltzmann viscoelastic model The integral and algebraic representations for the Boltzmann viscoelastic model are presented in expedite form, similarly, to perform in the last sub-section. The unknown values of tractions and displacements at the boundary of structures composed by materials governed by the Boltzmann viscoelastic model are obtained by the algebraic equation presented below: The total stress state components at internal points are evaluated by the following algebraic representation: As presented in ''Review over the viscoelastic models'', the total stress in this viscoelastic model is composed by Int J Adv Struct Eng (2017) 9:1-12 7 viscous and elastic portions. The elastic contribution is evaluated by the following equation: Consequently, the viscous contribution is determined as follows:

Examples
The presented viscoelastic BEM formulations are applied in the mechanical analysis of three cracked structures. In the first and second analyses, the results obtained by the BEM are compared with experimental and analytical responses provided by the literature. The last application concerns the structural analysis of a sandwich panel composed of nonhomogeneous viscoelastic materials including initial cracks. It is worth mentioning that nonhomogeneous viscoelastic cracked structures have not been properly analyzed in the literature. Therefore, the last application aims to be a benchmark for aiding the development of future viscoelastic numerical formulations.
In both analyses, the robustness and accuracy of the implemented viscoelastic BEM formulation are illustrated.

Griffith problem
The first example of this study concerns the structural analysis of the Griffith problem. Figure 2 illustrates the geometry and boundary conditions for this example, which involves a 2D structure in plane strain condition with a symmetric crack. The external load is composed by a remote tensile traction. The material was assumed as governed by the Boltzmann's model with the following values: E 1 = 22.5757 kgf/cm 2 , E 2 = 11 kgf/cm 2 , t = 0.0, c = 45.4545 days, P = 5.0 kgf/cm 2 , Dt = 1 day, total time = 300 days, and a = 1.5 cm.
The mechanical problem presented in Fig. 2 has an analytical solution. Such a solution associates the crack opening displacement (COD) values to the intensity of the external load by the relaxation Boltzmann's function, H(t), as follows: The numerical responses obtained by the BEM for COD are compared with the results provided by Eq. (48). The comparative results are presented in Fig. 3.
As illustrated in the results presented in Fig. 3, excellent agreement among the results is observed. The error between the response models is lower than 3% during the all time history analyzed.
The normal stress y was also determined by the viscoelastic BEM model, as presented in Fig. 4. According to this figure, the total normal stress y is constant along time, as expected, due to the equilibrium requirements. As a result, the elastic stress component grows along time, whereas the viscous stress part decreases.
The implemented viscoelastic BEM formulation enables for load phases conditions. Therefore, loading and unloading phases can be considered easily. To illustrate such a skill, the loading history presented in Fig. 5 was applied, in which one load and one unload phases are included.  The stress responses obtained by the BEM, accounting for the load presented in Fig. 5, are illustrated in Fig. 6.
As illustrated in Fig. 6, the total stress is constant along time, whereas its components have variation. Such behaviour is expected due to the equilibrium requirements. After 300 days of the loading process start, the structure is subjected to a complete unloading process. Thus, the structure recovers its undeformed configuration once the material is assumed as viscoelastic. Such a configuration is observed after 600 days of the loading process start. Then, the COD is nil after this time, as illustrated in Fig. 7.
It is worth mentioning that the unloading process causes an instantaneous recovering of strain. This behaviour is characteristic of the Boltzmann's model, once two spring elements are localized in a series sequence. In addition, the superposition of the Boltzmann's principle is observed in the results presented in the Figs. 6 and 7.
The structure presented in Fig. 2 may be analyzed considering another material mechanical behaviour model. For instance, the Maxwell's viscoelastic model may be applied. For this case, the problem has an analytical solution, which associates the COD to the load values as follows: The numerical viscoelastic BEM responses are compared with the results provided by the Eq. (49). The comparative results are illustrated in Fig. 8.  As observed in Fig. 8, extremely good agreement is achieved among the models applied, which indicates the relevance and accuracy of the implemented viscoelastic BEM formulation.

Central notched beam: three point bending test
The second example of the present study concerns the mechanical analysis of the notched beam illustrated in Fig. 9. This specimen was analyzed experimentally by (Zhou 1992), which utilized the following parameters values: a = 50 mm, b = 760 mm, L = 800 mm, e = 100 mm and h = 100 mm. The load is constant along time and applied at the middle span on the beam top. Its intensity is equal to 2.38 MPa.
The structural material is assumed as governed by the viscoelastic Boltzmann's model. The following parameters values are utilized: E 1 = 36 GPa, E 2 = 11 GPa, t = 0.2, c = 450.45 s, Dt = 1 s and total time = 500 s.
In the experimental study of (Zhou 1992) the values of crack mouth opening displacement (CMOD) along time were measured. The results obtained by (Zhou 1992) are compared with the responses achieved by BEM to validate the numerical BEM formulation. The experimental and numerical responses are illustrated in Fig. 10.
Excellent agreement is verified among the results. Particularly, until the structural fracture, which occurs in 500 s. During the structural fracture, the numerical BEM model loss convergence, as expected, once equilibrium is no longer verified.

Nonhomogeneous cracked panel
The last application of the present study refers to the structural analysis of the nonhomogeneous panel illustrated in Fig. 11. This sandwich panel is composed of three different materials, which have different mechanical behaviour. The mechanical coupling among each material is numerically performed by the sub-region BEM technique. The structure is clamped at its left boundary and at its right end a parabolic distributed load is applied. Three cracks F1, F2 and F3 are presented in the structure, which are positioned at the middle span of each material, as presented in Fig. 11.
The mechanical properties for each material (domain) that compose the structure are presented in Table 1.
To illustrate the mechanical time-effects, the crack mouth opening displacement (CMOD) was monitored along time for each crack simulated in the structure. The total time considered in this analysis is 300 days, which was simulated with time intervals of 1 day. The variation of CMOD along time is illustrated in Fig. 12.
According to Fig. 12, it is observed that the CMOD grows along time for all cracks studied. Such behaviour is observed independently of the mechanical material model adopted. Figure 13 presents the variation of the shear stress along time for a given point belonging to domain 3 (D3). The behaviour for the total shear stress and for its elastic and viscous parts is studied. The total shear stress is constant along time, as expected, due to the equilibrium requirements. The elastic   and viscous parts grow and decreases, respectively, to keep constant the total shear stress.
The results presented in this example may be utilized for validating future numerical formulations dedicated to the mechanical analysis of nonhomogeneous viscoelastic crack materials. Mechanical responses at the boundary and on the domain were presented.

Conclusion
The present study discussed a viscoelastic BEM formulation. Such formulation was applied for mechanical analyses of materials subjected to initial cracks. The numerical formulation is based on the dual BEM, in which singular and hyper-singular integral equations are applied to represent the mechanical behaviour of cracked bodies. This BEM approach is the most popular in the modelling of fracture problems and lead to accurate results. The mechanical material behaviour was modelled through the rheological models of Maxwell, Kelvin-Voigt and Boltzmann. The differential governing equations were presented and the algebraic BEM representations were achieved.
Fracture mechanics problems were analyzed considering bodies composed by either homogeneous or nonhomogeneous materials. The mechanical modelling of viscoelastic cracked bodies is still a challenge in engineering structures field. Especially, in the context of nonhomogeneous bodies where few advances were observed in the recent literature.
The presented BEM scheme was utilized to the mechanical analysis of fracture problems, where the numerical responses obtained were compared with analytical and experimental results available in the literature. Good agreement among the responses was observed for the first and second example presented in this study. Thus, the relevance, robustness and accuracy of the implemented formulation are demonstrated. Moreover, this formulation   is capable to represent loading phases, which enables the simulation of fatigue phenomenon considering time-dependent effects. Finally, the last example of this study presented the mechanical analysis of a nonhomogeneous viscoelastic structure, in which results from the boundary and the domain are illustrated. This application aimed to be a benchmark to aid the development and validation of future numerical formulations in this domain.