Modeling of viscoelastic structures with random material properties using time‐separated stochastic mechanics

Modeling and simulation of materials with stochastic properties is an emerging field in both mathematics and mechanics. The most important goal is to compute the stochastic characteristics of the random stress, such as the expectation value and the standard deviation. An accurate approach are Monte Carlo simulations; however, they consume drastic computational power due to the large number of stochastic realizations that have to be simulated before convergence is achieved.


INTRODUCTION
The incorporation of stochastic fluctuations has recently received increasing attention in modeling of materials. Microstructural irregularities or accidental variations in the production process result in material properties that do not always equal a fixed, idealized value, but exhibit unintended fluctuations. Therefore, the material parameters cannot be considered deterministic but can be modeled by random stochastic variables. With structures increasingly being constructed with a reduction in resources and with construction costs in mind, the buffer for structural integrity becomes smaller. The stochastic fluctuations of the material properties can then have a potentially drastic effect on the reliability of the construction. It is thus of great interest to estimate the effect of the random properties on effective quantities, such as the stress or the reaction force. This can be captured, for example, by calculating the expectation value of the stress together with the variance or standard deviation, which measures the order of resulting fluctuations of the stress. Models that account for stochastic aspects can be derived for different material classes and can be investigated at various length scales. In many cases, these models possess the advantage of an inherent inclusion of geometrical aspects. The starting point for various strategies is modeling of the elastic tensor as a random field. For instance, several approaches to modeling of stochastic materials at the material point level make use of a broken chaos polynomial expansion 1 or Karhunen-Loève expansion 2,3 ; see, eg, the works of Malyarenko and Ostoja-Starzewski 4,5 for elastic and Rosić and Matthies 6 for inelastic materials. Moreover, this mathematical concept can also be transferred to the spatial level constituting the so-called stochastic finite element method. [7][8][9] It represents a prominent tool for the computation of random fluctuations of the displacement field, from which random fluctuations of the strains and hence stresses can be derived. [10][11][12][13][14] The stochastic finite element method is based on a broken polynomial chaos expansion of the displacement field with the series terms as additional nodal unknowns. In this expansion, the influence of the random variables is determined by deterministic coefficients, which can be computed by solving an eigenvalue problem for the covariance operator of the random field. Corresponding coefficients of the resulting displacement fields can then be calculated in an enlarged system of equations. Although it can be proven mathematically that the stochastic series expansion converges to the exact stochastic field, quite a large number of series terms must be considered due to slow asymptotic behavior. 7,15 Therefore, the number of nodal unknowns drastically increases, turning the stochastic finite element method into a rather computationally expensive strategy.
Random fluctuations can also be incorporated into the modeling of materials by means of the stochastic perturbation method, which is based on a (direct) Taylor series expansion of the stochastic field variables. [16][17][18] Using a mixed perturbation Monte Carlo (MC) method, this strategy has proven advantageous for the computation of effective elastic constants for materials with stochastic inclusions. This approach operates on the scale of representative volume elements (RVEs) and allows for discretization of the microstructure. Unfortunately, implementation of the stochastic perturbation method into a finite element routine for macroscopic specimens requires solving a (nonlinear) finite element problem at the microscale for every integration point at the macroscale. This "inner" finite element method embedded into an "outer" finite element method is also utilized within the FE 2 approach. Although this strategy provides a detailed insight into the stochastic microstructure, 19,20 the computational costs are extreme, such that this method is limited to small-scale problems.
Recently, several other approaches have been developed, for instance, the generation of stochastic microstructures constitutes an important field to capture the multiscale effects on the random material behavior (cf other works [21][22][23] and the textbook of Sobczyk and Kirkner 24 ). Methods for the computation of effective material parameters with random input have been proposed, eg, for the elastic case in the works of Guilleminot and Soize 25 and Drygaś and Mityushev, 26 for bone material in the work of Basaruddin et al, 27 and for polymer composites employing a multiscale approach in the work of Shokrieh and Rafiee. 28 An approach for approximating the given random stresses for isotropic hyperelastic materials has been presented in the work of Mihai et al. 29 All approaches are converging approximations of the stochastic fields and thus constitute promising strategies for contributing to the fundamental problem of modeling materials with stochastic properties. Unfortunately, several unsolved drawbacks exist, depending on the individual concept, for instance, it is not a priori clear how to "extrapolate" an experimentally determined covariance in a testing device, which describes the spatial fluctuations of the microstructure, to a macroscopic construction part. This complicates direct usage of such methods. Furthermore, the application of some of the mentioned methods to large-scale macroscopic specimens and engineering structures is not yet possible due to extreme computational costs. Therefore, they are limited to rather small problems until now, which prohibits their usage for industrially relevant purposes. Lastly, stochastic approximation methods may not provide the full spectrum of requested information, ie, details on the (local) fluctuations might be lost and it is not clear how the output parameters (stresses) fluctuate when the input parameters are random (material parameters).
In a previous paper, 30 we developed a new approach to modeling random fluctuations of the stress in viscoelastic materials, which resulted in time-separated stochastic mechanics (TSM): an approximation of the resulting stress at the material point level with the main feature that it can be written in terms of functions depending, on the one hand, on the random fluctuations but which are constant in time (material parameters), and on the other hand, on the time step but which are deterministic (loads). To be more precise, the complete stochastic influence can be analysed in an a priori manner, independently of the individual boundary value problem. This enabled us to obtain a mathematically tractable approximation of the resulting stress and thus calculate its expectation and variance. In this new approach, there is no need to calculate series coefficients nor does the computation cost increase substantially. Instead, the expanded set of material parameters can be used to formulate a modified deterministic finite element problem. The numerical solution to this finite element problem then yields instantly the expectation and the variance of the random stress and of the reaction force.
We remark that there exists an essential difference between our approach and the other methods discussed above. In our approach, we focus on elastic constants, which are random but homogeneous in the macroscopic body. This means that our model takes account of the random fluctuations on the same scale on which the elastic constants are determined in an engineering fashion: macroscopic specimens are tested and, due to fluctuations in the microstructure (inside the specimen and on its surface), stochastic oscillations of the effective material parameters are observed. An experimental evidence for this modeling approach is given, eg, by Niyigena et al, 31 where exactly the mentioned measurement process has been used. The stochastic fluctuations of the Young's modulus in hemp concrete are reprinted in Figure 1, from which the wide range between 19.75 and 43.75 MPa is determined. In the very same manner as performed for receiving these results from the experimental investigations of specimens, our approach homogenizes spatial fluctuations but accounts for fluctuations among different samples. In view of practical applications, these fluctuations, which serve as input parameters, can be directly estimated by repeated measurements (cf the work of Niyigena et al 31 ). However, in general, we note that the method of TSM is not limited to a homogeneous distribution of elastic constants and could be generalized to a spatially varying random field. Furthermore, it is worth to mention that TSM is not restricted to Maxwell-type models but can also be applied to general Maxwell solid models and Kelvin-type models.
In this paper, we demonstrate that the TSM developed in the work of Junker and Nagel 30 for the material point level can also be applied to macroscopic full-body problems. In other words, the TSM allows to directly use the experimental findings based on testing specimens, eg, those by Niyigena et al, 31 to predict the stochastic behavior of macroscopic construction parts with arbitrary shape. Unfortunately, the original formulas in the work of Junker and Nagel 30 were introduced in a format that makes them very cumbersome for a finite element implementation. For instance, tensors of order eight have to be programmed and integral functions to be solved. We thus derive in the this paper an explicit representation of the equations that has a much simpler structure, enabling a numerical evaluation in higher dimensions than 1D. Furthermore, we derive formulas for the expectation and variance of the reaction force. After the derivation, we explain how to embed them in a finite element routine. We test our approach in two benchmark problems, for the brick with a centered hole and the spring, and for different classes of distribution of random elastic constants and varying loading velocities. In all cases, we can predict the true expectation and standard deviation of the stress and of the reaction force, as delivered by corresponding MC simulations, with very high accuracy. In addition, the computation cost compared to MC simulations is drastically reduced. In fact, the only additional computation cost is due to the estimation of the additional constants carrying the stochastic information, which has to be performed only once for each material and is then independent of the specific boundary value problem. Finally, we remark that this makes it also very easy to apply our method to other boundary value problems and to different distributions of elasticity tensors. The additional constants serve as a set of additional material parameters, which then capture not only the idealized deterministic behavior of the material but also incorporate the stochastic fluctuations. This paper is structured as follows. In Section 2, we recall the material model and derive the simpler formulas for the results in the work of Junker and Nagel. 30 In particular, Section 2.3 contains the equations for expectation and variance for the stress. Comparing with the aforementioned work, 30 these quantities are formulated in a simplified way that allows easy numerical implementation. Additionally, we also give in Section 2.4 the expectation and variance of the random reaction force. Detailed instructions for the finite element implementation of both the TSM approach and the MC simulations are provided in Section 3. The numerical study is then presented in Section 4. Finally, some technical derivations to obtain the simplified equations are deferred to Section 5.

