    # Numerical time integration of lumped parameter systems governed by implicit constitutive relations

Time-integration for lumped parameter systems obeying implicit Bingham-Kelvin constitutive models is studied. The governing system of equations describing the lumped parameter system is a non-linear differential-algebraic equation and needs to be solved numerically. The response of this system is non-smooth and the kinematic variables can not be written in terms of the dynamic variables, explicitly. To gain insight into numerical time-integration of this system, a new time-integration scheme based on the trapezoidal method is derived. This method relies on two independent parameters to adjust for damping and is stable. Numerical examples showcase the performance of the proposed time-integration method and compare it to a benchmark algorithm. Under this scheme, implicit-explicit integration of the governing equations is possible. Using this new method, the limitations of the trapezoidal time-integration methods when applied to a non-smooth differential-algebraic equation are highlighted.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Constitutive relations play a major role in modeling response of mechanical systems. These constitutive relations often can be written as a relation between kinematic and dynamic variables. This relation can be explicit, meaning that the value of a dynamic variable (such as force or stress) is given in terms of kinematic variables (such as displacement or strain). In this case, the governing differential equations of the mechanical system will be in the form of Ordinary or Partial Differential Equations (ODEs or PDEs). However, another possible case is the relation between kinematic and dynamic variables to be implicit, hence, the value of dynamic variable might not be a function (one-to-one) relation in terms of kinematic variables

(Rajagopal, 2003). For example, one might have kinematic variables in terms of dynamic variables. An example of these material models is the Bingham model (Bingham, 1922; Perzyna, 1966). This case, renders the governing equations of the mechanical system as Differential-Algebraic Equations (DAEs) or Partial Differential-Algebraic Equations (PDAEs) (Darbha et al., 2010). Obviously, these equations may not be solved exactly and one needs to resort to numerical approximations.

Historically, numerical methods have been developed for ordinary (partial) differential equations. Classical examples of time-integration algorithms are designed to perform well when applied to ODEs. The efforts by researchers in this case has been tremendously successful and forms the basis for most simulations in scientific computing. However, numerical integration of DAEs is still in a state of transition and development. The reason why numerical integration of DAEs can be more challenging than ODEs is that the presence of algebraic variables or relations can invoke instabilities (Brenan et al., 1995; Wanner and Hairer, 1991)

. In fact, a numerical time-integrator can be successful in approximating the solution to ODEs, but to fail in estimating the solution to a DAE. Considering the importance of DAEs to mathematical modeling of mechanical systems such as multi-body dynamics

(Schiehlen, 1997; Brüls and Cardona, 2010), domain decomposition for PDEs (Bursi et al., 2013; Karimi and Nakshatrala, 2014, 2015), as well as implicit constitutive relations in mechanical systems (Rajagopal, 2010; Darbha et al., 2010; Pražák and Rajagopal, 2012) and many others, many research endeavors have been centered around development of reliable numerical integration of DAEs for various applications. It is worthy to note that the algebraic constraints or relations in a DAE can be of different types. The algebraic constraint can be smooth (differentiable) in which case classical ODE solvers need minor adjustments to give accurate results. However, as the case in this paper, the algebraic constraint (the implicit constitutive relation) is continuous but non-smooth (not differentiable). Therefore, numerical time-integration needs to be designed with more care.

In this paper, time-integration of a lumped system (see figure 1) governed by the implicit Bingham-Kelvin constitutive relation is studied. The reason why we have chosen this rather simple system is to examine the basic performance of the derived time-integration methods while avoiding other numerical complexities. The system considered in this paper, has a non-smooth response to excitations. Based on findings in (Darbha et al., 2010), and generalized trapezoidal integration (Wood, 1990; Hughes, 2012), we derive an integration method that includes two independent parameters. A complete description of numerical integration under the proposed methods for cases of linear and nonlinear Bingham models is given. Finally, the performance of the derived methods is compared to the results from algorithm given in (Darbha et al., 2010) using several numerical examples.

### 1.1. Outline of the paper

