A temporally piecewise adaptive multiscale scaled boundary finite element 

method to solve two-dimensional heterogeneous viscoelastic problems 

ARTICLE INFO  

Keywords: 

Reduced order 

Heterogeneous viscoelastic materials 

Multiscale SBFEM 

Numerical base function 

Temporally piecewise adaptive algorithm 

ABSTRACT  

By combining the Multiscale Scale Boundary Finite Element Method (MsSBFEM) and the Temporally Piecewise 

Adaptive Algorithm (TPAA), a new numerical model is presented to reduce the solution scale of the two￾dimensional heterogeneous viscoelastic problems. Utilizing TPAA, a spatially and temporally related problem 

is transformed into a series of recursive spatial problems, which are solved by MsSBFEM. The solution scale can 

be effectively reduced by recourse of a bridge between small-scale and large-scale via Scaled Boundary Finite 

Element Method (SBFEM) based numerical base functions, and the solution accuracy can be improved only by 

increasing nodes of coarse elements without increasing any new node inside. By virtue of singular, polygon and 

Quadtree elements, SBFEM renders the proposed algorithm more efficient and convenient to tackle with the 

stress singularity, and to generate SBFE mesh. TPAA provides a measure to secure the temporal computational 

accuracy via an adaptive process when the step size varies. Numerical examples are provided to elucidate the 

effectiveness of proposed approaches, and satisfactory results are achieved at both the large and small scales. 

1. Introduction 

FEM based Direct Numerical Simulation (FE-DNS) for heterogeneous 

viscoelastic problems is often involved in large number of DOF required 

for discretization complying with heterogeneous materials distribution. 

For instance, in the DNS based effective evaluation on RVE for hetero￾geneous viscoelastic materials, the degree of freedom (DOF) of pixel/ 

voxel-based FE-DNS will scale 10^7 [1]. On the other hand, due to the 

materials heterogeneity, the FE mesh generation is challenged in dealing 

with irregular geometry of constitutes and transition between the areas 

with different mesh densities to meet the conforming conditions on the 

element/material interface [2]. In addition, when there exist cracks, the 

stress singularity becomes a difficulty in the heterogeneous analysis with 

a simple and efficient mean. Since viscoelastic analysis is time depen￾dent, and is usually conducted via step marching, the temporal solution 

accuracy needs to be secured for different step sizes which are usually 

hard to predict. Therefore, efficient numerical algorithms are definitely 

required with the comprehensive consideration of reduction of solution 

scale, mesh generation and stress singularity, and the temporal solution 

accuracy as well. However, in the context of FE-DNS based numerical 

analysis of heterogeneous viscoelastic problems, there seems few reports۔

directly related to this issue. 

The Multiscale Finite Element Method (MsFEM), proposed by Hou 

et al. [3,4], provides an effective mean to reduce the solution scale of FE 

analysis. The basic idea of the MsFEM is to construct a bridge between 

small scale and large-scale via numerical base functions, by which the 

DOFs at small-scale can be condensed via those at large scale, such that 

the solution scale can significantly be reduced. MsFEM has already been 

employed in various heterogeneous multiscale analysis related to heat 

transfer problems [5], elastic-plastic problems [6], and 

thermal-mechanical coupling problems [7]. Klimczak et al. [8] proposed 

a MsFEM based algorithm to reduce solution scale in modeling of het￾erogeneous viscoelastic materials, which was further enhanced via 

adaptive generation of coarse and fine meshes. However, this approach 

required an iterative computing at each time step. 

In recent years, by recourse of advantages of polygon element and 

quadtree mesh [9,10], the Scaled Boundary Finite Element Method 

(SBFEM) [11,12] frequently appears in FE form, instead of BEM, and has 

been efficiently used in the FE based multiscale heterogeneous analysis 

[13,14]. Praser et al. [15] presented a FE2 multiscale model for 

hyperelastic materials, in which quadtree SBFEM is employed in the 

numerical homogenization on RVE. Utilizing SBFEM with the idea of X. Wang at all.    

MsFEM, Egger [16] proposed a MsSBFEM to treat the two-dimensional 

linear elastic crack propagation problems, authors developed a poly￾gon and quadtree elements based MsSBFEM in the heterogeneous 

analysis for the steady state heat conduction problem [17]. The singular 

