Convergence Study of Variational Space-Time Coupled Least-Squares Frameworks in Simulation of Wave Propagation in Viscoelastic Medium

ere are many novel applications of space-time decoupled least squares and Galerkin formulations that simulate wave propagation through di erent types of media. Numerical simulation of stress wave propagation through viscoelastic medium is commonly carried out using the space-time decoupled Galerkin weak form in site response problem, etc. In a recent investigation into accuracy of this formulation in simulating elastic wave propagation, it was shown that the di usive and dispersive errors are greatly reduced when space-time coupled least squares formulation is used instead in variational form. is paper investigates convergence characteristics of both formulations. To this end, two test cases, which are site response and impact models for viscoelastic medium, are employed, one dominated by dispersive and the other by di usive numerical error. Convergence studies reveal that, compared to the commonly used space-time decoupled Galerkin and the coupled least squares formulation has much lower numerical errors, higher rates of convergence, and ability to take far larger time increments in the evolution process. In solving such models, coecient matrices would require updating after each time step, a process that can be very costly. However large time steps allowed by cLs are expected to be a signicant feature in reducing the overall computational cost.


Introduction
Vibration energy dissipation, damping, is involved in many elds of mechanics problems. Mostly, reducing response amplitudes of exible or solid bodies is important to be considered for engineering subjects. Employment of the elasticity theory to simplify the analysis proves to be inconsistent with the accurate behavior of materials, because most engineering materials have much dependency on their intrinsic properties due to internal friction. In mechanical problem, investigation of damping has a main role in smart mechanical tools, response free eldout and, structures such as tall building and high way bridges [1]. A number of researchers investigated damping property e ect of rubber-sand mixture material underlain by the structures to reduce peak of acceleration, displacement, and shear stress at their bases [2][3][4]. e study in [5] rst introduced a famous approach called viscous damping that had been obtained from rheology science. is was competent idealization to describe locally vague material properties which dissipated vibration energy by means of their internal friction. Other idealization is mentioned that the damping matrix assumed to be consisted of linear combination of mass and sti eness matrices. In this model, undamped systems can be used for the analysis of damped systems e ortlessly.
Overall, damping is divided into two main classes: (1) damping in discrete systems including SDOF and MDOF systems and (2) damping in continuous systems. Usually, the rst is in conjunction with ordinary di erential equations (ODEs), extracted from dynamic equilibrium equations, and the second is related to partial di erential equations (PDEs), such as wave di erential equation. e authors of [4,6,7] worked on viscoelastic problems in different damping models.
ere have been studies on one-dimensional viscoelastic analysis of axial wave propagation through rod with various damping models. e authors of [8][9][10] demonstrated realistic behavior of viscous damping for analyzing actual dynamic systems.
From before to recently, wave propagation simulation considering energy dissipation has been a desirable and attractive topic to research of viscoelastic problem all over the world. ere are two major computational methods to analyze viscoelastic dynamic problems, with discrete or continuous systems; (a) time-domain methods, in which the most general computational approaches are included, and (b) frequency-domain methods, which are favorable for linear or equivalent linear problems, extensively computed by mean of finite difference, boundary elements, and finite element methods along with the isogeometric and meshless variations. e study in [11,12] investigated the earthquake response of surcharged soil layers using Hybrid Frequency Time Domain (HFTD) approach. ey used and developed transfer function method, which is categorized into frequency domain analyses, for viscoelastic soil medium on which equivalent mass is surcharged. e study in [13] compared responses 1D viscoelastic horizontally layered soil model and 2D one and observed no significant difference between them.
In time-domain analyses, the solution is known for all points of the spatial domain at a given time t 0 and the objective is to determine the solution at time t 0 + Δt. e formulation needed for this evolution can be cast in a space-time decoupled or a space-time coupled computational framework. e decoupled Galerkin formulation is widely utilized in study of wave equations, as they appear in many engineering disciplines. Here, coefficient matrices representing discrete spatial distribution of medium's property form a set of ordinary timedependent equations, which are then solved over a given time span using difference-based techniques such as Newmark-β or Wilson-θ methods [14]. Many improvements have been done to the decoupled formulations [15][16][17]. Even though the semidiscrete finite elements approach has led to significant improvements in simulation of elastic wave propagation, problems with high frequency loading and sharp temporal gradients still present a significant challenge.
Space-time coupled formulations have been successfully applied in studying time-dependent problems [18]. In this framework, finite element interpolation is performed over an m + 1-dimensional computational domain with m being the number of spatial dimensions. e study in [19] was the first to carry out study over the coupled domain. As demonstrated in [20], in the study of elastic wave propagation due to impact, the decoupled Galerkin yields unsatisfactory results and a space-time coupled discontinuous Galerkin scheme was introduced in order to reduce the dispersive error in this problem. e author of [21] in PhD thesis proposes a discontinuous scheme with least squares stabilizers; this method requires auxiliary variables and, like other discontinuous schemes, upwinding parameters have to be tuned to a particular problem at hand.
Instead of using discontinuous formulations and dealing with complications of proper upwinding scheme, the study in [22] has suggested a comprehensive framework for casting the problem in a space-time coupled least squares framework with higher continuity elements. is framework presents means of convergence through a combination of element length, order of interpolation, and global smoothness. In an error analysis of the wave equation, a quadratic rate of convergence for the wave equation was predicted by [23]. e rate is with respect to element characteristic length, with approximation defined over the coupled domain. However, in context of finite element, it was established that higher smoothness is needed to obtain such rates.
Reference [24] proposed new nonsymmetric variational formulations which are employed to parallelize computations on MIMD-parallel computer for viscoelastic problems based on the continuous and discontinuous Galerkin method. In the study therein, the three-parameter Malvern model described viscoelastic pattern. Promising results were achieved for discontinuous Galerkin method with respect to continuous Galerkin method.
Reference [25] developed temporal decoupling of discontinuous Galerkian space-time finite element method, which is formulated by [26], which is applied only to parabolic differential equation, heat diffusion equation. Continuous Galerkin form in space and discontinuous Galerkin form in time were adopted for second-order wave equations including elastodynamic problem with and without Kelvin-Voigt and Maxwell-Zener, viscoelasticity. Acceptable results were extracted for moderately high-order (up to degree 7) temporal and spatial-temporal approximation. eir method, decoupling procedure, produced a set of boundary value problems that need to be solved for each time interval. e most popular methods in the engineering practice are the finite and the spectral element methods. ey present known advantages (deal with complex geometries, material nonlinearities, etc.) and drawbacks (numerical damping and dispersion, spurious reflections at artificial truncation boundaries). Although various numerical strategies exist to limit spurious reflections (e.g., absorbing boundary conditions or boundary layers), the boundary element method (BEM) remains a very effective approach for dynamic problems in spatially-extended regions (idealized as unbounded), especially so since the advent of fast BEMs such as the fast multiple method (FMM) used in this work. [27].
is macroelement allows one to model soil-footing geometric (uplift) and material (soil plasticity) nonlinearities that are coupled through a stiffness degradation model. Footing uplift is introduced by a simple non-linear elastic model based on the concept of effective foundation width, whereas soil plasticity is treated by means of a bounding surface approach in which a vertical load mapping rule is implemented [28]. Performance criteria are generally displacement-based in the performance-based design approach. Quality requirements in the approach to performance-based design are typically based on displacement. Firstly, keeping this displacement within an acceptable limit ensures the main purpose of planning to have sufficient strength and stiffness. In multistorey commercial and residential structures, coupled walls are a typical type of shear wall. rough a fusion of the coupling beam's frame action and the wall pier's flexural action, a coupled shear wall system resists lateral forces [29]. Carbon textile is considered an alternate material due to its corrosive resistance property, high tensile strength, and being perfectly elastic. Prestressing is also the only realistic way to utilize fully ultra-high tensile strength in carbon textile material [30].
In a recent paper [31], it was shown that numerical errors can be significantly reduced if the problem is cast in a spacetime coupled least squares finite element framework. As an extension to that work, in this paper, convergence characteristics of damped wave equation are studied for the decoupled Galerkin and coupled least squares formulations. In this investigation, since numerical errors are a combination of dispersive and diffusive types, two test cases were designed: one of them is site response model and the other is impact model, each dominated by one of the two error types. Using exact solution, convergence through spatial and temporal refinements and the effects of interpolation and global smoothness on the computational process are studied here.
e remainder of this manuscript is organized as follows: e mathematical model describing viscoelastic wave propagation is presented first followed by the decoupled Galerkin and coupled least squares weak formulations. e succeeding section presents two test cases; one simulates impact and the other base motion. Convergence studies are presented and overall patterns are deduced. is is followed by a discussion on findings and concluding remarks. is study tried to investigate about the relations of variational space-time coupled least squares frameworks using wave propagation in viscoelastic medium, while other research works concentrated on some aspect of this concept.

