\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{amssymb, color}
\usepackage{graphicx}
\usepackage{subfig}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2012 (2012), No. 194, pp. 1--42.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2012 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2012/194\hfil Dynamic contact of two Gao beams]
{Dynamic contact of two Gao beams}

\author[J. Ahn, K. L. Kuttler,  M. Shillor \hfil EJDE-2012/194\hfilneg]
{Jeongho Ahn, Kenneth L. Kuttler, Meir Shillor}  % in alphabetical order

\address{Jeongho Ahn \newline
Department of Mathematics and Statistics \\
Arkansas State University \\
State University,  AR 72467, USA}
\email{jahn@astate.edu}

\address{Kenneth L. Kuttler \newline
Department of Mathematics \\
Brigham Young University, Provo, UT 84602, USA}
\email{klkuttle@math.byu.edu}

\address{Meir Shillor \newline
 Department of Mathematics and Statistics \\
 Oakland University \\
 Rochester, MI 48309, USA}
\email{shillor@oakland.edu}

\thanks{Submitted June 7, 2012. Published November 6, 2012.}
\subjclass[2000]{74M15, 35L86,74K10, 74M15, 35L70}
\keywords{ Dynamic contact; Gao beam; mechanical joint; normal compliance;
\hfill\break\indent Signorini condition; weak solutions; numerical scheme}

\begin{abstract}
 The dynamic  contact of two nonlinear Gao beams that are connected
 with a joint  is modeled, analyzed, and numerically simulated.
 Contact is modeled with either (i) the normal compliance condition, or 
 (ii) the  unilateral Signorini condition. The model is in the form 
 of a variational  equality in case (i) and a variational  
 inequality in case (ii). The existence of the unique variational solution 
 is established for  the problem with normal compliance and the existence 
 of a weak solution  is proved in case (ii). 
 The solution in the second case is obtained as a limit
 of the solutions of the first case when the normal compliance stiffness
 tends to infinity.  A  numerical algorithm for the problem is constructed 
 using  finite elements and a mixed time discretization.  
 Simulation results, based on the  implementation of the algorithm,  
 of the two cases when the horizontal traction
 vanishes or when it is sufficiently large to cause buckling, are presented.
 The spectrum of the vibrations, using the FFT,  shows a rather noisy system.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\numberwithin{figure}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{algorithm}[theorem]{Algorithm}
\allowdisplaybreaks

\section{Introduction}\label{sec:intro}

This work models the vibrations of  two uniform elastic or viscoelastic nonlinear Gao beams
that are connected with a mechanical joint that has a gap or clearance. The Gao nonlinear
beam was introduced in \cite{Gao96, Gao00} to allow for the investigation of the vibrations of
beams about buckled  states. The model for the process  is in the form of a system of
nonlinear partial  differential equations that are coupled  across the joint. We establish
the existence of a weak solution to the system.  Then, we describe in detail a numerical
algorithm for the problem, based on mixed finite elements,  and present two typical
numerical simulations.  One deals with vibration when there is buckling, the other
one without buckling. This is a step in our investigation of the noise properties of
mechanical systems and the transmission of vibrations across joints.

Noise and vibration characteristics of mechanical assemblies, parts, and
components are very important in industry. Currently, the automotive industry
has considerable interest in the dynamic vibrations of mechanical systems, stemming
from the customer's perception of noise characteristics of cars. The topic is under
extensive research, mainly experimentally and in part computationally in the industry.
However, there exists no complete methods for such investigation, and full mathematical
analysis of typical automotive problems is still out of reach. Therefore, most of the
models use various ad hoc assumptions, many of which are too simplistic, and thus
unrealistic, and not very useful.

Here, we introduce a simple one-dimensional geometrical setting to avoid certain
mathematical complications related to two- or three-dimensions. On the other hand,
the physical content of the problem is not simplified. Our long-term goal is to
understand the transmission of vibrations across mechanical joints. This will help
industry in the design of parts, assemblies and systems with acceptable noise and
vibration characteristics. Eventually, theory will replace the current practice of
finding that the system needs to be modified at the testing or even production stages,
because of unacceptable noise generation.

Preliminary steps in this program have been taken in \cite{AS06, ASW96, KPSZ01, KS98ts}
where various aspects of vibrations and contact of the Euler-Bernoulli linear beams were
investigated. In particular the modeling, analysis, and simulations of the vibrations of a
beam between two stops can be found in \cite{AS06, KPSZ01, KS98ts}.  Here, we extend
this study to include nonlinear beams, a fact that makes the analysis and computations
much more complicated, since more specialized ideas and methods must be used.

Various mathematical and numerical results about the nonlinear Gao beam were reported
in  \cite{ahn:dfcgts, adbps:asndb, AMbS09, KPS11, MFM08, MbS08, MN11,SG12}.
The main interest there  was to study the buckling of the beam and its vibrations about a buckled state.
Contact conditions, even when the governing equations are linear, make the problems
nonlinear, which led to the currently developing  Mathematical Theory of Contact
Mechanics (\cite{KO, SST04} and references therein). Having a nonlinear equation
and nonlinear contact conditions, again, makes the problem much harder, both in terms
of its analysis and computationally. Computational issues related to contact can be found
in \cite{Du03, KO, MO, PS98, Wri02}

Recent numerical schemes  for dynamic contact  problems of linear beams have been
developed in \cite{sa:ebdcn, a:topdac, a:vtbclf}. In this paper, we propose a more advanced
numerical algorithm that is more technical since we have to handle properly
the nonlinear term in the equation. The theory of its convergence is investigated in
\cite{ahn:dfcgts}. Issues of control of the Gao beam can be found in \cite{Gao06, GaoR98}.

We assume, for the sake of generality, that the beams are viscoelastic, and describe the
viscosity by a  short-term memory of the Kelvin-Voigt type. Moreover,
we also deal with the the inviscid case by passing to the limit when the viscosity vanishes.
The dynamic contact is modeled by either the \emph{Signorini nonpenetration condition}
which describes \emph{unilateral contact} of two perfectly rigid stops, or by the
\emph{normal compliance condition,} which describes reactive stops.  This work
is an extension of \cite{AMbS09} where the problem of the vibrations of one Gao
beam between  two reactive stops was investigated. Some numerical simulations
of the vibrations of two Euler-Bernoulli linear beams with a reactive joint problem
can be found in \cite{KPSZ01}, where the method of lines was used.

Next, we present a  rather complex  numerical algorithm for the problems. The
complexity arises from the nonlinearity of the beam equation
coupled with the nonlinear contact condition. The main idea
of the numerical algorithm is to use time discretization  and the combination
of the central difference formula and the mixed finite element methods in space.
Two simulations for the vibrations of the system when buckling takes place and when it
is absent, are depicted, showing that the numerical scheme seem to be work
reasonably well. The finite elements naturally fit with the weak formulation of the
problem, since this is their usual setting. We note that, adding to the complexity, it
follows from the results of \cite{MoShow83} that even a much simplified
system of this type may exhibit chaotic behavior.

The rest of the work is organized as follows. In Section \ref{sec:model} we present the physical
setting, and the classical model. We describe either the
Signorini or the normal compliance contact conditions. The beams are assumed
to be either viscoelastic or elastic. Since we allow the left beam to be fixed to an
oscillating support, we introduce in Section \ref{sec:var} a change of variable that
allows to consider a problem in which the left end is stationary.
To present the variational formulation we introduce the relevant function spaces, and
state our main existence results, the existence of the unique weak solution
for the problem with normal compliance, Theorem \ref{t31}, and an existence
of a weak solution for the problem with the Signorini condition, Theorem \ref{t32}.

In Section \ref{sec:nc} we investigate the model with the
normal compliance condition, which in addition to having its own merit, may be
considered as regularization of the Signorini condition. We prove Theorem \ref{t31}.
In Section \ref{sec:Signorini} we establish Theorem \ref{t32}, in which a solution
is obtained as a limit of a sequence of solutions found in Section \ref{sec:nc}.
In Section \ref{sec:novisc} we show in Theorem \ref{t62}
that the problems with the normal compliance condition or the Signorini condition with
vanishing viscosity posses a weak solution, too.
We obtain the necessary a priori estimates (Theorem \ref{t61})  on the solutions
with viscosity, and pass to the limit when the viscosity vanishes and the normal compliance
stiffness coefficients tend to infinity.


We describe in Section \ref{sec:algorithm} a numerical algorithm for the model, based on the
mixed finite elements method. The central idea of the algorithm is to drive the
recursive formula which enables to compute next time step solutions of the
linear system per each step satisfying the contact conditions.
The results of our simulations are given in Section \ref{sec:Numerical_results}.
We present two sets of  simulations without buckling ($p=0$) and with the buckling of the right beam
($p=895$). We also  show the Fast Fourier Transform of the motion of the right end of the left beam,
and the noise, i.e., the spectrum of frequencies that it generates.

 The paper is summarized in Section \ref{sec:conclude}, where some
 unresolved issues and future work are  discussed, too.

\section{The Model }
\label{sec:model}

This section presents the classical formulation of the model, a change of variable,
the necessary function spaces, the variational form, and a statement of our theoretical results.

The physical setting, shown in Figure \ref{fig1} and  a detail of the joint in Figure \ref{fig2},
is as follows. Two uniform Gao beams are connected at at a joint with a
clearance.  The left beam is  attached rigidly at its left end to a supporting
device that may oscillate with time, possibly periodically.
The right end of the second beam is clamped to a rigid support, say a wall. The
motion of the  right end of the left beam is constrained by the motion of the left end of
beam 2, where it may move only within the given clearance. The
expanded view of the joint is depicted in Figure \ref{fig2}, where $g$ is the
\emph{clearance} or the \emph{gap}. For the sake of generality we assume that the joint may be
asymmetrical, and let $g=g_1+g_2$ ($g_1, g_2>0$), where $g_1$ is the upper
clearance and  $g_2$ is the lower one, when the system is at rest.

\begin{figure}[ht]
\begin{center}
%\begin{picture}(40,50)
\begin{picture}(280,50)(0,-20)
\thicklines
\qbezier(0,10)(30,10)(60,5)
\qbezier(60,5)(90,0)(120,0)
\qbezier(120,0)(150,0)(180,5)
\qbezier(180,5)(210,10)(240,10)
\qbezier(0,0)(10,0)(60,-5)
\qbezier(60,-5)(90,-10)(120,-10)
\qbezier(120,-10)(150,-10)(180,-5)
\qbezier(180,-5)(210,0)(240,0)
\put(120,-10){\line(0,1){10}}
\put(240,-20){\line(0,1){40}}
\put(0,5){\vector(1,0){285}}
\put(0,-20){\line(0,1){40}}
\put(20,10){\vector(0,1){25}}
\put(25,20){$u_1$}
\put(100,0){\vector(0,1){25}}
\put(85,15){$f_1$}
\put(120,3){\line(0,1){6}}
\put(120,-20){\color[rgb]{0,0,1} $l_*$}
\put(220,10){\vector(0,1){25}}
\put(258,8){\vector(-1,0){18}}
\put(225,20){$u_2$}
\put(272,-10){$x$}
\put(2,-20){$0$}
\put(242,-20){$1$}
\put(249,14){\color[rgb]{0,0,1} $p$}
\end{picture}
%}\end{picture}}
\end{center}
\caption{The two deflected Gao beams}
\label{fig1}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{picture}(180,60)(-30,-40)
\thicklines
\put(-30,0){\line(1,0){80}}
\put(-30,-30){\line(1,0){80}}
\put(50,-30){\line(0,1){30}}
\put(40,8){\line(0,1){15}}
\put(40,23){\line(1,0){20}}
\put(55,-35){\line(0,1){43}}
\put(40,8){\line(1,0){15}}
\put(40,-35){\line(1,0){15}}
\put(40,-40){\line(1,0){20}}
\put(60,-40){\line(0,1){63}}
\put(40,-40){\line(0,1){5}}
\put(60,-30){\line(1,0){90}}
\put(60,0){\line(1,0){90}}
\thinlines
\put(37,9){\line(-1,0){10}}
\put(37,-38){\line(-1,0){10}}
\put(35,-15){\line(-1,0){10}}
\put(30, 0){\vector(0,1){9}}
\put(30,-30){\vector(0,-1){8}}
\put(15,-38){$g_2$}
\put(15,6){$g_1$}
\end{picture}
\end{center}
\caption{The joint; $g=g_1+g_2$ is the gap}
\label{fig2}
\end{figure}


The area-centers of gravity of the beams in their (stress free) dimensionless
reference configurations coincide with the intervals $0\leq x\leq l_{*}$, and $
l_{*}\leq x\leq 1$, respectively.

A prescribed traction $p=p(t)$ acts at the right end $x=1$;  the beams' edges
in the joint are  assumed to be permanently in frictionless contact. Adding
friction may be of interest  in a future study.


We use dimensionless variables and denote by $u_1=u_1(x,t)$ and $
u_2=u_2(x,t)$ the vertical displacements of the beams, and by $\sigma
_i=\sigma _{i}(x,t)$  the shear stresses, for $i=1,2$, where here and below
$i=1$ indicates the left beam and $i=2$ the right one. The equations of motion are
\begin{gather}
\rho _1u_{1tt}- \sigma _{1x} =\rho _1 f_1, \quad 0<x<l_*, \label{21}\\
\rho _2u_{2tt}- \sigma _{2x} =\rho _2 f_2, \quad l_*<x<1,  \label{22}
\end{gather}
 where $f_i$ denote the density (per unit mass) of the vertical
applied forces, and $\rho _i$ are the beams' densities (per unit length).
Also, the subscripts $x$ and $t$ indicate partial derivatives. The time interval
of interest is $0\leq t \leq T$. The
beams are nonlinear and the \emph{constitutive relations} are,
\begin{equation}
\sigma _i(x,t)=-k_i  u_{ixxx}(x,t)-d_iu_{itxxx}( x,t)
+\frac 1 3 {a_i} u_{i x}^3-\overline{\nu_i}p  u_{i x},  \label{23}
\end{equation}
where  the $k_i$ are the coefficients of elasticity, $d_i$ are the coefficients of viscosity,
 $a_i$ are the Gao coefficients,
and $\overline{\nu_i}$ are positive coefficients related to the compressive or tensile traction,
 for $i=1,2$.  By assumption, the beams  are uniform and so the various coefficients are
 positive constants.

The initial conditions are
\begin{gather}
u_1(x,0)=u_{01}(x),\quad u_{1t}(x,0)=v_{01}(x)\quad 0\leq x\leq l_{*},
\label{24} \\
u_2(x,0)=u_{02}(x),\quad u_{2t}(x,0)=v_{02}(x)\quad l_{*}\leq x\leq 1,
\label{25}
\end{gather}
where $u_{0i}$ and $v_{0i}$ are prescribed functions
representing the initial deflections and velocities of the beams,
respectively.

We suppose that the left beam is rigidly attached at its left end to a device that may be
oscillating in time with vertical displacement $\phi$, thus, for $0\leq t \leq T$,
\begin{equation}
u_1(0,t)=\phi(t), \quad \text{and}\quad u_{1x}(0,t)=0.              \label{26}
\end{equation}
Here $\phi$ is known, and we have more to say about it below.
The right end of the right beam is clamped, so
\begin{equation}
u_2(1,t)=0,\quad \text{and}\quad u_{2x}(1,t)=0.                \label{27}
\end{equation}


We turn to describe the \emph{contact process} in the  joint. First, we model it with the
\emph{\it  normal compliance} contact condition (see, e.g., \cite{KO, Kl2, KMS, Ken97,
MO, SST04} and  the references therein), which describes reactive stops. The contact shears
satisfy
\begin{equation}
\sigma=\sigma(t)\equiv \sigma _1(l_{\ast },t)=-\sigma _2(l_{\ast },t).
 \label{28}
\end{equation}
When there is no contact  the shear vanishes, thus,
\[
u_2-g_2<u_1< u_2+g_1\; \Longrightarrow \sigma=0.
\]
The reaction of the stops takes place only when the displacement of the end
of the left beam exceeds the clearance,  is proportional to the interpenetration,
and when
\[
u_2( l_{\ast},t)  +g_1 \leq u_1(l_{\ast },t),
\]
the reaction is
\[
\sigma = -\lambda_1( u_1(l_{\ast },t) -u_2(l_{\ast},t),
-g_1)
\]
where $\lambda_1>0$ is the normal compliance stiffness of the top stop, and the reaction
(acting on the edge of the left beam) is downward. When
\[
u_1(l_{\ast },t) \leq u_2(l_{\ast},t)  -g_2,
\]
the reaction is
\[
\sigma = \lambda_2 ( u_2(l_{\ast },t) -u_1(l_{\ast },t) -g_2),
\]
where $\lambda_2>0$ is the normal compliance stiffness of the lower stop,
and the reaction
(acting on the edge of the left beam) is upward.
The discussion of the reaction of the stops can be summarized as follows,
\begin{equation} \label{29}
\begin{split}
\sigma(t)&=\Lambda(u_1(l_{*}, t)-u_2(l_{*}, t))\\
&\equiv -\lambda_1 ( u_1(l_{\ast },t) -u_2(l_{\ast},t)
-g_1) _{+}+ \lambda_2 ( u_2(l_{\ast },t) -u_1(l_{\ast },t) -g_2) _{+},
\end{split}
\end{equation}
and $(f)_+=\max(f, 0)$ is the positive part of the function.

We note that more general normal compliance power laws were introduced and used in
\cite{KO, Kl1, Kl2, KMS, Ken97, MO, SST04}, among many other references.

Finally, we assume that the
ends do not exert moments on each other, and thus,
\begin{equation}
w_{ixx}=0,\quad x=l_{*}.  \label{210}
\end{equation}

Given the problem data; i.e., $u_{0i}=u_{0i}(x)$,
$v_{0i}=v_{0i}(x)$, $(i=1,2)$,  $p(t)$,
and  $\phi=\phi(t)$, the physical parameters and the constitutive relations
\eqref{23}, the \emph{classical formulation} of \emph{the dynamic
contact at a joint of two Gao beams with the normal compliance contact condition} is:
\begin{problem}
\label{pbm21}
Find a pair of functions $(u_1,\,u_2)$ such that \eqref{21}--\eqref{210} hold.
\end{problem}
The problem, in its variational form, is studied in Section \ref{sec:nc} where we establish
the existence of its unique weak solution.



The second contact condition that is often used is the so-called \emph{
Signorini nonpenetration} or \emph{unilateral} condition in which the stops are
assumed to be \emph{perfectly rigid.} It is, essentially,  an idealization of the normal
compliance condition, and may be justified in certain settings. We discuss the
relationship between the conditions  below.

Since the stops are rigid right, the right  end of the left  beam must be within the
clearance of the left end of the right beam.  Thus,
\[
u_2(l_{\ast },t)-g_2\leq u_1(l_{\ast },t)\leq u_2(l_{\ast},t)+g_1.
\]
When strict inequalities hold, there is no contact, the
ends are free, and then $\sigma(t)=0.$
On the other hand, when the ends are in contact, the stresses are
equal but opposite, then \eqref{28} holds.

Next,  when contact takes place at the lower stop
$u_1(l_{\ast },t)=u_2(l_{\ast },t)-g_2$ and $\sigma \geq 0$, and when
the contact is at the upper stop, $u_1(l_{\ast },t)=u_2(l_{\ast},t)+g_1$
and $\sigma \leq 0.$ We may write it  concisely by using the indicator
function   $\chi$,
\begin{equation}
 \chi (r)=\begin{cases}
0 & \text{if }u_2-g_2\leq r \leq u_2+g_1, \\
+\infty & \text{ otherwise}.
\end{cases}  \label{211}
\end{equation}
Then, we let $\partial \chi $ denote the \emph{subdifferential} of $\chi $, i.e.,
\begin{equation}
\partial \chi (r)=\begin{cases}
0 & \text{if } u_2-g_2<r<u_2+g_1, \\
(-\infty ,0] & \text{if }r=u_2-g_2, \\
\lbrack 0,\infty ) & \text{if }r=u_2+g_1.
\end{cases}   \label{212}
\end{equation}
We note that  the effective domain of both $ \chi $
and $\partial \chi $ is $[ u_2-g_2, u_2+g_1]$.

The contact conditions may be summarized in the form of the
following set inclusion that represents the \textit{ complementary conditions}:
\begin{eqnarray}
u_2-g_2 \leq u_1\leq & u_2+g_1,  \label{213}\\
-\sigma \in &\partial \chi (u_1),  \label{214}
\end{eqnarray}
evaluated at $x=l_{*}$, together with \eqref{28}.