scaled boundary finite element inherits the advantage of SBFEM to 

tackle with the stress/heat flux singularity [18–21], and the polygon 

element with arbitrary number of edges and vertices provides conve￾nience to handle irregular geometric shape [9]. By combining with 

quadtree gridding techniques, the SBFE mesh generation becomes more 

efficient, and can even be directly conducted on the digital images of 

computing domains [22,23]. 

By integrating the advantages of the MsSBFEM and Temporally 

Piecewise Adaptive Algorithm (TPAA) that is able to secure a stable and 

high-fidelity temporal solution accuracy when the step size varies 

[24–28], a new reduced order approach, i.e., Multiscale Scaled Bound￾ary Finite Element Method (TPA-MsSBFEM), is put forward in the het￾erogeneous viscoelastic analysis. 

To some extent the presented work can be considered as an extension 

of our previous work [17], such as the extension of linear elastic MsFEM 

to elastic-plastic problems [29], thermal-mechanical coupling problems 

[30], and viscoelastic problems [8]. One of the key points or the major 

novelties of these extensions is to build a proper platform on which 

MsFEM in the linear elastic case can be implemented. In this sense, a 

main novelty of presented wok is to build a recursive SBFEM based 

framework, on which MsSBFEM can be implemented as in the linear 

elastic case, the evaluation of the coarse element stiffness matrices needs 

to generate only once via numerical base functions, and no iteration is 

required overall the time domain. Multiple benefits can be gained in 

terms of reduction of solution scale, convenience of dealing with the 

stress singularity and mesh generating with complex heterogeneous 

materials distribution, and the stable temporal solution accuracy as well. 

The paper is organized as follows. Section 2 gives the recursive 

governing and constitutive equations of viscoelastic problems; Section 3 

describes the construction process of numerical base functions and two 

kinds of SBFE gridding techniques; Section 4 addresses adaptive recur￾sive multiscale solutions at the small and large scales; Section 5 illus￾trates the effectiveness of proposed method via three numerical 

examples, and Section 6 provides conclusions and discussions. 

2. Recursive governing and constitutive equations 

2.1. Governing and recursive governing equations 

The governing equations of viscoelastic problem can be described by 

where {σ}, {ε} denote the vectors of stress and strain, respectively, {b} is 

the vector of the body force, {u} is the vector of displacement, and Ω 

stands for the domain of interest, and is enclosed by Γ =Γu∪Γσ. 

The boundary conditions are specified by

{u} = {u}on Γu (4)  

[n]{σ} = {p} on Γσ (5)  

where {u} and {p} are the vectors of prescribed displacements and 

tractions along the boundaries of Γu and Γσ, respectively, and [n] refers 

to the matrix of unit outward normal along Γσ. 

We divide time domain into a number of time intervals, the initial 

points and sizes of time intervals are defined by t1,t2,..., tk... andT1,T2,..., 

Tk..., respectively. At a discretized time interval, in order to describe the 

variation of variables more precisely, all variables are expanded in term 

of s 

where tk-1 and Tk represent the initial point and size of the k-th time 

interval, {σm}, {εm}, {bm}, {um}, {pm}, {um}, and {pm} denote the 

expanding coefficients of {σ(t)}, {ε(t)}, {b(t)}, {u(t)}, {u(t)} and 

{p(t)}, respectively. 

The conversion relationship between the differentiations respect to t 

and s is 

Substituting Eq. (6) - Eq. (12) into Eq. (1) - Eq. (5) and equating the 

power of two sides of the equation then yields 

[H]

T{

σm} + {bm} = 0 in Ω                                     (17)  

{εm} = [H]{um}                                             (18)  

[n]{σm} = {pm}on Γσ                                   (19)  

{um} = {um}on Γu                                        ( 20)

3.2. Polygon and quadtree gridding technique 

3.2.1. Polygon mesh 

The polygon mesh uses polygon elements with arbitrary number of 

sides and nodes on side [9], which brings convenience for characterizing 

complex geometries in the multiscale analysis, as the polygon mesh can 

be gridded by only one or very few SBFEs for the region with the same 

material property. Fig. 4 shows the polygon gridding for the heteroge￾neous plate which contains unit cells composed of two kinds materials 

represented by different colors, in which a SBFE mesh with only 5 