Stochastic modeling of material properties
The stochastic fluctuations of the material parameters is taken into account by stochastic elastic constants E. That is, E does not take some fixed value but is a random variable defined on some probability space (Ω pr ,  pr , P). * Any element of the sample space Ω pr can be interpreted as a specific realization, resulting in the elastic constant E( ). However, the modeling of stochastic fluctuations , and then also the value of E( ), cannot be assumed to be known. Instead, we compute for suitable real-valued functions f the expectation Once we specify the distribution of E, all expectations as in (1) are determined by this distribution and the precise definition of Ω pr is not relevant anymore. The fluctuations around the expectation can then be measured in terms of the variance and the standard deviation is given by Std( (E)) = √ Var( (E)). Of course, many quantities of interest are not scalar values, as for example, the stress tensor. For general tensor-valued f, we define the expectation ⟨f (E)⟩ componentwise, and the variance is Note that this increases the order of the tensor, as we also need to consider all possible covariances between entries. While this is important in the derivation of variances, eg, of the stress, we give in the following section equations for the variance of individual real-valued entries.
In this paper, we make the following assumptions regarding the distribution of the elastic constant E: • As described in the introduction, the value of E is random, but constant on the macroscopic body.
• The random realizations of E are positive definite and have finite exponential moments, that is, is finite for some positive constant c > 0.
These assumptions allow every general distributions for E, for example, entries which are normal distributed, gamma distributed, uniform on some interval, or parametrized by random Lamé parameters. In particular and in contrast to stochastic finite elements, we do not require the computation of any series expansion. Instead, all necessary coefficients can be directly estimated from experiments.
We remark that, in addition, the viscosity could also be modeled as a random variable. In the present work, we assume to be deterministic, but a possible approach to incorporate a stochastic viscosity based on an "emulated" microstructure has been presented in the work of Junker and Nagel. 5 *The sample space Ω pr is the set of all possible outcomes of an experiment and  pr a collection of subsets of Ω pr , to which the measure P can assign a probability (see the work of Billingsley 32 ).