Given the problem data, i.e., $u_{0i}=u_{0i}(x),\,v_{0i}=v_{0i}(x),\,\,(i=1,2)$, $p$, and
$\phi=\phi(x)$, the physical parameters and the constitutive relations \eqref{23},
 the `classical' formulation of \emph{the dynamic contact of two Gao
beams with the Signorini contact condition} is:

\begin{problem}
\label{pbm22}
Find a pair of functions $(u_1,\,u_2)$ such that \eqref{21}--\eqref{28},
\eqref{210} and \eqref{213}--\eqref{214} hold.
\end{problem}


However, the Signorini condition is an over-idealization when dynamic contact
takes place \cite{SST04}. In the literature, the normal compliance condition is often
used as a regularization of the Signorini condition, although physically, the Signorini
condition seems to be an idealization of the normal compliance condition,  as there are
no perfectly rigid objects. Moreover, the Signorini condition introduces
sever mathematical difficulties in multidimensional dynamic contact problems,
 that reflect the fact that there are no perfectly rigid objects. Moreover, there
 are issues of energy conservation when using it \cite{sa:ebdcn, PS98}. The Signorini condition
 can be obtained (formally) as a limit when the
stiffness in the normal compliance condition tends to infinity. 
This has been established rigorously in many works, see, e.g., 
\cite{Jiri, KS02nc} and the references therein.
Indeed, in this work and establish the solvability of the problem with the Signorini
condition in the limit $\lambda_1, \lambda_2 \to \infty$. However, we do it
mainly to show that the normal compliance problems behave well even when the
stiffnesses are arbitrarily large and, also,  for the sake of completeness.


\section{Variational formulation and results}
\label{sec:var}

We proceed to a weak or variational formulation of the two problems, since
the solutions may not be sufficiently regular. Indeed, this is the case in the
problem with the Signorini condition since the velocity is discontinuous upon impact.
Moreover, there is regularity ceiling for the problem with normal compliance (see
\cite{KS04nc} for the details). More importantly, we use it since the
mathematical tools for weak solutions are much more advanced. Moreover,
the variational formulation is the basis for the finite element algorithm
presented in Section \ref{sec:algorithm}.


We turn to introduce the function spaces that are needed. However, first we
change the variable $u_1$ so as to have a zero boundary condition at $x=0$,
so that the function spaces we use below that incorporate the boundary conditions
do not depend on time. To that end, we introduce a smooth function $\Phi=\Phi(x, t)$,
for $0\leq x \leq l_{*}$, such that
\begin{equation} \label{31}
\begin{split}
\Phi(0, t)= \phi(t),\quad \Phi_x(0,t)=0,\\
\Phi(l_{*}, t)=\Phi_x(l_{*},t)=\Phi_{xx}(l_{*}, t)=\Phi_{xxx}(l_{*},t)=0.
\end{split}
\end{equation}
We note that the function
\[
\Phi(x,t)=\phi(t) \cos^{4}(\frac{\pi x}{2l_{*}})
\]
satisfies all these conditions. The new dependent variable
\begin{equation}
\overline{u}_1(x,t)=u_1(x,t)-\Phi(x,t), \label{31b}
\end{equation}
satisfies,
\begin{gather*}
\overline{u}_1(0,t)={\overline{u}_1}_x(0,t)=0,\\
\overline{u}_1(l_{\ast },t)={\overline{u}_1}_x(l_{\ast },t)= {\overline{u}_1}_{xx}(l_{\ast },t)=
{\overline{u}_1}_{xxx}(l_{\ast },t)=0.
\end{gather*}
These properties guarantee that the contact conditions in terms
of $\overline{u}_1$ remain the same as above.

We perform the same change of the initial condition $u_{10}$ and $v_{10}$.

The shear stress in the first beam may be written as
\begin{equation} \label{33}
\begin{split}
\sigma _1(x,t)
&=-k_1{\overline{u}_1}_{xxx}-d_1{\overline{u}_1}_{txxx}+
\frac{1}{3}{a_1}( {\overline{u}_1}_x+\Phi_x ) ^3\\
&\quad  -\overline{\nu}_1p {\overline{u}_1}_x+\Psi (x,t,p),
\end{split}
\end{equation}
where we set
\begin{equation}
\Psi (x,t,p(t))=-k_1\Phi_{xxx}  -d_1\Phi_{txxx}
- \overline{\nu}_1 p \Phi _ x.  \label{34}
\end{equation}
We note that
$\Psi (l_{*},t,p(t))=0$,
so that the change of the variable does not change the contact shear stress.

Next, we rewrite the force as
\[
\overline{f}_1=f_1-\ddot{\phi}(t)\cos^{4} \big( \frac{
\pi x}{2l_{\ast }}\big).
\]


For the sake of simplicity, we do not use the bar below, and so $u_1$
denotes $\overline{u}_1(x,t)$ and $f_1$ stands for $\overline{f}_1$, so that the
equation of motion \eqref{21} remains the same.
We also rescale the variables so that $\rho _1=\rho _2=1$ for the sake
of simplicity of the presentation.


We now proceed to the variational formulation of the problem. We use the
following notation:
\begin{gather*}
V_1\equiv \{ w\in H^2(0,l_{\ast}) : w( 0)
=w_x( 0) =0\}, \\
V_2\equiv \{ w\in H^2( l_{\ast },1) : w( 1)
=w_x( 1) =0\}.
\end{gather*}
We note that the transformation of $u_1$ allows us to make $V_1$
independent of time. Also,
\begin{equation*}
H_1\equiv L^2(0,l_{\ast}) ,\quad H_2\equiv L^2(l_{\ast },1).
\end{equation*}
Let $V\equiv V_1\times V_2$ and $H\equiv H_1\times H_2.$ Then,
identifying $H$ and $H'$, we may write,
\begin{equation*}
V\subseteq H=H'\subseteq V'.
\end{equation*}


Let $\mathcal{V}_{i}\equiv L^2( 0,T;V_{i})$,
$\mathcal{H}_{i}\equiv L^2( 0,T;H_{i})$, for
$i=1,2 $, and $\mathcal{V}\equiv L^2( 0,T;V) $,
and $\mathcal{H}\equiv L^2( 0,T; H) $,
for $T>0$.  We denote by $(v,v)_{H_{i}}$ the inner product
on $H_{i}$ and the norm by $|v|_{H_{i}}^2=(v,v)_{H_{i}}$,
for $i=1,2$. The inner products and norms on the other spaces are defined similarly.
On $\mathcal{V}$ we may use the equivalent norm
\begin{equation*}
\| \mathbf{v} \| =| \mathbf{v}_{xx}| _H,
\end{equation*}
and we do so when it is convenient. We also define
\[
W_1 \equiv \{ w\in H^{1}(0,l_{\ast}) :w( 0)=0\},\quad
W_2 \equiv \{ w\in H^{1}( l_{\ast },1) :w( 1)=0\},
\]
with similar notation as above.

Let $\mathbf{w}=(w_1,w_2)\in \mathcal{V}$, and let $\phi \in
C^{\infty }( [0,T] ) $ with $\phi ( T) =0.$
We multiply \eqref{21} with $w_1(x)\phi(t) $ and integrate (formally) over
$(0,l_{\ast })  \times (0,T)$, and integrating
by parts we obtain the following variational equations
for $u_1$,
\begin{equation} \label{35}
\begin{aligned}
&-( v_{01},w_1) _{H_1}\phi ( 0)
-\int_0^T( u_{1t},w_1) _{H_1}\phi '(t) dt
\\
&-\int_0^T(\sigma _1( l_{\ast }) w_1( l_{\ast})  +(
k_1u_{1xx}+d_1u_{1txx},w_{1xx}) _{H_1})\phi ( t) dt
\\
&+\frac{1}{3}a_1\int_0^T\Big( ( u_{1x}+\Phi _x)
^3,w_{1x}\Big) _{H_1}\phi ( t)\,dt \\
&+ \int_0^T(-\overline{\nu }_1p( t) ( u_{1x},w_{1x})
_{H_1} - ( \Psi_x,w_1) _{H_1})\phi ( t)\,dt
\\
&=\int_0^T( f_1,w_1) _{H_1}\phi ( t) .
\end{aligned}
\end{equation}
Similarly,  we multiply \eqref{22} with $w_2(x)\phi(t) $,
 integrate over $(l_{\ast },1)\times ( 0,T)$ and find
\begin{equation} \label{36}
\begin{aligned}
&-( v_{02},w_2) _{H_2}\phi ( 0)
-\int_0^T( u_{2t},w_2) _{H_1}\phi '(t) dt \\
&+\int_0^T( \sigma _2( l_{\ast }) w_2( l_{\ast
})  +( k_2u_{2xx}+d_2u_{2txx},w_{2xx}) _{H_2}) \phi ( t) dt
\\
&+ \int_0^T( ( \frac{1}{3}a_2u_{2x}^3-\overline{\nu }
_2p( t) u_{2x}) ,w_{2x}) _{H_2}\phi (t) dt\\
&=\int_0^T( f_2,w_{2s}) _{H_2}\phi (t) dt
\end{aligned}
\end{equation}

Here, the shear stresses $\sigma _1( l_{\ast }) $ and $\sigma
_2( l_{\ast }) $ satisfy \eqref{28} and either the normal
compliance condition \eqref{29} or the Signorini condition \eqref{213}
and \eqref{214}.

The convex set of admissible pairs of displacements is
\begin{equation} \label{37}
\begin{split}
K=\big\{&\mathbf{w}=(w_1,w_2)\in V: w_1(0)=0, w_2(1)=0,  \\
&\quad w_2(l_{\ast })-g_2\leq w_1(l_{\ast})\leq w_2(l_{\ast })+g_1\big\}.
\end{split}
\end{equation}
and we set $\mathcal{K}=L^2(0, T; K)$.

We make the following assumptions on the problem data:
\begin{itemize}
\item[(A1)]  The initial displacements $\mathbf{u}_0=(u_{01}, u_{02})$  and velocities
$\mathbf{v}_0=(v_{01}, v_{02})$ satisfy
\begin{equation}\label{38}
\mathbf{u}_0\in K,\quad \mathbf{v}_0 \in H.
\end{equation}

\item[(A2)]  The traction $p=p(t)$ satisfies
\begin{equation}\label{39}
p\in C^{1}([0, T]), \quad |p|, |p'|\leq p_0.
\end{equation}

\item[(A3)]  The motion of the left support  $\phi=\phi(t)$ satisfies
\begin{equation} \label{310}
\phi\in C^2([0, T]),\quad u_{01}(0)=\phi(0).
\end{equation}

\item[(A4)]  The body forces  $\mathbf{f}=(f_1(x, t), f_2(x,t))$ satisfy
\begin{equation}\label{311}
\mathbf{f} \in H.
\end{equation}


\item[(A5)]  The coefficients
\begin{equation} \label{312}
g_1,  g_2,  k_{i},  d_{i}, a_{i},  \overline{\nu}_{i},  \lambda_{i}, \quad i=1,2,
\end{equation}
are positive constants.
\end{itemize}

The main theoretical result in this work is the following existence and
uniqueness result for the problem with the normal compliance contact condition.

\begin{theorem} \label{t31}
Assume that {\rm (A1)--(A5)} hold. Then,  there exists a unique solution
$\mathbf{u} =(u_1, u_2) \in \mathcal{V}$,
such that $\mathbf{u}(\cdot, 0)=\mathbf{u}_0$,  $\mathbf{u}_t(\cdot, 0)=\mathbf{v}_0$, and
the variational equations \eqref{35} and \eqref{36} hold. The stress $\sigma(t)$
satisfies \eqref{29} in a weak sense.
\end{theorem}

The proof of the theorem is presented in Section \ref{sec:nc}.
We conclude that Problem \ref{pbm21} has a unique weak solution.

 In the case of the Signorini condition, we consider the
following variational formulation of the problem  to find $\mathbf{u}$,
$\mathbf{v}\in \mathcal{V}$ such that if $\mathbf{w,w}'\in
\mathcal{V}$ is such that $\mathbf{w}( T) =\mathbf{u}(T) $ and
\begin{equation*}
-g_2\leq w_1(l_{\ast },t) -w_2(l_{\ast },t)
\leq g_1,
\end{equation*}
then,
\begin{align*}
&( \mathbf{v}_0,\mathbf{w}( 0) -\mathbf{u}_0)_H
-\int_0^T( \mathbf{v},\mathbf{v}-\mathbf{w}')_Hdt
\\
&+\int_0^T\Big( \Big( \int_0^{l_{\ast }}k_1u_{1xx}(
u_{1xx}-w_{1xx}) dx\Big) +\Big( \int_{l_{\ast
}}^{1}k_2u_{2xx}( u_{2xx}-w_{2xx}) dx\Big) \Big) ds
\\
&+\int_0^T\Big( \Big( \int_0^{l_{\ast }}d_1v_{1xx}(
u_{1xx}-w_{1xx}) dx\Big) +\Big( \int_{l_{\ast
}}^{1}d_2v_{2xx}( u_{2xx}-w_{2xx}) dx\Big) \Big) ds
\\
&+\int_0^T\frac{1}{3}a_1\Big( ( u_{1x}+\Phi _x)
^3,u_{1x}-w_{1x}\Big) _{H_1}-\overline{\nu }_1(
pu_{1x},u_{1x}-w_{1x}) _{H_1}dt
\\
&+\int_0^T-( \Psi _x,u_1-w_1) _{H_1}
\frac{1}{3}a_2\Big( ( ( u_{2x}) ^3-\nu
_2pu_{2x}) ,u_{2x}-w_{2x}\Big) _{H_2}dt
\\
&\leq \int_0^T(\mathbf{f},\mathbf{u}-\mathbf{w}) _Hds.
\end{align*}


The existence result for the problem with the Signorini condition is:
\begin{theorem}
\label{t32}
Assume that {\rm (A1)--(A5)} hold. Then,  there exists a function
 $\mathbf{u} =(u_1, u_2) \in \mathcal{K}$,
such that $\mathbf{u}(\cdot, 0)=\mathbf{u}_0$,  $\mathbf{u}_t(\cdot, 0)=\mathbf{v}_0$,
 and for a.a. $t\in(0, T)$ the variational equations \eqref{35} and \eqref{36} hold,
together with \eqref{213} and \eqref{214}.
\end{theorem}

The theorem is established in Section \ref{sec:Signorini}, and is based
on obtaining the necessary a priori estimates on the solutions of the
problems with normal compliance and passing to the limit
as the normal compliance stiffnesses become rigid; i.e.,
$\lambda_1, \lambda_2\to \infty$.

We conclude that Problem \ref{pbm22} has a weak solution and note that,
quite typically, the uniqueness of the solution is an unresolved question.

We also note that the corresponding problems without viscosity have variational
solutions, too, see Section \ref{sec:novisc}, however their uniqueness is open.

\section{The problem with normal compliance}
\label{sec:nc}

We establish Theorem \ref{t31} by transforming the variational equation into an abstract
evolution equation in $\mathcal{V}'$ to which we apply various results for such abstract equations.
To deal with the cubic terms in the constitutive relations \eqref{23} we introduce a truncation.
For detailed description of the various Sobolev Spaced used here, we refer to \cite{Ad}.

First, we depict in Figure \ref{fig3} the normal compliance condition \eqref{29} with the Lipschitz
continuous function $\Lambda$, and the contact shear is given by
\begin{equation}
\label{41}
\sigma(t)= \Lambda(u_1(l_{*}, t)-u_2(l_{*}, t)).
\end{equation}

\begin{figure}[ht]
\begin{center}
\begin{picture}(225,92)(30,-12)
\put(30,30){\vector(1,0){225}}
\put(140,-10){\vector(0,1){95}}
\put(245,18){$r$}
\put(142,18){$0$}
\thicklines
\put(130,30){\put(-30,0){\color[rgb]{0,0,1} \line(-1,-5){10}}
\put(50,0){\color[rgb]{0,0,1} \line(1,5){10}}
\put(-30,0){\color[rgb]{0,0,1} \line(1,0){80}}
 \put(-30,0){\line(1,0){80}}
\put(-54, -11){$-g_2$}
\put(44,-11){$g_1$}
\put(55,10){$\lambda_1$}
\put(-30,-25){$\lambda_2$}
\put(-14,42){$\Lambda(r)$}
}
\end{picture}
\end{center}
\caption{The normal compliance graph, the slopes are $\lambda_{i}$,
$i=1,2$, $r=0$ corresponds to the center at $u_2(l_{*}, t)$}
\label{fig3}
\end{figure}

When we study the problem with the Signorini condition, in the following section,
we let $\lambda_1, \lambda_2\to \infty$ leading the
 multifunction $\Lambda_{\infty}=\partial \chi$ given in \eqref{212} and
depicted in Figure \ref{fig4}.

To deal with the nonlinear cubic term,  we use truncation and replace the function
$Q( r) =r^3$ with
\begin{equation}
Q_{m}( r) \equiv
\begin{cases}
m^3 & \text{if }r>m \\
r^3 & \text{if }| r| \leq m \\
-m^3& \text{if }r<-m
\end{cases}   \label{42}
\end{equation}

Integration of \eqref{35} and \eqref{36} in time  over $(0, T)$ leads to
the approximate problem of finding, for  a fixed $m>0$,
a pair $\mathbf{u}_{m}=(u_{1m}, u_{2m})\in \mathcal{V}$ such that
\begin{equation}
\begin{split}
&\langle u_{1tt},w_1\rangle _{\mathcal{V}_1}
 -\int_0^T\sigma(t ) w_1( l_{\ast })\,dt
+( k_1u_{1xx}+d_1u_{1txx}, w_{1xx}) _{\mathcal{H}_1}
\\
&+\frac{1}{3}a_1( Q_{m}( u_{1x}+\Phi_x),w_{1x}) _{\mathcal{H}_1}
-\overline{\nu} _1 (  p  u_{1x},w_{1x}) _{\mathcal{H}_1}
-( \Psi _x,w_1) _{\mathcal{H}_1}
\\
&=( f_1,w_1) _{\mathcal{H}_1},
\end{split}\label{43}
\end{equation}
and
\begin{equation}
\begin{split}
&\langle u_{2tt},w_2\rangle _{\mathcal{V}_2}- \int_0^T\sigma
( t ) w_2( l_{\ast })\,dt
+( k_2u_{2xx}+d_2u_{2txx},w_{2xx}) _{\mathcal{H}_2}
\\
&+  \frac{1}{3}a_2 ( (Q_{m}( u_{2x}) -\overline{\nu}
_2pu_{2x}), w_{2x}) _{\mathcal{H}_2}\\
&=( f_2,w_2)_{\mathcal{H}_2},
\end{split}  \label{44}
\end{equation}
where we used condition \eqref{28} for $\sigma$, along with the initial conditions
$\mathbf{u}(\cdot, 0)=\mathbf{u}_0$,  $\mathbf{u}_t(\cdot, 0)=\mathbf{v}_0$.
 Here, $w_{i}\in \mathcal{V}_{i}$
and the time derivatives are understood in the sense of $V_{i}'$
 valued distributions.


Let $\mathbf{u} =(u_1, u_2) \in \mathcal{V}$, we use the notation
\begin{equation*}
\mathbf{v}=\mathbf{u}',\quad  \mathbf{u}( t) =\mathbf{u}_0+\int_0^t \mathbf{v} (s)\, ds,
\end{equation*}
where here and below the prime denotes the partial $t$ derivative. We use
the notation $\gamma $ for the trace map from $V$ on $x=l_{\ast }$
and the projection $\pi _{i}\mathbf{u}\equiv u_{i}$, for $i=1,2$, and $\gamma^{*}$
is the adjoint map of $\gamma$.  Next, we define the
operators $A,B,N_{m}( t,\cdot ) :V\to V'$ by
\begin{gather}
\langle A\mathbf{v},\mathbf{w}\rangle \equiv (
\int_0^{l_{\ast }}d_1v_{1xx}w_{1xx}dx) +( \int_{l_{\ast
}}^{1}d_2v_{2xx}w_{2xx}dx) ,  \label{45}
\\
\langle B\mathbf{u},\mathbf{w}\rangle \equiv (
\int_0^{l_{\ast }}k_1u_{1xx}w_{1xx}dx) +( \int_{l_{\ast
}}^{1}k_2u_{2xx}w_{2xx}dx) ,  \label{46}
\end{gather}
and
\begin{equation}
\begin{split}
\langle N_{m}( t,\mathbf{u}) ,\mathbf{w}\rangle
&=\frac{1}{3}a_1( Q_{m}( u_{1x}+\Phi _x) ,w_{1x}) _{H_1}-
\overline{\nu}_1( pu_{1x},w_{1x}) _{H_1}
\\
&\quad -( \Psi _x,w_1) _{H_1}+\frac{1}{3}a_2( (
Q_{m}( u_{2x}) -\overline{\nu }_2pu_{2x}) ,w_{2x})
_{H_2}.
\end{split} \label{47}
\end{equation}
Here and below, $\langle \cdot ,\cdot \rangle $ denotes the duality pairing
between $V$ and $V'$ and we consider these operators as also being
defined on $\mathcal{V}$ in the obvious way.


Then, \eqref{35} and \eqref{36}, along with the initial conditions, can be written
as an abstract equation in $\mathcal{V}'$ as follows:

\begin{problem} \label{pbm41} \rm
Find a pair $\mathbf{u}_{m}, \mathbf{v}_{m} \in \mathcal{V}$ with $\mathbf{v}'_{m} \in \mathcal{V}'$,
for integer $m>0$, such that,
\begin{gather}
\begin{aligned}
&\mathbf{v} '+A\mathbf{v} +B\mathbf{u}+N_{m}( \cdot ,\mathbf{u})
+\pi _1^{\ast }\gamma ^{\ast }\Lambda(u_1(l_{*}, \cdot)-u_2(l_{*}, \cdot))\\
& -\pi _2^{\ast}\gamma ^{\ast }\Lambda(u_1(l_{*}, \cdot)-u_2(l_{*}, \cdot)) =\mathbf{f},
\end{aligned}\label{48} \\
\mathbf{v} ( 0) =\mathbf{v} _0,  \label{49} \\
\mathbf{u}( t) = \mathbf{u}_0+\int_0^t \mathbf{v} (s) ds.  \label{410}
\end{gather}
\end{problem}

Here, and below, we omit the subscript $m$ from the solution, until we need it in
passing to the limit $m\to \infty$.

We note that the operator $N_{m}$ is Lipschitz continuous, and satisfies
\begin{equation*}
\| N_{m}( t,\mathbf{u}) -N_{m}( t,\mathbf{w} )
\| _{V'}\leq K_{m}\| \mathbf{u}-\mathbf{w} \| _{V},
\end{equation*}
where the constant $K_{m}$ which is independent of $t$.

Also, let $N( t,\mathbf{u}) :V\to V'$ be
defined by
\begin{equation}
\begin{split}
\langle N( \cdot ,\mathbf{u}) ,\mathbf{w}\rangle
&= \frac{1}{3}a_1( ( u_{1x}+\Phi _x) ^3,w_{1x})
_{H_1}-\overline{\nu }_1( pu_{1x},w_{1x}) _{H_1}
\\
&\quad -( \Psi _x,w_1) _{H_1}+\frac{1}{3}a_2( ( (
u_{2x}) ^3-\nu _2pu_{2x}) ,w_{2x}) _{H_2}.
\end{split} \label{411}
\end{equation}
Next, let
\begin{equation*}
R_{m}( r) =\int_0^{r}Q_{m}( z) dz,
\end{equation*}
Moreover, using assumption (A2), we obtain
\begin{align*}
\int_0^t p( s) ( u_{1x},u_{1tx}) _{H_1}\,ds
&=\frac{ 1}{2}\int_0^t p( s) \frac{d}{ds}| u_{1x}|_{H_1}^2
\\
&\geq -p_0| u_{1x}( t) | _{H_1}^2-C(
p( 0) ,u_{01}) -p_0\int_0^t | u_{1x}(
s) | _{H_1}^2ds.
\end{align*}
This implies that for suitable $\delta >0$, depending only on $a_1$ and
$a_2$, we can derive the estimate
\begin{align*}
&\int_0^t \langle N_{m}( s,\mathbf{u}( s) ),
\mathbf{v}( s) \rangle ds \\
&\geq \delta \int_0^{l_{\ast }}R_{m}( u_{1x}( t) +\Phi
_x( t) ) dx-\frac{a_1}{3}\int_0^t ( Q_{m}(
u_{1x}+\Phi _x) ,\Phi _x') _{H_1}dt
 \\
&\quad +\delta \int_{l_{\ast }}^{1}R_{m}( u_{2x}( t) )
dx-C( u_{10},\Phi _x( 0,\cdot ) ,\Phi _x) -(
2+p_0) \int_0^t \int_0^{l_{\ast }}| u_{1x}|
^2\,dx\,ds
\\
&\quad -C( \Psi _x,\Psi _{tx},\Psi _x( 0,\cdot ) )
-C( p) | u_{1x}( t) | _{H_1}^2.
\end{align*}
Consider that second term on the right of the inequality.
\begin{equation*}
\big| \int_0^t ( Q_{m}( u_{1x}+\Phi _x) ,\Phi _x') _Hdt\big|
\leq C( \Phi _x') \int_0^t \int_0^{l_{\ast }}Q_{m}( | u_{1x}+\Phi
_x| ) \,dx\,ds+C( \Phi _x',\Phi _x).
\end{equation*}
Then, adjusting the constants, we find
\begin{equation*}
\leq C( \Phi _x') \int_0^t \int_0^{l_{\ast
}}R_{m}( u_{1x}+\Phi _x) \,dx\,ds+C( \Phi _x',\Phi_x).
\end{equation*}
It follows that
\begin{align*}
&\int_0^t \langle N_{m}( s,\mathbf{u}( s) ),
\mathbf{v}( s) \rangle ds\\
&\geq \delta \int_0^{l_{\ast }}R_{m}( u_{1x}( t) +\Phi
_x( t) ) dx-C( \Phi _x')
\int_0^t \int_0^{l_{\ast }}R_{m}( u_{1x}+\Phi _x) \,dx\,ds
\\
&\quad -C( \Phi _x',\Phi _x) +\delta \int_{l_{\ast
}}^{1}R_{m}( u_{2x}( t) ) dx-C( u_{10},\Phi
_x( 0,\cdot ) ,\Phi _x)
\\
&\quad
-( 2+p_0) \int_0^t \int_0^{l_{\ast }}|
u_{1x}| ^2\,dx\,ds
-C( \Psi _x,\Psi _{tx},\Psi _x( 0,\cdot ) )
-C( p) | u_{1x}( t) | _{H_1}^2.
\end{align*}
Written in a simpler form, we have,
\begin{equation}
\begin{split}
&\int_0^t \langle N_{m}( s,\mathbf{u}( s) ),
\mathbf{v}( s) \rangle ds\\
&\geq \delta \int_0^{l_{\ast }}R_{m}( u_{1x}( t) +\Phi_x( t) ) dx
+\delta \int_{l_{\ast }}^{1}R_{m}(u_{2x}( t) ) dx\\
&\quad -C( p) \int_0^t \| \mathbf{u}\|_{W}^2ds
 -C( p) \| \mathbf{u}( t) \|_{W}^2
-C\int_0^t \int_0^{l_{\ast }}R_{m}( u_{1x}+\Phi
_x) \,dx\,ds -C,
\end{split} \label{412}
\end{equation}
where the constants depend on the initial data and the given functions, but
not on $a_{i},d_{i}$ or $k_{i}$, $i=1,2$.

The following existence and uniqueness result holds for each integer $m$.

\begin{lemma}\label{lemma42}
For each positive integer $m$,  there exists a unique solution $(\mathbf{u}_m, \mathbf{v}_m)$
to Problem   \ref{pbm41}.
\end{lemma}

\begin{proof} For the sake of convenience,
 let the operator $M(\cdot ,\mathbf{u})$ be defined as
\begin{equation*}
B\mathbf{u}+N_{m}( \cdot ,\mathbf{u}) +\pi _1^{\ast }\gamma
^{\ast }\Lambda ( u_1( l_{\ast }) -u_2( l_{\ast
}) ) -\pi _2^{\ast }\gamma ^{\ast }\Lambda ( u_1(
l_{\ast }) -u_2( l_{\ast }) ),
\end{equation*}
so that $M$ is Lipschitz as a map from $\mathcal{V}$ to $\mathcal{V}^{\prime}$
with a constant $K$ depending on $m$ and $\lambda _1,\lambda _2$.

For $\mathbf{u_1}\in \mathcal{V}$, let $\mathbf{v}_1$ be the solution of
\begin{equation*}
\mathbf{v}_1'+A\mathbf{v}_1+M( t,\mathbf{u}_1) =
\mathbf{f},\quad  {\mathbf{v}}_1( 0) =\mathbf{v}_0.
\end{equation*}
Then, let $\mathbf{w}_1$ be given by
\begin{equation*}
\mathbf{w}_1( t) =\mathbf{u}_0+\int_0^t \mathbf{v}
_1( s) ds.
\end{equation*}
Consider the map $\mathbf{u}_1\to \mathbf{v}_1\to
\mathbf{w}_1.$ Consider $\mathbf{u}_1,\mathbf{u}_2$ and the
corresponding  $\mathbf{w}_{i}$ that result in this way. First,
from the equation, we have
\begin{align*}
&\frac{1}{2}| \mathbf{v}_1( t) -\mathbf{v}_2(
t) | _H^2+\delta \int_0^t \| \mathbf{v}
_1( s) -\mathbf{v}_2( s) \| _{V}^2ds
\\
& \leq \int_0^t \langle M( s,\mathbf{u}_1) -M( s,
\mathbf{u}_1) ,\mathbf{v}_1( s) -\mathbf{v}_2(
s) \rangle ds
\\
&\leq C\int_0^t \| \mathbf{u}_1( s) -\mathbf{u}
_2( s) \| _{V}^2ds+\frac{\delta }{2}
\int_0^t \| \mathbf{v}_1( s) -\mathbf{v}_2(s) \| _{V}^2ds,
\end{align*}
and so
\begin{equation*}
| \mathbf{v}_1( t) -\mathbf{v}_2( t)| _H^2+\int_0^t \| \mathbf{v}_1( s)
-\mathbf{v}_2( s) \| _{V}^2ds
\leq C\int_0^t \| \mathbf{u}_1( s) -\mathbf{u}_2(s) \| _{V}^2ds.
\end{equation*}
Now, it follows  from the definition of $\mathbf{w}_{i}$ and the
above inequality, that
\begin{align*}
\| \mathbf{w}_1( t) -\mathbf{w}_2( t)
\| _{V}^2 &\leq C( T) \int_0^t \| \mathbf{v
}_1( s) -\mathbf{v}_2( s) \| _{V}^2ds \\
&\leq C\int_0^t \| \mathbf{u}_1( s) -\mathbf{u}
_2( s) \| _{V}^2ds.
\end{align*}
Iterating this inequality sufficient number of times shows that a high
enough power of this map is a contraction map on $\mathcal{V}$.
Therefore, it has a unique fixed point
$\mathbf{u}$ which yields the unique solution to Problem \ref{pbm41}.
\end{proof}


Next, we derive the necessary estimates that are independent of $m$.
We denote by $(\mathbf{u}_m, \mathbf{v}_m)$ the solution guaranteed in
Lemma \ref{lemma42},  however, below we omit the subscript $m$.

From \eqref{48} and the estimate
given for $N_{m}$ in \eqref{412}, it follows that there exists a constant
$\delta >0$, depending on the parameters of the problem, such that
\begin{align*}
&\frac{1}{2}| \mathbf{v}( t) | _H^2+\min
( d_1,d_2) \int_0^t \| \mathbf{v}( s)
\| _{V}^2ds+\min ( k_1,k_2) \| \mathbf{u}
( t) \| _{V}^2
\\
&+\delta \int_0^{l_{\ast }}R_{m}( u_{1x}( t) +\Phi
_x( t) ) dx+\delta \int_{l_{\ast }}^{1}R_{m}(
u_{2x}( t) ) dx
\\
&+\int_0^t \Lambda ( u_1( s,l_{\ast }) -u_2(
s,l_{\ast }) ) ( v_1( s,l_{\ast }) ) ds
\\
&-\int_0^t \Lambda ( u_1( s,l_{\ast }) -u_2(
s,l_{\ast }) ) ( v_2( s,l_{\ast }) ) ds
\\
&\leq \frac{1}{2}| \mathbf{v}_0| _H^2+
C( p') \int_0^t \| \mathbf{u}
\| _{W}^2ds+C( p) \| \mathbf{u}( t)\| _{W}^2
\\
&\quad +\int_0^t \int_0^{l_{\ast }}R_{m}( u_{1x}+\Phi _x) \,dx\,ds+ C.
\end{align*}
The two terms involving $\Lambda $ taken together are nonnegative. In fact,
they equal
\[
S( u_1( t, l_{\ast }) -u_2( t,l_{\ast }) ),
\]
where
\[
S( r) =\int_0^{r}\Lambda ( s)ds.
\]
Then, by using Gronwall's inequality to eliminate the term
\[
\int_0^t \int_0^{l_{
\ast }}R_{m}( u_{1x}+\Phi _x) \,dx\,ds,
\]
on the right-hand side, we can
adjust the constants and obtain
\begin{align*}
&| \mathbf{v}( t) | _H^2+\min (
d_1,d_2) \int_0^t \| \mathbf{v}( s)
\| _{V}^2ds+\min ( k_1,k_2) \| \mathbf{u}
( t) \| _{V}^2
\\
&+\delta \int_0^{l_{\ast }}R_{m}( u_{1x}( t) +\Phi
_x( t) ) dx+\delta \int_{l_{\ast }}^{1}R_{m}(u_{2x}( t) ) dx
+S( u_1( t, l_{\ast }) -u_2( t,l_{\ast }))
\\
&\leq C\Big( | \mathbf{v}_0|
_H^2+\int_0^t \| \mathbf{u}\| _{W}^2ds+\|
\mathbf{u}( t) \| _{W}^2+1\Big).
\end{align*}
The two terms involving $R_{m}$ are nonnegative, since $R_{m}$ is
nonnegative. Therefore, discarding these two terms yields
\begin{align*}
&| \mathbf{v}( t) | _H^2+\min (d_1,d_2) \int_0^t \| \mathbf{v}( s)
\| _{V}^2ds\\
&+\min ( k_1,k_2) \| \mathbf{u}( t) \| _{V}^2
+S( u_1( t,l_{\ast }) -u_2( t,l_{\ast })) \\
&\leq C\Big( | \mathbf{v}_0|
_H^2+\int_0^t \| \mathbf{u}\| _{W}^2ds+\|
\mathbf{u}( t) \| _{W}^2+1\Big).
\end{align*}
The compactness of the embedding of $V$ into $W$ implies
\begin{align*}
&| \mathbf{v}( t) | _H^2+\min (
d_1,d_2) \int_0^t \| \mathbf{v}( s)\| _{V}^2ds\\
&+\min ( k_1,k_2) \| \mathbf{u}( t) \| _{V}^2
 +S( u_1( t,l_{\ast }) -u_2( t,l_{\ast }))
\\
&\leq C\Big( | \mathbf{v}_0|
_H^2+\int_0^t \| \mathbf{u}\| _{W}^2ds+1\Big)
+\frac{1}{2}\min ( k_1,k_2) \| \mathbf{u}(
t) \| _{V}^2+C| \mathbf{u}( t)| _H^2.
\end{align*}
Also,
\[
\ | \mathbf{u}( t) | _H^2\leq C\Big(
| \mathbf{u}_0| ^2+\int_0^t  |\mathbf{v}( s) | ^2ds\Big),
\]
and so an application of the Gronwall inequality yields
\begin{equation}
\begin{aligned}
&| \mathbf{v}( t) | _H^2+\min (d_1,d_2) \int_0^t \| \mathbf{v}( s)
\| _{V}^2ds\\
&+S( u_1( t,l_{\ast }) -u_2(t,l_{\ast }) )
+\min ( k_1,k_2) \| \mathbf{u}( t) \| _{V}^2\\
&\leq C( \Phi ,u_0,v_0,T,p_0,k_1,k_2) .
\end{aligned} \label{413}
\end{equation}

 It follows from the continuity of the embedding of $H^{1}$ into
 $C([ 0,1])$ that this estimate provides an upper bound on
the values of $|u_{ix}( x,t) |$ that is independent of $m$.
Therefore, by taking $m$ sufficiently large, we find that the values of
$u_{ix}$ and $u_{1x}( t) +\Phi _x( t) $ remain in the
unmodified region in the definition of $Q_{m}$ in \eqref{42}. This
observation completes the proof of the following abstract version
of Theorem \ref{t31}.

\begin{theorem} \label{t43}
Let {\rm (A1)--(A5)} hold. Then, for each pair $(\lambda_1, \lambda_2)$
there exists a unique solution $(\mathbf{u}, \mathbf{v})$ to
\begin{gather}
\begin{aligned}
&\mathbf{v} '+A\mathbf{v} +B\mathbf{u}+N( \cdot ,\mathbf{u} )
+\pi _1^{\ast }\gamma ^{\ast }\Lambda( u_1(l_{\ast }) -u_2( l_{\ast }) ) \\
&-\pi _2^{\ast }\gamma ^{\ast }\Lambda ( u_1( l_{\ast })
-u_2( l_{\ast }) ) =\mathbf{f},
\end{aligned}\label{414}
\\
\mathbf{v} ( 0) =\mathbf{v} _0,\label{415}
\\
\mathbf{u}( t) = \mathbf{u}_0+\int_0^t \mathbf{v} (s) ds.  \label{416}
\end{gather}
This solution satisfies the estimate
\begin{equation}
\begin{aligned}
&| \mathbf{v} ( t) | _H^2+\min ( d_1,d_2)
\int_0^t \| \mathbf{v} ( s) \| _{V}^2ds
+\| \mathbf{u}( t) \| _{V}^2\\
&+\int_0^{l_{\ast}}u_{1x}^{4}( t) \,dx
+\int_{l_{\ast }}^{1}u_{2x}^{4}( t)\,dx
+S( u_1( t,l_{\ast }) -u_2( t,l_{\ast}) )
\\
&\leq C( p_0, u_0, v_0, C_0,k_1,k_2,T),
\end{aligned}\label{417}
\end{equation}
where the constant on the right depends on the indicated quantities but is
independent of $\lambda_1, \lambda_2, d_1$, and $ d_2$.
\end{theorem}

The estimate follows directly from \eqref{413} and the continuity of the
embedding of $H^{1}( I) $ into $C( I) $ when $I$ is an
interval. The continuity of this embedding implies that the solution to the
truncated problem is such that $u_{1x}$ and $u_{2x}$ remain in the
region where the truncation is inactive, thus yielding the unique solution
in the theorem. This proves Theorem \ref{t31} as well.

\section{The problem with the Signorini condition}
\label{sec:Signorini}

We turn to the idealized problem with perfectly rigid stops, and use the results of
the previous section. Indeed, we obtain a weak solution for the problem with the Signorini
contact condition as a limit of the solutions of the problem with normal compliance
in the limit when the normal compliance stiffness coefficients tend to infinity.
Therefore, in this section we replace the coefficients with
\begin{equation}
\label{51}
\lambda_1=\lambda_2= n,
\end{equation}
and obtain the necessary a priori estimates to pass to the limit as $n\to \infty$.

\begin{figure}[ht]
\begin{center}
\begin{picture}(225,95)(30,-10)
\put(30,30){\vector(1,0){225}}
\put(140,-10){\vector(0,1){95}}
\put(245,18){$r$}
\put(142,18){$0$}
\thicklines
\put(130,30){\put(-30,-30)
{\color[rgb]{0,0,1}  \line(0,1){30}}\put(-30,0) {\line(1,0){80}}
\put(50,0){\color[rgb]{0,0,1}  \line(0,1){30}}
\put(-30,0){\color[rgb]{0,0,1}  \line(1,0){80}}
%\put(0,0){\qbezier(0,-3)(0,0)(0,3)}
\put(-54, -11){$-g_2$}
\put(44,-11){$g_1$}
\put(-20,40){$\partial \chi(r)$}
  }
\end{picture}
\end{center}
\caption{The Signorini graph $\partial  \chi(r)$, $r=0$ corresponds
 to the center at $u_2(l_{*}, t)$}
\label{fig4}
\end{figure}

We have that the graph $\partial  \chi(r)$, depicted in Figure \ref{fig4}, can be
obtained from the
graph of $\Lambda$ in Figure \ref{fig3} in the limit $n \to \infty$.

For the sake of convenience we will assume in this section that
\begin{equation*}
u_{0ixx}\in L^{\infty }( I_{i}),\quad i=1,2.
\end{equation*}

Now, we recall the Signorini condition \eqref{213}-\eqref{214},
\begin{gather*}
-g_2 \leq u_1-u_2 \leq  g_1, \\
\sigma(t)=\sigma_1(l_*, t)=-\sigma_2(l_*, t) \in \partial \chi (u_1( l_{\ast }, t)
-u_2( l_{\ast }, t) ).
\end{gather*}
We note that assumption (A1) implies that
$\chi( u_{10}( l_{\ast }) -u_{20}( l_{\ast}) ) =0$.

We now add a subscript $n$ to the solution to the problem with the
normal compliance condition obtained in Section \ref{sec:nc} that corresponds
to the stiffnesses \eqref{51},
so we denote the solutions of \eqref{414}--\eqref{416} as $(\mathbf{u}_{n}, \mathbf{v} _{n})$.

We use below the following important theorem found in \cite{Sim87}.

\begin{theorem} \label{t51}
 Let $q>1$ and let $E\subseteq W\subseteq X$, where the
injection map from $W$ to $X$ is continuous and is compact from
$E$ to $W$. Let $S_R$ be defined by
\begin{equation*}
S_R=\big\{ u: \| u( t)
\| _{E}+\| u'\| _{L^{q}( [a,b] ;X) }\leq R, \;t\in [a,b] \big\}.
\end{equation*}
Then, $S_R\subseteq $ $C( [a,b] ;W) $ and if $\{
u_{n}\} \subseteq  S_R$, there exists a subsequence, $\{
u_{n_{k}}\} $ that converges uniformly to a function $u\in C(
[a,b] ;W) $.
\end{theorem}

It follows from the theorem and  estimate \eqref{417} that there exists a subsequence,
still denoted as $\mathbf{u}_{n}$, such that as $n\to \infty$ the
following convergences are obtained.
\begin{equation}
u_{inx}\to u_{ix}\quad \text{strongly in }C( [0,T];C( I_{i}) ),  \label{52}
\end{equation}
where here and below $I_1=[0,l_{\ast }] ,I_2=[l_{\ast },1]$, and $i=1,2$.
\begin{gather}
u_{ni}\to u_{i}\quad \text{strongly in }C( [0,T];C( I_{i}) ),  \label{53}
\\
\mathbf{u}_{n}\to \mathbf{u}\quad \text{weak}\ast \text{ in }L^{\infty}( 0,T;V),  \label{54}
\\
\mathbf{v} _{n}\to \mathbf{v}\quad\text{weakly in }\mathcal{V},\label{55}
\\
\mathbf{v} _{n}\to \mathbf{v} \quad \text{weak}\ast \text{ in }L^{\infty}( 0,T;H).  \label{56}
\end{gather}

Next, let $\Lambda_n$ be the function \eqref{29} with an obvious subscript,
then it follows from the bound in \eqref{417} that
\[
S_n( u_{1n}( l_{\ast }, t) -u_{2n}( l_{\ast}, t) ) \leq C,
\]
where $C$ is independent of $n$, and so in the limit it follows from
\eqref{53} that
\begin{equation*}
-g_2\leq ( u_1(l_{\ast },t) -u_2(l_{\ast},t) ) \leq g_1.
\end{equation*}%
Therefore, $\mathbf{u}\in \mathcal{K}$ and so it satisfies the necessary constraint.

It follows from the equations  estimate \eqref{52},
that $v_{i}'$ is bounded in
\[
L^2\Big( 0,T;( H_0^2( I_{i}) ) '\Big).
\]
More precisely, we mean $i^{\ast }v_{in}'$ is bounded, where
$i:H_0^2 \to V_0$ is the inclusion map. Therefore, we can also assume
that
\begin{equation}
v_{ni}'\to v'\quad \text{in }L^2( 0,T;(H_0^2( I_{i}) ) ').  \label{57}
\end{equation}

To continue, we need the following fundamental theorem in \cite{Lio69}.

\begin{theorem}\label{t52}
Let $E\subseteq W\subseteq X$, where the injection map
$i:W \to X$ is continuous and is compact from $E$ to $W$.
Let $p\geq 1$, let  $ q>1$, and define
\begin{align*}
S\equiv \big\{&u\in L^{p}( [a,b] ;E) :u'\in L^{q}( [a,b] ;X) \text{ and}\\
 &\| u\| _{L^{p}([a,b] ;E) }+\| u'\| _{L^{q}( [a,b] ;X) }\leq R\big\}.
\end{align*}
Then, $S$ is precompact in $L^{p}( [a,b] ;W) $, hence,
 if $\{ u_{n}\} _{n=1}^{\infty }\subseteq S$, it has a
subsequence $\{ u_{n_{k}}\} $ which converges in $L^{p}(
[a,b] ;W).$
\end{theorem}

We apply this theorem to the case where $W=H_{i}$ and $E=V_{i}$ and find that
addition to the above convergences, we also have
\begin{equation}
v_{ni}\to v_{i}\quad \text{strongly in }\mathcal{H}_{i},
\label{58}
\end{equation}
for a suitable subsequence.
We let
\[
V_0=\{\mathbf{u} =(u_1, u_2) \in V: u_{i}\in H_0^2( I_{i}) \; i=1,2\},
\]
and let $\mathcal{V}_0=L^2((0, T); V_0)$.
Then, it follows  from the definition of the operators in \eqref{414}
and the estimates above,  that the expression
\begin{equation}
u_{n1tt}+d_1v_{n1xxxx}+k_1u_{n1xxxx}
-\frac{1}{3}a_1( ( u_{n1x}+\Phi_x) ^3) _x
+\overline{\nu} _1p u_{n1xx}+\Psi _{xx}  \label{59}
\end{equation}
is bounded in $\mathcal{H}_1$, independently of $n$. Here, we use the
appropriate measurable representatives so that the result is product measurable.
It is obtained by applying the expression to $( \varphi ,0)$, where
$\varphi \in C_0^{\infty }( [0,T] \times [0,l_{\ast }] ) $,
and using  the density of the functions $\varphi $
in $\mathcal{H}_1$.  Similarly, we find that the expression
\begin{equation}
u_{n2tt}+d_2v_{n2xxxx}+k_2u_{n2xxxx}-  \frac{1}{3}a_2((
u_{n2x}) ^3-\nu _2pu_{n2x}) _x  \label{510}
\end{equation}
is bounded in $\mathcal{H}_2$, independently of $n$.

Because of the
strong convergence in  \eqref{52} - \eqref{56}, there exists a further subsequence,
still indexed by $n$, such that, in addition to \eqref{52} - \eqref{56}, it also follows that
the expression in \eqref{59} converges weakly in $\mathcal{H}_1$ to
\begin{equation}
m(x, t) \equiv u_{1tt}+d_1v_{1xxxx}+k_1u_{1xxxx}
-\frac{1}{3} a_1( ( u_{1x}+\Phi_x
) ^3) _x+\overline{\nu} _1pu_{1xx}+\Psi _{xx},  \label{511}
\end{equation}
and the expression in \eqref{510} converges weakly in $\mathcal{H}_1$ to
\begin{equation}
u_{2tt}+d_2v_{2xxxx}+k_2u_{2xxxx}-( \frac{a_2}{3}(
u_{2x}) ^3-\overline{\nu} _2pu_{2x}) _x.  \label{512}
\end{equation}
Now, let $\psi \in W_0^{2,\infty }( I_1) ,\varphi \in
W_0^{1,\infty }( 0,T)$ and consider
\begin{align}
&\int_{I_1}\int_0^Tm(x, t) u_1( x, t) \psi ( x) \varphi ( t) \,dt\,dx  \label{513}
\\
&= \int_{I_1}\int_0^T\Big( {-u_{1t}^2}
+d_1v_{1xx}u_{1xx}+k_1u_{1xx}^2 \nonumber \\
&\quad + \frac{a_1}{3}( u_{1x}+\Phi_x) ^3u_{1x}\Big) \psi(x) \varphi(t)s \,dt\,dx
 \label{514}
\\
&\quad +\int_{I_1}\int_0^T( -\overline{\nu} _1pu_{1x}^2-\Psi _xu_{1x})
\psi ( x) \varphi ( t) \,dt\,dx  \label{515}
\\
&\quad +\int_{I_1}\int_0^T( -u_{1t}\varphi '( t) \psi
( x) u_1+d_1v_{1xx}\psi _{xx}(x) \varphi ( t) u_1) \,dt\,dx
 \\
&\quad + \int_{I_1}\int_0^T \Big(k_1u_{1xx}\psi _{xx}(x)\phi(t) u_1+\frac{1}{3}a_1
( u_{1x}+\Phi_x ) ^3\psi _x(x)\phi(t) u_1\Big)\,dt\,dx
 \\
&\quad -\int_{I_1}\int_0^T(\overline{ \nu} _1pu_{1x}\psi _x(x)\phi(t) u_1
+\Phi _x\psi _x(x)\phi(t) u_1) \,dt\,dx.
 \label{516}
\end{align}

We proceed in a similar way and obtain the analogous expression for the
solutions that depend on $n$, so we replace $u_1$ with $u_{n1}$  and $v_1$
with $v_{n1}$  in \eqref{513}--\eqref{516}. We also replace $m( x, t) $
with  $m_{n}(x, t) $. We obtain,
\begin{align}
&\int_{I_1}\int_0^Tm_{n}( x, t) u_{n1}( x, t) \psi
( x) \varphi ( t) \,dt\,dx  \label{517}
\\
&=
\int_{I_1}\int_0^T\Big( {-u_{n1t}^2}
+d_1v_{n1xx}u_{n1xx}+k_1u_{n1xx}^2 \nonumber \\
&\quad + \frac{1}{3} a_1 ( u_{n1x}+\Phi_x
) ^3u_{n1x}\Big) \psi(x) \varphi(t) \,dt\,dx
\label{518}
\\
&\quad +\int_{I_1}\int_0^T( -\overline{\nu} _1pu_{n1x}^2-\Psi
_xu_{n1x}) \psi ( x) \varphi ( t) \,dt\,dx
\label{519}
\\
&\quad +\int_{I_1}\int_0^T( -u_{n1t}\varphi '( t) \psi
( x) u_{n1}+d_1v_{n1xx}\psi _{xx}(x)\varphi(t)
u_{n1}) \,dt\,dx  \nonumber \\
&\quad+
\int_{I_1}\int_0^T (k_1u_{n1xx}\psi _{xx}(x)\varphi(t) u_{n1}+\frac{1}{3}a_1
( u_{n1x}+\Phi_x ) ^3\psi _x(x)\varphi(t)
u_{n1}) \,dt\,dx
\nonumber \\
&\quad -\int_{I_1}\int_0^T( \overline{\nu} _1pu_{n1x}\psi _x(x)\varphi(t) u_{n1}+\Psi
_x\psi _x(x)\varphi(t) u_{n1}) \,dt\,dx.  \label{520}
\end{align}

Then, the convergences established above imply that $m_{n}$ converges weakly to $m$
in $\mathcal{H}_1$ as $n\to \infty$.  Indeed,
 \eqref{517} converges to  \eqref{513} and the part  between
\eqref{519} - \eqref{520}
converges to \eqref{515} - \eqref{516}, and the part between
 \eqref{518} - \eqref{519} converges
to \eqref{514} - \eqref{515}. This is summarized as the following lemma.

\begin{lemma}\label{lemma53}
Let $\psi \in W_0^{2,\infty }( I_1)
,\phi \in W_0^{1,\infty }( 0,T)$, then
\begin{align*}
&\lim_{n\to \infty }\int_{I_1}\int_0^T\Big(
-u_{n1t}^2 +d_1v_{n1xx}u_{n1xx}+k_1u_{n1xx}^2
\\
&+ \frac{1}{3} a_1( u_{n1x}+\Phi_x) ^3u_{n1x}\Big) \psi \phi \,dt\,dx
\\
&+\int_{I_1}\int_0^T( -\overline{\nu} _1pu_{n1x}^2-\Psi
_xu_{n1x}) \psi \phi \,dt\,dx
\\
&=\int_{I_1}\int_0^T\Big( {-u_{1t}^2}
+d_1v_{1xx}u_{1xx}+k_1u_{1xx}^2
+ \frac{1}{3}a_1 ( u_{1x}+\Phi_x
) ^3u_{1x}\Big) \psi \phi \,dt\,dx
\\
&\quad +\int_{I_1}\int_0^T( -\overline{\nu} _1pu_{1x}^2-\Psi _xu_{1x})
\psi \phi \,dt\,dx.
\end{align*}
\end{lemma}


We claim that the above estimate holds without the product $\psi \phi $.
To show this, we let $\varphi_{\delta }$ be a piecewise linear function
such that $\varphi _{\delta}(0)=\varphi _{\delta }(T)=0$, and
$\varphi _{\delta }(t)=1$ for $t\in [\delta ,T-\delta ]$, and
let  $\psi _{\delta }\in W_0^{2,\infty }( I_1) $
be a  function such that $\psi _{\delta }(0)=\psi
_{\delta }(1)=0$, and $\psi _{\delta }(x)=1$ for $x\in [\delta
^2,l_{\ast }-\delta ^2] $, with both $\psi _{\delta }$ and $\phi
_{\delta }$ having values in $[0,1] $. Then, we have the
following lemma.

\begin{lemma}\label{lemma54}
The following estimate holds,
\begin{align*}
&\limsup_{n\to \infty }\int_{I_1}\int_0^T\Big( {
u_{n1t}^2}-d_1v_{n1xx}u_{n1xx}-k_1u_{n1xx}^2
- \frac{1}{3}a_1 ( u_{n1x}+\Phi_x) ^3u_{n1x}\Big) \,dt\,dx
\\
&+\int_{I_1}\int_0^T( \overline{\nu} _1pu_{n1x}^2+\Psi _xu_{n1x})
\,dt\,dx
\\
&\leq \int_{I_1}\int_0^T\Big( {u_{1t}^2}
-d_1v_{1xx}u_{1xx}-k_1u_{1xx}^2 -
 \frac{1}{3}a_1 ( u_{1x}+\Phi_x) ^3u_{1x}\Big) \,dt\,dx
\\
&\quad +\int_{I_1}\int_0^T( \overline{\nu} _1pu_{1x}^2+\Psi _xu_{1x}) \,dt\,dx.
\end{align*}
\end{lemma}

\begin{proof}
Let $\phi _{\delta }$ and $\psi _{\delta }$ be the functions described
above. We denote by $Q_{\delta} = I_1\setminus [ \psi _{\delta }=1]$, and assume that
\[
\operatorname{meas}( Q_{\delta }) <\delta,
\]
and that $\delta $ is small enough so that  if $\eta >0$ is given then
the following estimate is  valid:
\begin{equation}
\frac{d_1}{2\delta }\int_{I_1}u_{01xx}^2( 1-\psi _{\delta
}) dx<\eta.  \label{521}
\end{equation}
This is easily obtained because of the assumption that $u_{01xx}\in
L^{\infty }$ and by the construction, $1-\psi _{\delta }=0$ off a set of measure
$2\delta ^2$. Also, let $\delta $ be small enough so that
\begin{equation}
\int_{I_1}\int_0^Tv_1^2( 1-\psi _{\delta }\phi _{\delta
}) \,dt\,dx<\eta.  \label{522}
\end{equation}
Next, consider the following integrals
\begin{gather}
J_{11}\equiv \int_{I_1}\int_0^Tv_{n1}^2( 1-\psi _{\delta }\phi _{\delta
}) \,dt\,dx,  \label{523}
\\
J_{12}\equiv \int_{I_1}\int_0^T-d_1v_{n1xx}u_{n1xx}( 1-\psi _{\delta }\phi
_{\delta }) \,dt\,dx,  \label{524}
\\
J_{13}\equiv\int_{I_1}\int_0^T-k_1u_{n1xx}^2( 1-\psi _{\delta }\phi
_{\delta }) \,dt\,dx.  \label{525}
\end{gather}
First, we note that it follows from \eqref{522}   and
 the strong convergence result in \eqref{58}, that if $n$ is
sufficiently large, then
\begin{equation*}
\int_{I_1}\int_0^Tv_{n1}^2( 1-\psi _{\delta }\phi _{\delta
}) \,dt\,dx<\eta.
\end{equation*}
Below, we only consider such $n$. Next, we consider  \eqref{525} and find
\begin{equation*}
\limsup_{n\to \infty
}\int_{I_1}\int_0^T-k_1u_{n1xx}^2( 1-\psi _{\delta }\phi
_{\delta }) \,dt\,dx<\eta
\end{equation*}
It remains to consider \eqref{524}, which is of the form
\begin{equation*}
-\frac{1}{2}d_1 \int_{I_1}\int_0^T\frac{d}{dt}(
u_{n1xx}^2) ( 1-\psi _{\delta }\phi _{\delta }) \,dt\,dx.
\end{equation*}
We integrate this by parts and obtain
\begin{align*}
&-\frac{1}{2}d_1\int_{I_1}\Big( u_{n1xx}^2( 1-\psi _{\delta }\phi
_{\delta }) |_0^T-\int_0^Tu_{n1xx}^2( -\psi _{\delta
}\phi _{\delta }') dt\Big) dx
\\
&=-\frac{1}{2}d_1 \Big(\int_{I_1}u_{n1xx}^2( T) dx-
\int_{I_1}u_{01xx}^2 + \int_{I_1}\int_0^Tu_{n1xx}^2
\psi _{\delta }\phi _{\delta }'\Big)\,dt\,dx
\\
&=-\frac{d_1}{2}\int_{I_1}u_{n1xx}^2( T) dx+\frac{d_1}{2}
\int_{I_1}u_{01xx}^2dx \\
&\quad-\frac{d_1}{2\delta }\int_{I_1}\int_0^{\delta }u_{n1xx}^2\psi
_{\delta }\,dt\,dx+\frac{d_1}{2\delta }\int_{I_1}\int_{T-\delta
}^Tu_{n1xx}^2\psi _{\delta }\,dt\,dx
\\
&\leq \frac{d_1}{2\delta }\int_{I_1}\int_{T-\delta }^T(
u_{n1xx}^2-u_{n1xx}^2( T) ) \,dt\,dx
+\frac{d_1}{2\delta }\int_{I_1}\int_0^{\delta }(
u_{01xx}^2-u_{n1xx}^2\psi _{\delta }) \,dt\,dx
\\
&=\frac{d_1}{2\delta }\int_{I_1}\int_{T-\delta }^T(
u_{n1xx}^2-u_{n1xx}^2( T) ) \,dt\,dx
+\frac{d_1}{2\delta }\int_{I_1}\int_0^{\delta }(
u_{01xx}^2-u_{n1xx}^2) \psi _{\delta }\,dt\,dx \\
&\quad +\frac{d_1}{2\delta }\int_{I_1}\int_0^{\delta }u_{01xx}^2(
1-\psi _{\delta }) \,dt\,dx
\end{align*}
and \eqref{521} implies
\begin{equation}
\leq \frac{d_1}{2\delta }\int_{I_1}\int_{T-\delta }^T(
u_{n1xx}^2-u_{n1xx}^2( T) ) \,dt\,dx
+\frac{d_1}{2\delta }\int_{I_1}\int_0^{\delta }(
u_{01xx}^2-u_{n1xx}^2) \psi _{\delta }\,dt\,dx+\eta l_{\ast }.
\label{526}
\end{equation}
Now, the second term in \eqref{526} is dominated by
\begin{align*}
&\frac{d_1}{2\delta }\int_{I_1}\int_0^{\delta }| (
u_{01xx}^2-u_{n1xx}^2) | \,dt\,dx\\
&\leq \frac{d_1}{\delta } \int_{I_1}\int_0^{\delta }\int_0^t |
v_{n1xx}u_{n1xx}| ds\,dt\,dx,
\\
&\leq \frac{d_1}{\delta }\int_0^{\delta }\int_{I_1}\int_0^{\delta
}| v_{n1xx}( x,s) u_{n1xx}( x,s) |\,ds\,dx\,dt
\\
&\leq d_1( \int_{I_1}\int_0^{\delta }| v_{n1xx}(
x,s) | ^2dsdx) ^{1/2}( \int_0^{\delta
}\int_{I_1}| u_{n1xx}( x,s) |
^2\,dx\,ds) ^{1/2} \\
&\leq Cd_1\sqrt{\delta }\equiv C\sqrt{\delta},
\end{align*}
thanks to estimate \eqref{417}. Similar considerations apply to the first
term of \eqref{526}. Therefore, if $\delta $ is small enough, $J_{12}$
in \eqref{524} is no larger than $\eta ( 1+l_{\ast }) $. We
assume below that $\delta $ is this small.

Next, recall the strong convergence results
on the lower order terms in \eqref{52}. It follows that
\begin{align*}
&\limsup_{n\to \infty }\int_{I_1}\int_0^T\Big({
u_{n1t}^2}-d_1v_{n1xx}u_{n1xx}-k_1u_{n1xx}^2
- \frac{1}{3}a_1 ( u_{n1x}+\Phi_x) ^3u_{n1x}\Big) \,dt\,dx
\\
&+\int_{I_1}\int_0^T( \overline{\nu} _1pu_{n1x}^2+\Psi _xu_{n1x})
\,dt\,dx
\\
&\leq \limsup_{n\to \infty }  \int_{I_1}\int_0^T(
u_{n1t}^2-d_1v_{n1xx}u_{n1xx}-k_1u_{n1xx}^2) \psi _{\delta
}\phi _{\delta }\,dt\,dx
\\
&\quad - \lim_{n\to \infty } \; \frac{1}{3}a_1 \int_{I_1}\int_0^T
( u_{n1x}+\Phi_x ) ^3u_{n1x}\, \,dt\,dx
\\
&\quad +\lim_{n\to \infty }\int_{I_1}\int_0^T(\overline{ \nu}
_1pu_{n1x}^2+\Psi _xu_{n1x}) \,dt\,dx
\\
&\quad + \limsup_{n\to \infty }  \int_{I_1}\int_0^T(
u_{n1t}^2-d_1v_{n1xx}u_{n1xx}-k_1u_{n1xx}^2) ( 1-\psi
_{\delta }\phi _{\delta }) \,dt\,dx.
\end{align*}
Now, by Lemma \ref{lemma53} and the above estimates on the integrals
\eqref{523} - \eqref{525}, we find
\begin{align*}
&\limsup_{n\to \infty }\int_{I_1}\int_0^T\Big({
u_{n1t}^2}-d_1v_{n1xx}u_{n1xx}-k_1u_{n1xx}^2
- \frac{1}{3}a_1( u_{n1x}+\Phi_x) ^3u_{n1x}\Big) \,dt\,dx\\
&+\int_{I_1}\int_0^T( \overline{\nu} _1pu_{n1x}^2+\Psi _xu_{n1x})
\,dt\,dx
\\
&\leq   \int_{I_1}\int_0^T(
u_{1t}^2-d_1v_{1xx}u_{1xx}-k_1u_{1xx}^2) \psi _{\delta }\phi
_{\delta }\,dt\,dx
\\
&\quad -\frac{1}{3}a_1 \int_{I_1}\int_0^T( ( u_{1x}+\Phi_x ) ^3u_{1x}
) \,dt\,dx
\\
&\quad +\int_{I_1}\int_0^T( \overline{\nu} _1pu_{1x}^2+\Psi _xu_{1x})
\,dt\,dx+( 3+l_{\ast }) \eta
\end{align*}
since $u_{1t}^2-d_1v_{1xx}u_{1xx}-k_1u_{1xx}^2$ is in $L^{1}$, this
implies that for $\delta $ possibly even smaller,
\begin{align*}
&\leq   \int_{I_1}\int_0^T(
u_{1t}^2-d_1v_{1xx}u_{1xx}-k_1u_{1xx}^2) \,dt\,dx
-\frac{1}{3}a_1 \int_{I_1}\int_0^T( ( u_{1x}+\Phi_x ) ^3
u_{1x}) \,dt\,dx
\\
&\quad +\int_{I_1}\int_0^T( \overline{\nu} _1pu_{1x}^2+\Psi _xu_{1x})
\,dt\,dx+( 4+l_{\ast }) \eta.
\end{align*}
Since $\eta $ is arbitrary, this establishes the lemma.
\end{proof}

The same arguments establish the following lemma.

\begin{lemma} \label{lemma55}
 The following estimate holds,
\begin{align*}
&\limsup_{n\to \infty }\int_{I_2}\int_0^T(
u_{n2t}^2-d_2v_{n2xx}u_{n2xx}-k_1u_{n2xx}^2-\frac{a_2}{3}(
u_{n2x}) ^{4}) \,dt\,dx
\\
&+\int_{I_2}\int_0^T\overline{\nu} _2pu_{n2x}^2\,dt\,dx
\\
&\leq \int_{I_2}\int_0^T(
u_{2t}^2-d_2v_{2xx}u_{2xx}-k_1u_{2xx}^2-\frac{a_2}{3}(
u_{2x}) ^{4}) \,dt\,dx
%\\&
+\int_{I_2}\int_0^T\overline{\nu} _2pu_{2x}^2\,dt\,dx.
\end{align*}
\end{lemma}


Lemmas \ref{lemma54} and \ref{lemma54} imply the following two inequalities.
\begin{equation} \label{527}
\begin{split}
&\liminf_{n\to \infty }\int_{I_2}\int_0^T(
-u_{n2t}^2+d_2v_{n2xx}u_{n2xx}+k_1u_{n2xx}^2+\frac{a_2}{3}(
u_{n2x}) ^{4}) \,dt\,dx
\\
&-\int_{I_2}\int_0^T\overline{\nu} _2pu_{n2x}^2\,dt\,dx
\\
&\geq \int_{I_2}\int_0^T(
-u_{2t}^2+d_2v_{2xx}u_{2xx}+k_1u_{2xx}^2+\frac{a_2}{3}(
u_{2x}) ^{4}) \,dt\,dx
\\
&\quad -\int_{I_2}\int_0^T\overline{\nu} _2pu_{2x}^2\,dt\,dx,
\end{split}
\end{equation}
and
\begin{equation}
\begin{split}
&\liminf_{n\to \infty }\int_{I_1}\int_0^T\Big({
-u_{n1t}^2}+d_1v_{n1xx}u_{n1xx}+k_1u_{n1xx}^2
+ \frac{1}{3}a_1 ( u_{n1x}+\Phi_x
) ^3u_{n1x}\Big) \,dt\,dx
\\
&-\int_{I_1}\int_0^T( \overline{\nu} _1pu_{n1x}^2+\Psi _xu_{n1x})\,dt\,dx
\\
&\geq \int_{I_1}\int_0^T\Big(-{u_{1t}^2+}
d_1v_{1xx}u_{1xx}+k_1u_{1xx}^2\
+\frac{a_1}{3}( u_{1x}+\Phi_x ) ^3u_{1x}) \,dt\,dx
\\
&\quad -\int_{I_1}\int_0^T( \overline{\nu} _1pu_{1x}^2+\Psi _xu_{1x})
\,dt\,dx.
\end{split} \label{528}
\end{equation}
Now, we define  $\Lambda_{n}$ as
\begin{equation} \label{529}
\begin{split}
&\Lambda _{n}(u_{n1}(l_{\ast },t)-u_{n2}(l_{\ast },t))  \\
&\equiv -n( u_{n1}(l_{\ast },t) -u_{n2}(l_{\ast },t) -g_1) _{+}+n( u_{n2}(l_{\ast },t)
-u_{n1}(l_{\ast },t) -g_2) _{+},
\end{split}
\end{equation}
and the two operators
\begin{gather*}
P_{n}( \mathbf{u}_{n}) \equiv \pi _1^{\ast }\gamma ^{\ast
}\Lambda _{n}( u_{n1}(l_{\ast },t) -u_{n2}(l_{\ast },t) ) -\pi _2^{\ast }\gamma ^{\ast }\Lambda _{n}(
u_1(l_{\ast },t) -u_2(l_{\ast },t) ) , \\
P( \mathbf{u}) \equiv \pi _1^{\ast }\gamma ^{\ast }\partial
\chi ( u_1(l_{\ast },t) -u_2(l_{\ast },t)
) -\pi _2^{\ast }\gamma ^{\ast }\partial \chi ( u_1(l_{\ast },t) -u_2(l_{\ast },t) ).
\end{gather*}
As was noted above, these operators are monotone because of the
graph of $\Lambda _{n}$. Also, the solution of the problem with
normal compliance satisfies
\begin{equation}
\mathbf{v} _{n}'+A\mathbf{v} _{n}+B\mathbf{u}_{n}+N( \cdot,
\mathbf{u}_{n}) +P_{n}( \mathbf{u}_{n}) =\mathbf{f},
\label{530}
\end{equation}
together with $\mathbf{v} _{n}( 0) =\mathbf{v} _0$ and
$\mathbf{u}_{n}( t) = \mathbf{u}_0+\int_0^t \mathbf{v} _{n}( s) ds$.  It was shown above that
\begin{equation*}
-g_2\leq u_1(l_{\ast },t) -u_2(l_{\ast },t) \leq g_1.
\end{equation*}

Inequalities \eqref{527} - \eqref{528}
can be written  simply in terms of the operators as
\begin{equation}
\begin{aligned}
&\liminf_{n\to \infty }\Big( -\int_0^T( \mathbf{v} _{n},
\mathbf{v} _{n}) _Hdt+\int_0^T\langle A\mathbf{v} _{n}, \mathbf{u}_{n}\rangle dt \\
&+\int_0^T\langle B\mathbf{u}_{n},\mathbf{u}_{n}\rangle
dt +\int_0^T\langle N( \cdot, \mathbf{u}_{n}),
\mathbf{u}_{n}\rangle dt\Big)
\\
&\geq \Big( -\int_0^T( \mathbf{v} ,\mathbf{v} )
_Hdt+\int_0^T\langle A\mathbf{v} ,\mathbf{u}\rangle
dt+\int_0^T\langle B\mathbf{u}, \mathbf{u}\rangle dt
 +\int_0^T\langle N( \cdot ,\mathbf{u}), \mathbf{u}
\rangle dt\Big).
\end{aligned} \label{531}
\end{equation}
Here, we used the notation $\langle\cdot, \cdot \rangle$ for the duality pairing
between $V$ and $V'$.

Now, let $\mathbf{w}, \mathbf{w}'\in \mathcal{V}$ such that
$\mathbf{w}( T) =\mathbf{u}( T) $ and for each $t$,
\begin{equation*}
-g_2\leq w_1(l_{\ast },t) -w_2(l_{\ast },t)
\leq g_1.
\end{equation*}
We apply \eqref{530} to $\mathbf{u}_{n}-\mathbf{w}$ and perform time
integration on both sides. Thus,
\begin{align*}
&-\int_0^T( \mathbf{v}_{n},\mathbf{v}_{n}-\mathbf{w}^{'}) dt+
( \mathbf{v}_0,\mathbf{w}( 0) -\mathbf{u}_0) _H
+\int_0^T\langle A\mathbf{v}_{n},\mathbf{u}_{n}-
\mathbf{w}\rangle ds
\\
&+\int_0^T\langle B\mathbf{u}_{n},\mathbf{u}_{n}-\mathbf{w}
\rangle ds+\int_0^T\langle N( \cdot ,\mathbf{u}
_{n}) ,\mathbf{u}_{n}-\mathbf{w}\rangle dt
+\int_0^T\langle P_{n}( \mathbf{u}_{n}) ,\mathbf{u}_{n}-
\mathbf{w}\rangle dt\\
&=\int_0^T( \mathbf{f},\mathbf{u}_{n}-\mathbf{w}) _Hds.
\end{align*}
Since $P_{n}( \mathbf{w}) =0$, it follows that
\begin{equation*}
\int_0^T\langle P_{n}( \mathbf{u}_{n}) ,\mathbf{u}_{n}-
\mathbf{w}\rangle dt=\int_0^T\langle P_{n}( \mathbf{u}
_{n}) -P_{n}( \mathbf{w}) ,\mathbf{u}_{n}-\mathbf{w}
\rangle dt\geq 0;
\end{equation*}
therefore,
\begin{align*}
&-\int_0^T( \mathbf{v} _{n},\mathbf{v} _{n}-\mathbf{w}') dt+( \mathbf{v} _0,\mathbf{w}( 0) -\mathbf{u}
_0) _H+\int_0^T\langle A\mathbf{v} _{n},\mathbf{u}_{n}-
\mathbf{w}\rangle ds
\\
&+\int_0^T\langle B\mathbf{u}_{n},\mathbf{u}_{n}-\mathbf{w}
\rangle ds+\int_0^T\langle N( \cdot ,\mathbf{u}
_{n}) ,\mathbf{u}_{n}-\mathbf{w}\rangle dt
\\
&\leq\int_0^T( \mathbf{f},\mathbf{u}_{n}-\mathbf{w}) _Hds.
\end{align*}
Passing to the $\liminf_{n\to \infty }$ of both sides, it follows from
\eqref{531} that
\begin{equation}
\begin{split}
&-\int_0^T( \mathbf{v} ,\mathbf{v} -\mathbf{w}')
dt+( \mathbf{v} _0,\mathbf{w}( 0) -\mathbf{u}_0)
_H+\int_0^T\langle A\mathbf{v} ,\mathbf{u}-\mathbf{w}\rangle ds
\\
&+\int_0^T\langle B\mathbf{u},\mathbf{u}-\mathbf{w}\rangle
ds+\int_0^T\langle N( \cdot ,\mathbf{u}) ,\mathbf{u}-
\mathbf{w}\rangle dt
\\
&\leq \int_0^T( \mathbf{f},\mathbf{u}-\mathbf{w}) _Hds.
\end{split} \label{532}
\end{equation}
This concludes the proof of Theorem \ref{t32}, which we restate as:

\begin{theorem}
Assume that {\rm (A1)--(A5)} hold along with the regularity assumption on
the initial data
\begin{equation*}
u_{0ixx}\in L^{\infty }( I_{i}) .
\end{equation*}
Then, there exists a pair $(\mathbf{u},\mathbf{v})$ such that
\begin{equation*}
\mathbf{u}( t) =\mathbf{u}_0+\int_0^t \mathbf{v}(s) ds,
\end{equation*}
$\mathbf{v}\in \mathcal{V}$, and $\mathbf{u}( t) \in \mathcal{\ K}$
such that whenever $\mathbf{w} ( t) \in \mathcal{K}$
with $\mathbf{w},\mathbf{w}'\in \mathcal{V}$, and
$\mathbf{w}(T) =\mathbf{u}( T)$, then the
variational inequality \eqref{532} holds.
\end{theorem}

\section{The inviscid case}
\label{sec:novisc}

In this section we study the case with the Signorini contact condition
when there is no viscosity, so, as above, we replace the normal compliance
stiffness coefficients $\lambda_1, \lambda_2$ and the viscosity coefficients
$d_1, d_2$ with the modified coefficients
\[
n \lambda_1,\quad  n \lambda_2, \quad  \frac{d_1}{n},\quad \frac{d_2}{n},
\]
and consider a sequence of solutions, one for each positive integer $n$.
Then, we obtain the necessary estimates and consider  the limit when $n \to \infty$.
The above discussion of the normal compliance case yields the following result.

\begin{theorem}\label{t61}
Assume that {\rm (A1)--(A5)} hold. Then, for each $n$ there
exists a unique solution to the abstract problem
\begin{equation}
\begin{split}
&\mathbf{v} '+\frac 1n A\mathbf{v} +B\mathbf{u}+N( \cdot ,\mathbf{u} )
+\pi _1^{\ast }\gamma ^{\ast }  n \Lambda ( u_1(
l_{\ast }) -u_2( l_{\ast }) ) \\
&-\pi _2^{\ast}\gamma ^{\ast }  n \Lambda ( u_1( l_{\ast })
-u_2( l_{\ast }) ) =\mathbf{f},
\end{split}\label{61}
\end{equation}
together with $\mathbf{v} ( 0) =\mathbf{v} _0$, and
$\mathbf{u}( t) = \mathbf{u}_0+\int_0^t \mathbf{v} ( s) ds$.

This solution satisfies for a.e. $t\in (0, T)$  the estimate
\begin{equation}
\begin{split}
&| \mathbf{v} ( t) | _H^2+
\frac1n \min(d_1, d_2) \int_0^t \| \mathbf{v} ( s) \| _{V}^2ds
+\int_0^{l_{\ast}}u_{1x}( t) ^{4}dx\\
&+\int_{l_{\ast }}^{1}u_{2x}( t)^{4}dx+\| \mathbf{u}( t) \| _{V}^2
+  nS ( u_1( t,l_{\ast }) -u_2( t,l_{\ast
}) )
\\
&\leq C( p_0,\mathbf{v}_0, \mathbf{u}_0,C_0,k_1,k_2,T)
\end{split}\label{62}
\end{equation}
where the constant $C$ depends on the indicated quantities but is
independent of $n$ and $d_{i}$.
\end{theorem}

Estimate \eqref{62} yields the same type of convergence results that were
obtained above, with the exception of \eqref{55}, which is lost but is replaced
with a related  convergence. Specifically, we have
\begin{equation} \label{63}
\begin{gathered}
u_{inx}\to u_{ix}\quad \text{strongly in }C( [0,T];C( I_{i}) ), \\
u_{ni}\to u_{i}\quad \text{strongly in }C( [0,T];C( I_{i}) ),  \\
\mathbf{u}_{n}\to \mathbf{u}\quad \text{weak}\ast \text{ in }L^{\infty}( 0,T;V) \\
\frac1{\sqrt{n}} \mathbf{v} _{nxx}\to 0 \text{ strongly in }\mathcal{H},\\
\mathbf{v} _{n}\to \mathbf{v}\quad  \text{weak}\ast \text{ in }L^{\infty}( 0,T;H).
\end{gathered}
\end{equation}
It follows from \eqref{62} and \eqref{63}, as above, that
\begin{equation*}
-g_2\leq ( u_1( t,l_{\ast }) -u_2( t,l_{\ast}) ) \leq g_1.
\end{equation*}
Therefore, $\mathbf{u}\in \mathcal{K}$ and so it satisfies the necessary constraint.

Estimates \eqref{62} and \eqref{61} imply that, just as above,
the function $v_{i}'$ is bounded in $L^2( 0,T;( H_0^2(I_{i}) ) ')$. Therefore, we can
also assume that
\begin{equation*}
v_{ni}'\to v'\text{ in }L^2\Big( 0,T;(
H_0^2( I_{i}) ) '\Big).
\end{equation*}
Now, the rest of the argument is similar to the above except for the issue
involving the viscous term. The formula \eqref{524} now takes the form
\begin{equation}
-\frac{d_1}n \int_{I_1}\int_0^T v_{n1xx}u_{n1xx}( 1-\psi _{\delta
}\phi _{\delta }) \,dt\,dx.
\end{equation}
Then, all the same considerations hold in the argument for approximation.
Consider next the integral
\begin{equation*}
\frac 1n \int_{I_1}\int_0^T v_{n1xx}u_{n1xx}\,dt\,dx,
\end{equation*}
which converges to 0 because of the estimate
\begin{align*}
& \big| \frac 1n \int_{I_1}\int_0^T v_{n1xx}u_{n1xx}\,dt\,dx\big| \\
&\leq \Big( \frac 1{n^2}\int_{I_1}\int_0^T v_{nxx}^2\,dt\,dx\Big)
^{1/2}\Big( \int_{I_1}\int_0^Tu_{n1xx}^2\,dt\,dx\Big) ^{1/2}
\\
&\leq \frac 1{\sqrt{n}}\Big( \int_{I_1}\int_0^T\frac 1n v_{nxx}^2\,dt\,dx\Big)
^{1/2}\Big( \int_{I_1}\int_0^Tu_{n1xx}^2\,dt\,dx\Big) ^{1/2},
\end{align*}
which converges to 0 by the estimates derived above.

It follows that
\begin{align*}
&\liminf_{n\to \infty }\Big( -\int_0^T( \mathbf{v} _{n},
\mathbf{v} _{n}) _Hdt+\int_0^T\langle \frac 1n A\mathbf{v} _{n},
\mathbf{u}_{n}\rangle dt \\
&+\int_0^T\langle B\mathbf{u}_{n},\mathbf{u}_{n}\rangle
dt +\int_0^T\langle N( \cdot ,\mathbf{u}_{n}),
\mathbf{u}_{n}\rangle dt\Big)
\\
&\geq \Big( -\int_0^T( \mathbf{v} ,\mathbf{v} )
_Hdt+\int_0^T\langle B\mathbf{u},\mathbf{u}\rangle
dt +\int_0^T\langle N( \cdot ,\mathbf{u}),
\mathbf{u}\rangle dt\Big),
\end{align*}
so that the viscous term disappears. Then, the same argument as above yields
the theorem for the inviscid case and with the Signorini contact condition.

\begin{theorem}\label{t62}
Assume that {\rm (A1)--(A5)} hold along with the regularity
assumption on the initial data
\begin{equation*}
u_{0ixx}\in L^{\infty }( I_{i}).
\end{equation*}
Then, there exists a pair $\mathbf{u},\mathbf{v}\in \mathcal{V}$,
$\mathbf{u}( t) =\mathbf{u}_0+\int_0^t \mathbf{v} (s) ds$,
$\mathbf{u}( t) \in \mathcal{K}$  such that whenever $\mathbf{w}
( t) \in \mathcal{K}$ with $\mathbf{w}'\in \mathcal{V}
$, $\mathbf{w}( T) =\mathbf{u}( T) $, it follows
\begin{align*}
&-\int_0^T( \mathbf{v} ,\mathbf{v} -\mathbf{w}')
_Hdt +\int_0^T\langle B\mathbf{u},\mathbf{u}-\mathbf{w}\rangle
ds+\int_0^T\langle N( \cdot ,\mathbf{u}) ,\mathbf{u}- \mathbf{w}\rangle dt
\\
&\leq  \int_0^T( \mathbf{f},\mathbf{u}-
\mathbf{w}) _Hds -( \mathbf{v} _0,\mathbf{w}( 0) -\mathbf{u}_0) _H.
\end{align*}
\end{theorem}

To obtain the existence  result for the inviscid problem with normal compliance we
replace the inequality with the following variational equation
\begin{align*}
&-\int_0^T( \mathbf{v} ,\mathbf{v} -\mathbf{w}')
_Hdt +\int_0^T\langle B\mathbf{u},\mathbf{u}-\mathbf{w}\rangle
ds+\int_0^T\langle N( \cdot ,\mathbf{u}) ,\mathbf{u}-
\mathbf{w}\rangle dt
\\
&+ \int_0^T\Lambda( u_1( t,l_{\ast })
-u_2( t,l_{\ast }) ) ( w_1( t,l_{\ast
}) ) dt \\
&-\int_0^T\Lambda ( u_1( t,l_{\ast })
-u_2( t,l_{\ast }) ) ( w_2( t,l_{\ast
}) ) dt
\\
&=  \int_0^T( \mathbf{f},\mathbf{u}-
\mathbf{w}) _Hds -( \mathbf{v} _0,\mathbf{w}( 0) -\mathbf{u}
_0) _H,
\end{align*}
and note that the requirement that $\mathbf{u}\in \mathcal{K}$ is not needed.


This concludes the analysis part of this work. The rest of the paper deals
with the numerical aspects  of the model.

\section{Computational method; numerical algorithm}
\label{sec:algorithm}

In this section we describe our numerical schemes for the model.
Unlike the Euler--Bernoulli or the Timoshenko beams, the two Gao beams
are nonlinear, which makes the approach more complex.

We choose a  fully discrete numerical method and use
a uniform discretization of the time domain $[0, T]$ and a mixture of the
Galerkin approximation and central difference formula to discretize the space
domain $[0, 1]$. The central difference formula is combined with the finite
element method (FEM) when we deal with the nonlinear Gao terms in the
two equations. Convergence results of similar numerical
schemes have been investigated in \cite{ahn:dfcgts}, using a
modified truncated operator in the discrete case. We note that a
different numerical algorithm for a Gao beam, but without contact, was presented
in \cite{adbps:asndb}, where the fully explicit finite difference method (FDM) was used.

Since the equation of a Gao beam contains second and
fourth order derivatives, we use a nonconforming mixed finite
element method. Indeed, nonconforming methods which have been studied in
many papers and books (e.g., \cite{sl:mtfem,p:femep} and the references therein)
provide more efficient and practical algorithms for the fourth order or even
higher order PDEs. The mixed method that yields a saddle-point problem
has been first proposed in \cite{pj:mfem} and its convergence analysis
has been studied extensively, see, e.g., \cite{ijj:ammumdn, db:fe} and
the references therein.

We now describe our numerical schemes. The interval
$[0,\, 1]$ is divided into $m$ subintervals $I_{i}=[x_{i},\, x_{i+1}]$,
for $0\leq i\leq m-1$, with grid points
\[
0=x_0<x_1<x_2<\dots<x_{i_{*}}<\dots<x_{m-2}<x_{m-1}<x_{m}=1.
\]
For the sake of simplicity, a uniform spacing is used and the mesh size is denoted
by $h_x=x_{i+1}-x_{i}$. We assume that the
two beams are joined at the node $x=x_{i_{*}}$. According to the
nonconforming finite element methods and Sobolev imbedding theorem,
we choose $V_{h_x}$ as the finite dimensional space
\[
V_{h_x}=\{ w_{h_x}\in C[0,\, 1]: w_{h_x}\big|_{I_{i}}\in\mathcal{L}
(I_{i})\cup\mathcal{B}(I_{i}),\,1\leq i\leq m-1\},
\]
where $\mathcal{L}$ is the family of piecewise linear functions and
$\mathcal{B}$ is the family of the piecewise cubic functions.
The  finite element space can be written in the equivalent way
\[
V_{h_x}=\text{span}\{ \psi_{i}\in C[0,\, x_{i_{*}}]\cup
C[x_{i_{*}},\, 1]: 1\leq i\leq m-1\} .
\]
The basis functions $\psi_{i}$ are constructed next.
We note that on the node $x_{i_{*}}$ two different basis functions  have
to be constructed separately, since one is associated with the right
end of the left beam and the other is with the left end of the right
beam. We construct the basis functions $\psi_{i}=\psi_{1,i}$ with
compact supports $[x_{i-1},\, x_{i+1}]$ and $\psi_{i}=\psi_{3,i}$
with compact supports $[x_{i-2},\, x_{i+2}]$; the first
is a typical piecewise linear function $L(s)$ for the basis functions
denoted by $\psi_{1,i}\in C[0,\, 1]$ and the other is a cubic spline
functions $B(s)$ for the basis functions denoted by
$\psi_{3,i}\in C^{1}[0,\, 1]$.
It is easy to construct the piecewise linear functions on the reference
interval $[-1,1]$ and the standard cubic spline functions $B(s)$
on the reference interval $[-2,\,2]$,
\[
L(s)=\begin{cases}
-|s|+1 &\text{if }|s|\leq1,\\
0&\text{if }|s|\geq1,
\end{cases}
\]
and
\[
B(s)=\frac{2}{3} \begin{cases}
1+\frac{3}{4}|s|^3-\frac{3}{2}|s|^2 &\text{for }|s|\leq1,\\
\frac{1}{4}(2-|s|)^3 &\text{for }1\leq|s|\leq2,\\
0 &\text{for }|s|\geq2.
\end{cases}
\]
For the piecewise linear functions with $1<i\leq m-1$ the basis functions
$\psi_{1,i}$ is defined by the shifted linear function
\[
\psi_{1,i}(x)=L\big(\frac{x-x_{i}}{h_x}\big),
\]
where each node is $x_{i}=i\, h_x$ and their supports
are $[x_{i}-h_x,\, x_{i}+h_x]$. For the B-spline functions with
$2<i\leq i_{*}-2$ and $i_{*}+2<i\leq m-1$ we use mostly the basis
functions $\psi_{3,i}$ defined by the shifted cubic splines
\[
\psi_{3,i}(x)=B\big(\frac{x-x_{i}}{h_x}\big),
\]
where each node is $x_{i}=i\, h_x$ and their supports
are $[x_{i}-2h_x,\, x_{i}+2h_x]$. While it is not hard to set
up the piecewise linear functions $\psi_{1,i}$ with the boundary
conditions, we need to be more careful about constructing the basis
functions $\psi_{3,i}$ so as to satisfy the boundary conditions. In
order to satisfy all the conditions, the cubic functions are modified
as follows:
\[
\psi_{3,1}(x)= \begin{cases}
-\frac{5}{3}(\frac{x}{h_x})^3+\frac{5}{2}(\frac{x}{h_x})^2
&\text{on }[0,\, h_x],\\
\frac{5}{6}(\frac{x}{h_x}-1)^3-\frac{3}{2}(\frac{x}{h_x}-1)^2+
\frac{5}{6} &\text{on }[h_x,\,2h_x],\\
\frac{1}{6}(2-(\frac{x}{h_x}-1))^3 &\text{on }[2h_x,\,3h_x],
\end{cases}
\]
and
\begin{align*}
&\psi_{3,m-1}(x)\\
&= \begin{cases}
\frac{5}{3}(\frac{x}{h_x}-m)^3+\frac{5}{2}(\frac{x}{h_x}-m)^2
&\text{on }[1-h_x,\, 1],\\
-\frac{5}{6}(\frac{x}{h_x}-(m-1))^3-\frac{3}{2}(\frac{x}{h_x}-(m-1))^2+\frac{5}{6}
&\text{on }[1-2h_x,\, 1-h_x],\\
\frac{1}{6}(2+(\frac{x}{h_x}-(m-1)))^3 &\text{on }[1-3h_x,\, 1-2h_x].
\end{cases}
\end{align*}
We need to construct the two basis functions
satisfying the natural boundary conditions imposed at $x=l_{*}$.
The basis functions $\psi_{i_{*}-1}$ and $\psi_{i_{*}}$ are
\begin{align*}
&\psi_{3,i_{*}-1}(x)\\
&= \begin{cases}
-\frac{5}{3}(\frac{x}{h_x}-(i_{*}-1))^3+\frac{5}{2}(\frac{x}{h_x}
-(i_{*}-1))^2 &\text{on } {  [x_{i_{*}}-h_x,\, x_{i_{*}}]},\\
\frac{5}{6}(\frac{x}{h_x}-(i_{*}-2))^3-\frac{3}{2}(\frac{x}{h_x}
-(i_{*}-2))^2+\frac{5}{6}
&\text{on }[{  x_{i_{*}}-\,2h_x,\, x_{i_{*}}-h_x}],\\
\frac{1}{6}(2-(\frac{x}{h_x}-(i_{*}-2)))^3 &\text{on }
[{  x_{i_{*}}-3h_x,\, x_{i_{*}}-2h_x}],
\end{cases}
\end{align*}
 and
\[
\psi_{3,i_{*}}(x)=-\frac{2}{h_x^3}(x-(i_{*}-1))^3-\frac{3}{h_x^2}
(x-(i_{*}-1))^2\quad\text{on }[{  x_{i_{*}}-h_x,\, x_{i_{*}}}].
\]
Similarly for the right Gao beam, we can obtain
\[
\psi_{3,i_{*}}(x)=\frac{2}{h_x^3}(x-(i_{*}+1))^3+\frac{3}{h_x^2}
(x-(i_{*}+1))^2\quad\text{on }[{  x_{i_{*}},\, x_{i_{*}}+h_x}].
\]
and
\begin{align*}
&\psi_{3,i_{*}+1}(x)\\
&= \begin{cases}
-\frac{5}{3}(\frac{x}{h_x}-i_{*})^3+\frac{5}{2}(\frac{x}{h_x}-i_{*})^2
&\text{on }[{  x_{i_{*}},\, x_{i_{*}}+h_x}],\\
\frac{5}{6}(\frac{x}{h_x}-(i_{*}+1))^3-\frac{3}{2}(\frac{x}{h_x}
-(i_{*}+1))^2+\frac{5}{6}
& \text{on }[{  x_{i_{*}}+h_x,\, x_{i_{*}}+2h_x}],\\
-\frac{1}{6}(2-(\frac{x}{h_x}-(i_{*}+2)))^3
&\text{on }[{  x_{i_{*}}+2h_x,\, x_{i_{*}}+3h_x}].
\end{cases}
\end{align*}


The time interval $[0, T]$ is also uniformly partitioned as
\[
0=t_0<t_1<t_2\dots<t_{k}<\dots<t_{n-1}<t_{n}=T,
\]
where $h_t=t_{k}-t_{k-1}$, for $1\leq k\leq n$. We suppose that
at each time step $t=t_{k}$ the fully discrete solutions $u_{1h_t,h_x}^k$
and $u_{2h_t,h_x}^k$ are represented by linear combinations
of B-splines over the domain $[0,\, 1]$ as
\[
u_{1h_t,h_x}^k=\sum_{i=1}^{i_{*}}u_{1i}^k\psi_{d,i}(x),\quad
u_{2h_t,h_x}^k=\sum_{i=i_{*}}^{m-1}u_{2i}^k\psi_{d,i}(x),
\]
where $d$ is the degree of the B-splines. When we deal with the $0$th
and $2$nd order partial derivative with respect to $x$, we write
the fully discrete solutions of the displacement and velocity as
\begin{equation} \label{eq:fully_discrete_linear_sols}
\begin{gathered}
u_{1h_t,h_x}^k=\sum_{i=1}^{i_{*}}u_{1i}^k\psi_{1,i}(x),  \quad
 u_{2h_t,h_x}^k
=\sum_{i=i_{*}}^{m-1}u_{2i}^k\psi_{1,i}(x), \\
(u_{1h_t,h_x}^k)_t=\sum_{i=1}^{i_{*}}v_{1i}^k\psi_{1,i}(x),  \quad
(u_{2h_t,h_x}^k)_t=\sum_{i=i_{*}}^{m-1}v_{2i}^k\psi_{1,i}(x).
\end{gathered}
\end{equation}
Similarly, the solutions of the fourth-order partial derivative are
\begin{equation}
u_{1h_t,h_x}^k=\sum_{i=1}^{i_{*}}u_{1i}^k\psi_{3,i}(x),
\quad u_{2h_t,h_x}^k=\sum_{i=i_{*}}^{m-1}u_{2i}^k\psi_{3,i}(x).
\label{eq:fully_discrete_bspline_sols}
\end{equation}
To derive the mixed method, we rewrite the original problem
as the following mixed form,
\[
\begin{bmatrix}
u_{1tt}\\
\omega_1\\
u_{2tt}\\
\omega_2
\end{bmatrix}
=\begin{bmatrix}
-k_1\omega_1-d_1\omega_{1t}+a_1(u_{1x})^2u_{1xx}
-\bar{\nu}_1\, p(t)\, u_{1xx}+f_1(t)\\
u_{1xxxx}\\
-k_2\varpi_2-d_2\varpi_{2t}+a_2(u_{2x})^2u_{2xx}-\bar{\nu}_2\,
p(t)\, u_{2xx}+f_2(t)\\
u_{2xxxx}
\end{bmatrix}.
\]
Multiplying each equation by a suitable test function and using
integration by parts, it is straightforward to obtain the following
weak formulation:
\\
Find $(u_1,\,\omega_1,u_2,\,\omega_2)\in H^2(0,\, l_{*})\times
L^2(0,\, l_{*})\times H^2(l_{*},\, 1)\times L^2(l_{*},\, 1)$
such that
\begin{align*}
0 & =  (u_{1tt}+k_1\omega_1+d_1\omega_{1t}-a_1(u_{1x})^2u_{1xx}
+\bar{\nu}_1\, p(t)\, u_{1xx}-f_1(t),\, w)_{L^2(0,\, l_{*})}\\
 & \quad +\sigma(t)\, w(l_{*})\quad\text{for all }w\in L^2(0,l_{*}),
\\
0 & =  (\omega_1,\, v)_{L^2(0,\, l_{*})}-(u_{1xx},\,
v_{xx})_{L^2(0,\, l_{*})}\quad\text{for all }v\in H^2(0,\, l_{*})
,\\
0 & =  \Big(u_{2tt}+k_2\varpi_2+d_2\varpi_{2t}-a_2(u_{2x})^2w_{2xx}
+\bar{\nu}_2\, p(t)\, w_{2xx}-f_2(t),\, u\Big)_{L^2(l_{*},\, 1)}
\\
 &\quad -\sigma(t)\, w(l_{*})\quad\text{for all }u\in L^2(l_{*},\, 1),\\
0 & =  (\omega_2,\,\mu)_{L^2(l_{*},\, 1)}-(w_{2xx},\,
\mu_{xx})_{L^2(l_{*},\, 1)}\quad\text{for all }\mu\in H^2(l_{*},\, 1),
\end{align*}
where at $x=l_{*}$ we have $\sigma(t)=\sigma_1(t)= - \sigma_2(t)$.
Then, the corresponding FEM
can be applied easily. We assume that the auxiliary variables $\omega_1=u_{1xxxx}$
and $\varpi_2=u_{2xxxx}$ and $\omega_{1t}=u_{1txxxx}$ and $\varpi_{2t}=u_{2txxxx}$
are  linear combinations of the piecewise constant functions
$\psi_{0,i}$ with $1\leq i\leq m-1$,
\begin{gather*}
  \omega_{1h_t,h_x}^k=\sum_{i=1}^{i_{*}}\omega_{1i}^k\psi_{0,i}(x),\quad
\omega_{2h_t,h_x}^k=\sum_{i=i_{*}}^{m-1}\omega_{2i}^k\psi_{0,i}(x),
\\
  (\omega_{1h_t,h_x}^k)_t=\sum_{i=1}^{i_{*}}\varpi_{1i}^k\psi_{0,i}(x),
 \quad (\omega_{2h_t,h_x}^k)_t=\sum_{i=i_{*}}^{m-1}\varpi_{2i}^k\psi_{0,i}(x)
 \quad\text{at }t=t_{k}.
\end{gather*}
When the weak formulation is written in the corresponding matrix
form, the supports of the $\psi_{0,i}$ change, since
when we deal with the B-spline of the degree $d=0$, we use $\psi_{0,i}(x)=1$
over $[x_{i-1}, x_{i+1}]$ and when we deal with cubic
B-spline, we use $\psi_{0,i}(x)=1$ over $[x_{i-2}, x_{i+2}]$.
In the discrete form, we obtain the following equation by using
the test functions. Let
\begin{gather*}
\boldsymbol{\omega}_{h_t,h_x}^k=(\omega_{1h_t,h_x}^k,\,
\omega_{2h_t,h_x}^k)^T,\quad
(\boldsymbol{\omega}_{h_t,h_x}^k)_t=(\varpi_{1h_t,h_x}^k,\,
\varpi_{2h_t,h_x}^k)^T,
\\
\mathbf{u}_{h_t,h_x}^k=(u_{1h_t,h_x}^k,\, u_{2h_t,h_x}^k)^T, \quad
\mathbf{v}_{h_t,h_x}^k=(v_{1h_t,h_x}^k,\, v_{2h_t,h_x}^k)^T,
\end{gather*}
at each time step $t_{k}$. For the sake of simplicity, we assume
that there are no body forces, i.e., $f_1=f_2=0$. In particular, this means that
the left end of the left beam is clamped.
Then, the following fully discrete numerical formulation holds in the distributional
sense; for $i=1,2$,
\begin{equation} \label{eq:fully1}
\begin{split}
&\frac{1}{h_t}(v_{i h_t,h_x}^{k+1}-v_{ih_t,h_x}^k)\\
&=  -\frac{k_{i}}{2}(\omega_{ih_t,h_x}^{k+1}+
\omega_{ih_t,h_x}^k) +(-1)^{i+1}\sigma_{ih_t,h_x}^k\\
&\quad -\frac{d_{i}}{2}
\Big((\omega_{ih_t,h_x}^k)_t+(\omega_{ih_t,h_x}^{k+1})_t\Big)
  +\Big(((u_{ih_t,h_x}^k)')^2-\overline{\nu}\,
 p^k\Big)(u_{ih_t,h_x}^k)'',
\end{split}
\end{equation}
 \begin{gather}
\omega_{ih_t,h_x}^k  =  (u_{ih_t,h_x}^k)'''',\quad
(\omega_{ih_t,h_x}^k)_t=(v_{ih_t,h_x}^k)'''',\label{eq:substituion}\\
\frac{1}{h_t}(u_{ih_t,h_x}^{k+1}-u_{ih_t,h_x}^k)
 = \frac{1}{2}(v_{ih_t,h_x}^{k+1}+v_{ih_t,h_x}^k).
\label{eq:fully2}
\end{gather}
By the contact condition \eqref{28} the numerical approximation of contact
forces becomes
$\sigma_{h_t,h_x}^k=\sigma_{1h_t,h_x}^k=-\sigma_{2h_t,h_x}^k$.

We turn to the contact of the two Gao beams at the joint.
The contact at  $x=x_{i_{*}}$ is described by
the normal compliance condition \eqref{214},
where we set $\lambda=\lambda_1=\lambda_2$,  thus at $t=t_{k}$,
\begin{equation} \label{eq:NC}
\begin{aligned}
\sigma_{h_t,h_x}^k:=\sigma_{h_t,h_x}(t_{k})
& =  \lambda (u_{2h_t,h_x}^{k+1}
(x_{i_{*}})-u_{1h_t,h_x}^{k+1}(x_{i_{*}})-g_2)_{+} \\
&\quad  -\lambda (u_{1h_t,h_x}^{k+1}(x_{i_{*}})-
 u_{2h_t,h_x}^{k+1}(x_{i_{*}})-g_1)_{+},
\end{aligned}
\end{equation}
where the normal compliance stiffness coefficient $\lambda>0$ is likely to be large.

As can be seen, over the time interval $[0,\, T]$, we use a hybrid of three
numerical schemes: for the elasticity and viscosity the midpoint rule is applied
and for the contact conditions the implicit Euler method is used,
and for the nonlinear term in the equation of motion the central difference
formula is applied. The fully discrete approximation of the displacement,
$\mathbf{u}_{h_t,h_x}$ is a linear interpolant with
\[
\mathbf{u}_{h_t,h_x}(\cdot,\, t_{k})
=\mathbf{u}_{h_x}^k, \quad  \text{and} \quad \mathbf{u}_{h_t,h_x}(\cdot,\, t_{k+1})
=\mathbf{u}_{h_x}^{k+1},
\]
and the fully discrete approximation of the velocity $\mathbf{v}_{h_t,h_x}$
is a constant interpolant with $\mathbf{v}_{h_t,h_x}(\cdot,\, t)=\mathbf{v}_{h_x}^{k+1}$
for $t\in(t_{k},\, t_{k+1}]$. We regard
the nonlinear term $((\mathbf{u}_{h_t,h_x}^k)')^2$
as the square of a strong derivative, so we apply the central difference
formula into it. Thus, at each time step $t_{k}=k\, h_t$, we deal with
the following numerical differentiations, at $x=x_{i}$,
\[
\Big((u_{1h_t,h_x}^k)'\Big)^2
=\Big(\frac{u_{1i+1}^k-u_{1i-1}^k} {2h_x}\Big)^2,\quad
\Big((u_{2h_t,h_x}^k)'\Big)^2
=\Big(\frac{u_{2i+1}^k-u_{2i-1}^k}{2h_x}\Big)^2.
\]
Next, we set up the linear system associated with
the nonlinear term. The coefficient matrix, at each
time step $t=t_{k}$, has the  form
\begin{gather*}
\mathbf{L}_1^k  \equiv  (L_1^k)_{ij}
=\int_0^{l_{*}}\psi_{1,i}'(x)\psi_{1,j}'(x)
\Big(a_1\Big(\frac{u_{1i+1}^k-u_{1i-1}^k}{2h_x}\Big)^2-\overline{\nu}\, p^k\big)dx,
\\
\mathbf{L}_2^k  \equiv  (L_2^k)_{ij}
=\int_{l_{*}}^{l}\psi_{1,i}'(x) \psi_{1,j}'(x)
\Big(a_2\Big(\frac{(u_{2i+1}^k-u_{2i-1}^k)}{2h_x}\Big)^2 -\overline{\nu}\, p^k\Big)dx.
\end{gather*}
Thus, we obtain the tridiagonal finite element coefficient matrix. The main
diagonal is assembled in the way that
if $i=j$, $(L_1^k)_{ij}=\varsigma_{i-1}^k+\varsigma_{i}^k$
for $1\leq i\leq i_{*}$. The two subdiagonals are assembled as follows;
if $j=i+1$ or $(L_1^k)_{ij}=-\varsigma_{i}^k$ for
$1\leq i\leq i_{*}-1$ and if $i=j+1$ or $(L_1^k)_{ij}=-\varsigma_{i}^k$
for $1\leq j\leq i_{*}-1$. If $|i-j|>1$, $(L_1^k)_{ij}=0$.
Here each element of the symmetric matrix $L_{ij}^k$ with $1\leq i,\, j\leq i_{*}$
is defined by
\begin{align*}
\varsigma_{i}^k
&=\int_{x_{i}}^{x_{i+1}}\Big(a_1(\Big(frac{u_{1i+1}^k-u_{1i-1}^k}{2h_x}\Big)^2
-\overline{\nu}\, p^k\Big) dx
\\
&=\frac{a_1(u_{1i+1}^k-u_{1i-1}^k)}{4h_x}^2-\overline{\nu}\, p^kh_x.
\end{align*}
Similarly, the coefficient matrix $\mathbf{L}_2^k$ can be assembled
at each time step $t=t_{k}$. When we compute the next step solutions
$(\mathbf{u}_{h_x}^{k+1},\,\mathbf{v}_{h_x}^{k+1})$
satisfying the fully discrete formulations \eqref{eq:fully1}--\eqref{eq:fully2},
we multiply both sides of \eqref{eq:fully1} by the basis functions
$\psi_{1,i}(x)$ and apply integration by parts, thus,
\begin{equation}
\begin{pmatrix}
\frac{1}{h_t}\sum_{j=1}^{i_{*}}(v_{1j}^{k+1}-v_{1j}^k)
 \int_0^{l_{*}}\psi_{1,j}(x)\psi_{1,i}(x)dx\\
\frac{1}{h_t}\sum_{j=i_{*}}^{m-1}(v_{2j}^{k+1}-v_{2j}^k)
 \int_{l_{*}}^{0}\psi_{1,j}(x)\psi_{1,i}(x)dx
\end{pmatrix}
=\begin{pmatrix}
\mathbf{D}_1\\
\mathbf{D}_2
\end{pmatrix},\label{eq:fully3}
\end{equation}
where
\begin{equation} \label{eq:d1}
\begin{split}
\mathbf{D}_1 & =  -\frac{k_1}{2}\sum_{j=1}^{i_{*}}(\omega_{1j}^{k+1}
+\omega_{1j}^k)\int_0^{l}\psi_{0,j}(x)\psi_{1,i}(x)dx \\
 & \quad -\frac{d_1}{2}\sum_{j=1}^{i_{*}}(\varpi_{1i}^k+
 \varpi_{1i}^k)\int_0^{l}\psi_{0,j}(x)\psi_{1,i}(x)dx+\boldsymbol{\sigma}_1^k
 \\
 & \quad -\sum_{j=1}^{i_{*}}\Big(a_1\Big(\frac{u_{1j+1}^k-u_{1j-1}^k}{2h_x}\Big)^2
-\overline{\nu}\, p^k\Big)u_{1j}^k\int_0^{l}\psi_{1,j}'(x)\psi_{1,i}'(x)dx,
\end{split}
\end{equation}
and
\begin{equation} \label{eq:d2}
\begin{split}
\mathbf{D}_2 & =  -\frac{k_2}{2}\sum_{j=i_{*}}^{m-1}(\omega_{2j}^{k+1}
+\omega_{2j}^k)\int_0^{l}\psi_{0,j}(x)\psi_{1,i}(x)dx \\
 &\quad -\frac{d_2}{2}\sum_{j=i_{*}}^{m-1}(\varpi_{2i}^k+\varpi_{2i}^k)
 \int_0^{l}\psi_{0,j}(x)\psi_{1,i}(x)dx-\boldsymbol{\sigma}_2^k \\
 &\quad -\sum_{j=i_{*}}^{m-1}\Big(a_2\Big(\frac{u_{2j+1}^k-u_{2j-1}^k}{2h_x}\Big)^2
 -\overline{\nu}\, p^k\Big)u_{2j}^k\int_0^{l}\psi_{1,j}'(x)\psi_{1,i}'(x)dx.
\end{split}
\end{equation}
The substitutions $\boldsymbol{\omega}_{h_t,h_x}^k=(\mathbf{u}_{h_t,h_x}^k)''''$,
and $(\boldsymbol{\varpi}_{h_t,h_x}^k)_t= (\mathbf{v}_{h_t,h_x}^k)''''$
in the equation \eqref{eq:substituion} allow us to obtain the following
intermediate equations:
\begin{equation}
\sum_{j=1}^{i_{*}}\omega_{1j}^k\int_0^{l}\psi_{0,j}(x)\psi_{3,i}(x)dx=
\sum_{j=1}^{i_{*}}u_{1j}^k\int_0^{l}\psi''_{3,j}(x)\psi''_{3,i}(x)dx,  \label{eq:sub1}
\end{equation}
and
\begin{equation}
\sum_{j=1}^{i_{*}}\omega_{2j}^k\int_0^{l}\psi_{0,j}(x)\psi_{3,i}(x)dx
=\sum_{j=1}^{i_{*}}u_{2j}^k\int_0^{l}\psi''_{3,j}(x)\psi''_{3,i}(x)dx.\label{eq:sub2}
\end{equation}
Since the volumes of the linear and the cubic B-splines near each
node are equal, it follows that, for $1\leq j\leq m-1$,
\begin{equation}
\int_0^{l}\psi_{1,j}(x)\, dx=\int_0^{l}\psi_{3,j}(x)\, dx. \label{eq:volume}
\end{equation}
Thus, it follows  from \eqref{eq:fully3}--\eqref{eq:volume} that the linear system
corresponding to equation \eqref{eq:fully3} can be written as follows,
\[
\begin{pmatrix}
\frac{1}{h_t}\mathbf{M}_1(\mathbf{v}_1^{k+1}-\mathbf{v}_1^k)\\[2pt]
\frac{1}{h_t}\mathbf{M}_2(\mathbf{v}_2^{k+1}-\mathbf{v}_2^k)
\end{pmatrix}
=\begin{pmatrix}
\mathbf{D}_1\\
\mathbf{D}_2
\end{pmatrix},
\]
where
\begin{gather*}
\mathbf{D}_1=-\frac{k_1}{2}\mathbf{K}_1(\mathbf{u}_1^{k+1}
+\mathbf{u}_1^k)-\frac{d_1}{2}\mathbf{K}_1(\mathbf{v}_1^{k+1}
+\mathbf{v}_1^k)-\mathbf{L}_1^k\mathbf{u}_1^k+\boldsymbol{\sigma}_1^k,
\\
\mathbf{D}_2=-\frac{k_2}{2}\mathbf{K}_2(\mathbf{u}_2^{k+1}
+\mathbf{u}_2^k)-\frac{d_2}{2}\mathbf{K}_2(\mathbf{v}_2^{k+1}
+\mathbf{v}_2^k)-\mathbf{L}_1\mathbf{u}_2^k-\boldsymbol{\sigma}_2^k.
\end{gather*}
Here,
\begin{gather*}
\mathbf{u}_1^k=(u_{11}^k,u_{12}^k,\dots,u_{1i-1}^k,u_{1i_{*}}^k)^T,\quad
\mathbf{u}_2^k=(u_{2i_{*}}^k,u_{2i_{*}+1}^k,\dots,u_{2m-2}^k,u_{2m-1}^k)^T,
\\
\mathbf{v}_1^k=(v_{11}^k,v_{12}^k,\dots,v_{1i_{*}-1}^k,v_{1i_{*}}^k)^T,\quad
\mathbf{v}_2^k=(v_{2i_{*}}^k,v_{2i_{*}+1}^k,\dots,v_{2m-2}^k,v_{2m-1}^k)^T,
\\
\boldsymbol{\sigma}_1^k=(0,0,\dots,0,\sigma_{h_t,h_x}^k)^T,\quad
\boldsymbol{\sigma}_2^k=(\sigma_{h_t,h_x}^k,0,\dots,0)^T.
\end{gather*}
The symmetric matrices $\mathbf{M}_1$, $\mathbf{M}_2$, $\mathbf{K}_1$, and
$\mathbf{K}_2$ are defined by
\begin{gather*}
\mathbf{M}_1\equiv (M_1)_{ij}=\int_0^{l_{*}}\psi_{1,i}(x)\psi_{1,j}(x),\\
\mathbf{M}_2\equiv (M_2)_{ij}=\int_{l_{*}}^{l}\psi_{1,i}(x)\psi_{1,j}(x)\, dx,\\
\mathbf{K}_1\equiv(K_1)_{ij}=\int_0^{l_{*}}\psi_{3,i}''(x)\psi_{3,j}''(x)\, dx,\\
\mathbf{K}_2\equiv(K_2)_{ij}=\int_{l_{*}}^{l}\psi_{3,i}''(x)\psi_{3,j}''(x)\, dx,
\end{gather*}
and they are banded with three subdiagonals and superdiagonals.
From the extra equation \eqref{eq:fully2}, we find
\begin{equation*}
\frac{1}{h_t}\sum_{j=1}^{i_{*}}(u_{1j}^{k+1}-u_{1j}^k)
 \int_0^{l}\psi_{1,j}(x)\psi_{1,i}(x)dx\\
=\frac{1}{2}\sum_{j=1}^{i_{*}}(v_{1j}^{k+1}+v_{1j}^k)
 \int_0^{l}\psi_{1,j}(x)\psi_{1,i}(x)dx,
%\label{eq:extra1}
\end{equation*}
and
\begin{equation*}
\frac{1}{h_t}\sum_{j=i_{*}}^{m-1}(u_{2j}^{k+1}-u_{2j}^k)
 \int_0^{l}\psi_{1,j}(x)\psi_{1,i}(x)dx
=\frac{1}{2}\sum_{j=i_{*}}^{m-1}(v_{2j}^{k+1}+v_{2j}^k)
\int_0^{l}\psi_{1,j}(x)\psi_{1,i}(x)dx.
%\label{eq:extra2}
\end{equation*}
Now,  we can use \eqref{eq:fully2} to set up the linear system that marches
the solution one time step,
for $\iota=1,\,2$
\begin{equation}
\begin{split}
\mathbf{u}_{\iota}^{k+1}
& =  \Big(\frac{2}{h_t^2}\mathbf{M}_{\iota}+(\frac{k_{\iota}}{2}
+\frac{d_{i}}{h_t})\mathbf{K}_{\iota}\Big)^{-1}\\
&\quad \times
 \Big(\Big(\frac{2}{h_t^2}\mathbf{M}_{\iota}-(\frac{k_{\iota}}{2}
 -\frac{d_{i}}{h_t})\mathbf{K}_{\iota}-\mathbf{L}_{\iota}^k\Big)\mathbf{u}_{\iota}^k
 +\frac{2}{h_t}\mathbf{M}_{\iota}\mathbf{v}_{\iota}^k+(-1)^{\iota+1}
\boldsymbol{\sigma}_{\iota}^k\Big).
 \end{split} \label{eq:linear_sys}
\end{equation}

We note that this linear system does not yet satisfy the normal
compliance contact condition \eqref{eq:NC}. So, we introduce additional notation
to deal with the condition and that will allow us to compute the next step solutions
$\mathbf{u}_{\iota}^{k+1}$ which satisfy \eqref{eq:NC} and \eqref{eq:linear_sys}.

Consider a matrix $\mathbf{A}\in\mathbb{R}^{i_{*}\times(m-1- i_{*})}$.
Then the notation of block matrices is as follows.
If $1\leq i_1\leq i_2\leq i_{*}$ and $1\leq j_1\leq j_2\leq m-1- i_{*}$, i.e.,
we extract $i_1$th row to $i_2$th row
and  $j_1$th column to $j_2$th column, the block matrix
is denoted by $A(i_1:i_2,\, j_1:j_2)\in\mathbb{R}^{(i_2-i_1-1)
\times(j_2- j_1-1)}$.
 We use the following substitutions in the linear system
\eqref{eq:linear_sys}  to simply the matrices involved,
\begin{gather*}
\mathbf{B}_{\iota}  =  \frac{2}{h_t^2}\mathbf{M}_{\iota}+(\frac{k{\iota}}{2}
+\frac{d_{\iota}}{h_t})\mathbf{K}_{\iota},\\
\mathbf{Q}_{\iota}^k  =  \frac{2}{h_t^2}\mathbf{M}_{\iota}-(\frac{k_{\iota}}{2}
-\frac{d_{\iota}}{h_t})\mathbf{K}_{\iota}-\mathbf{L}_{\iota}^k.
\end{gather*}
Note that the matrix $\mathbf{Q}_{\iota}^k$ changes at each
time step $t=t_{k}$. Next,  algebraic manipulations allow us to obtain
the normal compliance in terms of only $u_{1i_{*}}^{k+1}$ and $u_{2i_{*}}^{k+1}$,
indeed, we can write
\begin{gather}
q_1u_{1i_{*}}^{k+1}-b_1^k  =  \lambda (u_{1i_{*}}^{k+1}
-u_{2i_{*}}^{k+1}-g_1)_{+},\label{eq:w1^l+1}\\
q_2u_{2i_{*}}^{k+1}-b_2^k  =  \lambda (u_{2i_{*}}^{k+1}
-u_{1i_{*}}^{k+1}-g_2)_{+},\label{eq:u2^l+1}
\end{gather}
where
\begin{align*}
q_1 & =  (B_1)_{i_{*}i_{*}}-B_1(i_{*},\,1:i_{*}-1)[B_1
(1:i_{*}-1,\,1:i_{*}-1)]^{-1}\\
&\quad \times B_1(1:i_{*}-1,\, i_{*}),\\
b_1^k & =  -B_1(i_{*}\,,1:i_{*}-1)\times[B_1(1:i_{*}-1,\,1:i_{*}-1
)]^{-1}\\
 & \quad \times(Q_1^k(1:i_{*}-1,\,1:i_{*})\mathbf{u}_1^k+\frac{2}{h_t}
 M_1(1:i_{*}-1,\,1:i_{*})\mathbf{v}_1^k)\\
 & \quad +Q_1^k(i_{*},\,1:i_{*})\mathbf{u}_1^k+\frac{2}{h_t}M_1
 (i_{*}\,,1:i_{*})\mathbf{v}_1^k,\\
q_2 & =  (B_2)_{i_{*}i_{*}}-B_2(i_{*},\, i_{*}+1:m-1)\\
&\quad \times [B_1(i_{*}+1:m-1,\, i_{*}+1:m-1)]^{-1}
B_1(i_{*}+1:m-1,\, i_{*}),\\
b_2^k & =  -B_2(i_{*}\,,i_{*}:m-1)\times[B_1(i_{*}:m-1,\,
i_{*}:m-1)]^{-1}\\
&\quad \times \Big(Q_2^k(i_{*}:m-1,\, i_{*}:m-1)\mathbf{u}_2^k+\frac{2}
 {h_t}M_2(i_{*}:m-1,\, i_{*}:m-1)\mathbf{v}_1^k\Big)\\
 &\quad + Q_2^k(i_{*},\, i_{*}:m-1)\mathbf{u}_2^k+\frac{2}{h_t}
 M_2(i_{*}\,,i_{*}:m-1)\mathbf{v}_2^k.
\end{align*}
Once $u_{1i_{*}}^{k+1}$ and $u_{2i_{*}}^{k+1}$ have been computed from
\eqref{eq:w1^l+1} and \eqref{eq:u2^l+1}, the next step solutions
$\mathbf{u}_1^{k+1}$ and $\mathbf{u}_2^{k+1}$ can be computed.
Excluding the last component of $\mathbf{u}_1^{k+1}$ and the first
component of $\mathbf{u}_2^{k+1}$, we let
\begin{gather*}
\mathbf{x}_1^{k+1}=(u_{11}^{k+1},u_{12}^{k+1},\dots,
 u_{1i_{*}-2}^{k+1},u_{1i_{*}-1}^{k+1}) \in\mathbb{R}^{i_{*}-1},
\\
\mathbf{x}_2^{k+1}=(u_{2i_{*}+1}^{k+1},\dots,u_{2i_{*}+2}^{k+1},\dots,
 u_{2m-1}^{k+1}, u_{2m-2}^{k+1})\in\mathbb{R}^{m-i_{*}}.
\end{gather*}
Then
\begin{gather}
\begin{aligned}
\mathbf{x}_1^{k+1} & =  [B_1(1;i_{*}-1,1;i_{*}-1)]^{-1}\Big(Q_1^k
(1;i_{*}-1,1;i_{*})\mathbf{u}_1^k \\
 & \quad +\frac{2}{h_t}M_1(1;i_{*}-1,1;i_{*})\mathbf{v}_1^k
 -Q_1^k(1;i_{*}-1,i_{*})u_{1i_{*}}^{k+1}\Big),
\end{aligned}\label{eq:sol_next1}
\\
\begin{aligned}
\mathbf{x}_2^{k+1} & =  [B_1(i_{*};m,i_{*};m)]^{-1}\Big(Q_2^k
(i_{*}+1;m,i_{*};m)\mathbf{u}_2^k \\
 & \quad +\frac{2}{h_t}M_2(i_{*}+1;m,i_{*};m)\mathbf{v}_2^k-Q_2^k
 (i_{*}+1;m,i_{*})u_{2i_{*}}^{k+1}\Big).
\end{aligned}\label{eq:sol_next2}
\end{gather}
At the final step, the actual approximation of $u_1$ is computed.

The discussion above can be summarized as the following
algorithm.

\begin{algorithm}\label{alg:alg-1}\rm
 Assume that initial data $\mathbf{u}_1^{0}$
and $\mathbf{u}_2^{0}$ are given and $\lambda>0$ is sufficiently
large.

\noindent\textbf{for} $k=1:T/h_t$

 the previous solution $(\mathbf{u}_1^{k-1},\,\mathbf{u}_2^{k-1})^T$ is known

 Compute $(\mathbf{u}_1^k,\,\mathbf{u}_2^k)^T$ using the linear system \eqref{eq:linear_sys}

 Compute $(u_{1i_{*}}^{k+1},\, u_{2i_{*}}^{k+1})$ using \eqref{eq:w1^l+1} and \eqref{eq:u2^l+1}


\noindent\textbf{if } $ u_{1i_{*}}^k-u_{2i_{*}+1}^k<g_1$ \textbf{and}
$u_{2i_{*}+1}^k-u_{1i_{*}}^k<g_2$

 $\sigma^{k-1}\leftarrow0$, $u_{1i_{*}}^k\leftarrow b_1/q_1^{k-1}$,
$u_{2i_{*}+1}^k\leftarrow b_2/q_2^{k-1}$ from \eqref{eq:w1^l+1}--\eqref{eq:u2^l+1}

\noindent\textbf{else if $g_1\leq u_{1i_{*}}^k-u_{2i_{*}+1}^k \leq g_1+1/\lambda$}

 compute $(u_{1i_{*}}^k, u_{2i_{*}}^k)$
from \eqref{eq:w1^l+1} and \eqref{eq:u2^l+1} and then
\[
\sigma^k\leftarrow -\lambda(u_{1i_{*}}^{k+1}-u_{2i_{*}}^{k+1}-g_1)_{+}
\]

\noindent\textbf{else if $g_2\leq u_{2i_{*}+1}^k-u_{1i_{*}}^k\leq g_2+1/\lambda$}
 compute $(u_{1i_{*}}^k,\, u_{2i_{*}}^k)$
from \eqref{eq:w1^l+1} and \eqref{eq:u2^l+1} and then
\[
\sigma^k\leftarrow\lambda
(u_{2i_{*}}^{k+1}-u_{1i_{*}}^{k+1}-g_2)_{+}
\]
\textbf{end if}

 Compute $(\mathbf{u}_1^k,\,\mathbf{u}_2^k)$
by using \eqref{eq:sol_next1} and \eqref{eq:sol_next2}.

 Compute $(\mathbf{v}_1^k,\,\mathbf{v}_2^k)$
from the following identities
\[
\mathbf{v}_1^k=\frac{2}{h_t}(\mathbf{u}_1^k-\mathbf{u}_1^{k-1})
-\mathbf{v}_1^{k-1},\quad\mathbf{v}_2^k=\frac{2}{h_t}(\mathbf{u}_2^k
-\mathbf{u}_2^{k-1})-\mathbf{v}_2^{k-1}
\]
\textbf{end for}
\end{algorithm}

The algorithm was implemented and a typical run took on the
average 0.1274 seconds to compute one time step for the solutions with
$ p=0$ and an average of 0.1259 seconds to compute each time step in the
case of  $ p=895$.
The  numerical analysis of the algorithm is currently under study in \cite{ahn:dfcgts}.


\section{Numerical Results}
\label{sec:Numerical_results}

In this section we present  results of preliminary numerical  simulations obtained by
 using a code based on the algorithm described in the previous section. We show
 two sets of simulations that depict some of the behaviors of the system.
 Following the simulations, we make a few  remarks on the possible types of
 behavior and on the implementation. Our main interest  is in the simulations of  the
 dynamic behavior of the system in two cases: $p=0$, in which there is no buckling,
 and $p=895$, in which there is buckling.  This is a preliminary  study of the transmission
 of vibrations across joints between two Gao beams. As can be seen below,
even the two `simple' cases are not simple and lead to very  complicated behavior.
 Therefore,  a further in-depth study of the numerics has merit in and of itself.



