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 (onetoone) 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 DifferentialAlgebraic Equations (DAEs) or Partial DifferentialAlgebraic 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 timeintegration 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 timeintegrator 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 multibody 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 nonsmooth (not differentiable). Therefore, numerical timeintegration needs to be designed with more care.In this paper, timeintegration of a lumped system (see figure 1) governed by the implicit BinghamKelvin constitutive relation is studied. The reason why we have chosen this rather simple system is to examine the basic performance of the derived timeintegration methods while avoiding other numerical complexities. The system considered in this paper, has a nonsmooth 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
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:
(2.1a)  
(2.1b)  
(2.1c)  
(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:
(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 timeinterval . Since the initial conditions for both solutions and are the same, then for . From equation (2.1) we can derive:
(2.3) 
where and . Respectively, and correspond to the dashpot force for values of and . Multiplying equation (2.3) by gives
(2.4) 
Note that the nondecreasing assumption on of equation (2.2) gives
(2.5) 
Denoting and using relations (2.4) and (2.5) result in
(2.6) 
Note that based on definition of we have and . Integrating (2.6) over time interval gives
(2.7) 
The only way relation (2.7) can be satisfied is to have
(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:
(2.9a)  
(2.9b)  
(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 subintervals of equal size (timestep) . The value of variable at time level will be shown as . We will begin by integrating the linear momentum balance equation:
(2.10) 
Using trapezoidal integration, equation (2.10) results in:
(2.11) 
where is the integration parameter (). The equation (2.9)b can be integrated in a similar fashion to give:
(2.12) 
where is an integration parameter in . Equations (2.11) and (2.12) can be combined to give:
(2.13) 
where is defined as:
(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
(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:
(2.17) 
Hence a nonlinear solver is not required. Once the value of the is found, the value of can be found to be
(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:
(2.19) 
To complete the description of timestepping algorithm, we need to clarify two points:

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:
(2.21) 
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 timestepping 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 timesteps and timeintegration 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
(3.1) 
Note that because and have the same sign, the value of will be nonnegative at all timelevels. The error in kinematic quantities will be denoted by for , and is defined as
(3.2) 
where is the total number of timesteps and is the reference (benchmark) values. The numerical results are obtained from implementation in MATLAB (Mathworks, 2017).
3.1. The case of
In this problem, consider the SDOF system shown in figure 1 subject to external force given as:
(3.3) 
The initial conditions are and . The physical parameters for this problem are given in Table 1. To demonstrate the effect of timeintegration parameters and , we use different combinations of values of either 1 or 1/2 for each of them. Timeintegration 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 timesteps 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 timestep 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 nonsmooth nature of the constitutive relation considered in this paper, this property results in wrong values for other variables.
Parameter  

Value  1  100  1  1  1 
Case  

benchmark  1  1  
1  1  1/2  
2  1/2  1  
3  1/2  1/2 
3.2. The case of
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. Timeintegration 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 8–9. 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 timeinterval 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 nonsmooth 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.
Parameter  

Value  1  10  1  1  3 
Case  

benchmark  1  1  
1  1  1/2  
2  1/2  1  
3  1/2  1/2 
3.3. Implicitexplicit timeintegration
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 timesteps 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 timeintegration parameters. This method, therefore, does not suffer from excessive numerical damping for values of dashpot force .
4. Conclusion
In this paper the governing differentialalgebraic 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 timeintegration 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:

Because of the nonsmooth nature of the physical system considered in this paper, accuracy in estimation of dashpot force is of utmost importance. Choice of timeintegration 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.

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.

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.

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

Under the proposed method in this paper, one can use the and values for integration parameters. This means that one can have a convergent implicitexplicit timeintegration of the nonsmooth system considered here.
Future developments in this area can be extension of such methods to other implicit and nonsmooth material models Maxwell type (Christensen, 1982), extension to 2 and 3dimensional problems, and extension of Newmarktype integration methods (Newmark, 1959; Chung and Hulbert, 1993) for this class of problems.
References
 Bingham [1922] E. C. Bingham. Fluidity and Plasticity. McGrawHill, New York, 1922.
 Brenan et al. [1995] K. E. Brenan, S. L. Campbell, and L. R. Petzold. Numerical Solution of InitialValue Problems in DifferentialAlgebraic Equations. SIAM, 1995.
 Brüls and Cardona [2010] 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. [2013] O. S. Bursi, Z. Wang, C. Jia, and B. Wu. Monolithic and partitioned time integration methods for realtime heterogeneous simulations. Computational Mechanics, pages 1–21, 2013.
 Christensen [1982] R. M. Christensen. Theory of Viscoelasticity: An Introduction. Academic Press, New York, 1982.
 Chung and Hulbert [1993] 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. [2010] S. Darbha, K. B. Nakshatrala, and K. R. Rajagopal. On the vibrations of lumped parameter systems governed by differentialalgebraic equations. Journal of the Franklin Institute, 347(1):87–101, 2010.
 Hughes [2012] T. J. R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, Inc., Mineola, 2012.
 Karimi and Nakshatrala [2014] S. Karimi and K. B. Nakshatrala. On multitimestep monolithic coupling algorithms for elastodynamics. Journal of Computational Physics, 273:671–705, 2014.
 Karimi and Nakshatrala [2015] S. Karimi and K. B. Nakshatrala. A monolithic multitimestep computational framework for firstorder transient systems with disparate scales. Computer Methods in Applied Mechanics and Engineering, 283:419–453, 2015.
 Mathworks [2017] Mathworks. MATLAB User’s Guide. Natick, MA, 2017.
 Newmark [1959] N. M. Newmark. A method of computation for structural dynamics. ASCE Journal of Engineering Mechanics, 85(3):67–94, 1959.
 Perzyna [1966] P. Perzyna. Fundamental problems in viscoplasticity. Advances in Applied Mechanics, 9:243–377, 1966.
 Pražák and Rajagopal [2012] D. Pražák and K. R. Rajagopal. Mechanical oscillators described by a system of differentialalgebraic equations. Applications of Mathematics, 57(2):129–142, 2012.
 Rajagopal [2003] K. R. Rajagopal. On implicit constitutive theories. Applications of Mathematics, 48(4):279–319, 2003.
 Rajagopal [2010] K. R. Rajagopal. A generalized framework for studying the vibrations of lumped parameter systems. Mechanics Research Communications, 37(5):463–466, 2010.
 Schiehlen [1997] W. Schiehlen. Multibody system dynamics: roots and perspectives. Multibody System Dynamics, 1(2):149–188, 1997.

Wanner and Hairer [1991]
G. Wanner and E. Hairer.
Solving Ordinary Differential Equations: II. Stiff and DifferentialAlgebraic Problems
. SpringerVerlag, Berlin Heidelberg, 1991.  Wood [1990] W. L. Wood. Practical TimeStepping Schemes. Oxford University Press, Oxford, 1990.
Comments
There are no comments yet.