Stochastic solution for the stress
The application of the TSM for viscoelastic materials has been derived in the work of Junker and Nagel 30 and is briefly recalled here for convenience. For programming simplicity, we present all formulas here in the well-known Voigt notation. As usual, the viscous part v of the strain serves as the internal variable that is determined by the well-known evolution equation where denotes the viscosity, while the stress is given by with the elastic constants E. The deviator can be computed with the aid of the operator S, defined as Consequently, it holds that In the work of Junker and Nagel, 30 this material model was investigated with a random elastic constant E, which serves as first example of the TSM. The total strain is approximated by a deterministic function, but by (5) and (6), the viscous part v as well as the stress are random. It was then concluded that the solution of Equation (5) can be written as Here, the integrals I (1) and I (2) are given by and More details on the derivation can be found in the work of Junker and Nagel. 30 Inserting this formula into the stress in Equation (6) yields the stochastic stress The main feature of (9) is the TSM for v , and then also for , where the integrals I (1) and I (2) depend on time, ie, on the loads, but are deterministic. The random variable E appears in prefactors of these functions, which are constant in time, but fluctuate randomly. Using this decomposition, the expectation ⟨ ⟩ as well as the variance Var( ) can be analytically computed (see Section 4 in the work of Junker and Nagel 30 ). The formula for the variance in the original publication also included covariance terms that are of less interest, and tensors of order eight appeared in the formula which makes it quite cumbersome to program. Therefore, in the following section, we give simplified formulas for the diagonal entries of the variance tensor. Furthermore, we present formulas for the expectation and variance of the reaction force.