The two simulations are of  two long beams made of steel.
For the sake of simplicity, we assumed that the viscosity is negligible and  set $d_1=d_2=0$.
We assumed that no forces acted on the beams; i.e., $f_1=f_2=0$, and let $\Phi=0$.
Indeed, when the left end of the left beam is clamped to a stationary device, $\phi=0$
and the shift with $\Phi$ is not needed.  Therefore,  the equations of the two  Gao beams  are,
\begin{gather*}
 u_{1tt}+k_1u_{1xxxx}+(\overline{\nu}_1  p -a_1 u_{1x}^2)u_{1xx}=0,\\
 u_{2tt}+k_2u_{2xxxx}+(\overline{\nu}_2  p -a_2 u_{2x}^2)u_{2xx}=0.
\end{gather*}

\begin{table}[h]
\caption{Input data (dimensionless)} \label{tab:Input-data}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
$L_1$ & $L_2$ &  $k_1$ & $k_2$ & $\overline{\nu}_1=\overline{\nu}_2$
& $ \lambda $ & $a_1$ & $a_2$ & $g_1=g_2$ \\
\hline
$0.5$ & $0.5$  & $14.7$ & $7.35$ & $1.30$ & $10^{4}$ & $40$ & $110$ &
$4\times 10^{-7}$ \\
\hline
\end{tabular}
\end{table}
\smallskip