The rest of the paper is organized as follows. In section 2, the governing equations and the new numerical scheme is detailed and is followed by numerical examples in section 3. In section 4 the conclusions are drawn and suggestions for future research are given.

## 2. Method for numerical Integration

### 2.1. Model problem and governing equations

Consider the SDOF system shown in figure 1, where the block of mass is connected to a linear elastic spring of stiffness and a dashpot with Bingham constitutive relation. The governing equations for the motion of the block can be written as:

 ˙v(t)=1m(fext(t)−fs(t)−fd(t)) (2.1a) fs(t)=ku(t) (2.1b) v(t)=φ(fd(t)) (2.1c) ˙u(t)=v(t). (2.1d)

where and are the displacement and velocity of the block. The force of dashpot and spring are shown by and , respectively. The function is defined as follows:

 φ(f):={0∣∣f∣∣≤fyγ(∣∣f∣∣−fy)NSign[f]∣∣f∣∣≥fy, (2.2)

where , and are constitutive parameters. Equations (2.1) and (2.2) form a DAE system of equations. In this problem, the velocity is given as a function of the force of dashpot, hence, it is referred to as an implicit constitutive relation. Figure 2 shows the relation between the dashpot force and the velocity for both and . From figure 2 it can be observed that the function is a nondecreasing function. Herein, inspired by (Pražák and Rajagopal, 2012), we show that the solution to the equation (2.1) is unique.

Proposition. Consider initial conditions and and let be a nondecreasing function, then the solution to equation (2.1) is unique.

Proof. Suppose and are both solutions to equation (2.1), with the initial conditions and . Consider the arbitrary time-interval . Since the initial conditions for both solutions and are the same, then for . From equation (2.1) we can derive:

 ¨u−¨~u+1m(fs−~fs)+1m(fd−~fd)=0, (2.3)

where and . Respectively, and correspond to the dashpot force for values of and . Multiplying equation (2.3) by gives

 2(˙u−˙~u)(¨u−¨~u)+2m(˙u−˙~u)(fs−~fs)+2m(˙u−˙~u)(fd−~fd)=0. (2.4)

Note that the nondecreasing assumption on of equation (2.2) gives

 (˙u−˙~u)(fd−~fd)≥0. (2.5)

Denoting and using relations (2.4) and (2.5) result in

 d˙w2dt≤−k2dw2dt. (2.6)

Note that based on definition of we have and . Integrating (2.6) over time interval gives

 (˙w(t=δ))2+k2(w(t=δ))2≤0. (2.7)

The only way relation (2.7) can be satisfied is to have

 ˙w(t=δ)=0andw(t=δ)=0. (2.8)

This result shows that the value of the solution to (2.1) for a give initial condition is unique.

Note that the proposition given above applies to both linear and nonlinear Bingham models as in both cases is a nondecreasing function.

Here, we will differentiate equation (2.1)b with respect to time in order to reduce the index of this DAE system. This reduction in DAE index results in a more stable numerical integration. Doing so, results in the following system of equations:

 ˙v(t)=1m(fext(t)−fs(t)−fd(t)) (2.9a) ˙fs(t)=kv(t) (2.9b) v(t)=φ(fd(t)). (2.9c)

Arguments on existence and uniqueness of a solution to this system can be found in (Darbha et al., 2010; Pražák and Rajagopal, 2012), as well as using the uniqueness agument given above. In the following section, the efforts in reference (Darbha et al., 2010) will be extended and applied to solve the system given in (2.9).

### 2.2. Numerical time integration

Suppose that the time interval of interest is and is divided into sub-intervals of equal size (time-step) . The value of variable at time level will be shown as . We will begin by integrating the linear momentum balance equation:

 ∫t(n+1)t(n)˙v(t)dt=1m∫t(n+1)t(n)(fext(t)−fs(t)−fd(t))dt. (2.10)

Using trapezoidal integration, equation (2.10) results in:

 v(n+1)−v(n)=Δtm((1−α)(f(n)ext−f(n)s−f(n)d)+α(f(n+1)ext−f(n+1)s−f(n+1)d)), (2.11)

where is the integration parameter (). The equation (2.9)b can be integrated in a similar fashion to give:

 f(n+1)s−f(n)s=kΔt((1−β)v(n)+βv(n+1)), (2.12)

where is an integration parameter in . Equations (2.11) and (2.12) can be combined to give:

 (1+αβΔt2km)v(n+1)=αΔtm(^f(n+1)−f(n+1)d), (2.13)

where is defined as:

 ^f(n+1)=f(n+1)ext+(1α−1)f(n)ext−1αf(n)s−(1α−1)f(n)d+(mαΔt−kΔt(1−β))v(n). (2.14)

The constitutive relation for dashpot given in equation (2.2) shows us that at all time levels. Also, from equation (2.13) we can derive that , because

 ((1+αβΔt2km)∣∣v(n+1)∣∣+αΔtm∣∣f(n+1)d∣∣)≥0Sign[v(n+1)]=αΔtm∣∣^f(n+1)≥0∣∣Sign[^f(n+1)]. (2.15)

Equation (2.15) indicates that the value of can be used as a predictor to determine the relation between and . This equation implies that if , then we also have . Therefore, if then and from equation (2.13) we get . Otherwise, the force of dashpot needs to be solved for using the following equation:

 (2.16)

Equation (2.16) is obtained from equations (2.13) and (2.2). Obviously, this equation is in general a nonlinear equation. In the case of , equation (2.16) can be simplified as:

 f(n+1)d=11+αγm(1Δt+αβΔtkm)⎛⎜ ⎜⎝αγm(1Δt+αβΔtkm)^f(n+1)+Sign[^f(n+1)]fy⎞⎟ ⎟⎠. (2.17)

Hence a nonlinear solver is not required. Once the value of the is found, the value of can be found to be

 v(n+1)=γ(∣∣f(n+1)d∣∣−fy)NSign[f(n+1)d]. (2.18)

The value of the spring force can be updated using equation (2.12) and the displacement of the mass can be found from the following relation:

 u(n+1)=1kf(n+1)s. (2.19)

To complete the description of time-stepping algorithm, we need to clarify two points:

1. Initialization: In this paper, we will assume all the initial displacements to be purely elastic. In other words, given an initial displacement , the initial spring force will be taken to be and initial dashpot force should be found from

 (2.20)

if initial velocity is nonzero. Note that . Consequently, equation (2.20) for the case of gives:

 fd,0=(fy+∣∣v0∣∣γ)Sign[v0]. (2.21)
2. Initial guess for nonlinear solver: Equation (2.16) is a nonlinear equation with respect to . Hence, as with any nonlinear solver, an initial guess is needed. Based on our numerical experiments, the initial guess of seems to be a reasonable choice.

This concludes the description of the new time-stepping algorithm. A pseudocode for this method is given in algorithm 1. In the following section, numerical examples will be presented to demonstrate the effect of various time-steps and time-integration parameters and .

## 3. Numerical Results

In this section, we will showcase the performance of the derived method in comparison to the method proposed in (Darbha et al., 2010). The amount of damped (dissipated) energy will be denoted by and is defined as

 Ed(t):=∫t0v(t)fd(t)dt. (3.1)

Note that because and have the same sign, the value of will be nonnegative at all time-levels. The error in kinematic quantities will be denoted by for , and is defined as

 ep:=1M ⎷M∑i=1(p(i)−p(i)ref)2, (3.2)

where is the total number of time-steps and is the reference (benchmark) values. The numerical results are obtained from implementation in MATLAB (Mathworks, 2017).

### 3.1. The case of N=1

In this problem, consider the SDOF system shown in figure 1 subject to external force given as:

 fext(t)=2sin(2πt)e(−0.2t). (3.3)

The initial conditions are and . The physical parameters for this problem are given in Table 1. To demonstrate the effect of time-integration parameters and , we use different combinations of values of either 1 or 1/2 for each of them. Time-integration parameters used in all cases are given in Table 2. Note that the benchmark in Table 2 refers to the method proposed in (Darbha et al., 2010), which coincides with and under the method given in this paper. All the numerical solutions will be compared to that of benchmark method.

The numerical results for displacement and velocity, against time are shown in figures 3. The dashpot force and spring force is shown in figure 4. It can be seen that the numerical results from the benchmark solution and case 1 (i.e., and ) are very different from the cases 2 and 3. The numerical solution from case 1 is indeed very close to the solution from the benchmark method. Repeating the same simulation with cases 2 and 3 of Table 2 with smaller time-steps show that the damping properties in values of dashpot force do not improve (see figure 5). In fact the numerical values are very close to the ones with time-step of . An important physical quantity is the energy damping in equation (3.1), which shown in figure 6 for benchmark case and case 1 of Table 2. The amount of energy dissipated calculated from each of these cases are very similar, while cases 3 and 4 suggest zero energy damping. The solution from the benchmark method and case 1 converge to one another in order , shown in figure 7. However, cases 2 and 3 underestimate the dashpot force resulting in inaccurate values for displacement and velocity. In fact, the value of in cases 2 and 3 seems to dampen excessively.

From this numerical experiment, we find that among the choices for parameters and , the values for case 1 and the benchmark method are the only viable ones. It seems that choosing , introduces excessive damping to the value of . Due to the non-smooth nature of the constitutive relation considered in this paper, this property results in wrong values for other variables.

### 3.2. The case of N≠1

In this case, the value of dashpot force is the solution of the nonlinear equation (2.16). For this case, the physical parameters are given in Table 3. Time-integration parameters are given in Table (4). Here, we take the value of and the benchmark numerical solution is from the method in (Darbha et al., 2010). The external force is the same as in equation (3.3), and the initial displacement and velocity are both zero.

The numerical results are given in figures 89. Similar to the case of , the benchmark numerical solution and case 1 are close and compatible with one another. However, the parameters selected for cases 2 and 3 are very different and predict a value of zero throughout the time-interval of interest for velocity and displacement. In fact, in cases where the smoothing results in damping in values of . It is worth noting that unlike the case of , the displacement seems to converge to a point other than as time proceeds. In the case of , the amount of energy damped shown in figure 10 suggests a similar conclusion to the case of . In this example, the energy dissipation predicted by case 1 is very similar to the benchmark solution. Other cases of and predict zero energy dissipation in this system.

This experiment concludes that the only viable value for under the proposed method is 1. The numerical solution to this nonlinear and non-smooth system is sensitive to smoothing effects of trapezoidal integration in linear momentum balance equation. The numerical results seem to be much less sensitive to the choice of between 1/2 and 1.

### 3.3. Implicit-explicit time-integration

Thus far, the numerical experiments suggest that the only viable value for is 1. However, for , the numerical results for various values of seem to be very close. In this section, we showcase the numerical results for the case of and . For , the equation (2.12) will become an explicit Euler integration scheme while equation (2.11) remains implicit Euler integration. In this case, the constitutive relation of the elastic spring is integrated using an explicit scheme and the momentum balance is integrated via an implicit relation.

The problem is solved for the external force in equation (3.3) using time-steps and . The error (compared to the benchmark numerical result) is shown in figure 11. It can be seen that the numerical solutions from the two methods converge as , which implies that the choice of and is a viable choice of time-integration parameters. This method, therefore, does not suffer from excessive numerical damping for values of dashpot force .

## 4. Conclusion

In this paper the governing differential-algebraic equations for a SDOF system obeying Bingham constitutive model were reviewed and a proof for uniqueness of the solution was presented. Using the generalized trapezoidal integration a method for numerical integration of the given equations was derived with two independent time-integration parameters. A pseudocode for this numerical algorithm along with initialization and initial guess for the nonlinear equations was provided. This derivation was followed by numerical examples showcasing its performance compared to a benchmark method. The findings from these numerical experiments can be summarized as below:

1. Because of the non-smooth nature of the physical system considered in this paper, accuracy in estimation of dashpot force is of utmost importance. Choice of time-integration parameters can lead to excessive numerical damping in value of dashpot force , which in turn gives completely different numerical values for other kinematic and dynamic variables.

2. Under the presented numerical scheme in this paper, the only viable choice of integration parameter is 1. Other values of result in additional (unphysical) smoothing in the numerical solution.

3. Using the method of numerical time integration in this paper the value of integration parameter can be anything in range . It seems that the values of do not introduce any more significant damping to the estimated value of dashpot force.

4. Given and under the numerical scheme here, the numerical solution converges to the solution from the benchmark method in order .

5. Under the proposed method in this paper, one can use the and values for integration parameters. This means that one can have a convergent implicit-explicit time-integration of the non-smooth system considered here.

Future developments in this area can be extension of such methods to other implicit and non-smooth material models Maxwell type (Christensen, 1982), extension to 2 and 3-dimensional problems, and extension of Newmark-type integration methods (Newmark, 1959; Chung and Hulbert, 1993) for this class of problems.

## References

• Bingham  E. C. Bingham. Fluidity and Plasticity. McGraw-Hill, New York, 1922.
• Brenan et al.  K. E. Brenan, S. L. Campbell, and L. R. Petzold. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. SIAM, 1995.
• Brüls and Cardona  O. Brüls and A. Cardona. On the use of Lie group time integrators in multibody dynamics. Journal of Computational and Nonlinear Dynamics, 5(3):031002, 2010.
• Bursi et al.  O. S. Bursi, Z. Wang, C. Jia, and B. Wu. Monolithic and partitioned time integration methods for real-time heterogeneous simulations. Computational Mechanics, pages 1–21, 2013.
• Christensen  R. M. Christensen. Theory of Viscoelasticity: An Introduction. Academic Press, New York, 1982.
• Chung and Hulbert  J. Chung and G. M. Hulbert. A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized- Method. Journal of Applied Mechanics, 60:371–375, 1993.
• Darbha et al.  S. Darbha, K. B. Nakshatrala, and K. R. Rajagopal. On the vibrations of lumped parameter systems governed by differential-algebraic equations. Journal of the Franklin Institute, 347(1):87–101, 2010.
• Hughes  T. J. R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, Inc., Mineola, 2012.
• Karimi and Nakshatrala  S. Karimi and K. B. Nakshatrala. On multi-time-step monolithic coupling algorithms for elastodynamics. Journal of Computational Physics, 273:671–705, 2014.
• Karimi and Nakshatrala  S. Karimi and K. B. Nakshatrala. A monolithic multi-time-step computational framework for first-order transient systems with disparate scales. Computer Methods in Applied Mechanics and Engineering, 283:419–453, 2015.
• Mathworks  Mathworks. MATLAB User’s Guide. Natick, MA, 2017.
• Newmark  N. M. Newmark. A method of computation for structural dynamics. ASCE Journal of Engineering Mechanics, 85(3):67–94, 1959.
• Perzyna  P. Perzyna. Fundamental problems in viscoplasticity. Advances in Applied Mechanics, 9:243–377, 1966.
• Pražák and Rajagopal  D. Pražák and K. R. Rajagopal. Mechanical oscillators described by a system of differential-algebraic equations. Applications of Mathematics, 57(2):129–142, 2012.
• Rajagopal  K. R. Rajagopal. On implicit constitutive theories. Applications of Mathematics, 48(4):279–319, 2003.
• Rajagopal  K. R. Rajagopal. A generalized framework for studying the vibrations of lumped parameter systems. Mechanics Research Communications, 37(5):463–466, 2010.
• Schiehlen  W. Schiehlen. Multibody system dynamics: roots and perspectives. Multibody System Dynamics, 1(2):149–188, 1997.
• Wanner and Hairer  G. Wanner and E. Hairer.

Solving Ordinary Differential Equations: II. Stiff and Differential-Algebraic Problems

.
Springer-Verlag, Berlin Heidelberg, 1991.
• Wood  W. L. Wood. Practical Time-Stepping Schemes. Oxford University Press, Oxford, 1990.