Expectation and variance of the stress
For the sake of a clearer presentation, we introduce the following matrices, which appear as stochastic coefficients in formula (10), Then, the stress has the compact representation and for each component By linearity of expectation, the expectation of the stress components is then given by The variance of i is more intricate and computed in Section 5. The resulting formula is The coefficients E pq ii , p, q ∈ {0, 1, 2} for the variance are obtained from the matrices in (11) by Equation (15) contains the expectation of these coefficients, which, again by linearity, is determined by These expectations can be computed by MC simulations once the stochastic properties of the elastic constants of a material are known. It is important to notice that these parameters are computed only once for each material. To be more precise, they increase the set of elastic constants such that it includes sufficient stochastic information. Consequently, they do not have to be computed for each individual boundary value problem but are fixed for a given material with specific distribution. In other words, the entire stochastic analysis can be performed independently of any boundary value problem in an a priori manner.

Expectation and variance of the reaction force
An additional important quantity is the reaction force that reflects the global response of the construction part, which has not been considered so far in the work of Junker and Nagel. 30 For the case of stochastic material properties, this reaction force is also stochastic. We derive closed formulas for the computation of both expectation and variance of the reaction force, based on our new approach. The calculations are postponed to Section 5.
The reaction force for each component i ∈ {1, 2, 3} is given by wherẽ=̃i e i e denotes the stress in the usual matrix notation. The quantity t is the traction vector on the Neumann boundary Γ ⊂ Ω. Let us introduce the nabla operator in Voigt notation by B, defined by = B · u, such that We note that the product in the last integral is actually the scalar product of two vectors B i · when we write B i = B ji e j , j ∈ {1, … , 6}. Again due to its linearity, the expectation is fairly simple to compute and reads where ⟨ j ⟩ is as in (14). For the variance, we need to calculate Var The first moment may be calculated according to (20), and we show in Section 5 that the second moment is given by where E (pq) for p, q ∈ {0, 1, 2} is the fourth-order tensor defined by and the matrix coefficients  (p) i are defined as As may be expected, covariances are present in formula (21), which, unfortunately, increase the complexity. However, the evaluation effort in a finite element context still remains moderate. It is a remarkable consequence of the structure of as in (10) that the second moment of F i , which involves a double integral according to its definition, may actually be calculated only by evaluating the single integrals (23).

Evolution equations for I (1) and I (2)
The integral representation of the deterministic functions I (1) and I (2) are cumbersome for numerical implementation, in particular due to the matrix-valued exponential function. Hence, we present here an alternative representation in terms of a differential equation, which avoids calculation of the exponential. The derivation is again postponed to Section 5. The first integral function can be computed according to Computation of the rate for the second integral function yields .
For both I (1) and I (2) , the initial value at t = 0 is 0 for a virgin material with no viscous strain contribution. The differential Equations (24) and (25) have a mathematical structure that is also well known for modeling of deterministic viscoelastic materials and can then be easily evaluated in a numerical computation. We remark that we have to solve two (deterministic) differential equations instead of the single ODE for the deterministic case. However, computation of this solution is basically for free in terms of numerical costs; the solution of the system of algebraic equations provided by the finite element method requires a computational effort that is several orders of magnitude larger than the numerical solution of Equations (24) and (25). Of course, evaluation of the variance in Equation (15) requests the evaluation of several scalar products. The impact on the global computation time is still very small such that our simulations, which cover the full stochastic information, never need more than 1.05 of the computation of one deterministic simulation (whereas a computation in the MC setting requires at least 400 repetitions of one deterministic simulation for convergence).