The input data (in dimensionless form) is provided in Table \ref{tab:Input-data}.
We note that the stops are relatively very rigid, the left beam is twice as  stiff
as the right one, the Gao coefficients are large, and the gap is very small.

 In the first simulation the horizontal traction is $p=0$, and so we do not expect
 bucking of the two beams. In the second simulation  the horizontal traction is
 $p=895$ which allows for buckling. Moreover, we  assumed that $l_{*}=0.5$,
 so the two beams are of equal length. The fact that the left beam is stiffer than
 the right one causes asymmetrical behavior that is clearly seen in the simulations.
 We used the normal compliance contact condition with very large stiffness
 ($\lambda=10^{4}$). The case of the Signorini contact condition was deemed
 not realistic and was not simulated.


The initial displacements and velocities of the  beams were, respectively,
\begin{gather*}
u_{01}(x)=2.5\cdot10^{-5}\,x^3,\quad v_{01}(x)=0, \quad 0\leq x \leq 0.5,\\
u_{02}(x)=-2.5\cdot10^{-5}(x-1)^3,\quad v_{02}(x)=0,\quad 0.5 \leq x \leq 1.
\end{gather*}
Both displacements were very small and there was no contact initially,
also, the the vibrations started from rest.

We now describe the numerical experiments.  The size of the subinterval,
i.e., the space step, was $h_x=2\times 10^{-4}$ and the time step was
$h_t=1\times 10^{-4}$.  We computed the solution over $1\times 10^{4}$
time steps, that is over the time interval $0\leq t \leq 1$.
At each time step the $2500\times2500$ linear system, composed of
band matrices,  was computed by using the sparse matrix function in MATLAB,
which is based on the Gaussian elimination.  This allowed us to save
memory by compressing the zero elements in the band matrices.