polygon SBFEs can be observed at the small-scale. In addition, it is more 

important that this kind of mesh provides a flexible way to increase 

nodes of coarse element without adding any fine mesh nodes. As shown 

in Fig. 4, on the basis of the commonly used four-node coarse element, a 

multi-node coarse element can be constructed by adding any number of 

coarse grid nodes (denoted as red color) at any position along the side of 

coarse grid. 

3.2.2. Quadtree mesh 

The quadtree gridding provides an effective tool of the mesh gen￾eration for the image-based analysis. In Fig. 5, each pixel in the image is 

represented as a square domain which is divided into a number of ele￾ments, and all the color densities of the pixels are recorded. If the dif￾ference between the maximum and minimum color intensity in an 

element is larger than the color threshold, the element is recursively 

divided further into four equal-sized elements until all the elements 

satisfy the criterion of homogeneity or reaches the minimum edge length 

[25]. When a 2:1 rule is used in the process of quadtree gridding in 

two-dimensional case, 16 possible element patterns have been provided 

for a convenient generation of SBFE stiffness matrix [25]. Besides, the 

‘hanging nodes’ produced in the match of adjacent elements in the 

quadtree mesh generation can be easily treated as nodes of a new 

polygon element, instead of regular ones, as shown in Fig. 5. 

Fig. 16. The curves ofstress with distance to crack tip in feature coarse element I 

color) will be used to test the influence of the type of coarse element. 

Firstly, the computational accuracy is verified by the reference so￾lutions provided by ABAQUS based on FE solution (6669 nodes) for 

entire structure shown in Fig. 8. Table 2 shows the variations of 

displacement at large-scale feature point 2 with different sizes of time 

steps, in which the error tolerance is β= 10-6, the material parameter is 

Case 2 and Model C is selected as the coarse element. Also, the 

displacement of feature point iii (shown in Fig. 7) in the coarse element 

II (shown in Fig. 6(b)) at small-scale is illustrated in Table 3. It is shown 

that the maximum eu

i is 0.95% in large-scale and 0.44% in small-scale, 

respectively, and the proposed method has a stable accuracy both in 

large and small scales when the time step increases from 0.001 s to 0.1 s. 

Fig. 9 shows the variation of displacement with time for the feature 

point 2 at large-scale and the feature point iii in the coarse element II at 

small-scale. It is demonstrated that the proposed algorithm can ensure 

the computational accuracy with different time steps both in large and 

small scales, while the solutions from ABAQUS change relatively larger 

as the time step increase from 0.001 s to 0.1 s. 

The variations of the number of recursion steps in time domain is 

shown in Fig. 10. As the β changes from 10− 6 to 10− 12 and the size of 

time step changes from 0.5 s to 0.1 s, the proposed algorithm can 

adaptively adjust the number of recursive steps according to the pre￾scribed error tolerance and the size of time step to maintain computa￾tional accuracy. 

Secondly, the influence of material heterogeneity for the result is 

investigated. Table 4 provides the comparison of proposed algorithm 

and reference solution on displacement ux at t= 10s and the corre￾sponding error indicator ψ for all the feature nodes at the large-scale in 

Fig. 6(b), respectively. The Model A coarse element is used. The results 

for displacements in large-scale with different cases of materials are 

shown in Table 1. It is found that as the ρ = EMat.1 2 /EMat.2 2 increases from 2 

to 100, the maximum eu

i of a single feature point increase from 2.10% to 

4.70%, and ψ increases from 0.0011 to 0.0081. With the increase of the 

material heterogeneity, in term of the increase of material ratio, the 

solution accuracy of the proposed algorithm decreases gradually. 

Thirdly, the influence of the coarse element with different number of 

coarse node is tested. Table 5 provides the comparison of the proposed 

algorithm and the reference solutions on displacement ux at t= 10s and 

the corresponding error indicator ψ in small-scale with corresponding 

feature points shown in Fig. 7 in the feature coarse elements I-III shown 

in Fig. 6(b). The Case 3 material parameter and all the three models for 

coarse element are used. It is shown that as the coarse nodes increase 

from 4 to 12, the maximum eu

i decreases from 5.28% to 2.48%, and the ψ 

(a) Influence of models for coarse element (b) Influence of size of time steps.