An Improved Dynamic Contact Model for Mass–Spring and Finite Element Systems Based on Parametric Quadratic Programming Method

Abstract

An improved dynamic contact model for mass-spring and finite element systems is proposed in this paper. The proposed model avoids the numerical troubles of spurious high-frequency oscillations for mass-spring and finite element systems in dynamic contact problems by using the parametric quadratic programming technique. The iterative process for determination of contact states are not required for each time step in the proposed method, as the contact states are transformed into the base exchanges in the solution of a standard quadratic programming problem. The proposed methodology improves stability and has good convergence behavior for dynamic contact problems. Numerical results demonstrate the validity of the proposed method.

Keywords
Dynamic Contact Problem; Finite Element Method; Mass-Spring Model; Parametric Quadratic Programming; Linear Complementarity

1 INTRODUCTION

Various types of impact problems are encountered in engineering structural analysis ( Balestrino et al., 2010 ; Rabb and Fahrenthold, 2010 ; and Auersch, 2005 ). It will lead to a heavy computational burden for dynamic contact analysis, as the solution of a discretization non-linear system is required at each integration time step ( Liu and Wang, 2007 ; Padmanabhan and Laursen, 2001 ). Also, the analysis of dynamic contact problems present strong nonlinear characteristic, thus, efficiency and instability of the algorithm for dynamic contact problems are challenges. In view of its great significance of shock or impact analysis in civil and aerospace engineering, it is necessary to develop an efficient and stable numerical method for dynamic contact or impact problems.

Various numerical methods are used to solve contact problems including the Lagrange multiplier methods ( Chen and Yeh, 1990 ; Hu, 1997 ; Laursen and Chawla, 1997 ; Jiang and Rogers, 1988 ; Simo et al., 1985 ), penalty methods ( Kanto and Akawa, 1990 ; Peric and Owen, 1992 ; Chen et al., 2015 ), contact element methods ( Dowling and Hall, 1989 ), contact force methods ( Wang et al., 2007 ), etc. Most conventional numerical methods for dynamic contact first assume the contact status, then solve the problem, and finally check if the assumption is correct. The methods iteratively solve the problem by trial-and-error, thus requiring significant computation. These iterative methods are vulnerable to convergence behavior. In a contact analysis involving dynamics, shock loading is often encountered if two moving bodies hit one another. Depending on the velocity magnitudes, material properties, etc., such shock loading can trigger spurious high-frequency oscillations, which may cause numerical troubles. The Newmark method is the most popular implicit scheme for time-integration, but it results in excessive numerical oscillations for dynamic contact problems ( Kanto and Akawa, 1990 ). Although the Newmark method is unconditionally stable, it may still result in spurious high-frequency modes for small time steps. When Newmark method is used for time-integration in dynamic contact problems, the specific damping properties must be known in order to get rid of the undesired oscillations ( Kim and Kwak, 1996 ). However, the choice of damping properties is not trivial, usually also causing the damping-out of relevant oscillations.

The objective of this paper is to develop a stable and efficient numerical dynamic contact model for mass–spring and finite element systems. The proposed method does not depend on displacement iteration, but on the base exchanges in the solution of a parametric quadratic programming (PQP) problem. PQP algorithms can significantly weaken the non-linear effects on the numerical computation, and the computational process can thus be simplified ( Zhang et al., 2012 ). Thus, the dynamic contact model and method developed present good convergence behavior in computation. This paper is organized as follows: Section 2 describes the basic theory of PQP-based dynamic contact model. It then presents the improved PQP-based technique for dynamic contact model. Section 3 provides numerical examples that demonstrate the validity of the proposed method. Finally, important conclusions are summarized in Section 4.

2 Numerical model for dynamic contact problem

2.1 Basic formulation of PQP-based dynamic contact model

For the sake of simplicity, the discussion is started with a simple mass-spring system with two degree-of-freedoms (DOF) as shown in Fig. 1 . The dynamic equations for the discretized system are as follow

$$m_1{\ddot{x}}_1+c_1{\dot{x}}_1+k_1{x}_1=f_1+{\lambda }_p$$ (1)

$$m_2{\ddot{x}}_2+c_2{\dot{x}}_2+k_2{x}_2=f_2-{\lambda }_p$$ (2)