In the actual computations, about $1275$ seconds of real time
elapsed for the results in the case $p=0$ and over $1260$ seconds in the case
of  $p=895$. Therefore, the average time to compute one step of the numerical
solutions was approximately $0.126$ second.

The system started without contact, and in the case  $p=0$
the first contact was at the upper stop at about $t=0.4546$, and for $p=895$
the first contact was also at the upper stop at about $t=0.4490$. Then,
the oscillations exhibited a slower wave that moved from the joint to the ends,
upon which oscillations with higher  frequency were superimposed.


\begin{figure}[ht]
\subfloat[$p=0$]{
\includegraphics[width=6cm,height=3.5cm]{fig81a} %contatct_p0.eps
}
\subfloat[$p=895$ ]{
\includegraphics[width=6cm,height=3.5cm]{fig81b} % contact_p895.eps
}
\caption{Contact forces for $0\leq t \leq 1$}
\label{cap:contact_case1}
\end{figure}


The contact force, computed from the normal compliance condition
\eqref{eq:w1^l+1} and \eqref{eq:u2^l+1}, is displayed in  Figure \ref{cap:contact_case1}.
The sampling was done at every $41$ time steps. 
It reflects a very interesting feature of both
simulations, namely, up to about $t=0.45$ the amplitude of the oscillations was comparable
to the initial displacements, and then, once first contact was established,
the amplitude switched to about $.04$ and then decreased to about $0.02$,
where it stabilized. We see that once contact started, the left beam's end oscillated
and contacted both stops consecutively, and that the contact stress was negative
at the top stop and positive at the bottom stop. The contact stress amplitude in both
cases shows a very similar trend and magnitude.
However, the oscillations in the case $p=895$ were more orderly.

