\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2013 (2013), No. 85, pp. 1--23.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2013/85\hfil Dynamic contact of viscoelastic bodies] {Dynamic contact of viscoelastic bodies \\ with two obstacles: mathematical and \\ numerical approaches} \author[J. Ahn, J. Calhoun \hfil EJDE-2013/85\hfilneg] {Jeongho Ahn, Jon Calhoun} % in alphabetical order \address{Jeongho Ahn \newline Department of Mathematics and Statistics\\ Arkansas State University, PO Box 600\\ State University, AR 72846, USA} \email{jahn@astate.edu} \address{Jon Calhoun \newline Department of Computer Science \\ University of Illinois, Urbana-Champaign \\ 201 North Goodwin Avenue, Urbana, IL 61801-2302, USA} \email{jccalho2@illinois.edu} \thanks{Submitted September 5, 2012. Published April 5, 2013.} \subjclass[2000]{74M20, 74K10, 35L85} \keywords{Dynamic contact; variational inequality; signorini contact conditions; \hfill\break\indent complementarity conditions; strong pointedness; numerical scheme} \begin{abstract} The motion of viscoelastic (Kelvin-Voigt model) bodies between an upper and a lower obstacle is studied both mathematically and numerically. The two obstacles are assumed to be stationary perfect rigid, therefore, Signorini contact conditions are imposed at each obstacle, which can be interpreted as a couple of complementarity conditions (CCs). The convergence of numerical trajectories for general dimensional problems is shown based on the box constrained variational inequality (VI) which is equivalent to the two CCs. A one-dimensional example is provided. Unlike higher dimensional cases, different perspectives are used to prove the results of its existence. Numerical results are also presented and discussed, showing a typical behavior of the system \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \newtheorem{algorithm}[theorem]{Algorithm} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} This article considers a dynamic model for frictionless contact of an elastic or a viscoelastic body with two flat rigid obstacles situated above and below it, while most articles (e.g., see \cite{as:fncv,ks:dcscsrdf,ps:vmcs,ps:trcfaler}) deal with the related problem with only one obstacle. Thus, we impose Signorini contact conditions on the two obstacles, since the body may bounce off each of them. We note that mathematical and numerical approaches for more generalized obstacles seem to be very challenging; for instance, a more general dynamic case is when a viscoelastic body moves around inside a room enclosed by a rigid obstacle. Therefore, we commence by considering dynamic frictionless contact problems with two flat obstacles in order to extend into more general types of obstacles in our future research. We remark that many questions on the class of \emph{dynamic} problems with \emph{purely} elastic bodies $\Omega\subset\mathbb{R}^{d}$ with $d\geq3$, for example, the existence of solutions and their regularity, remain still open. Adding a \emph{viscosity} term into the equation of motion allows us to avoid some of these mathematical difficulties. Petrov and Schatzman \cite{ps:vmcs} study a one-dimensional problem where a viscoelastic rod is considered in the half space and the original partial differential equations (PDEs) with unilateral boundary conditions are switched to a VI on the end of the rod, using both convolution and CCs. In the paper \cite{as:fncv}, Ahn and Stewart extend the viscoelastic rod problem into a higher dimensional case. They use the trace theorem to set up the VI on the boundary $\partial\Omega$ which is equivalent to the complementarity problem (CP). In addition to that, numerical schemes are developed, convergence of numerical trajectories are shown, and some numerical results are presented, too. One of their major concerns is to show energy balance because the viscosity in the system causes energy loss. A significant issue is addressed in the two papers on how the viscosity affects the magnitude of contact forces, i.e., a viscoelastic body produces higher singularity than a purely elastic body. The paper \cite{ps:vmcs} proves the issue theoretically, and the paper \cite{as:fncv} provides numerical evidences for it. In Kuttler and Shillor's paper~\cite{ks:dcscsrdf}, \emph{frictional} contact is considered which is described by a slip rate dependent coefficient. They also regularize the nonlocal stress by using an averaging operator. In the paper~\cite{ps:trcfaler}, Shi studies a contact model with a \emph{purely} elastic rod and derives explicit formulas for the rebounding time period and the dependence of the coefficient of restitution on the initial condition. Indeed, imposing the coefficient of restitution (see the paper~\cite{aw:dip}) may be a good idea to show the uniqueness of solutions. \begin{figure}[ht] \centering \setlength{\unitlength}{0.8cm} \begin{picture}(10,7.5)(2.3,0) \put(2.3,0){\includegraphics[width = 8cm, height = 6cm]{fig1.png}} \put(7.5,3.2){$ {\Omega}$} \put(5.5,3.3){$ {\Gamma_1}$} \put(9.3,3.3){$ {\Gamma_2}$} \put(8.2,2.2){$ {\Gamma_{c}^{B}}$} \put(8.2,4.3){$ {\Gamma_{c}^{T}}$} \put(6.5,4.5){$ {N}$} \put(7.2,1.7){$ {N}$} \end{picture} \caption{Dynamic contact model with two obstacles} \label{fig:dynamic_contact} \end{figure} Among the recent works on the contact of beams, the Gao beam \cite{gao:netaxva}, which is highly nonlinear, has received considerable attention. In the recent papers \cite{ap:dcgt,aks:dcgb}, numerical algorithms are proposed and numerical results (simulations) are presented. In particular, the convergence theory for the Gao beam with dynamic contact has been initially studied in the paper \cite{ap:dcgt}. Even though the existence of solutions for the purely elastic case ($d\geq3$) is still an open question, many numerical schemes which are based on the Newmark schemes \cite{n:mcsd} have been proposed recently (see the paper~\cite{rm:fsdsecp} and the references therein). Applying Newmark schemes into the viscoelastic cases can provide a great idea to obtain higher stability of numerical solutions. This will be our future work when we develop numerical schemes on higher dimensional problems. The remaining sections of this paper is structured as follows. Section~\ref{sec:The-general-dimensional} which deals with higher dimensional problems ($d\geq2$) consists of two Subsections. In Subsection~\ref{sec:Mathematical-Preliminaries}, some mathematical background and notations are introduced. In Subsection~\ref{sub:The-existence-results}, the dynamic contact problem is formulated and its existence results are shown, using a box constrained VI on the boundary $\partial\Omega$. In Section~\ref{sec:An-one-dimensional}, a one dimensional problem is considered with different methods, since it satisfies the strong pointedness unlike the general dimensional problems. The convergence of numerical trajectories is shown through the first two Subsections~\ref{sub:Continuous-formulations} and \ref{sub:Numerical-formulations}. Subsection~\ref{sub:Fully-discrete-numerical} proposes the fully numerical schemes, and numerical results are presented and discussed in Subsection~\ref{sub:Numerical-Results-and}. This work concludes with some remarks in the last Section~\ref{sec:Conclusion}. \section{General dimensional problems} \label{sec:The-general-dimensional} This dynamic contact problem is considered with the Lipschitz domain $\Omega\subset\mathbb{R}^{d}$ (with $d=2, 3$ in applications) over the bounded time interval $[0, T]$. The deformation field and velocity of viscoelastic bodies are denoted by the functions $\mathbf{u}=\mathbf{u}(t, \mathbf{x})$ and $\mathbf{v}=\mathbf{v}(t, \mathbf{x})$, where $\mathbf{x}=(x_1,x_2,\dots,x_{d})\in\Omega$ is a material particle and time $t\in[0, T]$, respectively. The bodies move between two fixed rigid foundations which are expressed by functions $x_{d}=\varphi_{B}(x_1,x_2,\dots,x_{d-1})$ and $x_{d}=\varphi_T(x_1,x_2,\dots,x_{d-1})$ and $\varphi_{B}<\varphi_T$ for all points $(x_1,x_2,\dots,x_{d-1})$ on the boundary $\partial\Omega\subset\mathbb{R}^{d-1}$. By the proper transformation, we shall assume that $\varphi_{B}<0<\varphi_T$, which may assist to avoid geometrical complexities. The stress and the strain are denoted by the tensors $\boldsymbol{\sigma}=(\sigma_{ij})$ and $\boldsymbol{\varepsilon}=(\varepsilon_{ij})$ for $1\leq i, j\leq d$. In this paper, all vectors or tensors are written by boldface characters or component forms. Since our viscoelastic bodies are of type of the Kelvin-Voigt, we introduce its constitutive relation \begin{equation} \sigma[\mathbf{u},\mathbf{v}]=A \varepsilon(\mathbf{u}) +B\varepsilon(\mathbf{v}),\label{eq:stress} \end{equation} where $A$ is is called a linear elasticity operator and $B$ is called a linear viscosity operator and the linearized strain tensor is expressed by \[ \varepsilon(\mathbf{u})=\frac{\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}}{2}. \] Those tensors are the symmetric; i.e., $(\sigma_{ij})=(\sigma_{ji})$ and $(\varepsilon_{ij})=(\varepsilon_{ji})$. In the tensor form \eqref{eq:stress}, $A$ and $B$ can be reformed by the fourth order tensor in the component form; \[ \sigma_{ij}=a_{ijkl}\varepsilon_{ki}(\mathbf{u})+b_{ijkl}\varepsilon_{ki}(\mathbf{v})=a_{ijkl}u_{k,l}+b_{ijkl}v_{k,l}, \] where $A=(a_{ijkl})$ and $B=(b_{ijkl})$ with $1\leq k$, $l\leq d$ are an elasticity tensor and a viscosity tensor, respectively, and $u_{k,j}=\partial u_{k}/\partial x_{j}$ and $v_{k,j}=\partial v_{k}/\partial x_{j}$ are partial derivatives. Since there are the upper and lower rigid obstacles, we assume that the boundary of the body consists of four disjoint subsets; $\partial\Omega=\Gamma_1\cup\Gamma_{c}^{T}\cup\Gamma_2\cup\Gamma_{c}^{B}$ and $\Gamma_1\cap\Gamma_{c}^{T}\cap\Gamma_2\cap\Gamma_{c}^{B}=\emptyset$ with $\operatorname{meas}(\Gamma_1),\operatorname{meas}(\Gamma_2)\geq0$ such that contact forces do not take a place at $\Gamma_1,\Gamma_2\subset\partial\Omega$. Thus, on $\Gamma_{c}^{T}$ and $\Gamma_{c}^{B}$ the body may come in contact with the upper obstacle $x_{d}=\varphi_T$ and the lower obstacle $x_{d}=\varphi_{B}$. Those stationary obstacles are perfectly rigid and thus do not allow penetration of the viscoelastic body. The boundary $\partial\Omega$ of the body is assumed to have the unit outward normal vector $\mathbf{n}(\overline{\mathbf{x}})=(n_1(\overline{\mathbf{x}}), n_2(\overline{\mathbf{x}}),\dots,n_{d}(\overline{\mathbf{x}}))$ for almost all $\overline{\mathbf{x}}\in\partial\Omega$. This physical situation is illustrated in Figure~\ref{fig:dynamic_contact}. Thus, we are led to formulate the following PDEs: \begin{gather} \dot{\mathbf{v}} = \nabla\cdot\sigma[\mathbf{u}, \mathbf{v}] +\mathbf{f}\quad\text{in }(0, T]\times\Omega,\label{eq:eq_motion1}\\ \sigma[\mathbf{u}, \mathbf{v}] = A \varepsilon(\mathbf{u})+B \varepsilon(\mathbf{v}),\label{eq:stress_visco1}\\ \sigma[\mathbf{u}, \mathbf{v}]\cdot\mathbf{n} = N(t)\,\mathbf{n}\quad\text{in }(0, T]\times(\Gamma_{c}^{T}\cup\Gamma_{c}^{B}),\label{eq:contact_top_botttom}\\ \sigma[\mathbf{u}, \mathbf{v}]\cdot\mathbf{n} = \mathbf{0}\quad\text{in }(0, T]\times\partial\Omega \backslash(\Gamma_{c}^{T}\cup\Gamma_{c}^{B}),\label{eq:contact_zero}\\ 0\geq N(t) \perp \mathbf{u}\cdot\mathbf{n}-\varphi_T\leq0\quad\text{in }(0, T]\times\Gamma_{c}^{B},\label{eq:cc_top}\\ 0\leq N(t) \perp \mathbf{u}\cdot\mathbf{n}-\varphi_{B}\geq0\quad\text{in }(0, T]\times\Gamma_{c}^{T},\label{eq:cc_bottom}\\ \mathbf{u}^{0}=\mathbf{u}(0, \mathbf{x}), \quad \mathbf{v}^{0}=\mathbf{v}(0, \mathbf{x})\quad\text{in }\Omega,\label{eq:initial_data} \end{gather} where $\mathbf{u}^{0}$ and $\mathbf{v}^{0}$ are initial data and $\mathbf{f}$ is a body force and the normal contact forces $N$ are occurring at the top and the bottom obstacles. Defining the normal contact forces $\sigma_{ij}n_{j}=Nn_{i}$ in \eqref{eq:contact_top_botttom} enables them to be frictionless, since the normal stress $\sigma_{n}=\sigma_{ij}n_{i}n_{j}$ and the tangential stress $(\sigma_T)_{i}=\sigma_{ij}n_{j}-\sigma_{n}n_{i}=N\, n_{i} -\sigma_{ij}n_{i}n_{j}n_{i}=0$. A couple of CCs \eqref{eq:cc_top}--\eqref{eq:cc_bottom} which is equivalent to Signorini's contact conditions will be explained in the next Subsection~\ref{sec:Mathematical-Preliminaries}. The two natural boundary conditions \eqref{eq:contact_top_botttom}--\eqref{eq:contact_zero} do not guarantee the coercivity of operators for the static or even quasistatic problems. Nevertheless, we are still able to verify the existence of solutions for this type of dynamic frictionless contact problems. \subsection{Mathematical preliminaries} \label{sec:Mathematical-Preliminaries} In this Subsection, mathematical background and notation are introduced to present the existence results. The solution spaces that we mostly deal with are based on the Gelfand triples; $V\subset H=H'=V'$ (see the book~\cite{Wloka_j}), where all spaces are separable Hilbert spaces and all inclusions are densely compact. Here $(')$ denotes the dual space. For Banach space $X$, the duality pairing between $X'$ and $X$ is denoted by $\langle \cdot, \cdot\rangle _{X'\times X}$. When a duality pairing is defined on a known space, the simpler notation $\langle \cdot, \cdot\rangle $ may be used. In our problem, the Hilbert space $V$ will be Sobolev spaces typically. Now we define the linear operators $\mathcal{A}$ and $\mathcal{B}$ from $V$ to $V'$ with appropriate boundary conditions by \begin{gather*} \langle \mathcal{A}\mathbf{u}, \mathbf{w}\rangle := \int_{\Omega}A\varepsilon(\mathbf{u}):\varepsilon(\mathbf{u})dx,\\ \langle \mathcal{B}\mathbf{u}, \mathbf{w}\rangle := \int_{\Omega}B\varepsilon(\mathbf{u}):\varepsilon(\mathbf{u})dx, \end{gather*} where the notation $(:)$ means the product of tensors. In general, the CCs $0\leq a\,\bot\, b\geq0$ implies that $a, b\geq0$ and $a\cdot b=0$, and similarly, $0\ge a\,\bot\, b\leq0$ means that $a, b\leq0$ and $a\cdot b=0$. Due to the nonpenetration of the rigid obstacles, the CCs \eqref{eq:cc_top}--\eqref{eq:cc_bottom} can be understood in the following way. If viscoelastic bodies are not in contact with either of the obstacles; i.e., $N(t)=0$, then $\mathbf{u}(t, \overline{\mathbf{x}})\cdot\mathbf{n}(\overline{\mathbf{x}}) -\varphi_{B}>0$ for $\overline{\mathbf{x}}\in\Gamma_{c}^{B}$ and $\mathbf{u}(t, \overline{\mathbf{x}}) \cdot\mathbf{n}(\overline{\mathbf{x}})-\varphi_T<0$ for $\overline{\mathbf{x}}\in\Gamma_{c}^{T}$. However, if the body is in contact with the the bottom obstacle, then $N(t)\geq0$ and $\mathbf{u}(t, \overline{\mathbf{x}}) \cdot\mathbf{n}(\overline{\mathbf{x}})-\varphi_{B}=0$ for all $\overline{\mathbf{x}}\in\Gamma_{c}^{B}$. If the body is in contact with the top obstacle, then $N(t)\leq0$ and $\mathbf{u}(t, \overline{\mathbf{x}})\cdot\mathbf{n} (\overline{\mathbf{x}})-\varphi_T=0$ for all $\overline{\mathbf{x}}\in\Gamma_{c}^{T}$. Interpreting Signorini's contact conditions as CCs is an easier way to develop finite dimensional approaches and numerical algorithms, while using the indicator function and its subdifferential for Signorini's contact conditions can be a more theoretical way to prove the existence of solutions in the infinite dimension (see the paper~\cite{km:dcncwdfc} and references therein). In this dynamic contact problem, applying trace theorems (see the books~\cite[pp.83-89]{ko:cpevifem}) plays an essential role in deriving a box constrained (VI) on the boundary $\partial\Omega$ which is equivalent to the PDEs and CCs \eqref{eq:eq_motion1}--\eqref{eq:initial_data}. Since our dynamic contact problem is frictionless, we are able to apply the trace operator $\gamma$ from $\mathbf{H}^{1}(\Omega)=(H^{1}(\Omega))^{d}$onto $\mathbf{H}^{1/2}(\partial\Omega)=(H^{1/2}(\partial\Omega))^{d}$ such that \[ \gamma_{n}(\mathbf{w})=\gamma(\mathbf{w})\cdot\mathbf{n} \quad\text{for }\mathbf{w}\in\mathbf{H}^{1}(\Omega), \] where $\gamma_{n}$ from $\mathbf{H}^{1}(\Omega)$ onto $H^{1/2}(\partial\Omega)$ is the normal trace operator. Let the Sobolev spaces $\mathbf{H}^{1}(\Omega)$ and $H^{1/2}(\partial\Omega)$ be $V$ and $W$, respectively. Now the formal adjoint operator $\gamma_{n}^{*}:W'\to V'$ can be defined to be \begin{equation} \langle \varpi, \gamma_{n}\mathbf{w}\rangle _{W'\times W} =\langle \gamma_{n}^{*}\varpi, \mathbf{w}\rangle _{V'\times V}\quad \text{for all }\varpi\in W'\text{ and }\mathbf{w}\in V.\label{eq:self_adjoint} \end{equation} Since $\gamma_{n}$ is subjective, there is an extension operator $\mathcal{E}:W\to V$ so that $I_{W}=\gamma_{n}\mathcal{E}:W\to W$ is the identity operator. \subsection{Existence results for the general dimensional problems} \label{sub:The-existence-results} Let $K=[a, b]$ be a closed interval with $a0$ is the coefficient of elasticity. Since there are not simultaneously occurring contact forces at the two obstacles, we could impose the orthogonal condition $\int_{[0,T]}N_{B}(t) N_T(t)dt=0$ for all $t\in[0, T]$. However, it seems not necessary for improving the existence results. In addition, another possible contact model will be the compression of the rod applied by contact forces at each obstacle, without taking into consideration moving of the rod. For now, $H, V$ will become the spaces $L^{2}(0, l), H^{1}(0, l)$, respectively. We also define the self adjoint operator $\mathcal{A}:V\to V'$ by \[ \langle \mathcal{A}u, w\rangle :=\int_0^{l}u_{x}w_{x}dx. \] While we apply the trace operator for higher dimensional problems, we use two restriction operators $\beta_{B}(u)=u(0)$ and $\beta_T(u)=u(l)$ which are bounded operators $\beta_{B}, \beta_T:V\to W$. Moreover, the bounded operators provide the adjoint operators $\beta_{B}^{*}, \beta_T^{*}:W'\to V'$, where $\beta_{B}^{*}(N_{B})=N_{B}(t)\,\delta(x)$ and $\beta_T^{*}(N_T)=N_T(t)\,\delta(x-l)$. Here $\delta$ is the Dirac delta function which is the identity in the set of all distributions. Using integration by parts with the two natural boundary conditions \eqref{eq:BC1}--\eqref{eq:BC2} and applying the convergence of distributions, we can get to the following. \begin{lemma} \label{lem:distributions} If $N_T(t), N_{B}(t)\in W'$ with all $t\in[0, T]$ are given, then \eqref{eq:PDE}--\eqref{eq:BC2} are equivalent to the variational formulation; for all $t\in[0, T]$ \begin{align*} u(t)\in V:\langle v_{t}, w\rangle _{V'\times V} & = -\langle c\,\mathcal{A}u+\alpha\,\mathcal{A}v, w\rangle _{V'\times V}+(f, w)_{H}\\ &\quad +\langle N_T(t)\delta(x-l)-N_{B}(t)\delta(x), w\rangle _{V'\times V}\quad\text{for all }w\in V. \end{align*} \end{lemma} Referring to \cite[Lemma~1.1]{as:vtbdfi} will be helpful for understanding the proof. The variational formulation in Lemma~\ref{lem:distributions} will be used to set up the semi discrete numerical formulations in the time interval $[0, T]$. A set $K$ is called a cone if $\lambda x\in K$ when $x\in K$ and $\lambda\geq0$. Then, its dual cone denoted by $K'$ is defined as \[ K'=\{ y\in W'\mid\langle y, x\rangle _{W'\times W}\geq0 \text{ for all }x\in K\} . \] To prove the existence of solutions, the CCs will be interpreted with closed convex cones rather than with the original CCs \eqref{eq:CC1} and \eqref{eq:CC2}. Let $K_1$ and $K_2$ be closed convex cones such that $K_1\cap K_2=\{0\}$. Finally, the equation of motion with two contact conditions is reformulated in a more abstract setting: \begin{gather} v_{t}=c\,\mathcal{A}u+\alpha\,\mathcal{A}v+\beta_T^{*}(N_T) -\beta_{B}^{*}(N_{B})+f(t, x)\quad\text{in }(0, T]\times(0, l),\label{eq:abs1}\\ W'\supset K_1'\ni N_{B}(t)\,\perp\,\beta_{B}u-\varphi_{B}\in K_1 \subset W\quad\text{on }(0, T],\label{eq:abs_CC1}\\ W'\supset K_2'\ni N_T(t)\,\perp\,\beta_Tu+l-\varphi_T\in K_2 \subset W\quad\text{on }(0, T],\label{eq:abs_CC2}\\ u^{0}(x)=u(0, x),\quad v^{0}(x)=v(0, x)\quad\text{in }(0, l). \label{eq:abs_initials} \end{gather} Unlike the higher dimensional problems, we are able to show the boundedness of contact forces $N_T, N_{B}$ on a suitable space. In order to do so, we require the following definition. \begin{definition}\label{def:strong_point} \rm A dual cone $K'$ is said to be strongly pointed if there exist $\kappa\in K$ and $\eta>0$ such that for any $\xi\in K'$, \[ \eta\| \xi\| _{W'}\leq\langle \xi, \kappa\rangle _{W'\times W}. \] \end{definition} This definition can be founded in the papers~\cite{as:esciv,s:fmdiid}. Indeed, definition~\ref{def:strong_point} is not a perfect sufficient condition for purely elastic problems, since they require a gap in the scale of interpolation spaces. See the paper~\cite{as:esciv} for the detailed illustration. However, in this viscoelastic problem, being strongly pointed without considering a gap between interpolations spaces is sufficient to show the existence results. For now, we set $W=\mathbb{R}$ and thus $W'=\mathbb{R}$. Then we can have $K_1=\mathbb{R}_{+}\subset W$ and $K_2=\mathbb{R}_{-}\subset W$ and thus their dual cones become $K_1'=\mathbb{R}_{+}\subset W'$ and $K_2'=\mathbb{R}_{-}\subset W'$. Thus it is easy to see that the two dual cones $K_1'$ and $K_2'$ satisfy the strong pointedness. When we show the compactness of solutions, the H\"older space $C^{\theta}(0,T;V)$ with the norm is needed: \[ \| u\| _{C^{\theta}(0,T;V)}=\| u\| _{C(0,T;V)} +\sup_{s\neq t}\frac{\| u(t)-u(s)\| _{V}}{|t-s|^{\theta}}\quad\text{for } 0\leq s0$. \item The linear operator $\mathcal{A}:V\to V'$ is self adjoint. \item The space of solutions over the spacial domain is based on Gelfand triples; $V\subset H=H'\subset V'$. \item There exists a linear operator $\beta_{B}, \beta_T:V\to W$ such that its adjoint operator is $\beta_{B}^{*}, \beta_T^{*}:W'\to V'$ and $\beta_{B}, \beta_T$ are subjective and $\beta_T^{-1}, \beta_{B}^{-1}$ are bounded right inverse, i.e., $\beta_T^{-1}\beta_T=\beta_{B}^{-1}\beta_{B}=I_{W}$. \item The convex cones $K_1, K_2\subset W$ are closed and its dual cones $K_1', K_2'$ satisfy the strong pointedness. \end{itemize} Here is an important remark; for the convergence of numerical trajectories, two CCs \eqref{eq:abs_CC1} and \eqref{eq:abs_CC2} have to be understood in the sense of measures; i.e., \[ \int_{B}\langle N_{B}(t),\beta_{B}u(t)-\varphi_{B}\rangle _{W'\times W} =\int_{B}\langle N_T(t),\beta_Tu(t)-\varphi_T\rangle _{W'\times W}=0, \] where any Borel set $B\subseteq[0, T]$. Under those assumptions, we can arrive at the following main result. \begin{theorem} \label{thm:one_dim_result} There are solutions $(u, v)$ with $u\in C(0, T;C[0, l])\cap C^{1/2}(0,T;V)$ and $v\in L^{2}(0, T;V)\cap L^{\infty}(0,T;H)$ satisfy the equations \eqref{eq:abs1}--\eqref{eq:abs_initials}. \end{theorem} The proof of this theorem will be shown through several steps in the following Subsection~\ref{sub:Numerical-formulations}. \subsection{Numerical formulations with time discretization and its convergence} \label{sub:Numerical-formulations} Time discretization is carried out by partitioning the time interval with a small time step size $h_{t}>0$; \[ t_0=00$ \[ a(u, w)=\langle (I+h_{t}(c\, h_{t}+\alpha)\mathcal{A})u, w\rangle , \] where $I$ is the identity operator. Then it is easy to check that the linear operator $I+h(ch+\alpha)\mathcal{A}$ with any $h_{t}>0$ is bounded and elliptic on the space $V$. Substituting \eqref{eq:num_extra} in the left side of the equation \eqref{eq:num_motion}, we can easily obtain \begin{equation} a(u^{k+1}, w)=\langle \Phi^k, w\rangle \quad \text{for all }w\in V,\label{eq:bilinear} \end{equation} where $\Phi^k=u^k+h_{t}^{2}v^k+\alpha h_{t}\mathcal{A}u^k +h_{t}^{2}(N_T^k-N_{B}^k+f)\in V'$ can be found at the previous step. Therefore, by the Lax Milgram Lemma we can conclude that there is a unique next time stepping solution $u^{k+1}\in V$. Unlike the general dimensional problems, we have to impose an additional condition in order to prove boundedness of numerical trajectories in appropriate spaces. The condition \eqref{eq:condition_for_discrete} in Lemma~\ref{lem:energy_bd} is required only for numerical trajectories. The reason can be justified from the CCs \eqref{eq:CC1}. If $N_{B}=0$, then $N_{B} v(t, 0)=0$, but if $N_{B}>0$, then $u(t, 0)=\varphi_{B}$ and thus $N_{B} v(t, 0)=0$. From the numerical point of view, we can guess that the left side will be very close to zero, which implies that we can choose any suitable viscosity quantity $\alpha>0$. This argument will be support by numerical results in Subsection~\ref{sub:Numerical-Results-and}. For the rest of this article, a quantity $C>0$ does not depend on any parameters but it may be different in each occurrence. \begin{lemma} \label{lem:energy_bd} Suppose that the numerical solutions $(u^{k-1},v^{k-1},N_T^{k-1},N_{B}^{k-1})$ satisfy the discrete formulations \eqref{eq:num_motion}--\eqref{eq:num_cp_t}. If there exists an $\alpha>0$ such that \begin{equation} -\int_0^{t_{k}}\langle (N_{B})_{h_{t}}(\tau+h_{t}), \beta_{B}(v_{h_{t}}(\tau))\rangle d\tau\leq\alpha\int_0^{t_{k}}\langle \mathcal{A}v_{h_{t}}(\tau), v_{h_{t}}(\tau)\rangle d\tau,\label{eq:condition_for_discrete} \end{equation} then each time step solutions $u^k$ and $v^k$ with $k\geq1$ are uniformly bounded, independent of $h_{t}>0$. Furthermore, we can obtain the following estimates; \begin{equation} \begin{gathered} \max_{k\geq1}\| v^k\| _{H}\leq\sqrt{E^{0}(1+C T e^{T})}<\infty\quad \text{and} \\ \max_{k\geq1}\| u^k\| _{V}\leq\sqrt{(T\, E^{0}(2T+C)(1+C T e^{T}) +\| u^{0}\| _{H}^{2})}<\infty. \end{gathered} \label{eq:energy_est} \end{equation} \end{lemma} \begin{proof} Multiplying the left side of \eqref{eq:num_motion} by the left of \eqref{eq:num_extra} we have \begin{equation} \begin{aligned} \frac{1}{h_{t}}(v^{k+1}-v^k, v^{k+1}) & = \frac{1}{2h_{t}}(v^{k+1}-v^k, v^{k+1}+v^{k+1}+v^k-v^k)_{H} \\ & = \frac{1}{2h_{t}}(\| v^{k+1}\| _{H}^{2}-\| v^k\| _{H}^{2}+(v^{k+1}-v^k, v^{k+1}-v^k)_{H}). \end{aligned}\label{eq:left_num} \end{equation} Similarly, we use \eqref{eq:num_extra} to modify the right side of \eqref{eq:num_motion}; \begin{equation} \begin{aligned} & -\frac{1}{h_{t}}\langle \mathcal{A}u^{k+1}, u^{k+1}-u^k\rangle -\alpha\langle \mathcal{A}v^{k+1}, v^{k+1}\rangle +\frac{1}{h_{t}}\langle N_T^k, \beta_T(u^{k+1}-u^k)\rangle \\ & -\frac{1}{h_{t}}\langle N_{B}^k, \beta_{B}(u^{k+1}-u^k)\rangle -(f, v^{k+1}) \\ & =-\frac{1}{2h_{t}}\langle \mathcal{A}u^{k+1}, u^{k+1}\rangle +\frac{1}{2h_{t}}\langle \mathcal{A}u^k, u^k\rangle -\frac{1}{2h_{t}}\langle \mathcal{A}(u^{k+1}-u^k), u^{k+1}-u^k\rangle \\ & \quad -\alpha\langle \mathcal{A}v^{k+1}, v^{k+1}\rangle +\frac{1}{h_{t}}\langle N_T^k, \beta_T(u^{k+1})+l-l-\beta_T(u^k)+\varphi_T-\varphi_T\rangle \\ & \quad -\frac{1}{h_{t}}\langle N_{B}^k, \beta_{B}(u^{k+1})-\beta_{B}(u^k)+\varphi_{B}-\varphi_{B}\rangle -(f, v^{k+1}).\label{eq:right_num} \end{aligned} \end{equation} From \eqref{eq:left_num}--\eqref{eq:right_num}, it follows that \begin{align*} E^{k+1} & \leq E^k-\alpha h_{t}\langle \mathcal{A}v^{k+1}, v^{k+1}\rangle +\langle N_T^k, \beta_T(u^{k+1})+l-\varphi_T\rangle \\ & \quad -\langle N_T^k, \beta_T(u^k)+l-\varphi_T\rangle -\langle N_{B}^k, \beta_{B}(u^{k+1})-\varphi_{B}\rangle \\ & \quad +\langle N_{B}^k, \beta_{B}(u^k)-\varphi_{B}\rangle +h_{t}(f, v^{k+1}). \end{align*} It follows from the two numerical CCs \eqref{eq:num_cp_b}--\eqref{eq:num_cp_t} that \begin{equation} E^{k+1}\leq E^k-\alpha h_{t}\langle \mathcal{A}v^{k+1}, v^{k+1}\rangle +\langle N_{B}^k, \beta_{B}(u^k)-\varphi_{B}\rangle +h_{t}(f, v^{k+1}). \label{eq:energy_dec} \end{equation} By using the CCs \eqref{eq:num_cp_t} and the extra equation \eqref{eq:num_extra}, the second term on the right side of \eqref{eq:energy_dec} becomes \[ \langle N_{B}^k, \beta{}_{B}(u^k)-\varphi_{B}\rangle =\langle N_{B}^k, \beta_{B}(u^{k+1})-h_{t}\beta_{B}(v^{k+1}) -\varphi_{B}\rangle =-h_{t}\langle N_{B}^k, \beta_{B}(v^{k+1})\rangle . \] Let $N_{h_{t}}(t)=0$ for all $t\in[-h_{t}, 0]$. The telescoping series and the assumption enable us to get to the following; for any integers $k\geq1$ \begin{equation} \begin{aligned} E^k &\leq E^{0}-\int_0^{t_{k}}\langle (N_{B})_{h_{t}}(t+h_{t}), \beta_{B}(v_{h_{t}}(t))\rangle dt +\int_0^{t_{k}}(f(t), v_{h_{t}}(t))dt \\ &\quad -\alpha\int_0^{t_{k}}\langle \mathcal{A}v_{h_{t}}(t), v_{h_{t}}(t)\rangle dt\leq E^{0}+\int_0^{t_{k}} \| f(t)\| _{H}\| v_{h_{t}}(t)\| _{H}dt. \end{aligned}\label{eq:energy_bd} \end{equation} Now, we can use H\"older's inequality to see that \[ \| v^k\| _{H}^{2}\leq E^{0}+C\int_0^{t_{k}}\| v_{h_{t}}(t)\| _{H}^{2}dt. \] It follows from Grownall's inequality that \[ \| v^k\| _{H}^{2}\leq E^{0}(1+C\, T\, e^{T})\quad\text{for any }k\geq1. \] Since $\langle \mathcal{A}(\cdot), \cdot\rangle $ is not equivalent to the $V$ norm $\| \cdot\| _{V}$, we claim that $\| u^k\| _{H}$ is uniformly bounded for any $k\geq0$. Since $u^k(\cdot)=\int_0^{t_{k}}v_{h_{t}}(t,\cdot)\, dt+u^{0}(\cdot)$, we can see that for any $k\geq1$ \[ \| u^k\| _{H}^{2}\leq2(T\int_0^{t_{k}}\| v_{h_{t}}\| _{H}^{2}dt +\| u^{0}\| _{H}^{2}) \leq 2(T^{2}E^{0}(1+C\, T\, e^{T})+\| u^{0}\| _{H}^{2}). \] Therefore, we can obtain the two estimates \eqref{eq:energy_est}. \end{proof} As we observe Lemma~\ref{lem:energy_bd}, the interpolants $u_{h_{t}}$ is uniformly bounded in $C(0, T;V)$ and $v_{h_{t}}$ is uniformly bounded in $L^{\infty}(0, T;H)$. \begin{lemma} \label{lem:velocity_bd} Under the same assumption of Lemma~\ref{lem:energy_bd}, we have the following estimate, independent of $h_{t}>0$: \[ \int_0^{T}\| v_{h_{t}}\| _{V}^{2}dt\leq M<\infty. \] \end{lemma} \begin{proof} It follows from the two estimates \eqref{eq:right_num} and \eqref{eq:energy_bd} that for any $h_{t}>0$, \begin{align*} \infty & > E^{0}(1+C\, T(1+C\, T\, e^{T}))\\ & \geq \alpha h_{t}\sum_{k=0}^{m-1}\langle \mathcal{A}v^{k+1}, v^{k+1}\rangle =\alpha\int_0^{T}\| v_{h_{t}}\| _{V}^{2}dt, \end{align*} as required. \end{proof} Using energy boundedness and Lemma~\ref{lem:velocity_bd}, it is easy to prove that $u_{h_{t}}$ is uniformly bounded in the H\"older space $C^{1/2}(0,T;V)$. Thus, we can apply the Arzela-Ascoli Theorem~\cite[pp.114]{mr:ipde} and Sobolev imbedding Theorem~\cite[pp.215]{mr:ipde} to show that there is a subsequence, denoted by $u_{h_{t}}$ such that $u_{h_{t}}\to u$ in $C[0, T]\times C[0, l]$, as $h_{t}\downarrow0$. Next, we want to show that $(N_{B})_{h_{t}}, (N_T)_{h_{t}}$ are bounded in the measure senses. It follows from \eqref{eq:num_contact_bot}--\eqref{eq:num_contact_top} that \[ \int_0^{T}\| (N_{B})_{h_{t}}\| _{W'}dt=h_{t} \sum_{k=0}^{m-1}\| N_{B}^k\| _{W'},\quad\int_0^{T}\| (N_T)_{h_{t}}\| _{W'}dt=h_{t}\sum_{k=0}^{m-1}\| N_T^k\| _{W'}. \] Since our one dimensional problem satisfies strong pointedness, there are $\eta_1, \eta_2>0$ and $\kappa_1\in K_1, \kappa_2\in K_2$ such that \begin{align*} & \int_0^{T}\| (N_{B})_{h_{t}}\| _{W'}dt +\int_0^{T}\| (N_T)_{h_{t}}\| _{W'}dt\\ & \leq\eta_1\, h_{t}\sum_{k=0}^{m-1}\langle N_{B}^k,\kappa_1\rangle +\eta_2\, h_{t}\sum_{k=0}^{m-1}\langle N_T^k, \kappa_2\rangle . \end{align*} We can choose $w$ as $w(x):=-x+l/2$ and thus $\kappa_1=l/2\in K_1$ and $\kappa_2=-l/2\in K_2$. Then by using the discrete equations \eqref{eq:num_motion}--\eqref{eq:num_extra}, we can show easily that two contact forces are bounded as $W'$-measures, independent of $h_{t}>0$. Therefore, there are subsequences, denoted by $(N_{B})_{h_{t}}$ and $(N_T)_{h_{t}}$ such that $(N_{B})_{h_{t}}\rightharpoonup^{*}N_{B}$ and $(N_T)_{h_{t}}\rightharpoonup^{*}N_T$ as $h_{t}\downarrow0$ in the sense of measures. Finally, we need to prove that the solutions which are convergent by the subsequences satisfy the CCs \eqref{eq:CC1} and \eqref{eq:CC2}. Since $(N_T)_{h_{t}}\leq0$ and $(N_{B})_{h_{t}}\geq0$, it turns out that $N_T\leq0$ and $N_{B}\geq0$. Similarly, since $u_{h_{t}}-\varphi_{B}\geq0$ and $u_{h_{t}}+l-\varphi_T\leq0$, $u-\varphi_{B}\geq0$ and $u+l-\varphi_T\leq0$. We also need to show that the limits of subsequences satisfy the CCs in a measure sense. It is easy to see from \eqref{eq:num_contact_bot}--\eqref{eq:num_contact_top} that \[ \int_0^{T}\int_0^{l}((N_{B})_{h_{t}}, (N_T)_{h_{t}})^{T}\cdot(u_{h_{t}} (t, x)-\varphi_{B}, u_{h_{t}}(t, x)+l-\varphi_T)dx\, dt=0. \] By the convergence of subsequences $(N_T)_{h_{t}}, (N_{B})_{h_{t}}$ and $u_{h_{t}}$ as $h_{t}\downarrow0$, we can obtain \begin{align*} 0 & = \int_0^{T}\int_0^{l}(N_{B})_{h_{t}}(u_{h_{t}}-\varphi_{B})\,dx\,dt +\int_0^{T}\int_0^{l}(N_T)_{h_{t}}(u_{h_{t}}+l-\varphi_T)\,dx\,dt\\ & \to\int_0^{T}\int_0^{l}N_{B}(u-\varphi_{B})\,dx\,dt +\int_0^{T}\int_0^{l}N_T(u-\varphi_T)\,dx\,dt. \end{align*} From Lemmas~\ref{lem:energy_bd} and \ref{lem:velocity_bd}, $v_{h_{t}}$ is bounded in $L^{\infty}(0,T;H)\cap L^{2}(0,T;V)$. Therefore, by applying the Alaoglu's theorem, we can say that there exists a subsequence, denoted by $v_{h_{t}}$ such that $v_{h_{t}}\rightharpoonup^{*}v$ in $L^{\infty}(0,T;H)\cap L^{2}(0,T;V)$ as $h_{t}\downarrow0$, since $L^{\infty}(0,T;H)\simeq L^{1}(0,T;H)'$ and $L^{2}(0,T;V)\simeq L^{2}(0,T;V)'$. \subsection{Fully discrete numerical schemes} \label{sub:Fully-discrete-numerical} Finite Element Methods (FEMs) are popular numerical schemes which are used to find approximations of solutions for elliptic PDEs. Our FEM is accomplished by the uniform partitioning: \[ x_0=00$ for non-negative integers $j\geq0$. Then, the finite dimensional space is chosen by \[ V_{h_{x}}=\{ w_{h_{x}}\in H^{1}(0, l)\mid w_{h_{x}}\in\mathcal{L}([x_{j}, x_{j+1}]),0\leq j\leq n-1\} , \] where $\mathcal{L}$ is a family of piecewise linear functions. Thus, the basis functions, associated with the each node $x_{j}$ with $1\leq j\leq n-1$ are set up as follows; \[ \Psi_{j}(x)= \begin{cases} (x-x_{j-1})/h_{x} &\text{on }[x_{j-1}, x_{j}],\\ (x_{j+1}-x)/h_{x} &\text{on }[x_{j}, x_{j+1}],\\ 0 &\text{on }[0, l]\backslash[x_{j-1}, x_{j+1}], \end{cases} \] and the first and last basis functions are \begin{gather*} \Psi_0(x)= \begin{cases} (x_1-x)/h_{x}&\text{on }[0, h_{x}],\\ 0 &\text{on }[0, l]\backslash[0, h_{x}], \end{cases}\\ \Psi_{n}(x)= \begin{cases} (x-x_{n-1})/h_{x} &\text{on }[l-h_{x}, l],\\ 0 &\text{on }[0, l]\backslash[l-h_{x}, l]. \end{cases} \end{gather*} Having applied time discretization into the time interval $[0, T]$ and partitioned the spacial domain $[0, l]$ into small sub-intervals, we now assume that all full approximations of deformation and velocity, denoted respectively by $u_{h_{t},h_{x}}^k$ and $v_{h_{t},h_{x}}^k$ are written at each time step $t_{k}$ as follows; \[ u_{h_{t},h_{x}}^k(t_{k}, x):=u_{h_{x}}^k(x) =\sum_{j=0}^{n}u_{j}^k\Psi_{j}(x),\quad v_{h_{t},h_{x}}^k(t_{k}, x):=v_{h_{x}}^k(x) =\sum_{j=0}^{n}v_{j}^k\Psi_{j}(x). \] Recalling the semi-discrete formulations \eqref{eq:num_motion}--\eqref{eq:num_cp_t}, we establish fully discrete formulations with two natural boundary conditions \eqref{eq:BC1} and \eqref{eq:BC2}: \begin{gather} \frac{v_{h_{x}}^{k+1}-v_{h_{x}}^k}{h_{t}} = c(u_{h_{x}}^{k+1})''+\alpha(v_{h_{x}}^{k+1})''+f,\label{eq:fully_motion-1}\\ \frac{u_{h_{x}}^{k+1}-u_{h_{x}}^k}{h_{t}} = v_{h_{x}}^{k+1},\label{eq:fully_extra-1}\\ 0\leq N_{B}^k \perp u_{h_{x}}^{k+1}(0)-\varphi_{B}\geq0,\label{eq:fully_cc_b-1}\\ 0\geq N_T^k \perp u_{h_{x}}^{k+1}(l)+l-\varphi_T\leq0,\label{eq:fully_cc_t-1} \end{gather} where $('')$ is the second derivative with respect to $x\in(0, l)$. We notice that the contact forces $N_{B}^k, N_T^k$ are approximated only over the time space $[0, T]$, since those are computed only at two end points $x=0, l$. Using the extra equation \eqref{eq:fully_extra-1} we can set up the following recursive formula; \begin{equation} \frac{1}{h_{t}^{2}}u_{h_{x}}^{k+1}-\frac{\alpha}{h_{t}}u_{h_{x}}^{k+1} -c(u_{h_{x}}^{k+1})''=\frac{1}{h_{t}^{2}}u_{h_{x}}^k -\frac{\alpha}{h_{t}}(u_{h_{x}}^k)''+\frac{1}{h_{t}}v_{h_{x}}^k+f. \label{eq:recursive-1} \end{equation} For our actual simulations, we shall assume that the body force $f$ is expressed by $f=\sum_{j=0}^{n}g_{j}\,\Psi_{j}(x)$. Multiplying both sides in \eqref{eq:recursive-1} by the basis function $\Psi_{i}(x)$ with $0\leq i\leq n$ and then using integration by parts with the boundary conditions, we obtain \begin{equation} \begin{aligned} &\frac{1}{h_{t}^{2}}\Big(\sum_{j=0}^{n}u_{j}^{k+1} \int_0^{l}\Psi_{j}(x)\Psi_{i}(x)\, dx -\sum_{j=0}^{n}u_{j}^k\int_0^{l}\Psi_{j}(x)\Psi_{i}(x)\, dx\Big) \\ &-\frac{1}{h_{t}}\sum_{j=0}^{n}v_{j}^k\int_0^{l}\Psi_{j}(x)\Psi_{i}(x)\, dx\\ &=c\sum_{j=0}^{n}u_{j}^{k+1}\int_0^{l}\Psi'_{j}(x)\Psi'_{i}(x)\, dx +\frac{\alpha}{h_{t}}\sum_{j=0}^{n}u_{j}^{k+1} \int_0^{l}\Psi'_{j}(x)\Psi'_{i}(x)\, dx +N_T^k\Psi_{i}(1)\\ &\quad -\frac{\alpha}{h_{t}}\sum_{j=0}^{n}u_{j}^k \int_0^{l}\Psi'_{j}(x)\Psi'_{i}(x)\, dx +\sum_{j=0}^{n}g\int_0^{l}\Psi_{j}(x)\Psi_{i}(x)\, dx -N_{B}^k\Psi_{i}(0). \end{aligned} \label{eq:Integrated_GA-1} \end{equation} Before switching the integrations in \eqref{eq:Integrated_GA-1} into a linear system, we introduce two matrices, mass $\mathbf{M}$ and stiffness $\mathbf{K}$, which are defined in \eqref{eq:matrices-1}, respectively; \begin{equation} \mathbf{M}=M_{ij}=\int_0^{l}\Psi_{i}(x)\Psi_{j}(x)\, dx,\quad \mathbf{K}=K_{ij}=\int_0^{l}\Psi'_{i}(x)\Psi'_{j}(x)\, dx.\label{eq:matrices-1} \end{equation} The next step deformation vector $\widetilde{u}^{k+1}\in\mathbb{R}^{n+1}$ can be computed by the following linear system at each time step; \begin{equation} \Big(\frac{1}{h_{t}^{2}}\mathbf{M}+(c+\frac{\alpha}{h_{t}})\mathbf{K}\Big) \widetilde{u}^{k+1} =\big[\big(\frac{1}{h_{t}^{2}}\mathbf{M}+\frac{\alpha}{h_{t}}\mathbf{K}\big) \widetilde{u}^k +\frac{1}{h_{t}}\mathbf{M}\widetilde{v}^k+\widetilde{N}^k +\mathbf{M}\widetilde{f}\big],\label{eq:linear system-1} \end{equation} where the previous fully discrete approximations are given by the following vector forms; \begin{gather*} \widetilde{u}^k=(u_0^k,u_1^k,\dots,u_{n}^k)^{T}, \quad\widetilde{v}^k=(v_0^k,v_1^k,\dots,v_{n}^k)^{T},\\ \widetilde{N}^k=(N_{B}^k,0,\dots,0,-N_T^k)^{T}, \quad\widetilde{f}=(g_0,g_2,\dots,g_{n})^{T}. \end{gather*} Here each component $g_{i}$ will be taken to be $g_{i}=-9.81$ for all $0\leq i\leq n$ in our actual simulations. The linear system \eqref{eq:linear system-1} is incomplete, because the next step solution $\widetilde{u}^{k+1}$ needs to satisfy the two numerical CCs \eqref{eq:fully_cc_b-1} and \eqref{eq:fully_cc_t-1}. Now, we consider the linear system \eqref{eq:linear system-1} as the simplified form $\mathbf{A}\,\widetilde{u}^{k+1}=\widetilde{b}^k$ with \begin{gather*} \mathbf{A} = (\frac{1}{h_{t}^{2}}\mathbf{M}+(c+\frac{\alpha}{h_{t}}) \mathbf{K})\in\mathbb{R}^{(n+1)\times(n+1)}\quad \text{and}\\ \widetilde{b}^k = (\frac{1}{h_{t}^{2}}\mathbf{M} +\frac{\alpha}{h_{t}}\mathbf{K})\widetilde{u}^k +\frac{1}{h_{t}}\mathbf{M}\widetilde{v}^k+\widetilde{N}^k +\mathbf{M}\widetilde{f}\in\mathbb{R}^{(n+1)}. \end{gather*} Next, we break apart the matrix $\mathbf{A}$ into the submatrices $\mathbf{A}_1,\mathbf{A}_{4}\in\mathbb{R}^{n\times n}$, column vectors $\widetilde{a}_2=a_{n-1n}\mathbf{e}_2$, $\widetilde{a}_{5}=a_{21}\mathbf{e}_{5}$, row vectors $\widetilde{a}_3=a_{nn-1}\mathbf{e}_3$, $\widetilde{a}_6=a_{12}\mathbf{e}_6$, and entries $a_{nn}$, $a_{00}$ as shown below; \[ \mathbf{A=}\begin{bmatrix} & & & |\\ & \mathbf{A}_1 & & | & \widetilde{a}_2\\ & & & |\\ - & - & - & | & -\\ & \widetilde{a}_3 & & | & a_{nn} \end{bmatrix}]\quad\text{or } \mathbf{A}=\begin{bmatrix} a_{00} & | & & \widetilde{a}_6\\ - & | & - & - & -\\ & |\\ \widetilde{a}_{5} & | & & \mathbf{A}_{4}\\ & | \end{bmatrix}, \] where $\mathbf{e}_2=(0,\dots,0,1)^{T}\in\mathbb{R}^{n}$, $\mathbf{e}_{5}=(1,\dots,0,0)^{T}\in\mathbb{R}^{n}$, $\mathbf{e}_3=\mathbf{e}_2^{T}$, and $\mathbf{e}_{5}=\mathbf{e}_6^{T}$. From those submatirces we split the linear system $\mathbf{A}\,\widetilde{u}^{k+1}=\widetilde{b}^k$ into a vector equation \eqref{eq:Vector Equation-2} and a scalar equation \eqref{eq:Scalar Equation-2}; \begin{gather} \mathbf{A}_1\widetilde{x}^{k+1}+u_{n}^{k+1}\widetilde{a}_2 =\widetilde{y}^k,\label{eq:Vector Equation-2}\\ \widetilde{a}_3\widetilde{x}^{k+1}+a_{nn}u_{n}^{k+1} =b_{n}^k,\label{eq:Scalar Equation-2} \\ 0\geq q\, u_{n}^{k+1}+z^k\perp u_{n}^{k+1}+l-\varphi_T\leq0,\label{eq:CC-linear system-2} \end{gather} where $\widetilde{x}^{k+1}=(u_0^k,u_1^k, \dots,u_{n-1}^k)^{T}\in\mathbb{R}^{n}$, $\widetilde{y}^k=(b_0^k,b_1^k,\dots,b_{n-1}^k)^{T}\in\mathbb{R}^{n}$, $q=(-\widetilde{a}_3\mathbf{A}_1^{-1}\widetilde{a}_2+a_{nn})$, and $z^k$ contains the quantities coming from the previous data. One can notice that $N_T^k$ is replaced by $q\, u_{n}^{k+1}+z^k$ in the CCs \eqref{eq:CC-linear system-2}. Thus we can compute $u_{n}^{k+1}$ from the CCs \eqref{eq:CC-linear system-2} and use the vector equation \eqref{eq:Vector Equation-2} to compute the rest of components in the next step solution $\widetilde{u}^{k+1}$ by finding $\widetilde{x}^{k+1}$. Similarly, using the submatrices $\mathbf{A}_{4}$, we can arrive at the following equations: \begin{gather} \mathbf{A}_{4}\widetilde{r}^{k+1}+u_0^{k+1}\widetilde{a}_{5} =\widetilde{s}^k,\label{eq:Vector Equation-1-1} \\ \widetilde{a}_6\widetilde{r}^{k+1}+a_{00}u_0^{k+1}=b_0^k, \label{eq:Scalar Equation-1-1} \\ 0\leq p\, u_0^{k+1}+d^k\perp u_0^{k+1}-\varphi_{B}\geq0,\label{eq:CC-linear system-1-1} \end{gather} where $\widetilde{r}^{k+1}=(u_1^k,u_2^k,\dots,u_{n}^k)^{T}\in\mathbb{R}^{n}$, $\widetilde{s}^k=(b_1^k,b_2^k,\dots,b_{n}^k)^{T}\in\mathbb{R}^{n}$, and $p=(-\widetilde{a}_6\mathbf{A}_{4}^{-1}\widetilde{a}_{5}+a_{00})$. Therefore, $\widetilde{u}^{k+1}$ can be computed through \eqref{eq:Vector Equation-1-1}--\eqref{eq:CC-linear system-1-1}. Finally, we present the numerical algorithm which summarizes our numerical schemes proposed above. Additionally, the initial contact forces are assumed to be zero, since the rod moves down initially without any contact. \subsection*{Algorithm} Suppose that the initial data $\widetilde{u}^{0}$, $\widetilde{v}^{0}$, and $\widetilde{N}^{0}=\mathbf{0}$ are given. \begin{tabbing} \textbf{for}\= $k=1:T/h_{t}$\\ \> \textbf{if} \= $N_T^{k-1}=0$ \% Assume that a rod drops down\\ \> \> \textbf{if} \= $u_0^k=\varphi_{B}$\\ \> \> \> $N_{B}^{k-1}\leftarrow\varphi_{B}\, p+d^k$ \% use \eqref{eq:CC-linear system-1-1}\\ \>\> \textbf{elseif} \= $u_0^k>\varphi_{B}$ \\ \> \> \> $N_{B}^{k-1}\leftarrow0$\\ \> \> \> $u_0^k\leftarrow-d^k/p$ \% use \eqref{eq:Vector Equation-1-1}--\eqref{eq:Scalar Equation-1-1}\\ \> \> \textbf{endif}\\ \> \> Compute $\widetilde{r}^k$ and then obtain $\widetilde{u}^k$ from \eqref{eq:Vector Equation-1-1}\\ \> \textbf{endif} \\[4pt] \> \textbf{if} \= $N_T^{k-1}<0$\\ \> \> \textbf{if} \= $u_{n}^k=\varphi_T-l$\\ \> \> \> $N_T^{k-1}\leftarrow(u_{n}^k-l)\, q+z^k$ \% use \eqref{eq:CC-linear system-2}\\ \> \> \textbf{elseif} \= $u_{n}^k<\varphi_T-l$\\ \> \> \> $N_T^{k-1}\leftarrow0$\\ \> \> \> $u_{n}^k\leftarrow-z^k/q+l$ \% use \eqref{eq:Vector Equation-2}--\eqref{eq:Scalar Equation-2}\\ \> \>\textbf{endif}\\ \> \> Compute $\widetilde{x}^k$ and then obtain $\widetilde{u}^k$ from \eqref{eq:Vector Equation-2}\\ \> \textbf{endif} \\[4pt] \> Compute $\widetilde{v}^k$ from \eqref{eq:fully_extra-1}\\ \> Compute the actual displacement using the change of a variable\\ \textbf{endfor} \end{tabbing} In the fully discrete case the energy function at each time step $t=t_{k}$ can be defined as \begin{equation} E^k:=E(t_{k})=\frac{1}{2}[(\widetilde{v}^k)^{T} \mathbf{M}\widetilde{v}^k+c\,(\widetilde{u}^k)^{T} \mathbf{K}\widetilde{u}^k]-(\widetilde{f})^{T}\mathbf{M} \widetilde{u}^k.\label{eq:Energy-1} \end{equation} The energy function will be evaluated at each time step $t_{k}$. We monitor the energy throughout the simulation because our numerical scheme does not permit the energy to be unbounded. The numerical evidence of bounded energy function is supported in Lemma~\ref{lem:fully_energy_bd-1}. \begin{lemma} \label{lem:fully_energy_bd-1} Suppose that numerical solutions satisfy the fully discrete formulation \eqref{eq:fully_motion-1}--\eqref{eq:fully_cc_t-1}. If there is an $\alpha>0$ such that \begin{equation} \alpha\sum_{\iota=1}^k\widetilde{v}^{\iota}\mathbf{K}\widetilde{v}^{\iota} +\sum_{\iota=1}^kN_{B}^{\iota-1}v_0^{\iota}\geq0\label{eq:energy_cond-1} \end{equation} for any integer $k\ge1$. Then $E^k\leq E^{0}$. \end{lemma} \begin{proof} The next solutions are decomposed as follows; the displacement can be $(u_{h_{x}}^{k+1})''=\frac{1}{2}[((u_{h_{x}}^{k+1})''-(u_{h_{x}}^k)'') +((u_{h_{x}}^{k+1})''+(u_{h_{x}}^k)'')]$ and the velocity can be $v_{h_{x}}^{k+1}=\frac{1}{2}[((v_{h_{x}}^{k+1})-(v_{h_{x}}^k)) +((v_{h_{x}}^{k+1})+(v_{h_{x}}^k))]$. We integrate over the length of the original rod, and by recalling the extra equation \eqref{eq:fully_extra-1}, the CCs \eqref{eq:fully_cc_b-1}--\eqref{eq:fully_cc_t-1} and the boundary conditions \eqref{eq:BC1}--\eqref{eq:BC2} allow us to cancel some terms; \begin{align*} &\frac{1}{2h_{t}}\sum_{i,j=0}^{n}(v_{j}^{k+1}-v_{j}^k) \int_0^{l}\Psi_{j}(x)\Psi_{i}(x)dx\,(v_{i}^{k+1}-v_{i}^k)\\ & +\sum_{i,j=0}^{n}v_{j}^{k+1}\int_0^{l}\Psi_{j}(x)\Psi_{i}(x)dx\, v_{i}^{k+1} -\sum_{i,j=0}^{n}v_{j}^k\int_0^{l}\Psi_{j}(x)\Psi_{i}(x)dx\, v_{i}^k\\ &=-\frac{c}{2h_{t}}\sum_{i,j=0}^{n}(u_{j}^{k+1}-u_{j}^k) \int_0^{l}\Psi'_{j}(x)\Psi'_{i}(x)dx\,(u_{i}^{k+1}-u_{i}^k)\\ &\quad -\frac{c}{2h_{t}}\sum_{i,j=0}^{n}u_{j}^{k+1} \int_0^{l}\Psi'_{j}(x)\Psi'_{i}(x)dx\, u_{i}^{k+1}+\frac{c}{2h_{t}} \sum_{i,j=0}^{n}u_{j}^k\int_0^{l}\Psi'_{j}(x)\Psi'_{i}(x)dx\, u_{i}^k\\ &\quad +\frac{1}{h_{t}}\sum_{i,j=0}^{n}g_{j}\int_0^{l}\Psi_{j}(x)\Psi_{i}(x) dx(u_{i}^{k+1}-u_{i}^k)-\alpha\sum_{i,j=0}^{n}v_{j}^{k+1}\int_0^{l}\Psi'_{j}(x)\Psi'_{i}(x)dx\, v_{i}^{k+1}\\ &\quad +\frac{1}{h_{t}}[N_T^k(u_{n}^{k+1}-u_{n}^k)-N_{B}^k(u_0^{k+1}-u_0^k)]. \end{align*} Thus, we can obtain the following equations in terms of matrices and vectors; \begin{align*} & \frac{1}{2h_{t}}[(\widetilde{v}^{k+1}-\widetilde{v}^k)^{T} \mathbf{M}(\widetilde{v}^{k+1}-\widetilde{v}^k) +(\widetilde{v}^{k+1}){}^{T}\mathbf{M}\widetilde{v}^{k+1} -(\widetilde{v}^k){}^{T}\mathbf{M}\widetilde{v}^k]\\ & =-\frac{c}{2h_{t}}[(\widetilde{u}^{k+1}-\widetilde{u}^k)^{T} \mathbf{K}(\widetilde{u}^{k+1}-\widetilde{u}^k)+(\widetilde{u}^{k+1}){}^{T} \mathbf{K}\widetilde{u}^{k+1}-(\widetilde{u}^k){}^{T}\mathbf{K}\widetilde{u}^k]\\ & \quad+\frac{1}{h_{t}}[(\widetilde{f})^{T}\mathbf{M}\widetilde{u}^{k+1} -(\widetilde{f})^{T}\mathbf{M}\widetilde{u}^k]-\alpha(\widetilde{v}^{k+1})^{T} \mathbf{K}\widetilde{v}^{k+1}\\ & \quad+\frac{N_T^k}{h_{t}}[(u_{n}^{k+1}+l-\varphi_T)-(u_{n}^k+l-\varphi_T)] -\frac{N_{B}^k}{h_{t}}[(u_0^{k+1}-\varphi_{B})-(u_0^k-\varphi_{B})]\\ & =-\frac{c}{2h_{t}}[(\widetilde{u}^{k+1}-\widetilde{u}^k)^{T} \mathbf{K}(\widetilde{u}^{k+1}-\widetilde{u}^k)+(\widetilde{u}^{k+1}){}^{T} \mathbf{K}\widetilde{u}^{k+1}-(\widetilde{u}^k){}^{T}\mathbf{K}\widetilde{u}^k]\\ & \quad+\frac{1}{h_{t}}[(\widetilde{f})^{T}\mathbf{M}\widetilde{u}^{k+1} -(\widetilde{f})^{T}\mathbf{M}\widetilde{u}^k]-\alpha(\widetilde{v}^{k+1})^{T} \mathbf{K}\widetilde{v}^{k+1}\\ & \quad-\frac{1}{h_{t}}[N_T^k(u_{n}^k+l-\varphi_T)-N_{B}^k(u_0^k-\varphi_{B})]. \end{align*} Since $\mathbf{M}$ is a positive definite matrix and $\mathbf{K}$ is a semipositive definite matrix, recalling the energy function in the fully discrete case, we can obtain the following inequality; \begin{equation} \begin{aligned} E^k & \ge E^{k+1}+\frac{1}{2}(\widetilde{v}^{k+1} -\widetilde{v}^k){}^{T}\mathbf{M}(\widetilde{v}^{k+1} -\widetilde{v}^k)+\frac{c}{2}(\widetilde{u}^{k+1} -\widetilde{u}^k)^{T}\mathbf{K}(\widetilde{u}^{k+1}-\widetilde{u}^k) \\ &\quad +[N_T^k(u_{n}^k+l-\varphi_T)-N_{B}^k(u_0^k-\varphi_{B})] +\alpha h_{t}(\widetilde{v}^{k+1})^{T}\mathbf{K}\widetilde{v}^{k+1} \\ &\quad \geq E^{k+1}-N_{B}^k(u_0^k-\varphi_{B}) +\alpha h_{t}(\widetilde{v}^{k+1})^{T}\mathbf{K}\widetilde{v}^{k+1}. \end{aligned}\label{eq:ineq_fully-1} \end{equation} Now, we use the extra equation \eqref{eq:fully_extra-1} to switch the second term in \eqref{eq:ineq_fully-1} into the following; $N_{B}^k(u_0^k-\varphi_{B})=N_{B}^k(u_0^{k+1}-\varphi_{B}-h_{t}v_0^{k+1}) =-h_{t}N_{B}^kv_0^{k+1}$. By using the telescoping sum from time $t=t_0$ to $t=t_{k}$, it follows from \eqref{eq:energy_cond-1} that \[ E^{0}\ge E^k+\alpha\sum_{\iota=1}^k(\widetilde{v}^{\iota})^{T} \mathbf{K}\widetilde{v}^{\iota}+\sum_{\iota=1}^kN_{B} ^{\iota-1}v_0^{\iota}\geq E^k, \] as desired. \end{proof} \subsection{Numerical results and discussion} \label{sub:Numerical-Results-and} In this Subsection, the numerical results (simulations) are presented and discussed. We simulate the almost elastic case ($\alpha=0.0001$) and the viscoelastic case ($\alpha=0.01$). For both cases we use the data displayed in Table~\ref{tab:Data-1}. Note that $g=g_{i}$ with $0\leq i\leq n$ and the unit of measure shall not be considered in our simulations. We also provide numerical evidences for a pure elastic case ($\alpha=0$) in support for the conclusions in the paper~\cite{ps:trcfaler}, although the existence results for the purely elastic case are not mentioned in this paper. \begin{table} \caption{Data with $\alpha=0.0001$ and $\alpha=0.01$} \label{tab:Data-1} \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|} \hline $l$ & $u^{0}$ & $u_{t}^{0}$ & $T$ & $\varphi_T$ & $\varphi_{B}$ & $h_{x}$ & $h_{t}$ & $c$ & $g$\tabularnewline \hline $1$ & $[0, 1]$ & $-20$ & $1.0$ & $3$ & $0$ & $0.0002$ & $10^{-5} \text{ or }10^{-6}$ & $100$ & $-9.81$\tabularnewline \hline \end{tabular} \end{table} Before we show any numerical results, we consider the Courant-Friedrichs-Lewy (CFL) condition $(\sqrt{c}\, h_{t})/h_{x}\leq1$. The CFL condition is necessary for the convergence of a finite difference scheme for hyperbolic PDEs. Therefore our selection of $h_{t}, h_{x}$ for both simulations conform to the CFL condition and may be helpful to obtain numerically stable results. Numerical simulations were preformed using Matlab R2010b on a Windows 7 workstation computer with an Intel Core i5 650 processor running at 3.20 GHz. Since the construction of the mass and stiffness matrices have nonzero entries only on the main, upper, and lower diagonals we chose to use the sparse matrices in order to save system memory. From Table \ref{tab:Data-1} we can see that each matrix in our linear system is of size $5000\times5000$ and each vector is therefore of size $5000$. The linear system is solved by performing a direct method, Gaussian elimination. Thanks to the use of a sparse matrix which greatly reduces computation time dealing with matrix operations we are able to complete each time step in about $0.0085$ seconds. Over the time interval $[0, T]$, an adaptive method is used to maintain reasonable computation time. When the rod is not in contact with either obstacle we take a larger time step size, $h_{t}=10^{-5}$, because it is more important to accurately show what happens to the rod as it contacts the obstacles (the time step size $h_{t}=10^{-6}$ used) than how it travels between them. \begin{figure}[ht] \begin{center} \includegraphics[width=0.8\textwidth]{fig3a} \\ % elastic_small \includegraphics[width=0.8\textwidth]{fig3b} % visco_small \end{center} \caption{Motion of two rods: $\alpha=0.0001$ (top) and $\alpha=0.01$ (bottom)} \label{fig:sim1_Position-1} \end{figure} Shown in Figure~\ref{fig:sim1_Position-1} is the actual displacement, $w$ for both cases over the whole simulation. The almost elastic rod $(\alpha=0.0001)$ falls as expected toward the bottom obstacle, and when contact occurs, the rod starts to undergo deformation by being compressed. As time passes we see the rod expand and lift off the bottom obstacle and regain its original length. A similar effect can be seen during the contact on the top obstacle and subsequent contacts on either obstacle. As we see the viscoelastic case $(\alpha=0.01)$, similarities to the almost elastic simulation can be drawn, but it is of greater interest to investigate the changes in the model caused by moving from an almost elastic to a viscoelastic rod. When the first contact occurs, we notice that the degree of deformation is nowhere near that of the almost elastic case. Over the duration of the simulation we can see that the viscosity of the rod acts as a damper that retards velocity causing it contact the obstacles less frequently than in the almost elastic case. Each subsequent contact occurs later in the simulation than its corresponding contact in the almost elastic simulation. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig4a} \quad % Paper_elastic/energy \includegraphics[width=0.45\textwidth]{fig4b} \\ % Paper_elastic/energy (A) $\alpha=0.0001$ \hfil (B) $\alpha=0.01$ \end{center} \caption{Energy function} \label{fig:Sim1_EnergyGraph-1} \end{figure} Figure~\ref{fig:Sim1_EnergyGraph-1} displays a plot of the energy function \eqref{eq:Energy-1} for the two simulations which supports numerical stability. From the numerical observation, we would not have to impose the condition \eqref{eq:energy_cond-1}, due to the fact that the contact forces and the velocity are orthogonal. This has been already justified theoretically in Subsection~\ref{sub:Numerical-formulations}. Lemma~\ref{lem:fully_energy_bd-1} provides validity to our simulation as evidence of its numerical stability. While the rod is not in contact with either obstacle and regaining its original length, we see the energy function decreasing less and less returning to a more stable constant state. Comparing the two energy functions, the graph (A) is more flat than (B), because of a bigger viscous quantity used in (B). \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig5a} \quad % force \includegraphics[width=0.45\textwidth]{fig5b} \\ % Paper_elastic/force_zoom (A) \hfil (B) \end{center} \caption{Contact force for $\alpha=0.0001$} \label{fig:Sim1_ContactGraph-1} \end{figure} In Figure~\ref{fig:Sim1_ContactGraph-1} (A) shows the graph of the contact force for the almost elastic case. (B) is only a zoomed in version of (A). It is zoomed in to show the contact force on the bottom obstacle in more detail. In Figure~\ref{fig:Sim2_ContactGraph-1} (A) and (B) show the graph of the contact force for the viscoelastic case. As when using an almost elastic rod there is a non-zero contact force only when there is contact, and the magnitude of the contact force on the upper obstacle is much higher than the contact force on the lower obstacle. This may happen due to the fact that our physical configuration is not symmetric with respect to the origin. In Figure~\ref{fig:Sim2_ContactGraph-1} (B) shows in more detail the contact force on the lower obstacle. This difference in magnitude is only magnified during the viscoelastic simulation with the contact force's magnitude for each obstacle being larger than the corresponding magnitude from the almost elastic case. We observe that contact force on the lower obstacle is not as uniform as it was in the almost elastic case. \begin{figure}[ht] \begin{center} \includegraphics[width=0.46\textwidth]{fig6a} \quad % Paper_visco/force \includegraphics[width=0.45\textwidth]{fig6b} \\ %Paper_visco/force_zoom (A) \hfil (B) \end{center} \caption{Contact force for $\alpha=0.01$} \label{fig:Sim2_ContactGraph-1} \end{figure} Longitudinal waves of an almost elastic rod is depicted in Figure~\ref{fig:sim1_velocity-1}. While the rod is not in contact with a obstacle, in (B) and (D) of Figure~\ref{fig:sim1_velocity-1}, the waves are sampled every $600$ time steps, but because the duration of contact is short in comparison, the waves in (A) and (C) of Figure~\ref{fig:sim1_velocity-1} are sampled every 186 time steps. Note that before the initial velocity is uniform across the rod and is not presented in a graph. We see in (A) of Figure~\ref{fig:sim1_velocity-1} over the local time interval $[0.05, 0.05637)$, that as the rod contacts the bottom obstacle waves propagate across the compressed rod with the end that is in contact with the obstacle having a zero velocity. As velocity becomes all positive we see (B) of Figure~\ref{fig:sim1_velocity-1}, the local time interval $[0.05637, 0.1621)$, which details the rod's ascent to the top obstacle. When the rod contacts the upper obstacle, in (C) of Figure~\ref{fig:sim1_velocity-1} over the local time interval $[0.1621, 0.1687]$, we can see that velocity becomes zero near the top of the rod. In (D) of Figure~\ref{fig:sim1_velocity-1} during the local time period $[0.1687, 0.277]$, the waves of velocity traverse the length of the rod as it falls toward the bottom obstacle. As shown in (A) of Figure~\ref{fig:sim1_Position-1} there is more than one contact on either obstacle, but since each subsequent contact is similar to the corresponding previous contact, the graphs of wave propagation for the later contacts are omitted. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig7a} \quad %Paper_elastic/vel1 \includegraphics[width=0.45\textwidth]{fig7b} \\ %Paper_elastic/vel2 (A)\hfil (B)\\ \includegraphics[width=0.45\textwidth]{fig7c} \quad %Paper_elastic/vel3 \includegraphics[width=0.45\textwidth]{fig7d} \\ %Paper_elastic/vel4 (C) \hfil (D) \end{center} \caption{Longitudinal waves of the rod for $\alpha=0.0001$} \label{fig:sim1_velocity-1} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig8a} \quad %Paper_visco/vel1 \includegraphics[width=0.45\textwidth]{fig8b} \\ %Paper_visco/vel2 (A)\hfil (B)\\ \includegraphics[width=0.45\textwidth]{fig8c} \quad %Paper_visco/vel3 \includegraphics[width=0.45\textwidth]{fig8d} \\ %Paper_visco/vel4 (C) \hfil (D) \end{center} \caption{Longitudinal waves of the rod for $\alpha=0.01$} \label{fig:sim2_velocity-1} \end{figure} (A)--(D) in Figure~\ref{fig:sim2_velocity-1} contain velocity along the rod for the viscoelastic case. The waves in (B) and (D) of Figure~\ref{fig:sim2_velocity-1} are sampled every 625 time steps, and the waves in (A) and (C) of Figure~\ref{fig:sim2_velocity-1} are sampled every 200 time steps after a contact occurs. In Figure~\ref{fig:sim2_velocity-1} (A) over the time period $[0.05, 0.05625]$, as the rod contacts the bottom obstacle we see waves of velocity appear, but they are not as extreme in magnitude as they were in (A) of Figure~\ref{fig:sim1_velocity-1}. We also can see that the wave itself is much smoother as it propagates along the rod. This is evident in the near flat lines that appear in (B) and (D) of Figure~\ref{fig:sim2_velocity-1}, over the two time intervals $[0.05625, 0.1866)$ and $[0.1929, 0.3523]$ respectively, as the rod is not undergoing deformation when traveling between the obstacles. We can see in (C) of Figure~\ref{fig:sim2_velocity-1} over time interval $[0.1866, 0.1929)$ when contact occurs on the top obstacle the wave behaves similarly to the previous contact on the bottom obstacle. Evidence shows that viscosity acts as a damper by decreasing the magnitude of velocity dramatically after each contact more so than in Figure~\ref{fig:sim1_velocity-1}. \begin{table}[ht] \caption{Numerical results for contact duration} \label{tab:PS_case1} \begin{tabular}{|c|c|c|c|} \hline $c$ & Expected contact duration & Observed contact duration & Rel. error\tabularnewline \hline $1$ & $2$ & $0.0633$ & $0.96835$\tabularnewline \hline $100$ & $0.02$ & $0.00234$ & $0.883$\tabularnewline \hline $200$ & $0.01$ & $0.00448$ & $0.552$\tabularnewline \hline $500$ & $0.004$ & $0.00285$ & $0.2875$\tabularnewline \hline $1000$ & $0.002$ & $0.00202$ & $0.01$\tabularnewline \hline \end{tabular} \end{table} Another interesting numerical experiment is performed to support important theoretical results in the paper~\cite{ps:trcfaler}. The results for contact characteristics can be summarized in the notation of our model with $h$ being the distance from the bottom of the rod to the lower obstacle: the first case states that if $2gl