where $m_1$ , $c_1$ , $k_1$ and $f_1$ is the mass, dump, stiffness and external force of the first DOF; $m_2$ , $c_2$ , $k_2$ and $f_2$ is the mass, dump, stiffness and external force of the second DOF. ${\lambda }_p$ is the contact force.

Figure 1
Oscillations results of PQP-based dynamic contact model for mass-spring system with two DOFs without external force. The mass and stiffness are 1kg and 1N/m for both two DOFs. The initial conditions are $x_1\left(0\right)=1$ , ${\dot{x}}_1\left(0\right)=1$ for DOF of 1, and $x_2\left(0\right)=0$ , ${\dot{x}}_2\left(0\right)=0$ for DOF of 2.

Several different time-stepping methods can be used for Eqs. (1) and (2) , For example, New mark method can be utilized as such:

$$x_1^{n+1}=A_1{\lambda }_p^{n+1}+B_1^{n+1}$$ (3)

$$x_2^{n+1}=-A_2{\lambda }_p^{n+1}+B_2^{n+1}$$ (4)

$$A_1={\left(k_1+\frac{1}{\beta \Delta t^2}{m}_1+\frac{\gamma }{\beta \Delta t}{c}_1\right)}^{-1}$$ (5)

$$A_2={\left(k_2+\frac{1}{\beta \Delta t^2}{m}_2+\frac{\gamma }{\beta \Delta t}{c}_2\right)}^{-1}$$ (6)

$$B_1^{n+1}=\frac{f_1^{n+1}+m_1\left(\frac{1}{\beta \Delta t^2}{x}_1^n+\frac{\gamma }{\beta \Delta t}{\dot{x}}_1^n+\frac{1-2\beta }{2\beta }{\ddot{x}}_1^n\right)+c_1\left(\frac{\gamma }{\beta \Delta t}{x}_1^n+\frac{\gamma -\beta }{\beta }{\dot{x}}_1^n+\Delta t\frac{\gamma -2\beta }{2\beta }{\ddot{x}}_1^n\right)}{k_1+\frac{1}{\beta \Delta t^2}{m}_1+\frac{\gamma }{\beta \Delta t}{c}_1}$$ (7)

$$B_2^{n+1}=\frac{f_2^{n+1}+m_2\left(\frac{1}{\beta \Delta t^2}{x}_2^n+\frac{\gamma }{\beta \Delta t}{\dot{x}}_2^n+\frac{1-2\beta }{2\beta }{\ddot{x}}_2^n\right)+c_2\left(\frac{\gamma }{\beta \Delta t}{x}_2^n+\frac{\gamma -\beta }{\beta }{\dot{x}}_2^n+\Delta t\frac{\gamma -2\beta }{2\beta }{\ddot{x}}_2^n\right)}{k_2+\frac{1}{\beta \Delta t^2}{m}_2+\frac{\gamma }{\beta \Delta t}{c}_2}$$ (8)

$$\gamma \ge \mathrm{0.5,}\beta \ge 0.25{\left(0.5+\gamma \right)}^2$$ (9)

where $\Delta t$ is the time step, and the superscript n means the n-th time step.

Then the contact status constraint should be

$$x_1^{n+1}-x_2^{n+1}+{\delta }_N\ge 0$$ (10)

where δ_Ν is the distance between two moving bodies.

Substituting Eqs. (3) and (4) into Eq.(10)

$$x_1^{n+1}-x_2^{n+1}+{\delta }_N=\left(A_1+A_2\right){\lambda }_p^{n+1}+\left(B_1^{n+1}-B_2^{n+1}+{\delta }_N\right)\ge 0$$ (11)

Let

$$\Delta x^{n+1}=x_1^{n+1}-x_2^{n+1}+{\delta }_N$$ (12)

$$M=A_1+A_2$$ (13)

$q^{n+1}=B_1^{n+1}-B_2^{n+1}+{\delta }_N$ (14)

The equations become

$$\begin{array}{l}\Delta x^{n+1}=M\cdot {\lambda }_p^{n+1}+q^{n+1}\\ \Delta x^{n+1}\ge \mathrm{0,}{\lambda }_p^{n+1}\ge \mathrm{0,}\Delta x^{n+1}\cdot {\lambda }_p^{n+1}=0\end{array}$$ (15)