ALGORITHMIC TREATMENT
All results obtained with our novel approach will be compared to MC simulations that serve as reference. The respective algorithmic treatment is presented in this section. Both the new model and the MC simulations are evaluated within a finite element context where the weak form of the balance of linear momentum ∇ · + b = 0 for the quasi-static case is solved. To this end, the residual R n+1 for the current time step n + 1 is defined in usual manner as with the shape functions S and the external body forces b.
The roots of R n+1 are usually found by employing iterative gradient solvers. Therefore, the tangent operator needs to be computed. Here, the property n+1 = B · u n+1 has been used. The implementation of a material model via a common material model interface thus demands the computation of the stresses as well as its derivative with respect to the strains, both at the current time-step n + 1, ie, n+1 and d n+1 ∕d n+1 , respectively. We present the individual formulas for the novel approach and the MC simulations in the following subsections.

Novel approach
In the novel approach, the viscous strain is approximated by the two internal variables I (1) and I (2) , which can be expressed in terms of the differential equations (24) and (25). These differential equations are approximated via implicit Euler approaches, yielding .
where the subscript n + 1 is related to the current time step and n denotes the previous time step. The time increment is given by Δt. The values of I (1) and I (2) at the current time step, after rearranging Equations (28) and (29), thus read and It is worth to mention that they can be evaluated in a staggered scheme: the value for I (1) is updated first by usage of Equation (30); inserting this value into Equation (31) enables us to update I (2) .
For the finite element implementation, the expectation of stress is used to compute the residual R n+1 in Equation (26). For the current time step, it is given by For application of the usual Newton-Raphson iteration procedure to compute the roots of the residual, the algorithmic tangent operator has to be calculated, with the derivatives and which is inserted into Equation (27). The variance of each stress component i at the current time step n + 1 is given by Equation (15) Var with the expectation of the matrices E (pq) ii according to (17). The values for each component of the variance are computed within the material model interface and stored for later visualization. Obviously, they serve as postprocessing output variables that do not need to be calculated for the finite element routine to work. However, the variance provides important information regarding the probability estimation for the stress values to exceed some critical value ⟨ i ⟩ + crit : even when no information regarding the probability distribution function of i is known (which is usually the case), Chebyshev's inequality may serve as (rough) estimator. The reaction force is an additional postprocessing output variable whose expectation according to Equation (20) can be computed in a straightforward manner using Equation (32), yielding with n e being the number of finite elements with nonzero prescribed displacements. The variance of the reaction force for the current time step is given by where the first term is computed by and, following (23), the matrices  (p) i,n+1 are given as We mention that the reaction forces are identically zero when no boundary values are applied on the specific node (balance of linear momentum). Consequently, the result is deterministic for almost all nodes for which the variance is zero. Equations (38) and (41) are thus only evaluated for elements that contain nodes with nonzero boundary conditions.

MC simulation
For the MC simulations, deterministic finite element simulations are performed, each with a different realization for the stochastic elastic constants. The expectation and variance of the individual quantities are computed in a postprocessing step after the last realization has been computed for the entire time interval being investigated.
Of course, the material model is viscoelastic, reading in Voigt notation which is approximated by an implicit Euler approach as The viscous strains of the current time step thus read v n+1 = For the finite element implementation, the current stress is given by with the consistent tangent operator which is necessary for the evaluation of the residual and its derivative according to Equations (26) and (27), respectively. To compute the expectation and variance of the stress, the above procedure is repeated for a series of N realizations E MC 1 , … , E MC N of the elastic constant generated according to its distribution. The values MC 1 , … , MC N of the corresponding stress and its square at each integration point within each finite element are stored for each time point and each MC iteration. After all values have been computed and collected, the expectation is estimated by The variance of the stress is estimated by Analogously, the estimated mean of the reaction force is and its variance As for the TSM formulas, Equations (49) and (50) are only evaluated for elements with nonzero boundary conditions for the displacements.