The displacements behaved in a way that is very different from that of
the related  linear beams  ($a_1=a_2=0$). Indeed, in addition to the
buckling of the right beam,
the oscillations were quite complex, as can be seen in Figure \ref{cap:disp}.

\begin{figure}[ht]
\subfloat[$0<t\leq0.45$]{
\includegraphics[width=6cm,height=3.5cm]{fig82a} % first_left_p895.eps
}
\subfloat[ $0<t\leq0.45$]{
\includegraphics[width=6cm,height=3.5cm]{fig82b} % first_right_p895.eps
}
\subfloat[$0.45<t\leq0.75$]{
\includegraphics[width=6cm,height=3.5cm]{fig82c} % second_left_p895.eps
}
\subfloat[$0.45<t\leq0.75$]{
\includegraphics[width=6cm,height=3.5cm]{fig82d} % second_right_p895.eps
}
\subfloat[ $0.75<t\leq1$]{
\includegraphics[width=6cm,height=3.5cm]{fig82e} % third_left_p895.eps
}
\subfloat[ $0.75<t\leq1$]{
\includegraphics[width=6cm,height=3.5cm]{fig82f} % third_right_p895.eps
}
\caption{Displacements of the beams (left and right)
about three different times, $p=895$}
\label{cap:disp}
\end{figure}