Eq. (15) is a standard linear complementary problem (LCP) ( Zhang et al., 2005 ). Its physical significance is that, when the contact force is not equal to zero, the distance between the two moving bodies must be equal to zero, and vice versa. In mathematics, a quadratic programming problem can be equivalent to an LCP according to the Karush-Kuhn-Tucker (KKT) theorem. Several quadratic programming methods are available for this well-studied problem ( Zhang et al., 2005 ; Ferris and Pang, 1997 ; Fischer, 1992 ). Details of these methods will not be elaborated upon here due to space limitations.

Once the contact force ${\lambda }_p^{n+1}$ is solved through Eq. (15) , the displacements $x_1^{n+1}$ and $x_2^{n+1}$ can be obtained by Eqs. (3) and (4) . PQP-based dynamic contact model is similar to the numerical technique ( Kim and Kwak, 1996 ; Zhang et al., 2012 ), the method is potentially good algorithms for dynamic contact problems because it requires neither the assumptions of iterative process or contact states in each time step. These assumptions are not required because the method is based on the quadratic programming method in which the contact conditions are treated as a series of LCPs.

However, as shown in Fig. 1 , the PQP-based dynamic contact model suffers numerical problems of oscillations ( Kim and Kwak, 1996 ; Solberg and Papadopoulos, 1998 ), even the time step is very small for time integration. The dynamic contact condition of Eq. (16) is too strict for the time domain, leading to numerical troubles associated with spurious high-frequency oscillations. We will discuss a strategy to address this issue in next section.

2.2 Improved PQP-based technique for dynamic contact model

In performing dynamic contact analysis, a constrained minimization problem arises with the contact status constraint ( Taylor et al., 1991 ). The penalty method is one of the most elegant procedures to apply mathematical constraints to a system. The penalty function method can be utilized to introduce constraint conditions to the modified variational principle

$$\Pi ={\Pi }_u+{\Pi }_{CP}$$ (16)

where

$${\Pi }_u={\displaystyle {\int }_{\Omega }\frac{1}{2}}{u}^T{B}^{T}DBudv-{\displaystyle {\int }_{\Omega }{u}^T{N}^{T}f}dv$$ (17)

$${\Pi }_{CP}={\displaystyle {\int }_{S_c}\left({\alpha }_N{\left(u_N^A-u_N^B-{\delta }_N\right)}^2+{\alpha }_1{\left(u_1^A-u_1^B\right)}^2+{\alpha }_2{\left(u_2^A-u_2^B\right)}^2\right)}ds$$ (18)

where $\Omega$ is volume, $S_c$ is area of the contact surface, $D$ is elasticity matrix, $N$ is the shape function, $B$ is the strain matrix, $u$ is displacement, $f$ is the volume force. As shown in Fig. 2 , $u_N^A$ is the normal displacement of contact surface A, $u_N^B$ is the normal displacement of contact surface B. $u_1^A$ and $u_2^A$ are the tangential displacements of contact surface A, $u_1^B$ and $u_2^B$ are the tangential displacements of contact surface B. α_N, α1 and α2 are the penalties. Without considering friction, α1 and α2 are equal to zero. Usually, α_N is a very large positive number that is related to the dynamic contact condition: when α_N gets larger, the dynamic contact condition becomes stricter, and spurious, high-frequency oscillations occur.

Figure 2
The contact surface between body A and body B.

For the sake of simplicity, friction is not considered. Introducing the parametric variable, the contact states (stick and separation) in Eq. (19) can be unified as follows:

$${\Pi }_{CP}={\displaystyle {\int }_{S_c}{\alpha }_N{\left(u_N^A-u_N^B+{\delta }_N-\lambda \right)}^2}ds$$ (19)

$$\begin{array}{l}\lambda =u_N^A-u_N^B+{\delta }_N,if{u}_N^A-u_N^B>-{\delta }_N\\ \lambda =\mathrm{0,}if{u}_N^A-u_N^B\le -{\delta }_N\end{array}$$ (20)

Introducing the slack variable $v$ , the Eq.(21) can be rewritten as