SIMULATION RESULTS
We consider two different boundary value problems, the brick with centered hole and the spring, for which we compare our equations from Section 2 with the results of MC simulations. The number of iterations in the MC simulation is set so high that the resulting estimates for expectation and standard deviation obtained from the MC simulation can be treated as the reference values (N = 400), to which we compare the expectation and standard deviation calculated according to Equations (14) and (15) for the stresses and Equations (20) and (21) for the reaction force, respectively.
We use trilinear shape functions with eight Gauß points per element. The stress plots over curves are obtained by extrapolating the Gauß point quantities to the nodes. We use ParaView for visualization.
The probability distribution of the random elastic constant E is parametrized in terms of the Lamé parameters and . The two parameters are gamma distributed in Section 4.1 and stochastically independent. The mean of is 1000 MPa, and the mean of is 800 MPa. The standard deviation for both parameters is 20% of their mean. To demonstrate the robustness of the method, we also consider in Section 4.2 a uniform distribution with the same mean and standard deviation.

Main simulation results
For improved presentation purposes, we begin with the discussion of the main results for both boundary value problems. Here, we restrict ourselves to the gamma distribution and one loading velocity: the maximum prescribed displacement is applied within 20 seconds. Each second is discretized with five time increments, ie, Δt = 0.2 seconds. Consequently, 100 time increments are required in total for each boundary value problem.

Brick with a centered hole
The bottom right diagram in Figure 2 shows the boundary value problem of the brick with a centered hole. Only one quarter is discretized by a finite element consisting of 8606 nodes that form 4183 elements. One element is used in the thickness direction, such that a plane stress state is approximated. Symmetry boundary conditions are applied along the left vertical line and the bottom horizontal line, whereas all nodes at the top horizontal line are subject to a prescribed displacement that is a linear function in time.
In the bottom right part of Figure 2, the distribution of the norm of the standard deviation of the stress is plotted for maximum prescribed displacement. For the same point in time, the diagrams (A) to (E) show the norm of the stress plus/minus the standard deviation, ie, | ± Std( )|. Here, we let | · | be the standard Euclidian norm, ie, |a| = √ a 2 1 + … + a 2 6 . This is plotted along the five radial directions that are drawn in the bottom right diagram with relative positions indicated by r. The corresponding angles are = {0, 1, 2, 3, 4} ∕8. The red curves in (A) to (E) represent the results obtained with our model; the reference MC results are represented by dotted black lines.
It can be clearly be seen that we obtain extraordinarily good agreement between the results with our equations and the MC simulations. Furthermore, the results demonstrate that the standard deviation of the stress is not simply proportional to the standard deviation of the elastic constants nor to the expectation of the stress: it varies strongly. For instance, in the radial direction with = 3 ∕8 (plot (D)), almost no standard deviation is present at 18% in the radial direction. In contrast, the relative standard deviation is the highest at the top right corner of the brick (plot (C) at 100% in the radial direction).
The global reaction of the brick with a centered hole is given in terms of the force/time and force/displacement diagrams, respectively (see Figure 3). On the left-hand side, the reaction force in the vertical direction is plotted versus time; on the right-hand side, the force is plotted versus the displacement. Again, both the expectation (solid curve) plus/minus the standard deviation (dashed) are given. The red curves show the results with our model, and the black curves show the results of the MC simulations. As for the stress plots in Figure 2, we see remarkable agreement between our results and the reference ones.
The characteristic of increasing and decreasing standard deviation is also present in the force diagrams. Interestingly, the reduction of the standard deviation takes place over time and space, as observed in Figure 2.

The spring
The second boundary value problem is the spring, as depicted in the center of Figure 4. Here, we used a finite element mesh with 8607 nodes and 7200 elements. In total, five windings with an escaping arc of 30 • result in a spring length of 16.82 mm. The inner radius is 2 mm, and the radius of the bar is 1 mm. Consequently, the outer radius of the spring is 4 mm. The nodes at the left-hand end face are fixed; those at the right-hand end face are subjected to a displacement in the longitudinal direction (y-direction). All other displacements at this end face are zero.
In the top and the bottom diagrams of Figure 4 show the norm of the stress plus/minus the standard deviation along two different curves: the stress along the magenta spiral is plotted at the top; the stress along the blue spiral is plotted at the bottom. The magenta curve runs along the bottom of each cross-sectional area; the blue one runs along the inner radius of the spring.
The results with our model are again plotted in red and the MC result in black. The agreement between our results and the reference solution is once more very good, although larger deviations are present, which result from the much more complex stress state that is present for the spring: the dominant torsion load state is superposed by remarkable normal stresses due to bending and tension in the spring. As for the brick with a centered hole, spatial positions with vanishing standard deviation can be identified.
The global reaction of the spring is given in Figure 5 in terms of a force/time diagram on the left-hand side and a force/displacement diagram on the right-hand side. Both the expectation and the expectation plus/minus the standard deviation are present. Red curves indicate our results, and the black curves indicate the reference MC results.  As for the brick with a centered hole, we observe excellent agreement between the reference solution, displaying the real answer, and our approach. Due to the complex fully three-dimensional stress state, this example serves as impressive demonstration of the validity of our model.