The displacements of the two beams at six different times in the interval  $0<t\leq0.45$ are depicted
in parts (a) and (b) of Figure \ref{cap:disp}. The left beam (which is stiffer)  oscillates about
$u_1=0$, that is its zero equilibrium  state, while  the right beam (which is softer)
oscillates about a buckled state. In both cases small amplitude high frequency
oscillations are superimposed on the main ones. We note that the scale is $10^{-6}$,
which is the same as the scale of the initial displacements.
Next, i(c) and (d) show the displacements of the two beams  at six  different times in the time
interval $0.45<t\leq0.75$. Here, the scale is $10^{-3}$ for the left beam and $10^{-5}$ for the right one.
The left beam is oscillating about the zero equilibrium state with a much higher amplitude, while
the right beam is oscillating about a buckled state.  Again, small amplitude high frequency
oscillations are superimposed on the main ones.
Finally, (e) and (f) depict the displacements of the two beams  at 6 different times
in the interval $0.75<t\leq1$.  The scale is $10^{-3}$ for the left beam that oscillates about
the zero state., while the scale is $10^{-4}$ for the right beam, and it  oscillates about its buckled
state. However,  there are no small amplitude high frequency oscillations that were observed
earlier.

As was noted above, the difference in the beams' behavior is related to the difference
in their stiffness and the Gao coefficient ($a_1$ and $a_2$ in Table \ref{tab:Input-data}).
Concerning the small amplitude high frequency oscillations, at this stage it is
not clear if these are just  `noise,' since they seem to be very regular, therefore,
additional investigation is needed to determine if these are `real' or just numerical noise.