$$\begin{array}{l}v-\lambda =-{\delta }_N-{\theta }^T{N}_{c}u\\ v\ge \mathrm{0,}\lambda \ge \mathrm{0,}v\cdot \lambda =0\end{array}$$ (21)

where $\theta$ is the transformation matrix of normal direction of the contact surface and global coordinate system, and $N_c$ is the shape function of the contact point in global coordinate system.

Eq. (22) is a standard LCP. Several highly efficient computational methods ( Zhang et al., 2005 ; Ferris and Pang, 1997 ; Fischer, 1992 ) from computational mathematics can be used to solve the final quadratic programming model. Therefore the trial-and-error and iterative process in each time step is not required in the proposed method.

The Lagrangian function of the spatial discretized system is

$$L\left(\dot{u},u\right)=T-\Pi =\frac{1}{2}{\dot{u}}^{T}M\dot{u}-\frac{1}{2}{u}^{T}Ku+u^{T}F+u^{T}h\lambda$$ (22)

where

$$M={\displaystyle {\int }_{\Omega }\rho N^T}Ndv$$ (23)

$$K={\displaystyle {\int }_{\Omega }{B}^{T}DBdv}+{\displaystyle {\int }_{S_c}{\alpha }_N}{N}_c^T\theta {\theta }^T{N}_{c}ds$$ (24)

$$F={\displaystyle {\int }_{\Omega }{N}^{T}f}dv-{\displaystyle {\int }_{S_c}{\alpha }_N{N}_c^T\theta {\delta }_{N}ds}$$ (25)

$$h={\displaystyle {\int }_{S_c}{\alpha }_N}{N}_c^T\theta ds$$ (26)

where $\rho$ is the material density.

The equation of motion can be obtained by using the Lagrangian equation

$$M\ddot{u}+Ku=F+h\lambda$$ (27)

The formulation for analysis of impact problems can be established. The formulations are summarized as