Detailed analysis of the simulation results
We present a more detailed analysis of our results in this subsection. Here, we also investigate the uniform distribution and a different loading velocity for the spring when the maximum load is then applied within 40 seconds. We use 100 time increments also for this simulation, demanding Δt = 0.4 seconds here. This, of course, results in large viscous effects.

The brick with a centered hole Individual stress components
The expectation and variance of the individual stress components over the radial directions (A) to (E) in Figure 2 are plotted in Figures 6 to 10. Again, our results are plotted as red curves, those of the reference MC simulation in black.   The agreement between our result and the MC simulation is very good, particularly for the stress components that are of a larger amount, ie, for the yy components. Also for other rather complex stress distributions, we receive a very good agreement, eg, for shear stresses. However, larger deviations are present for the xx and zz components in some directions. These components are not predominant since other stress components are of a larger amount, which is the reason for the excellent agreement of the norm in Figure 2

The spring Individual stress components
For completeness, we present the expectation and standard deviation of each stress component along the magenta and blue spirals in Figure 11 and Figure 12, respectively, for which the norms have been presented in Figure 4.
The individual stress components are very complex curves with varying "wave lengths". As already mentioned, there are some stress components that show a vanishing standard deviation. Our results (red curves) are in very good agreement with the reference MC solutions (black curves), particularly for those components that are of a larger amount.

Halved loading velocity
The force/time and force/displacement diagrams for the spring loaded within 40 seconds is given in Figure 13. For improved comparability, we also plot the result for the reference loading velocity, as given in Figure 5, as transparent curves.
Again, our results are almost identical to those from the reference MC simulation. Only at larger time steps does our solution start to diverge. This effect is expected since we use a broken series expansion for the exponential function (cf the work of Junker and Nagel 30 ). However, even for the much more prominent viscous behavior, very convincing agreement can be achieved by our approach.

Uniform distribution and relaxation
As a final example, we once more investigate the spring and apply the load with halved velocity. However, we skip the unloading but keep the prescribed displacement at its maximum value (10 mm, which will be applied within 20 seconds) to analyze the relaxation behavior. We, furthermore, use a uniform probability distribution function to demonstrate that our equations are valid also for this case. We present the global reaction in terms of the force/displacement diagram in Figure 14.
As for the other examples presented, we also observe for this boundary value almost no deviation between our results and those from the reference MC simulations. Only at larger time steps do the curves start to diverge. This could be improved by approximating the solution of the evolution equation for the viscous strains with higher terms (see     (9) and the derivation in the work of Junker and Nagel 30 ). Nevertheless, the results are more than convincing, particularly when the numerical costs are compared as investigated in the following section. Moreover, we note that the uniform distribution evaluated here is quite different from the gamma distribution, thus demonstrating the robustness of the method.