Mathematical Model
Propagation of viscoelastic waves in a one dimensional homogeneous medium can be modeled with where E and η respectively are the elastic and viscous modulus; ρ denotes mass density; and is model (1) together with suitable boundary and initial conditions constitute the strong formulation of viscoelastic wave propagation in one dimension. e exact solution to (1) can be approximated by u h (x, t) as where superscript h signifies characteristic length in the finite element mesh over the spatial domain X. N is the row vector of the approximating shape functions and Φ represents the degrees of freedom vector 4. is approximation leads to the residual function which is subsequently utilized to develop variation-based integral statements that can be used to determine Φin (2). Next, two different integral-based computational frameworks, namely the space-time coupled least squares and the space-time decoupled Galerkin are generated for (3). Numerical experiments are then conducted to compare computational errors and convergence rates in each of the frameworks.
Viscoacoustic seismic modeling is much cheaper, but at the tradeoff of using incomplete physics. e algorithm contains two viscoacoustic forward modeling steps; the first is the same as the traditional viscoacoustic modeling, while the second propagation is generated using a residual error source, which is derived by comparing the viscoacoustic and viscoelastic wave equations in the form of stress-particle velocity formulations. e corrected P-wave particle velocities can be obtained by adding the wavefield from the second simulation step to the original (the first simulation step) viscoacoustic wavefield. Only P waves are modeled. e overall cost is about twice that of viscoacoustic modeling, but it is significantly less than a viscoelastic propagation, because there are fewer calculations, and we can use a coarser grid and larger time steps for the same accuracy [32].