$$\{\begin{array}{c}M\ddot{u}+Ku=F+h\lambda \\ v-\lambda =-{\delta }_N-{\theta }^T{N}_{c}u\\ v\ge \mathrm{0,}\lambda \ge \mathrm{0,}v\cdot \lambda =0\end{array}$$ (28)

The proposed formulation (29) combines Eq. (28) with the linear complementary problem of Eq. (21) , yielding the displacement for each time step from a unified formula. Eq. (28) can be solved by the follow formulation:

$$u^{n+1}=Ah{\lambda }^{n+1}+B^{n+1}$$ (29)

where

$$A={\left(K+\frac{1}{\beta \Delta t^2}M\right)}^{-1}$$ (30)

$$B^{n+1}=\frac{F^{n+1}+M\left(\frac{1}{\beta \Delta t^2}{u}^n+\frac{\gamma }{\beta \Delta t}{\dot{u}}^n+\frac{1-2\beta }{2\beta }{\ddot{u}}^n\right)}{K+\frac{1}{\beta \Delta t^2}M}$$ (31)

Furthermore, substituting Eq. (30) into Eq. (22) yields

$$\begin{array}{l}{v}^{n+1}+\left({\left({\theta }^{n+1}\right)}^T{N}_c^{n+1}Ah-I\right){\lambda }^{n+1}=-{\delta }_N^{n+1}-{\left({\theta }^{n+1}\right)}^T{N}_c^{n+1}{B}^{n+1}\\ v^{n+1}\ge \mathrm{0,}{\lambda }^{n+1}\ge \mathrm{0,}{v}^{n+1}\cdot {\lambda }^{n+1}=0\end{array}$$ (32)

And, eigenvalue analysis is performed to estimate a proper time step.

3 RESULTS

3.1 Impact of circular plate simulated by finite element system

A circular plate with a thickness of 2.5×10^-3 m and radius of 5×10 ^-2 m is clamped around its circumference and hit by a ball with a radius of 2×10 ^-2 m. The initial velocity of the ball is 2.5 m/s. Both the plate and the ball are considered to be elastic. The plate has a Young’s modulus of 7×10⁹ N/m², a Poisson’s ratio of 0.3 and a density of 2.5×10³ kg/m3. The ball has a Young’s modulus of 2×10¹¹ N/m² , a Poisson’s ratio of 0.3 and a density of 7.8×10³ kg/m³ . The plate and the ball are modeled using 3-node axis-symmetric elements with full integration. The mesh of the plate is refined around the contact area. A total time of 6×10 ^-3 s are analyzed in 4000 equally sized steps, corresponding to a time step of 1.5×10 ^-6 s.

Fig. 3 shows the model of the ball and plate. The finite element mesh and boundary conditions for the model are shown in Fig. 3(a) , and the equivalent Mises stress in deformed configuration at increment 1600 is plotted in Fig. 3(b) . In order to verify the validity of the proposed method, the problem was solved using the MSC.Marc software. Fig. 4 shows the displacement of the center of the circular plate versus time, obtained by the proposed method and benchmark MSC.Marc software with refined mesh. The results agree well. Table 1 shows the computational cost for the proposed method and MSC.Marc software. It indicates that the proposed method is more efficient than MSC.Marc software.

Figure 3
The model of a Ball on a Plate (a) Finite element mesh and boundary conditions for impact of a ball on a plate problem. (b) Equivalent Mises stress in deformed configuration at increment 1600.

Figure 4
The displacements of the center of circular plate versus time obtained by the proposed method and benchmark MSC.Marc software.

Table 1
Computational cost of the proposed method and MSC.Marc software

Method	Time step(ns)	Computational cost(s)
The proposed method	15.703	2534
MSC.Marc	15.703	3225

3.2 Simply supported beam crossed by moving mass-spring system

It is of great interest in vehicle and bridge interaction problem ( Stăncioiu et al., 2008 ; Baeza and Ouyang, 2008 ; Nguyena et al., 2009 ) and it can be modeled as a simply supported beam crossed by moving oscillators. For determining the dynamic response of a structure, the two most important factors are inertia and stiffness. The complexity of dynamic contact/impact problems is usually due to large deformations of the structure as well as nonlinear properties of the material. Therefore, in order to reveal the fundamental features of structural dynamic contact/impact with reduced complexity, various mass-spring systems have been widely employed in models ( Ruan and Yu, 2005 ; Ouyang, 2011 ).

Fig. 5 shows a simply supported beam carrying a moving mass-spring system. This model is analyzed by the PQP-based dynamic contact model and the proposed method in the paper. In this model, the gravity acceleration is considered as 9.8 m/s². The mass of the first DOF is m₁=1.8×10⁴ kg and the mass of the second DOF is m₂ =7.2×10³ kg. The stiffness and the damping of the mass-spring system are k=1.012×10⁸ N/m and c=1.3496×10⁶ Ns/m, respectively. The mass-spring system is moving along the beam at a constant velocity V= 114.5617 m/s. The bending stiffness of the beam is EJ=7.48×10¹⁰ Nm², the mass per unit length of the beam is A=10⁴ kg/m, the length of the beam is L=60 m. Fig.6 shows the time history of the displacement of middle point in the beam obtained by PQP-based dynamic contact model and the proposed method. The results agree well with the solution of the proposed method and the PQP-based dynamic contact model. Fig.7 shows the time history of beam-oscillator interaction forces obtained by PQP-based dynamic contact model and the proposed method. The high-frequency oscillation is observed for the PQP-based dynamic contact model, while the proposed method can solve this problem without numerical oscillations. Although the artificial damping is required to get rid of the undesired oscillations for PQP-based dynamic contact model, the choice of artificial damping properties is not trivial, and usually also relevant oscillations may be damped out, and it often leads to a large computational error. Table 2 shows the computational cost for the proposed method and PQP-based dynamic contact model. It indicates that the proposed method is more stable and efficient than PQP-based dynamic contact model with larger time step.

Figure 5
Simply supported beam crossed by moving mass-spring system with two DOFs

Figure 6
The time history of the displacement of middle point in the beam obtained by PQP-based dynamic contact model and the proposed method

Figure 7
The time history of the beam-oscillator interaction forces obtained by PQP-based dynamic contact model and the proposed method

Table 2
Computational cost of the proposed method and PQP-based dynamic contact model

Method	Time step(ms)	Computational cost(s)
The proposed method	5.2370	0.0449
PQP-based model	0.5237	0.1480

3.3 Dynamic behaviors of safety valve disc–seat impact at low velocity

Industrial safety valve is an important safety system component, which aims at preventing the system from over pressure and reducing further risk.The sealing performance is an important parameter for safety valve. However, instant impact is generated between the valve seat and the valve disc when the safety valve is shut down quickly. This condition seriously damages the sealing performance of components, which significantly decreases the life span of the safety valve.

In this study, the dynamic behaviors on a valve seat subjected to the impact of the valve disc at low velocities were investigated by using numerical analysis with the proposed method. The impact of the valve disc against seat was evaluated through an axis-symmetric model. Isotropic linear elastic response is assumed for valve disc and seat. The mechanical properties of valve disc and seat are listed in Table 3 . Accordingly, the simplified finite element model of the valve established, as illustrated in Fig. 8 . The mesh density of disc and seat is also displayed in Fig. 9 .

Table 3
Material parameters for valve disc and seat in the simulations

	E (GPa)	υ	ρ (kg/m3)
valve disc	193	0.3	7900
valve seat	214	0.3	7900

Figure 8
The simplified axis-symmetric finite element model of safety valve.

Figure 9
Finite element mesh for valve disc and seat.

The impact velocity on the elastic impact behavior of valve disc and seat were investigated. The simulation results could be served as a guide for impact velocity design of safety valve. Accordingly, Fig. 10 and Fig. 11 show the equivalent Mises stress on contact surface of valve disc and seat after the onset of the impact for different velocities at different velocities. Fig. 12 shows the contact pressure on the surface of valve disc and seat under different impact velocities. The contact stress at the impact velocity 3 m/s increased significantly corresponding to the results at the impact velocity 0.48 m/s. Fig. 11 indicates that the highest peaks of the contact stresses of valve disc at the impact velocities 3 m/s and 0.48 m/s were 2258 MPa and 360 MPa, respectively. The maximum equivalent stress occurs at the outer edge of valve disc.

Figure 10
The equivalent Mises stress of contact surface of valve disc and seat after the onset of the impact at velocity 0.48 m/s. (a) valve disc (b) valve seat

Figure 11
The equivalent Mises stress of contact surface of valve disc and seat after the onset of the impact at velocity 3.0 m/s. (a) valve disc (b) valve seat

Figure 12
The contact pressure on the surface of valve disc and seat under different impact velocities. (a) 0.48 m/s (b) 3.0 m/s

4 CONCLUSIONS

(1) This paper develops an improved dynamic contact model for mass-spring and finite element systems. In the proposed model, the contact states (stick and separation) are transformed into a simple unified formula with parametric variables. Consequently the complex trial-and-error and iterative processes in conventional dynamic contact numerical methods are avoided.

(2) The improved PQP-based dynamic contact model avoids the numerical troubles of spurious high-frequency oscillations. Artificial damping is not required to get rid of undesired oscillations or instability problems in the proposed method.

(3) The proposed model is convenient to apply to the mass-spring model and finite element model. Several engineering applications, such as analysis of wheel-rail contact problem and valve disc–seat impact problem. The numerical examples indicate that the proposed method improves stability and has good convergence behavior for dynamic contact problems.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant No. 51401045 and the Fundamental Research Funds for the Central Universities under Grant No. DUT16RC(4)36.

• http://dx.doi.org/10.1590/1679-78254420

References

• Auersch L., (2005). The excitation of ground vibration by rail traffic: theory of vehicle-track-soil interaction and measurements on high-speed lines. Journal of Sound and Vibration 284 (1-2): 103-132.
• Baeza L., Ouyang H., (2008). Dynamics of a truss structure and its moving-oscillator exciter with separation and impact-reattachment. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 464: 2517-2533.
• Balestrino A., Bruno O., Landi A., Sani L., (2010). Innovative solutions for overhead Catenary-Pantograph system: wire actuated control and observed contact force, International Journal of Vehicle Mechanics and Mobility, 33(2): 69-89.
• Chen J., Jiang D., Wang N., An S., (2015). Dynamic response of rub caused by a shedding annular component happening in a steam turbine, Shock and Vibration, 607674.
• Chen W., Yeh J., (1990). Three-dimensional finite element analysis of static and dynamic contact problems with friction, Computers & Structures, 35 (5): 541-552.
• Dowling M.J., Hall J. F., (1989). Nonlinear seismic analysis of arch dams, ASCE Journal of Engineering Mechanics, 115 (4): 768-789.
• Ferris M. C., Pang J. S., (1997). Engineering and economic applications of complementarity problems, SIAM Review, 39 (4): 669-713.
• Fischer A., (1992). A special newton-type optimization method, Optimization, 24 (3-4) 269-284.
• Hu N., (1997). A solution method for dynamic contact problems, Computers & Structures, 63 (6): 1053-1063.
• Jiang L., Rogers R. J., (1988). Combined Lagrangian Multiplier and Penalty Function Finite Element Technique for Elastic Impact Analysis, Computers & Structures, 30 (6): 1219-1229.
• Kanto Y., Akawa G. Y., (1990). A dynamic contact buckling analysis by the penalty finite element method, International Journal for Numerical Methods in Engineering, 29 (4): 755-774.
• Kim J. O., Kwak B. M., (1996). Dynamic analysis of two-dimensional frictional contact by linear complementarity problem formulation, International Journal of Solids Structures, 33 (30): 4605-4624.
• Laursen T.A., Chawla V., (1997). Design of energy conserving algorithms for frictionless dynamic contact problems, International Journal for Numerical Methods in Engineering, 40: 863-886.
• Liu X., Wang G. G., (2007). Non-linear dimensional variation analysis for sheet metal assemblies by contact modeling, Finite Elements in Analysis and Design, 44 (1-2): 34-44.
• Nguyena D. V., Kim K. D., Warnitchai P., (2009). Simulation procedure for vehicle–substructure dynamic interactions and wheel movements using linearized wheel–rail interfaces, Finite Elements in Analysis and Design, 45 (5): 341-356.
• Ouyang H., (2011). Moving-load dynamic problems: A tutorial (with a brief overview). Mechanical Systems and Signal Processing, 25 (6): 2039-2060.
• Padmanabhan V., Laursen T. A., (2001). A framework for development of surface smoothing procedures in large deformation frictional contact analysis, Finite Elements in Analysis and Design, 37 (3): 173-198.
• Peric D., Owen D.R.J., (1992). Computational model for 3-D contact problems with friction based on the penalty method, International Journal for Numerical Methods in Engineering, 35 (6):1289-1309.
• Rabb R., Fahrenthold E., (2010). Impact dynamics simulation for multilayer fabrics, International Journal for Numerical Methods in Engineering, 83(5): 537-557.
• Ruan H. H., Yu T. X., (2005). Collision between mass-spring systems, International Journal of Impact Engineering, 31 (3): 267-288.
• Simo J. C., Wriggers P., Taylor R. L., (1985). A perturbed Lagrangian formulation for the finite element solution of contact problems, Computer Methods in Applied Mechanics and Engineering, 50: 163-180.
• Solberg J. M., Papadopoulos P., (1998). A finite element method for contact/impact, Finite Elements in Analysis and Design, 30 (4): 297-311.
• Stăncioiu D., Ouyang H., Mottershead J. E., (2008). Vibration of a beam excited by a moving oscillator considering separation and reattachment. Journal of Sound and Vibration, 310 (4): 1128-1140.
• Taylor R. L., Carpenter N. J., Katona M. G., (1991). Lagrange constraints for transient finite element surface contact, International Journal for Numerical Methods in Engineering, 32 (1): 103-128.
• Wang F.J., Wang L. P., Cheng J. G., Yao Z. H., (2007). Contact force algorithm in explicit transient analysis using finite-element method, Finite Elements in Analysis and Design, 43 (6-7):580-587.
• Zhang H., He S., Li X., (2005). Two aggregate-function-based algorithms for analysis of 3D frictional contact by linear complementarity problem formulation, Computer Methods Applied Mechanics and Engineering, 194 (50-52): 5139-5158.
• Zhang H., Zhang L., Gao Q., (2012). Numerical method for dynamic analysis of two-dimensional bimodular structure, AIAA Journal, 50 (9): 1933-1945.

Publication Dates

• Publication in this collection
  2018

History

• Received
  22 Aug 2017
• Reviewed
  12 Dec 2017
• Accepted
  18 Dec 2017