Numerical cost
The TSM and the MC model were implemented into the finite element program FEAP. All simulations were executed on a desktop PC with an i7 CPU at 2.93 GHz and 8 GB memory. The respective computation times are given in Table 1. The simulations are performed by discretizing the time into several time steps and thus load steps. For each load step, one finite element iteration is required to solve the updated displacement field and viscous strains that balance the linear momentum. One finite element iteration for one stochastic realization in the MC simulations needs the same computational effort as a simple deterministic simulation. To achieve reliable estimates for expectation and variance, this needs to be repeated for a larger number of iterations. That is, to obtain the stochastic information, the computation cost using the MC method amounts to the computation cost of a single finite element computation, multiplied by the amount of iterations, often several hundred.
In the novel approach using the TSM, the computational effort consists of computing the coefficients in the equations for the mean and the variance and then one single finite element realization using these deterministic coefficients. The former step of computing the coefficients is also realized by the MC method, for which, however, no finite element iteration needs to be performed. Consequently, this MC simulation is very fast (in the range of few seconds). As mentioned, this has to be performed only once for a particular material with the respective distribution function. It thus does not affect the finite element evaluation of the TSM equations. This evaluation only demands the solution of two ODEs and several dot products, which amounts such that the novel approach needs about 5% more computation time than one deterministic simulation for the entire time interval considered. This increase is very moderate since we are able to compute the basically identical result that is found using the MC simulation. In contrast to our moderate time consumption, the MC simulations needed for the presented boundary value problems 400 stochastic realizations, which means that we consume less than 0.3% of the computation time of the MC simulations. To give an example, our simulations finished in less than 5 minutes, whereas the MC simulations needed more than 24 hours.
In conclusion, the numerical extra costs required by our approach are negligible compared to those for a deterministic simulation but are capable of accurately computing the expectation and the variance of both the stresses and the reaction forces. This makes our approach very promising for computing the stochastic behavior of viscoelastic materials at reasonable numerical effort.

DERIVATION OF THE VARIANCES
In this section, we present the detailed derivation for all equations presented above.

Derivation of the variance of the stress
We recall that all calculations are in Voigt notation. We concentrate on the variance of individual entries of the stress tensor, that is, of i , i ∈ {1, … , 6}, given by For these scalar random variables, the variance can be computed via We calculate first the square of the stress component as Here and in the following calculations, i is a fixed index and there is no summation over i. For each component, we then obtain as products of random variables the nine matrices E pq ii , p, q ∈ {0, 1, 2}, as defined in (16). We compute the individual terms in Equation (52). Multiplying out the first one leads to i · I (2) .
The second term yields ii · I (2) .
To obtain the variance, we have to take the expectations in (58) and subtract ⟨ i ⟩ 2 , which proves (15).

Derivation of the variance of the reaction force
The second moment of the reaction force can be calculated for each component i ∈ {1, 2, 3} as where we indicated by (·)| z that the respective functions have to be evaluated at the point z. The respective stresses read Consequently, the product j | x k | y is given by Recall the definition of the fourth-order tensor E (i ) abcd = E (i) ab E ( ) cd . Taking the expectation in (60) and multiplying by the matrices B ji , we get This may be written in tensor notation as  (1) )| x ∶ ⟨E (11) ⟩ ∶ (B i ⊗ I (1) )| + (B i ⊗ I (1) )| x ∶ ⟨E (12) ⟩ ∶ (B i ⊗ I (2) )| − (B i ⊗ I (2) )| x ∶ ⟨E (20) ⟩ ∶ (B i ⊗ )| + (B i ⊗ I (2) )| x ∶ ⟨E (21) ⟩ ∶ (B i ⊗ I (1) )| Since all these terms are a product of functions evaluated at x and functions evaluated at y, the double integral in (59) simplifies to After integrating, we identify again several symmetry relations in (63), eg, Therefore, the integral reduces to with the matrices  (i) i as in (23). The structure of the second moment of the reaction forces is thus quite similar to the second moment of the stresses (see Equation (58)).
Then, we have answer because they reflect the experimental procedure. For all results, individual components of the stress, the Euclidean norm of the stress, and components of the reaction forces, we receive strikingly good agreement between our approach and the MC results. Our model, however, demands only negligibly more computational effort than one deterministic simulation. Since the MC procedure requires the repetition of hundreds of simulations, our approach is faster by order of magnitudes, simultaneously being in excellent agreement with the "true" stochastic material behavior. Due to the very simple implementation of our model, a broad applicability is given, and large-scale simulations can be performed in future investigations.
Our method, referred to as time-separated stochastic mechanics (TSM), opens a new perspective on how to deal accurately and effectively with modeling of materials with stochastic properties. All stochastic information is included in deterministic additional parameters that can be easily found once the stochastic behavior of the material parameters has been determined experimentally, eg, as done in the work of Niyigena et al. 31 Then, the entire stochastic analysis can be performed in an a priori manner.