A new nonlinear viscoelastic model and mathematical solution of solids for improving prediction accuracy

We developed an innovative material nonlinear viscoelastic model with physical mechanism and mathematical solution to improve existing ones. The relaxation modulus transits from the glassy stage to the rubbery stage through a time-dependent viscosity in a continuous spectrum considering the nonlinear strain hardening. Experimental results of differential solid materials including asphalt concrete, agarose gel, vaginal tissue, polymer, agar, bone, spider silk, and hydrogel demonstrate that the developed model is superior to generalized Maxwell model or Prony series for more accurate prediction outside of the range for data fitting while using much less model parameters. Numerical simulation results indicate that the new model has improved accuracy. It is stable numerically, and does not reduce computation speed. Therefore, the model may be used to simulate a broad range of viscoelastic solids for predicting experimental data and responses with improved accuracy.

The deformation of viscoelastic material is temperature dependent with thermal transition. Free volume changes or relaxation time are used to describe this transition behavior. Viscoelastic deformation has been interpreted by two atomistic mechanisms 1 : (1) distortion of chemical bonds' lengths and angles connecting atoms in a small and quick motion; (2) large-scale rearrangements of atoms and molecules. The thermal transitions with elevated temperatures may involve a few stages (e.g. for polymers) including the γ (local motion of molecular), β(bend/stretch of molecular), glass (from the glassy to the rubbery stage), and terminal transitions (melts into liquids) 2 . Different materials involve different transitions, but glass transition is the majority format of deformation for viscoelasticity as the focus of this study. Material viscoelasticity can be characterized by relaxation modulus E t ( ) or dynamic modulus ). E t ( ) and ω ⁎ E ( ) are dependent on temperature which can be converted to an equivalent time or frequency according to the time-temperature superposition principal. Researchers have developed physical models to describe the (linear) viscoelastic behavior. Based on the molecular dynamics theory, the Rouse model 3 uses the Brownian motion theory to simulate the single chain diffusion of beads that are connected via harmonic springs. The Kremer-Grest model 4 uses up to hundreds of chains 5 for simulating polymer elements. The single-chain theories such as the tube theory 6 and arm retraction model with arm-starts 7 (see Fig. 1a) have also been proposed to describe the linear viscoelastic behavior of entangled polymers. Among these models, the generalized Maxwell (GM) model appears to be mostly used for describing the glass transition of linear viscoelastic solids. For the unentangled polymers the spring-dashpot series of GM model physically represent different molecular chains with variable lengths under time distributions 1 (see Fig. 1b). E t ( ) of GM model is expressed as: where ∞ E is the modulus at infinite time, E i is the elastic modulus of the spring, η i is viscosity of the linear dashpot in series, and n is the number of spring-dashpot terms. This summarized exponential terms in a discrete spectrum is also named as Prony series (PS). It becomes a standard solid model when = n 1. The Maxwell or GM-based model and its PS formula have been widely adopted to fit modulus of linear viscoelastic materials 8 . These materials include polymers 9 ,dielectric elastomers 10,11 , glasses 12 , silicon 13 , tissue 14 , brain 15 , ligament 16 , blood vessels 17 , worm 18 , and asphalt concrete (AC) 19 . GM model has also been used as a base to be extended to plasticity model including the viscoplastic Bingham-Maxwell model 20 . Other formulas similar to GM model includes the (generalized) Kelvin model and Burgers model 21 . PS is computationally efficient with its exponential formula for time integration 22 , and thus it has been implemented as a standard material model in most numerical software such as the ANSYS and ABAQUS for simulating material and structural responses. The nonlinear viscoelastic models have also been proposed for simulating materials with large deformation or properties that change with deformation or time 23,24 . For example, Schapery's model has considered the spring's elastic modulus as a nonlinear function of time 25 , and some other models have presented the dashpot's stress in a nonlinear function dependent on the strain rate 26,27 or relaxation time 28 . Please note that our literature review doesn't cover fluid models.
However, a few critical questions arise. As discussed above, the models based on the theory of multiple molecular chains were derived and validated primarily for polymers, and they are less suitable for describing the physics of other materials with different morphology (e.g. asphalt with amorphous structure). For the GM model or PS: (1) the PS formula can produce an instability in fitting experimental data 19,29 as will be discussed later; (2) it is relatively difficult to determine a great number of model parameters using experimental data, (3) with a great number of model parameters the accuracy for fitting experimental data has improved mathematically. However, the physical interpretation of this large spring and dashpot system becomes less clear and more challenging.
These shortcomings of existing models are our motivations to seek development of a new material viscoelastic model that represents a broad range of materials more accurately. We have developed a theoretical model and mathematical solution. We conducted experimental validations for a number of materials that range from the inorganics to biomaterials. We have demonstrated that the model has improved the accuracy in both fitting experimental data as well as predicting modulus outside of the experimental range. The model is stable numerically, and it doesn't reduce computation speed. It uses much less model parameters than the GM model or PS. We also account for the nonlinear strain-hardening behavior in the proposed model.

