Abstract
There are many novel applications of space-time decoupled least squares and Galerkin formulations that simulate wave propagation through different 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 diffusive and dispersive errors are greatly reduced when space-time coupled least squares formulation is used instead in variational form. This 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 diffusive 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, coefficient 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 significant feature in reducing the overall computational cost.
1. Introduction
Vibration energy dissipation, damping, is involved in many fields of mechanics problems. Mostly, reducing response amplitudes of flexible 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 fieldout and, structures such as tall building and high way bridges [1]. A number of researchers investigated damping property effect of rubber-sand mixture material underlain by the structures to reduce peak of acceleration, displacement, and shear stress at their bases [2–4].
The study in [5] first introduced a famous approach called viscous damping that had been obtained from rheology science. This 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 stiffeness matrices. In this model, undamped systems can be used for the analysis of damped systems effortlessly.
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 first is in conjunction with ordinary differential equations (ODEs), extracted from dynamic equilibrium equations, and the second is related to partial differential equations (PDEs), such as wave differential equation. The authors of [4, 6, 7] worked on viscoelastic problems in different damping models.
There have been studies on one-dimensional viscoelastic analysis of axial wave propagation through rod with various damping models. The authors of [8–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. There 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. The study in [11, 12] investigated the earthquake response of surcharged soil layers using Hybrid Frequency Time Domain (HFTD) approach. They used and developed transfer function method, which is categorized into frequency domain analyses, for viscoelastic soil medium on which equivalent mass is surcharged. The 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 and the objective is to determine the solution at time . The formulation needed for this evolution can be cast in a space-time decoupled or a space-time coupled computational framework. The 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 time-dependent 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–17]. Even though the semi-discrete 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 -dimensional computational domain with being the number of spatial dimensions. The 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. The 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. This 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]. The 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. Their method, decoupling procedure, produced a set of boundary value problems that need to be solved for each time interval.
The most popular methods in the engineering practice are the finite and the spectral element methods. They 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].
This 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. Through 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 space-time 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.
The remainder of this manuscript is organized as follows: The mathematical model describing viscoelastic wave propagation is presented first followed by the decoupled Galerkin and coupled least squares weak formulations. The succeeding section presents two test cases; one simulates impact and the other base motion. Convergence studies are presented and overall patterns are deduced. This is followed by a discussion on findings and concluding remarks. This 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.
2. Mathematical Model
Propagation of viscoelastic waves in a one dimensional homogeneous medium can be modeled withwhere and respectively are the elastic and viscous modulus; denotes mass density; and is the displacement field over space-time domain defined by intervals and .
This model (1) together with suitable boundary and initial conditions constitute the strong formulation of viscoelastic wave propagation in one dimension.
The exact solution to (1) can be approximated by aswhere superscript signifies characteristic length in the finite element mesh over the spatial domain . is the row vector of the approximating shape functions and represents the degrees of freedom vector 4. This approximation leads to the residual functionwhich 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. The 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. The 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. The 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].
2.1. 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 space; i.e. . The least square error functional is defined aswhich is bounded, nonnegative, and quadratic by construction. First variation of yields the least squares minimization statement.where
Solution vector that satisfies the minimization statement given in (5) is the best approximation to that may be found in solution space spanned by basis .
Note that the least squares weak form, (5), requires global continuity of 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 elements [33].
The 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 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].
2.2. Space-Time Decoupled Galerkin Formulation
Space-time decoupled Galerkin formulation (dGa) is the most commonly used computational framework in study of viscoelastic vibrations. The space-time decoupled Galerkin weak form of (1) can be derived by letting , leading to a set of time dependent ordinary differential equations, i.e.,where is the mass matrix, is the damping matrix, and is the stiffness matrix.
The highest order of derivative appearing in the definition of each coefficient matrix is unity; hence, utilization of linear finite elements provides the minimum global continuity requirement in dGa computations.
The set of linear ordinary differential equations in (7) is evolved through time using the unconditionally-stable Newmark- integration. The method assumes linear acceleration profile; i.e., the displacement is cubic in time. This 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.
3. 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.
3.1. Evaluation Tools
For evaluating the merits and drawbacks of each formulation, means of convergence and measurement of numerical error are defined as follows.
3.1.1. Error Measurement
Since the exact solutions are available in all case studies, relative error is measured using where is the weighted average difference and is the weighted average exact value. Here, , is the approximate value, and is the exact value at the time point of time segment .
3.1.2. 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 elements is defined over the spatial domain of length with . In all dGa studies, is taken from the set and in cLs studies from the set , where
Note that the sequence uses a growth factor of 2. Hence, the characteristic length h is halved in each refinement.
3.1.3. R-Convergence
It is also possible to reduce numerical error using smaller time increments, or R-refinement. Here, time increment is kept constant and it is set equal to a fraction of the reference time increment ; i.e., ; where, and 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., , R may be regarded as the ratio of mesh speed to the stress wave speed, i.e., .
In dGa studies, R-convergence is applied by halving time increment’s duration, starting with ; which in turn doubles the number of time steps required for evolving the solution over the specified time span, . In cLs studies values of are also considered for the sake of speeding the computational process.
Time increment factor, R, is selected from sets and in dGaand cLs convergence studies, respectively.
Note that sequence has a growth factor of 2; i.e., in dGa R-convergence studies time increment is halved in each refinement, with being the smallest time increment. Also note that set contains values less than unity as well. For mesh speed is less than stress wave’s speed, i.e., . This allows for larger time increments to be taken; which translates to fewer number of evolution steps needed to reach a given final time.
3.1.4. p-Convergence
Minimum polynomial order required by dGa weak form (7) is unity with global continuity; and for cLs weak form (5), a bicubic polynomial with 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. Finite elements used here are Lagrange-based hierarchical elements [33].
3.1.5. -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 to divisions.
3.2. 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. The 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 for Imp. and it is for base motion case. Errors are measured using the computed displacement or stress where is the strain field.
3.2.1. Impact: Model Specifications
The medium considered is a rod of unit area and unit length. Mass density and elastic modulus are and , respectively; hence, the wave speed in this model . The rod is fixed at and loaded with point load at for .
Closed form solution to (1) based on conditions stated here will bewhere is the applied strain at the loaded end and
3.2.2. 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, , elastic modulus , and weight density , where . The model is at rest prior to imposition of harmonic displacement at . Boundary and initial conditions can be stated as, and . The displacement condition has amplitude of and frequency of , with Furthermore, since finite elements are employed in studying (5), the stress free boundary condition is imposed explicitly; a feature that is not available in computations based on the -dGa formulation. A closed form solution of (1) subjected to boundary and initial conditions stated here, may be found through separation of variables technique, yieldingwhere and the remaining parameters are as defined in (13).
3.2.3. Damping Ratio
All studies are carried out by considering damping ratio where critical damping for a homogeneous axial rod is
Values of considered here are .
3.2.4. 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. The 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. The 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 .
(a)
(b)
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.
3.3. 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.
(a)
(b)
(c)
(d)
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).
(a)
(b)
(c)
(d)
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. The 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.
3.3.1. 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 R-values. 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 most double the accuracy. R-convergence in cLs framework will be discussed separately. For all Imp. studies presented in Figure 2 except the critical damping case (2d) the convergence rate varies between 3.3 and 4.4 for cLs and between 1.3 and 2 for dGa cases. Errors in Bms. studies are high and difficult to reduce. Uniform mesh refinements 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).
(a)
(b)
As can be seen from Figure 4(a), the error in cLs computations is confined to the starting cycle, which disappears quickly as time passes. This 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 at . Since the error is spatially and temporally confined, uniform space-time mesh refinement is not the best approach for removing this error; a simple mesh grading that refines 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-refinement at 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 ; which increases with mesh refinement. Error plots also indicate that using higher values of , i.e., smaller time increments along with sufficient spatial meshing reduce the errors significantly, as shown in Figure 3. However, the fictitious strain profiles shown in Figure 4(b) can be suppressed if similar to cLs computations; finite elements are employed in dGa study of this ill posed problem. Figure 5(b) shows the h- and R-convergence curves from dGa: computations and h-convergence curves from dGa: and cLs: 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. Furthermore, using elements, the errors are pushed down significantly through R-refinement.
(a)
(b)
3.3.2. CPU Usage
The two formulations of cLs and dGa studied here pass through different 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 R-convergence study and the cLs h-convergence study.
(a)
(b)
It can be noted that cLs produces least error for least amount of CPU used, regardless of the error type.
3.4. 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 affect R-convergence of stress away from the reflecting boundary in dGa case, which could also be realized at lower number of degrees of freedom if elements are used.
Comparison of Figures 7(a) and 7(b) to Figure 2(a) shows that error in stresses computed away from the reflecting boundary cannot be reduced through R-refinement; 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 finite elements at p = 3, which requires less Dofs for the same accuracy, as shown in Figure 8(b).
(a)
(b)
(a)
(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. The issue is best demonstrated by Figure 9. Here, for , i.e., large time increment of , the cLs:h-convergence curve shows hyperconvergence, whereas increasing p in either direction ( or ) has no benefit to the process and in fact can cause divergence in case of crude spatial meshing when polynomial’s order is increased in time.
3.5. R- and -Convergence in cLs Framework
In the space-time 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 is controlled by R, i.e., . When mesh speed equals the undamped wave speed; and in dGa studies, solutions were found to improve for , i.e., finer meshing in time. Obviously when , the error increases and h-convergence in dGa would have higher error compared to 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 space-time slab can be meshed in the temporal direction, a series of studies into benefits of temporal meshing are also presented. Here, is meshed uniformly to study -convergence in cLs formulation.
3.5.1. R-Convergence
Based on R-values listed in (11), cLs:R-convergence for both test cases was investigated. Figure 10 indicates that for values of , i.e., smaller time increments, little if any improvement is seen for the extra computational effort spent in stepping through the extra number of time increments. This is in contrast to R-convergence in the dGa-framework, which shows linear rate of convergence in most cases.
(a)
(b)
(c)
(d)
An important point to note here is that errors in computations with values of close to one are almost the same as those computed using . Also, note that for significant error exists when crude spatial-mesh is employed, however, h-convergence shows hyper-rates that reduce the error quickly. Convergent process for is a significant characteristic of cLs framework since it reduces the computational cost significantly by taking large time increments. Convergence for is not necessarily monotonic, as shown in Figures 10(a) and 10(c), but for the significantly low value of h-convergence is monotonic and hyper and requires ten times less number of evolution steps as compared to Dof-equivalent computations.
3.5.2. -Convergence
Over the space-time domain, h-refinement means increasing the number of divisions in spatial direction. Similarly, -refinement implies an increase in the number of divisions in temporal direction over a given time slab with duration. Three uniform temporal divisions of were considered for -convergence study.
As can be observed from Figure 11(a), -refinement does not improve the solution if . For , the -refinement shows hyper-convergence in improving a crude solution, as shown in Figure 11(b). However, note that -refinement applied to spatial meshes with cannot improve the solution beyond the accuracy offered by for the same spatial meshing. Therefore, -refinement does not offer 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.
(a)
(b)
4. 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 diffusive error types, were employed to investigate convergence characteristics of dGa and cLs weak formulations of damped wave (1). To this end, standard refinements were applied and numerical errors were plotted against the degrees of freedom (Dofs) used by refined 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 diffusive, mechanical damping present in the system affects computational accuracy differently. Dispersive errors are higher for lower damping values and diffusive errors are higher for higher damping values. This characteristic is true for both dGa and cLs.
The Bms. problem is ill posed in dGa- framework as shown in Figures 3 and 4(b); employment of finite elements removes the issue and increases the accuracy in dGa-R convergence significantly 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 in cLs framework.
5. 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. The focus of this research has been to investigate the computational characteristics of cLs and compare them to those of dGa for viscoelastic wave (1). The first is the space-time coupled least squares (cLs) formulation, which was recently employed for solving the undamped wave equation [31]. The 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 ; however, for small values, e.g., for for which the time increments are ten times larger, h-convergence is of high rate and stable Figure 10. It was established that for 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., , yields comparable accuracy to cLs at . 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 . 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 is meshed was shown to be counterproductive as shown in Figure 11. Also contrary to dGa, where reducing (or increasing R) can reduce the error at a linear rate, cLs shows negligible gain in accuracy for , 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. This rate does not change with higher p-level or higher global continuity; however, employment of 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, h-convergence rates in cLs studies were at most quadratic for . While taking smaller time increments, i.e., was found to be without merit, large time steps, in particular was found to be stable with hyper-h-convergence allowing for reasonably refined spatial mesh evolving at to have comparable accuracy to the 10 times slower evolution of the same spatial mesh at . However, it should be noted that while evolution of the solution over large time increments of 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 ; the resulting system of equations is then evolved in time according to . The cLs coefficient matrix is formed over the two-dimensional domain 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 , 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 [email protected].
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.