Space-Time Coupled Least Squares Formulation.
Solutions from space-time coupled least squares formulation (cLs) are a linear combination of two dimensional shape functions formed over ( e least square error functional is defined as which is bounded, nonnegative, and quadratic by construction. First variation of J yields the least squares minimization statement. where Solution vector Φthat satisfies the minimization statement given in (5) is the best approximation to u that may be found in solution space spanned by basis N(x, t).
Note that the least squares weak form, (5), requires global continuity of C 1 over both space and time computational domains. In all computations carried out here, space-time elements are constructed from tensor product of two cubic ordered C 1 elements [33]. e representation of viscoelastic media in the time domain becomes more challenging with greater bandwidth of the propagating waves and number of travelled wavelengths. With the continuously increasing computational Advances in Civil Engineering power, more extreme parameter regimes become accessible, which requires the reassessment and improvement of the standard memory variable methods to implement attenuation in time-domain seismic wave propagation methods [34].

Space-Time Decoupled Galerkin Formulation.
Space-time decoupled Galerkin formulation (dGa) is the most commonly used computational framework in study of viscoelastic vibrations. e space-time decoupled Galerkin weak form of (1) can be derived by letting u(x, t) � N(x).Φ(t), leading to a set of time dependent ordinary differential equations, i.e., where M � X ρN T NA dx is the mass matrix, C � X ηdN T / dxdN/dxA dx is the damping matrix, and K � X EdN T / dxdN/dxA dx is the stiffness matrix. e highest order of derivative appearing in the definition of each coefficient matrix is unity; hence, utilization of linear C 0 finite elements provides the minimum global continuity requirement in dGa computations. e set of linear ordinary differential equations in (7) is evolved through time using the unconditionally-stable Newmark-β integration. e method assumes linear acceleration profile; i.e., the displacement is cubic in time. is order of interpolation in time is equivalent to the cubic distribution employed in approximation of displacement in temporal direction. However, in the cLs framework, continuity of displacement and velocity are strictly enforced, whereas in the Newmark method continuity of displacement and velocity are controlled by a shift in continuity of velocity defined as averaged acceleration across the time-increment.

Numerical Experiments
In studying the computational characteristics of dGa and cLs, dynamic response predictions made by the weak forms (7) and (5) of (1) are compared to exact solution. Test cases are designed to allow for separate studies in connection with dispersive and diffusive numerical errors. Errors in displacement and stress distributions over the first few cycles of wave reflection and propagation are measured; and are utilized to identify convergence characteristics of each computational framework.

Evaluation Tools.
For evaluating the merits and drawbacks of each formulation, means of convergence and measurement of numerical error are defined as follows.

Error Measurement.
Since the exact solutions are available in all case studies, relative error is measured using

h-Convergence.
In any finite element computational process, errors can be reduced by increasing the number of approximating elements (mesh refinement), i.e., h-Convergence. In studies conducted here uniform spatial meshing is employed, i.e., a uniform mesh of n elements is defined over the spatial domain of length L with h � L/n. In all dGa studies, n is taken from the set N Ga and in cLs studies from the set N Ls , where Note that the sequence uses a growth factor of 2. Hence, the characteristic length h is halved in each refinement.

R-Convergence.
It is also possible to reduce numerical error using smaller time increments, or R-refinement. Here, time increment Δt is kept constant and it is set equal to a fraction of the reference time increment Δt; i.e., Δt � Δt/R; where, Δt � h/c and c � ��� E/ρ is the axial stress wave propagation speed in domain Ω with zero damping. Defining mesh speed as the ratio of characteristic length to time-increment's size, i.e., V � h/Δt, R may be regarded as the ratio of mesh speed to the stress wave speed, i.e., R � V/c.
In dGa studies, R-convergence is applied by halving time increment's duration, starting with R � 1; which in turn doubles the number of time steps required for evolving the solution over the specified time span, T. In cLs studies values of R < 1 are also considered for the sake of speeding the computational process.
Time increment factor, R, is selected from sets R Ga and R Ls in dGaand cLs convergence studies, respectively.
Note that R Ga sequence has a growth factor of 2; i.e., in dGa R-convergence studies time increment is halved in each refinement, with Δt � Δt/32 being the smallest time increment. Also note that R Ls set contains values less than unity as well. For R < 1 mesh speed is less than stress wave's speed, i.e., V < c. is allows for larger time increments to be taken; which translates to fewer number of evolution steps needed to reach a given final time.

p-Convergence.
Minimum polynomial order required by dGa weak form (7) is unity with C 0 global continuity; and for cLs weak form (5), a bicubic polynomial with C 1 continuity in both space and time is the minimum requirement. Possibility of reducing the approximation error through polynomial-order increase or p-convergence is also investigated here. 4 Advances in Civil Engineering Finite elements used here are Lagrange-based hierarchical elements [33].

τ-Convergence.
In particular, since the coupled formulation allows for temporal discretization as well as the usual spatial discretization over each time slab (space-time domain over a time increment), possibility of reducing the error through temporal meshing or τ-convergence is also investigated for the least squares formulation. Here, the computational domain is meshed over the time increment, where three cases are considered by dividing Δt to 1, 2, 4 { } divisions.

Test Models.
In general, numerical errors have dispersive and diffusive (dissipative) characteristics [31]. Each of the two test models employed here exhibits dominance of one type of numerical error. e first test case, Impact (Imp.), simulates the response of a relatively stiff medium to a constant load, where dispersive numerical error is dominant. Diffusive error is dominant in the second test case, Base-motion (Bms.), which simulates response of a relatively flexible medium to an imposed harmonic displacement boundary condition. Relative stiffness in each case refers to stress wave speed; which is 1000m/s for Imp. and it is 60m/s for base motion case. Errors are measured using the computed displacement u or stress σ � Eϵ � Ezu/zx where ϵ is the strain field.

Impact: Model Specifications.
e medium considered is a rod of unit area and unit length. Mass density and elastic modulus are ρ � 1000 kg/m 3 and E � 1 GPa, respectively; hence, the wave speed in this model c � 1000 m/s. e rod is fixed at x � 0 and loaded with point load P � 1 at Closed form solution to (1) based on conditions stated here will be where ϵ is the applied strain at the loaded end and

Site Response Model: Model Specifications.
In seismic numerical simulations of wave propagation, it is very important for us to consider surface topography and attenuation, which both have large effects (e.g., wave diffractions, conversion, amplitude/phase change) on seismic imaging and inversion. An irregular free surface provides significant information for interpreting the characteristics of seismic wave propagation in areas with rugged or rapidly varying topography, and viscoelastic media are a better representation of the earth's properties than acoustic/elastic media [35].
In this model, L � 10m, elastic modulus E � 7.5 MPa, and weight density ρg � 20kN/m 3 , where g � 9.81m/s 2 . e model is at rest prior to imposition of harmonic displacement at x � 0. Boundary and initial conditions can be stated as, u(0, t) � A sin (ωt) and zu/zx(L, t) � u(x, 0) � zu/zt (x, 0) � 0. e displacement condition u(0, t) has amplitude of A � 0.01 m and frequency of ω � 2πf, with f � 10 Hz. Furthermore, since C 1 finite elements are employed in studying (5), the stress free boundary condition zu/zx| (x�L,t) � 0 is imposed explicitly; a feature that is not available in computations based on the C 0 -dGa formulation. A closed form solution of (1) subjected to boundary and initial conditions stated here, may be found through separation of variables technique, yielding where α � L 2 ρω and the remaining parameters are as defined in (13).

Damping Ratio.
All studies are carried out by considering damping ratio ξ � η/η c , where critical damping for a homogeneous axial rod is Advances in Civil Engineering Values of ξ considered here are ξ ∈ 0.01, 0.1, 0.5, 1.0 { }.

Solution Characteristics and Profiles.
In definition of the impact problem (Imp.), a constant force is applied at the free end of a relatively stiff rod, producing a constant-stress wave which reflects from the fixed boundary and interacts with the incoming imposed wave. e result is an oscillating rectangular wave with sharp temporal gradient, the sharpness of which depends on the amount physical damping present in the system. Numerical simulation of undamped impact case [31] shows that the dominant error is dispersive in nature and it is most significant at the fixed boundary, as shown in Figure 1(a). Fictitious oscillations observed here decrease naturally with increase in system damping; i.e., depending on amount of damping present in the system, the square profile of the stress wave will assume smoother corners, thereby reducing the oscillatory overshoots. In definition of the base motion (Bms.) test case, medium is relatively flexible and low frequency harmonic displacement is imposed at the base of a long rod. e overall motion is harmonic with smooth gradients in time. In this setup, the dominant numerical error is of diffusive nature; and it increases with increase in system's damping. Figure 1 shows the diffused profile of computed displacement from a study using 384-element mesh with R � 1.
In general, numerical error from any computational model is a combination of both diffusive and dispersive types. Convergence characteristics and computational time of both formulations for each dominant error type are examined next.

h-and R-Convergence.
As noted earlier, the dominant numerical error is of dispersive type in case of impact and of diffusive type in case of base-motion disturbance. Numerical studies indicate that regardless of the framework employed, increasing damping has contradictory effects on dispersive and diffusive errors. For impact loading, note from Figure 2 that, given the same computational resources, errors are much smaller in cases that have higher damping values.
Note that physical viscosity reduces dispersive error; this is because the sharp corners of the rectangular pulse in undamped case, shown in Figure 1(a), which are responsible for solution dispersion, are smoothened in direct proportion to physical damping present in the system, hence reducing the dispersion. On the other hand, since viscosity is diffusive in nature, it increases the numerical error of diffusive type; i.e., given the same computational resources, in the case of base-motion disturbance, where diffusive-type numerical error is dominant, increase in damping increases the computational error ( Figure 3).
Comparing cLs and dGa h-curves in sub-figures of Figures 2 and 3 shows a difference of at least an order of magnitude in numerical error. e R-convergence curves demonstrate that dGa solutions can be improved upon when smaller time increments are employed; however, as can be noted from Figures 2(a) and 3(a) for dGa:n12 curves, insufficient spatial meshing, and, as will be shown later, insufficient global continuity inhibit the R-convergence process.

h-and R-Convergence Rates.
Effectiveness of h and R refinements can be assessed by measuring the rate at which such refinements reduce the computational error. h-convergence rates based on mesh refinements defined in (8) and (9) are computed by comparing the ratio of errors from two consecutive refinements to the ratio of their corresponding degrees of freedom. Similarly, in computing R-convergence rates, in accordance to (10) and (11), the ratio of numerical errors from two R-refinements is compared to the ratio of corresponding Rvalues. It is found that except for the computed stresses in Bms. case, rate of h-convergence for all other cases is at most linear in dGa and at most quadratic in cLs frameworks; i.e., by halving h the accuracy will at most double in dGa and quadruple in cLs computations. Furthermore, R-convergence rate in dGa studies is also at most linear; i.e., halving time increment's size will at  Uniform mesh re nements and time increment reduction seem to yield at most a linear convergence rate in both dGa and cLs studies. A closer look at the solution reveals the error's source in each case ( Figure 4).
As can be seen from Figure 4(a), the error in cLs computations is con ned to the starting cycle, which disappears quickly as time passes. is error is due to inconsistency between the zero initial velocity and the nonzero velocity condition at time zero coming from the imposition of displacement function A sin (ωt) at x 0. Since the error is spatially and temporally con ned, uniform space-time mesh re nement is not the best approach for removing this error; a simple mesh grading that re nes the region close to the disturbed boundary resolves the issue, as shown in Figure 5(a). Considering dGa test cases presented in Figure 3, the nonconvergent h-re nement at R 1 would have actually diverged if the unconditionally stable Newmark-β method was not employed for evolving the solution. Source of this error in dGa simulation is due to lack of boundary conditions at the free end. Figure 4(b) displays the violation of free-stress condition at x L; which increases with mesh re nement. Error plots also indicate that using higher values of R, i.e., smaller time increments along with su cient spatial meshing reduce the errors signi cantly, as shown in Figure 3. However, the ctitious strain pro les shown in Figure 4(b) can be suppressed if similar to cLs computations; C 1 nite elements are employed in dGa study of this ill posed problem. Figure 5(b) shows the h-and Rconvergence curves from dGa:C 1 computations and hconvergence curves from dGa:C 0 and cLs:C 1 computations. Among test cases considered here, Bms. has shown the slowest rates of convergence, especially in cLs computations. However, comparing Figures 3(b) and 5(b) clearly shows that imposition of stress-free condition yields convergence rates that are consistent with other dGa studies.         Advances in Civil Engineering Furthermore, using C 1 elements, the errors are pushed down signi cantly through R-re nement.

CPU Usage.
e two formulations of cLs and dGa studied here pass through di erent sets of computational steps; a meaningful comparison of the two would be accuracy versus cost; with cost being the CPU time used in each study. Furthermore, since highest errors in Imp. studies occur at lower damping and in Bms. studies at higher damping values, two representative graphs, one for Imp. and one for Bms., are presented in Figure 6. Each dGa curve shows the relative computational time spent in each Rconvergence study and the cLs h-convergence study.
It can be noted that cLs produces least error for least amount of CPU used, regardless of the error type.

p-Convergence.
Minimum polynomial order required by dGa formulation is one and for cLs, it is three. Increasing p is known to increase computational accuracy in boundary value problems. In studies conducted here, increasing p in spatial interpolation is found to a ect R-convergence of stress away from the re ecting boundary in dGa case, which could also be realized at lower number of degrees of freedom if C 1 elements are used.
Comparison of Figures 7(a) and 7(b) to Figure 2(a) shows that error in stresses computed away from the re ecting boundary cannot be reduced through R-re nement; i.e., reducing time increment's size does not reduce the error. However, this issue can be remedied by increasing the order of polynomial, p, from one to three, as shown in Figure 8(a) or simply using C 1 nite elements at p 3, which requires less Dofs for the same accuracy, as shown in Figure 8(b).  In cLs studies, increasing p in space does not change the error present in computations, regardless of space-time meshing used. However, increasing p in time's direction can make the evolution divergent. e issue is best demonstrated by Figure 9. Here, for R 0.1, i.e., large time increment of Δt 10 Δt, the cLs:h-convergence curve shows hyperconvergence, whereas increasing p in either direction (x or t) has no bene t to the process and in fact can cause divergence in case of crude spatial meshing when polynomial's order is increased in time.

R-and τ-Convergence in cLs Framework.
In the spacetime least squares process duration of each time increment is equal to the space-time slab's size in temporal direction. Duration of the time increment Δt is controlled by R, i.e., Δt Δt/R. When R 1 mesh speed equals the undamped wave speed; and in dGa studies, solutions were found to improve for R > 1, i.e., ner meshing in time. Obviously when R < 1, the error increases and h-convergence in dGa would have higher error compared to R 1 curve. However, the coupled nature of space and time in cLs formulation allows for hyper-h-convergence to take place for large time steps.
Furthermore, since over each time increment the spacetime slab can be meshed in the temporal direction, a series of studies into bene ts of temporal meshing are also presented. Here, Δt is meshed uniformly to study τ-convergence in cLs formulation.

R-Convergence.
Based on R-values listed in (11), cLs:Rconvergence for both test cases was investigated. Figure 10

10
Advances in Civil Engineering indicates that for values of R > 1, i.e., smaller time increments, little if any improvement is seen for the extra computational e ort spent in stepping through the extra number of time increments. is is in contrast to R-convergence in the dGa-framework, which shows linear rate of convergence in most cases. An important point to note here is that errors in computations with values of R close to one are almost the same as those computed using R 1. Also, note that for R 0.1 signi cant error exists when crude spatial-mesh is employed, however, h-convergence shows hyper-rates that reduce the error quickly. Convergent process for R < 1 is a signi cant characteristic of cLs framework since it reduces the computational cost signi cantly by taking large time increments. Convergence for R < 1 is not necessarily monotonic, as shown in Figures 10(a) and 10(c), but for the signi cantly low value of R 0.1h-convergence is monotonic and hyper and requires ten times less number of evolution steps as compared to Dof-equivalent R 1 computations.

τ-Convergence.
Over the space-time domain, h-renement means increasing the number of divisions in spatial direction. Similarly, τ-re nement implies an increase in the number of divisions in temporal direction over a given time slab with Δt duration. ree uniform temporal divisions of τ ∈ Δt/ 1, 2, 4 { } were considered for τ-convergence study. As can be observed from Figure 11(a), τ-re nement does not improve the solution if R 1. For R < 1, the τ-re nement shows hyper-convergence in improving a crude solution, as shown in Figure 11(b). However, note that τ-re nement applied to spatial meshes with R 0.1 cannot improve the solution beyond the accuracy o ered by R 1 for the same spatial meshing. erefore, τ-re nement does not o er any advantage in cLs framework; in fact, far more resources would be needed when more than one division is used in temporal direction of space-time slab.

Discussion
Accurate simulation wave propagation through media is an important issue. Test cases simulating impact (Imp.) and base motion (Bms.), with respective dominance in dispersive and di usive error types, were employed to investigate convergence characteristics of dGa and cLs weak formulations of damped wave (1). To this end, standard re nements were applied and numerical errors were plotted against the degrees of freedom (Dofs) used by re ned models. Furthermore, since dGa and cLs pass through completely different sets of computational steps, graphs of error versus total CPU time were also generated. Based on these studies, several general remarks can be made with regard to accuracy and convergence characteristics of cLs and dGa computational frameworks.
Depending on the numerical error type, being dispersive or di usive, mechanical damping present in the system a ects computational accuracy di erently. Dispersive errors are higher for lower damping values and di usive errors are higher for higher damping values. is characteristic is true for both dGa and cLs.
e Bms. problem is ill posed in dGa-C 0 framework as shown in Figures 3 and 4(b); employment of C 1 nite elements removes the issue and increases the accuracy in dGa-R convergence signi cantly as shown in Figures 5(b) and 8. Rate of h-convergence is at most linear for dGa and at most quadratic for cLs. A quadratic rate predicted by [23] can not be realized in dGa because of its decoupled constitution. Rate of R-convergence is at most linear for dGa and is hyper h-convergence for small values of R in cLs framework.

Conclusions
Even though space-time coupled framework has been promoted by different researchers, the extra dimension of the computational domain makes it unappealing and hence not commonly utilized in study of real world problems. e focus of this research has been to investigate the computational characteristics of cLs and compare them to those of dGa for viscoelastic wave (1). e first is the space-time coupled least squares (cLs) formulation, which was recently employed for solving the undamped wave equation [31]. e second is the widely utilized space-time decoupled Galerkin (dGa) formulation along with the unconditionally stable Newmark-β method. Comparisons were made based on accuracy versus Dofs as well as accuracy versus CPU time spent. Two test cases, each susceptible to a particular type of numerical error, i.e., dispersive or diffusive, were studied and some general conclusions on computational characteristics of the methods were made.
R-convergence in cLs has little or no effect for R values around one and R > 1; however, for small values, e.g., for R � 0.1 for which the time increments are ten times larger, h-convergence is of high rate and stable Figure 10. It was established that for R � 1 numerical errors from dGa computations are in general one to two orders of magnitude larger than error in cLs results for the same number of Dofs, as shown in Figures 2 and 3. Results from dGa studies require significantly smaller time increments, e.g., R � 32, yields comparable accuracy to cLs at R � 1. Convergence rate is not the same for all spatial points; this was seen in Imp. problem where increasing R has no effect in lowering error in stress profile predicted by dGa at points away from the reflecting boundary, as shown in Figure 7. p-convergence in dGa and cLs has little to no effect on numerical accuracy, as shown in Figure 8(a). However, it has significant effect in dGa framework where R-convergence, for points away from the boundaries, was shown to be inhibited when p � 1. In the cLs framework, increasing p from its minimum value of three in the spatial direction has no effect on accuracy. However, increasing p in temporal direction makes the evolutionary process divergent as shown in Figure 9. τ-convergence in cLs, where space-time slabs are refined in temporal direction, shows no improvement on computational accuracy as shown in Figures 11(a) and 11(b). CPU time used by each framework while improving computational accuracy was constructed for cases with highest numerical error. Comparison shows cLs requires less time than dGa to reach a certain level of accuracy.
It was found that mechanisms which allow for refinements over the coupled domain in cLs are not necessarily beneficial to the process; e.g., τ refinement where time increment Δt is meshed was shown to be counterproductive as shown in Figure 11. Also contrary to dGa, where reducing Δt (or increasing R) can reduce the error at a linear rate, cLs shows negligible gain in accuracy for R > 1, as shown in Figure 10, hence rendering the process costly.
Convergence rates were found to be at most linear in dGa formulation from the h-and R-convergence studies.
is rate does not change with higher p-level or higher global continuity; however, employment of C 1 elements was found to be essential in applying free end condition in base-motion test case. Except for a meshing issue in study of Bms. case, hconvergence rates in cLs studies were at most quadratic for R � 1. While taking smaller time increments, i.e., R > 1 was found to be without merit, large time steps, in particular R � 0.1 was found to be stable with hyper-h-convergence allowing for reasonably refined spatial mesh evolving at R � 0.1 to have comparable accuracy to the 10 times slower evolution of the same spatial mesh at R � 1. However, it should be noted that while evolution of the solution over large time increments of R � 0.1 appears to be stable for cases considered here, it would most likely fail if discontinuities exist, e.g., multimaterial interface. Coefficient matrices of dGa formulation are constructed over the one-dimensional domain x; the resulting system of equations is then evolved in time according to Δt. e cLs coefficient matrix is formed over the two-dimensional domain (x, t) which requires relatively more computational resources than dGa. However, it was observed that the extra demand by cLs is offset by the need of dGa for very small time steps; i.e., for comparable accuracy Δt cLs ≫ Δt dG a , i.e., dGa requires time steps with more orders of magnitude compared to cLs to yield the same accuracy. It should be noted that, for very small time increments, round-off errors can become a significant source of error. In continuation of this work, nonlinear models are considered.

Data Availability
Requests for access to these data should be made to the corresponding author at amin.saffarian@eng.uk.ac.ir.