Results and Analysis
proposed new material model. Micro-mechanisms and network-based models have been proposed for simulating material viscoelasticity in polymer science. For example, the coarse-grain-polymer network model represents the stretched mesh-like networks in a viscosity medium 30 , and the micro-network model symmetrically grows like tree 31 . We proposed a network model in a relatively simple format to interpret the physical mechanism of general materials as shown in Fig. 1c. The model domain consists of an elastic network and viscous media or "fluid" filled within the network branches. The elastic network carries a modulus ranging from the minimum ∞ E ( = ∞ t ) at rubbery stage to the maximum E 0 ( = t 0) at glass stage. Viscosity may be dependent on time and strain rate. E.g., it may follow a power-law function dependent on strain rate 32 or a linear function of time 33,34 . Here we proposed the viscosity model following a power function of time to consider its nonlinear behavior, as follows: where t 0 is a reference time (s), α is model parameter, μ is a viscosity parameter (MPa.s), and τ (s) is a retardation time (s) that accounts for temperature change as well as vanishes the unit of absolute time − t t ( 0 ). τ = 1 s when = T T 0 as constant, where T 0 is the reference or constant temperature. When ≠ T T 0 τ can be estimated using the WLF rule for time-temperature superposition as . Energy dissipation is produced by internal friction according to the interfacial friction theory-based viscous mechanism 35 . Viscosity is essentially a "fluid" friction 35 that transforms kinetic energy to heat through the interfacial friction. Therefore, we may understand μ as a friction coefficient between the elastic network and viscous media in this model. A higher μ value induces a higher force of friction, resulting in relatively higher E value ( Supplementary Fig. 1a). The fractional number α indicates the time or thermal sensitivity (here "thermal" has the equivalent effect of relaxation time according to the time-temperature superposition rule. A higher α value indicates a higher modulus variation rate ∂ ∂ E t t ( )/ (Supplementary Fig. 1b). The fractional number and derivatives are related to molecular theories describing the macroscopic behavior of viscoelastic materials 36 . www.nature.com/scientificreports www.nature.com/scientificreports/ The model can also be interpreted mathematically by a nonlinear spring (elastic network) and dashpot (viscosity) system as shown in Fig. 1d.
Applying a constant strain ε H t ( ) 0 on the model shown in Fig. 1d (H t ( ) is the Heaviside step function that = H t ( ) 1 when ≥ t t 0 and = H t ( ) 0 when < t t 0 ), the stress equilibrium satisfies the followings: where ε 1 and ε 2 are strain of the spring and dashpot, respectively, and is the modulus variation range.
Substitute Eq. (3) into (5), ε 1 can be solved as follows: can be derived as follows: starts at the plateau modulus E 0 at the glassy stage ( = t t 0 ), relaxes with time via the viscosity reaction rate and time sensitivity parameter α, and then transits to ∞ E at the rubbery stage until the infinite time. At the frequency (ω) domain, the viscosity can be expressed as follows: , under the sinusoidal loading signal the stress-strain relationship can be expressed as: According to the stress equilibrium the complex modulus ω ⁎ E ( ) of the model can be derived as follows: can be re-derived as follows (Supplementary information: Model derivation for complex modulus): In comparison, the complex modulus of GM model is expressed as 19 : The model becomes a standard solid model ( = n 1) when α = 0 and φ π = /2. To consider more parallel series like PS does, our model can also be extended as follows: where E i is the modulus of the i th parallel series of the system. The extended model may not be necessarily warranted for improving the accuracy of fitting experimental data, but it may offer a more flexible structure to fit and predict modulus in a wide range of relaxation time as if needed. We didn't derive the model function from experimental data, instead the model parameters can be determined by fitting the model formula on the measured experimental data of modulus. This approach is similar to some other models including the popular Prony series, for which the model parameters including viscosities are not directly measured instead back-calculated by fitting on the measured experimental data of E t ( ) or ω ⁎ E ( ). Strain hardening may occur for viscoelastic material. The model can be extended to consider the nonlinear strain hardening behavior of the elastic network E R as follows: (2020) 10:2202 | https://doi.org/10.1038/s41598-020-58240-y www.nature.com/scientificreports www.nature.com/scientificreports/ β is a coefficient to account for the nonlinear strain hardening. All model parameters are determined by fitting on experimental data using optimization method. The proposed model may satisfy the thermodynamic consistency (Supplementary information: Thermodynamic consistency) to ensure nonnegative energy dissipation. The viscosity is more often dominated by shear modulus 37,38 , and thus the shear and bulk relaxation modulus G t ( ) and K t ( ) are often not equivalent. However, in our numerical implementation we considered Poisson's ratio v as constant which is the often case for engineering analysis. Therefore, G t ( ) and K t ( ) can be derived following the same formula as E t ( ) according to their linear relationship 1,28, numerical solution method. Consequently, we developed a numerical solution to implement the proposed model for simulating responses of materials and structures. Numerical solution methods for the PS were presented in existing literature. To implement the model in numerical solution, different methods exist with their own advantages or disadvantages. Either the time domain or frequency domain can be used in the numerical solution 39,40 . Our objective for the numerical implementation is to obtain similar computation speed and numerical stability versus the PS as the basis. We developed a Galerkin-based finite element (FE) scheme within a time domain, which is robust and fast for computation without requiring transformation as the frequency-domain method does.
For a solid under a dynamic loading, we can express the stress equilibrium based on a strong form as follows: d that t d is the end of loading time, ρ is the density, b is the force of body weight, and Ω is the space domain. σ t ( ) can be expressed as follows 39 : where R is relaxation modulus matrix as a function of E t ( ), G t ( ) and K t ( ). Following the Galerkin method, multiply a test function p t ( ) and then apply space integrations on Eq. (15): Substitute Eqs. (16) to (17) and apply Green's function, the following weak form can be attained: is the loading vector, and f t ( ) is the external loading applied on the surafce ∂Ω as a natural bounary conditon.
For numerical solution purpose the time integration and differential are discretized to finite time steps. Different methods exist for nuermical approximation of differentials, and here we used the Eurler's rule and backward method to discretize Eq. (18) as follows:    . Follwing the deformation-strain-stress constitutive relationship, the responses of strain and stress can also be solved. experimental evaluation. We conducted experimental validations on different materials including inorganics, biomaterials and tissues. We conducted dynamic modulus tests on the AC material, at four temperatures (5, 10, 25, and 40 °C) and six frequencies (0.1, 0.5, 1, 5, 10, and 25 Hz). The measured modulus values of ⁎ E f ( ) at different temperatures can be converted to that at a reference temperature (25 °C) in a reduced frequency based on the WLF temperature-time/frequency superposition rule. Figure 2a plots the experimental results of ⁎ E f ( ) and its model fittings using the generalized reduced gradient optimization algorithm. We compared the optimization results of the proposed model and GM model or PS (see Supplementary A). Asphalt material is amorphous and has a low molecular weight (see the molecular structure example of a crude asphalt in Fig. 2b). Thus, we may not represent its structure with individual segment components or polymer-chain-like fractions 41 , and therefore its physical behaviors may not be properly interpreted by polymer models and the PS. Results have shown that the new model is able to fit experimental data accurately with a smooth curve. It is not surprising that the standard solid model (n = 1) can only fit a small frequency range of the experimental data. With n ≥ 7 (total ≥ 15 model parameters), GM model can capture almost a full frequency range of experimental data. However, it produces local "oscillations", and at the high frequencies ⁎ E f ( ) converges to E 0 very sharply beyond the experimental range, which seems less reasonable than expected. The average norm value of the 2nd order differential MP/Hz 2 ). The higher the norm differential, the lower fitting accuracy is. Fig. 2b plots the loss moduli ′′ E f ( ) within the frequency domain, in which again GM model produces local "oscillations" even with a fairly large number of terms (n = 20). Existing literature of this subject matter does not provide a proper interpretation of this phenomenon. In comparison, the proposed model produces smooth ⁎ E f ( ) fit on experimental data. The ⁎ E f ( ) value also gradually and smoothly converged when transmitting from the glassy stage to the rubbery stages outside of experimental data range. This result clearly illustrates that the proposed model has improved both the fitting and prediction accuracy as compared to GM model.
We also evaluate the proposed model versus GM model for both fitting and predicting E t ( ) of other materials including polymer, agar, bone, and tissue. Model parameters are determined by fitting on experimental data using generalized reduced gradient optimization algorithm. It is known that optimization may achieve different model parameters to satisfy the objective function. Our results have shown that the proposed model achieves more unique and stable optimization results than GM model even using the same number of model parameters, i.e. = n 2 for GM model with five model parameters (Supplementary Fig. 3). Generally, the fitting accuracy of GM model increases with a higher term number, but it may produce more variability of output parameters without a unique solution. We used a term number of 14 (total 29 model parameters) of GM model for the following validations since it can achieve satisfied accuracy on fitting experimental data. When only part of the experimental data are used for fitting model parameters, the rest of the data is used for predictions as validation. Figure 3a presents model fitting and predictions on experimental E t ( ) of agarose gel (data was reproduced from literature 42 , agarose is a polysaccharide polymer derived from seaweeds 43 ). Figure 3b presents model fitting and prediction for vaginal tissue (data was reproduced from literature 44 ). Again, the PS (corresponding to GM model in frequency domain) produces instability with many local "oscillations", as well as sharp transitions from the glassy stage to the rubbery stage beyond the experimental range which is likely false prediction. In comparison, the proposed model www.nature.com/scientificreports www.nature.com/scientificreports/ achieves very smooth curve fitting on experimental data. It also more accurately predicts the modulus value which gradually converges.
Secondly, the middle time-range of experimental data is used for fitting model parameters which are then used to predict E t ( ) at the lower and higher time ranges of experimental data. Results indicate the proposed model achieves fairly accurate predictions. In comparison, PS under-predicts E t ( ) at low t while over-predicts E t ( ) at high t as it tends to converge to E 0 and ∞ E more sharply than expected as shown in Fig. 3. Based on the thermodynamics theory, modulus of the viscoelastic materials is expected to grow and transmit smoothly when no sudden triggering occurs to change molecular structures 45 , which conflicts with the PS prediction results. The ) of the proposed model is 7.0% and 0.8% for the agarose gel and vaginal tissue, respectively, which are apparently lower than that of PS (12.0% and 1.2%), indicating its improved accuracy of predictions (see Fig. 3b). By using the extended model in Eq. (13) to fit experimental data, the accuracy has very slight improvement in fitting experimental data (e.g. for the agarose gel, the R 2 is 4.51%, 4.34%, and 2.10%, for n = 1, 2, and 3, respectively), indicating that a single series with only a few model parameters is capable of affording satisfied accuracy for both data fitting and prediction.
Validations on agar and bone materials have shown similar findings ( Supplementary Figs. 4 and 5). The model parameters of PS for all these materials are presented in the supporting information (Supplementary Table 2). numerical modeling results. After the theoretical demonstration and experimental validation, we implemented the proposed model for simulating responses of materials. We developed a computer coding for simulating responses in one-dimensional domain. We implemented the three-dimensional (3-D) model with user-defined material model in ANSYS software, in which the modulus and viscosity are defined as functions of relaxation time. We considered Poisson's ratio as constant, and thus calculated G t ( ) and K t ( ) from E t ( ) values based on their linear relationships for the 3-D model. Figure 4a,b present the simulated axial deformations of AC material under a constant loading pressure (1 MPa) and a sinusoidal loading with an amplitude of 1 MPa, respectively, for ∈ t [0, 10 s] and time step length ∈ . . dt [0 1, 0 001 s]. This case study indicates that at some times the simulated deformations using the proposed model can be significantly different from that of PS. For example, a maximun deformation difference of 0.06 mm (27.4%) and 0.022 mm (12.2%) of these two models is obsered at the  For PS model, the simulation accuracy may increase with a larger number of model parameters. However, it has been usual to use a term number of ≤ n 7 for PS in numerical practice for simplicity purose 46,47 , which may cause potential numerical errors. Resutls also indicate that the developed model is always numerically stable. Simulation accuracy increaes with reducing dt, but longer computation time. We have tested variable dt and find that the proposed model converges at = .
dt 0 01 s, the same as that of PS. This indicates that the proposed model does not require longer computation time than the PS. Simulated results clearly show the time lag of deformation to loading due to viscoelastic effect (see Fig. 4b).
We performed a 3-D dynamic viscoelastic simulation of a hexagon honeycomb structure made of agarose gel (see Fig. 5a), with inputs of E t ( ) exhibited in Fig. 3 and a density of 924 kg/m 3 . A sinusoidal loading pulse = +°− P t 5 sin(7200 90 ) 5(Pa) for ∈ t s [0, 2 ] with 20 cycles was applied on the top. All degrees of freedom are fixed at the bottom. We also simulated a solid model without hexagon holes with the same model size for comparison. Simulated deformation increases with time and remains positive at the end of time with zero loading (see Fig. 5b), illustrating the viscoelastic behavior with energy dissipation. It is known that honeycomb structure is excellent for sustaining moment loading. Our results show that the honeycomb structure is also a good option for sustaining compression loading, as it only cause slightly larger (i.e. 1.9%, see Fig. 5b) or ignorable deformation than the full solid model without honeycombs, while it has 26.1% less weight and two-times free surface area.
We also simulated the nonlinear stress-strain relationship of a spider silk thread and a polyampholyte-formed hydrogel material, taking into account the nonlinear strain hardening. In the lab test the spider silk was subjected to a stretching loading (data reproduced from literature 48 with permission, see Fig. 6a). The hydrogel displays viscoelastic behavior with stress recovery even at very large strains (laboratory data reproduced from literature 49 with permission). Results have shown that the proposed model is able to capture the nonlinear strain hardening www.nature.com/scientificreports www.nature.com/scientificreports/ behavior, and can accurately simulate the stress-strain curves, as illustrated in Fig. 6. It shall be noted that these materials may perform plastic deformations and fractures at larger strain ranges, which are not studied in this paper.