We now describe the frequency spectrum or the `noise' characteristics of the system.
To that end we present the Fast Fourier Transform (FFT)
of the displacements $u_1(0.25, t)$ and $u_2(0.75, t)$ in 
Figs.\,\ref{cap:FFTp0}
and \ref{cap:FFTp895}.  These provide the system behavior in  the frequency domain.

\begin{figure}[ht]
\subfloat[left beam]{
\includegraphics[width=6cm,height=3.5cm]{fig83a} % FFT-L-p0.eps
}
\subfloat[right beam]{
\includegraphics[width=6cm,height=3.5cm]{fig83b} % FFT-R-p0.eps
}
\caption{The frequency spectrum (FFT) of the beams, $p=0$}
\label{cap:FFTp0}  
\end{figure}


\begin{figure}[ht]
\subfloat[left beam]{
\includegraphics[width=6cm,height=3.5cm]{fig84a} % FFT-L-p895.eps
}
\subfloat[ right beam ]{
\includegraphics[width=6cm,height=3.5cm]{fig84b} % FFT-R-p895.eps
}
\caption{The frequency spectrum (FFT) of the beams, $p=895$}
\label{cap:FFTp895}    
\end{figure}

The number of points in the FFT is $f_{s}=512$,
and  the scale has been normalized so the frequency is in between $0$ and $1$,
also the sample was of $200$ signals.  We depict the frequencies in $[-0.5, 0.5]$
to better capture the signal near the origin.

It is seen in the case $p=0$, Figure \ref{cap:FFTp0}, that the main two frequencies lie
near zero, but there are many activated frequencies, and also the spectrum near $0.5$ is quite
dense. It is noted that the system is rather noisy, with many small amplitude frequencies present.

Next, Figure \ref{cap:FFTp895} depicts the FFT of the motion of the two beams
 in the case of
$p=895$. The FFT is similar the one in the first case, but the central peak and
the second
peak are much more pronounced in the right beam, indeed, they are an order of magnitude larger.
On the other hand the other frequencies are distributed somewhat differently.

We note that the presence of contact  between the beams makes the results much more noisy.
Indeed, as is indicated in the FFT, a whole spectrum of frequencies is present, which means
plenty of noise. A further study is warranted to determine how much of the noise is numerical.

As was pointed above, these numerical experiments are preliminary and they open a whole
new direction for research, and more will be done
in our future studies. An interesting set would be the study of the motion when one beam is vibrating
about an upper buckled state and the other one about a lower state.

\section{Conclusions}
\label{sec:conclude}

This work studies the vibrations of two nonlinear Gao beams that are connected via a joint
in which there are two stops at the left end of the right beam with
a clearance between them (see Figs.\,\ref{fig1} and \ref{fig2}). The right end of the left beam is constrained to
move within this clearance. The motivation for the study is our long-term goal of studying
the transmission of vibrations across joints. The interest in using the Gao beam model
arises from the fact that it can naturally predict oscillations about a buckled steady state.


The process was modeled using the system
\eqref{21}--\eqref{210} when contact was modeled with the \emph{normal compliance contact condition,}
which means that the stops were reactive. When the stops were perfectly rigid, contact
was modeled with the \emph{Signorini nonpenetration condition,} and the system consisted
of  \eqref{21}--\eqref{28}, \eqref{210} and \eqref{213}--\eqref{214}.

The existence of the unique weak solution of the problem with the normal compliance
condition was stated in Theorem \ref{t31} and established in Section \ref{sec:nc}. The main issue was
to deal with the cubic term $(u_x)^3$  in the variational formulation. The existence
of the weak solution to the problem with the Signorini condition was stated 
in Theorem \ref{t32}, and
shown in Section \ref{sec:Signorini}. It was based on obtaining a priori estimates on
the solutions of the problem with the normal compliance condition and passing to the limit
when the normal compliance stiffness tends to infinity. As usual, the question of the uniqueness
of the solutions remains unresolved. The existence of weak solutions to the respective problems
without viscosity were obtained in Section \ref{sec:novisc}, again, as the limits of
convergent sequences of solutions with viscosity when the viscosity coefficients vanish.

The mathematical and numerical difficulties associated with the Signorini contact condition
are very challenging, and it seems to us of very little applied relevance, since there are no
perfectly rigid materials or stops. In passing we just mention that, in addition to the weak regularity
of the solutions with this condition, the non-uniqueness seems to be related to the fact
that following contact there is no clear rule about the rebound, which often is masked
by supplementing the so-called \emph{restitution coefficient.}


The numerical algorithm for the computer simulations of the model with the normal
compliance condition was developed in Section \ref{sec:algorithm}. It is a
fully discrete numerical method and uses
a uniform discretization of the time interval $[0, T]$ and a mixture of the
Galerkin approximation and central difference formula to discretize the space
interval $[0, 1]$. The central difference formula is combined with the finite
element method (FEM) when we deal with the nonlinear Gao terms in the
two equations.  The basis functions for the FEM were piecewise linear and piecewise
cubic splines, with special attention to the two elements on both sides of the joint.
The normal compliance contact condition was implemented by using an implicit
discretization. The weak formulation led to a linear system that was sparse and solved using
\eqref{eq:linear_sys}. The convergence of the algorithm is under study in \cite{ahn:dfcgts}.

Finally, the results of the numerical simulations obtained by implementing the algorithm,
were presented in Section \ref{sec:Numerical_results}. We used two sets of simulations for
beams with different Gao coefficients and stiffnesses. In the first set the applied horizontal
traction was $p=0$, so the beams oscillated about the zero equilibrium. The oscillations were
quite interesting and there was superposition of fast frequencies on top of the basic frequency.
The vibrational spectrum, or the noise characteristics in this case, depicted in
Figure \ref{cap:FFTp0},  is seen to be quite noisy.
In the second case of $p=895$,  the right beam was in a buckled state and the left beam was not.
As was found in \cite{adbps:asndb, KPS11}  in the case of one Gao beam, the oscillations
were about the buckled state, as can be seen clearly in Figure \ref{cap:disp} (b, d, f).
The system was also noisy, but the basic frequency was more pronounced.


These results open the way to a thorough investigation of the vibrations transmission across joints,
which, as stated in the introduction, is the long-term goal of this work. However, on the way
we still need to make the implementation faster, and to run numerical experiments to
get a more comprehensive understanding of the vibrations of Gao beams.
Then, we plan to use the driving function $\phi=\phi(t)$ for the investigation.

Another issue of interest is to add friction to the contact condition in the
joint. That is likely to make the model substantially more complex both
mathematically and computationally.


\begin{thebibliography}{99}

\bibitem{Ad}  R. Adams; \emph{Sobolev Spaces}, Academic Press, 1975.

\bibitem{a:topdac} J. Ahn;
\emph{Thick obstacle problems with dynamic adhesive contact},
	Math. Model. Numer. Anal. (M2AN) \textbf{42}(6) (2008),  1021--1045.

\bibitem{a:vtbclf} J. Ahn;
\emph{A viscoelastic Timoshenko beam with Coulomb law of friction},
	 Appl. Math. \& Computation \textbf{218} (2012),  7078--7099.
		
\bibitem{ahn:dfcgts} J. Ahn and Eun-jae Park;
\emph{Dynamic frictionless contact of a	Gao beam with two stops}, preprint.

\bibitem{sa:ebdcn} J. Ahn, and D. E. Stewart;
\emph{Euler--Bernoulli beam with dynamic contact: Discretization, convergence,
	and numerical results}, SIAM J. Numer. Anal. \textbf{43}(4) (2005),
	1455--1480 (electronic).

\bibitem {Kl1} L.-E. Andersson;
\emph{A quasistatic frictional problem with a normal compliance penalization term},
 Nonlinear Analysis  \textbf{37} (1999),  689--705.

\bibitem{adbps:asndb}  K. T. Andrews, M. F. M'Bemgue,  Y. Dumont, J. Purcell,
	and M. Shillor;
\emph{Analysis and Simulations of a Nonlinear Dynamic Beam},
	ZAMP, to appear.

\bibitem{AMbS09}  K. T. Andrews, M. F. M'Bengue, and M. Shillor;
\emph{Vibrations of a nonlinear dynamic beam between two stops},
	Discrete and Continuous Dynamical System (DCDS-B ) \textbf{12}(1)
	(2009),  23--38.

\bibitem{AS06}  K. T. Andrews and M. Shillor;
\emph{Thermomechanical
	behaviour of a damageable beam in contact with two stops}, Applicable
	Analysis \textbf{85}(6-7) (2006), 845--865.

\bibitem{ASW96}  K. T. Andrews, M. Shillor, and S. Wright;
\emph{On the dynamic Vibrations of an elastic beam in frictional contact
with a rigid obstacle},	J. Elasticity, \textbf{42}(1996), 1--30.

\bibitem{ijj:ammumdn} I. Bab{\'{s}}ka, J. Osborn, and J. Pitk\"{a}ranta;
\emph{Anaylsis of mixed methods using mesh dependent norm},
	Math. Computation \textbf{35}(152)(1980), 1039--1062

\bibitem{db:fe} D. Braess;
 \emph{Finite Elements}, Cambridge University Press, 2001.

\bibitem{sl:mtfem} S. C. Brenner and L. R. Scott;
\emph{The Mathematical Theory of Finite
	Element Methods}, Texts in Applied Mathematics \textbf{15}, Springer,  1994.

\bibitem{p:femep}  P. G. Ciarlet;
\emph{The Finite Element Method for Elliptic Problems},
	Applied Mathematics \textbf{40}, SIAM,   2002.

\bibitem{Du03} Y. Dumont;
\emph{Some remarks on a vibro-impact scheme},
	Numerical Algorithms \textbf{33}(2003), 227--240.

\bibitem{Jiri}  C. Eck, J. Jaru\v sek, and M. Krbec;
\emph{Unilateral Contact Problems,}
	Pure and Applied Mathematics \textbf{270}, Chapman \& Hall/CRC Press,  2005.

\bibitem{Gao96}  D. Y. Gao;
\emph{Nonlinear elastic beam theory with application in
	contact problems and variational approaches},
 Mechanics Research	Communications \textbf{23}(1)  (1996),  11--17.

\bibitem{Gao00}  D.Y. Gao;
\emph{Finite deformation beam models and triality theory
	in dynamical post-buckling analysis},
 Intl. J. Non-Linear Mechanics \textbf{35} (2000),	103--131.

\bibitem{Gao06}  D. Y. Gao;
\emph{Duality in distributed-parameter control of nonconvex
	and nonconservative dynamical systems with applications},
 Nonlinear Dynamics and Systems Theory \textbf{6}(3)(2006), 257--279.

\bibitem{GaoR98}  D. Y. Gao  and D. L. Russell;
\emph{An extended beam theory for	smart materials applications:
II Static formation problems}, Appl. Math. Optim.
	\textbf{38}(1)(1998), 69--94.

\bibitem{KO}  N. Kikuchi and J. T. Oden;
\emph{Contact Problems in 	Elasticity}, SIAM,  1988.

\bibitem {Kl2}  A. Klarbring, A. Mikelic, and M. Shillor;
\emph{On friction	problems with normal compliance},
  Nonlinear Analysis  \textbf{13}(8) (1989),	935--955.

\bibitem{KMS}  A. Klarbring, A. Mikelic, and M. Shillor;
\emph{Frictional contact problems with normal compliance},
 Intl. J. Engng. Sci.	\textbf{26}(8) (1988), 811--832.
	
\bibitem{Ken97}  K. L. Kuttler;
\emph{Dynamic friction contact problems for general normal
	and friction laws}, Nonlinear Analysis TMA  \textbf{28}(3)(1997), 559-- 575.

\bibitem{KPSZ01}  K. L. Kuttler, A. Park, M. Shillor, and W. Zhang;
\emph{Unilateral dynamic contact of two beams},
	Math. Comput. Modelling \textbf{34}(2001), 365--384.

\bibitem{KPS11} K. L. Kuttler, J. Purcell, and M. Shillor;
\emph{Analysis and simulations of a contact problem for a nonlinear
	dynamic beam with a crack}, Quarterly J. Mech. Appl. Math.  (2011);
	doi: 10.1093/qjmam/hbr018

\bibitem{KS98ts}  K. L. Kuttler and M. Shillor;
\emph{Vibrations of a beam	between two stops},
  Dynamics of Continuous, Discrete and	Impulsive Systems \textbf{8}(1)(2001), 93--110.

\bibitem{KS02nc}  K. L. Kuttler and M. Shillor;
\emph{Dynamic contact with	normal compliance, wear, and discontinuous
friction coefficient},	SIAM J. Math. Anal.  \textbf{34}(1)(2002), 1--27.	
	
\bibitem{KS04nc} K. L. Kuttler and M. Shillor;
\emph{Regularity of solutions to a dynamic
	frictionless contact problem  with normal compliance}, Nonlinear Analysis
	\textbf{59} (2004), 1063--1075. 	

\bibitem{Lio69}  J. L. Lions;
 \emph{Quelques M\'ethodes de R\'esolution des
	Probl\`emes aux Limites Non Lin\'eaires,} Dunod, Paris, 1969.

\bibitem{MN11}  J. Machalova and H. Netuka;
\emph{Bending of a nonlinear beamreposing on an unilateral foundation},
 Appl. Comput. Mechanics \textbf{5}(2011),  45--54.


\bibitem{MO}  J. A. C. Martins and J. T. Oden;
\emph{A numerical analysis of a class of problems in elastodynamics with friction},
 Comput. Meth. Appl.Mech. Engnr., \textbf{40}(1983), 327--360.
	
\bibitem{MFM08}  M. F. M'Bengue;
\emph{Analysis of a Nonlinear Dynamic Beam
	with Material Damage or Contact,} Ph.D. Dissertation, Oakland University,
	March 2008.
	
\bibitem{MbS08}  M. F. M'Bengue and M. Shillor;
\emph{Regularity result For the problem of vibrations of a nonlinear beam},
	Electron. J. Diff. Eqns., Vol. 2008(2008), No. 27, 1--12.

\bibitem{MoShow83}  F. C. Moon and S. W. Shaw;
\emph{Chaotic vibration of a beam with	nonlinear boundary conditions},
 J. Nonlinear Mech. \textbf{18}(1983),	465--477.

\bibitem{pj:mfem}  P. A. Raviart and M. F. Wheeler;
\emph{A mixed finite elementmethod for 2nd order elliptic problems},
 in \emph{Mathematical Aspects of	the Finite Element Method},
 Lecture Notes in Math., Springer,	Berlin, 1977.

\bibitem{PS98} L. Paoli and M. Schatzman;
\emph{Resonance in impact 	problems},'
 Math. Comput. Modelling, \textbf{28}(4-8)(1998), 385--406.

\bibitem{SG12}  H. A .F. A. Santos  and D. Y.  Gao;
\emph{Canonical dual finite	element method for solving post-buckling problems
of a large deformation elastic beam}, Int. J. Nonlinear Mechanics \textbf{47}(2)(2012),
240--247,doi:10.1016/j.ijnonlinmec.2011.05.012

\bibitem{SST04} M. Shillor,  M. Sofonea, and J. J. Telega;
 \emph{Models and Analysis  of	Quasistatic Contact},
 Lecture Notes in Physics,   \textbf{655}, Springer, Berlin, 2004.

\bibitem{Sim87}  J. Simon;
\emph{Compact sets in the space $L^{p}(	0,T;B) $},
 Ann. Mat. Pura. Appl. \textbf{146}(1987), 65--96.

\bibitem{Wri02}  P. Wriggers;
\emph{Computational Contact Mechanics,} Wiley,  2002.

\end{thebibliography}
\end{document}