\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 $a**0$ 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**