Discussion
We developed an innovative material nonlinear viscoelastic model to describe glass transition of solid materials for overcoming shortcomings of existing models. The model describes modulus with only five to six (w/o one additional parameter to consider nonlinear strain hardening) model parameters in a continuous spectrum. Experimental validations on different material types have demonstrated that the new model has improved accuracy in both fitting experimental data as well as predicting relaxation modulus beyond the experimental range as compared to GM model or PS which is the widely-used model for solid materials. Accurate prediction of modulus could be very useful since laboratory tests can only capture a limited range of reduced frequency or time.
The new model is able to predict E t ( ) outside of the experimental range more accurately with a smoother curve than the PS. The predicted curve captures the glass transition more naturally and smoothly than the PS. In comparison, the PS tries to achieve an accurate fitting only on given experimental data, but it may over-or underestimate modulus values outside of the experimental range, resulting in likely false E 0 or ∞ E values. The comparison between the new model and the PS (or GM model) is fair as we have used the same optimization scheme to determine/fit model parameters. When using the same amount of model parameters, the proposed model has achieved more unique solution than the PS as illustrated in the Supplementary A. Given a relatively very large term number of the GM model ( ≥ n 30), its accuracy may be improved for the data fitting (very likely) and data prediction outside of experimental range (less likely). Some researchers used a pre-smooth method to improve the fitting accuracy of the PS 29 . However, the new model could be superior than the GM model or PS considering multiple aspects as follows: (1) a model like the new one with less model parameters that satisfies the fitting accuracy is often preferred for the simplicity purpose; (2) a large number of model parameters of the PS (GM model) can produce more variability with non-unique solutions for determining the model parameters using mathematical optimization schemes; (3) it is more difficult to explain the physical mechanism for the model with a relatively much larger number of model parameters other than improving the accuracy for fitting the experimental data; (4) the proposed model offers more flexibility, e.g. it is a standard solid-model when α = 1, and (5) we also include the nonlinear strain-hardening behavior in the new model.
Our numerical solution have also indicated that the new model has improved the prediction accuracy. It is stable numerically. It does not reduce computation speed for convergence as compared to the PS.
Therefore, the new model may be used as an alternative for describing viscoelastic behavior of a broad range of solid materials to improve the data prediction accuracy. However, as GM model does the model is unable to simulate the special or abnormal physical behaviors such as the stress overshot of some materials including structural glasses 50 and amorphous solids 51 .

Materials and Methods
We validated the model based on laboratory testing data and provided physical interpretation. We developed numerical solution -a FE method to implement the model for simulating responses of a number of different materials.
We conducted relaxation and frequency-sweep (dynamic modulus) tests for part of these materials to determine their relaxation and complex modulus, while the left experimental data were obtained from existing literature. For the dynamic modulus test on AC material, four temperatures and six frequencies were used. The spider silk has a very small size (i.e. 20 μm) 48 . To measure the stress vs. strain behavior of the silk 48 , the material was adhered on a glass substrate under a unidirectional stretch (see Fig. 6(a)).
To fit model parameters from experimental data, we used the generalized reduced-gradient optimization algorithm which is one of the most frequently used optimization methods in practice. We also evaluated the model for predicting modulus values outside of the frequency or reduced-time range of experimental data, which was not (often) conducted by existing literature. By using the major middle part of the experimental data in time or www.nature.com/scientificreports www.nature.com/scientificreports/ frequency domain, we fitted the model parameters. Consequently, we used the fitted model to predict the modulus values outside of the time/frequency domain used for model fitting, and compared prediction values with the actual experimental measurements.
We developed a Galerkin-based FE method to implement the proposed E t ( ) formula and derived its mathematical solution. The method is based on a time-domain scheme that integrates E t ( ) with time and space for a robust and fast computation [52][53][54] .