\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{algorithm, algorithmic}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2017 (2017), No. 301, pp. 1--42.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2017/301\hfil Debonding of bonded rod-beam system]
{Modeling, analysis and simulations of debonding of bonded rod-beam
system caused by humidity and thermal effects}
\author[K. L. Kuttler, S. Kruk, P. Marcinek, M. Shillor \hfil EJDE-2017/301\hfilneg]
{Ken L. Kuttler, Serge Kruk, Pawel Marcinek, Meir Shillor }
\address{Ken L. Kuttler \newline
Brigham Young University,
Department of Mathematics,
Provo, UT 84602, USA}
\email{klkuttler@gmail.com}
\address{Serge Kruk \newline
Oakland University,
Department of Mathematics and Statistics,
Rochester, MI 48309, USA}
\email{kruk@oakland.edu}
\address{Pawel Marcinek \newline
Oakland University,
Department of Mathematics and Statistics,
Rochester, MI 48309, USA}
\email{pbmarcin@oakland.edu}
\address{Meir Shillor \newline
Oakland University,
Department of Mathematics and Statistics,
Rochester, MI 48309, USA}
\email{shillor@oakland.edu}
\thanks{Submitted October 16, 2017. Published December 6, 2017.}
\subjclass[2010]{74K10, 74F25, 74M99, 35L86, 74H15}
\keywords{Rod-beam system; debonding; single lap joint; differential inclusions;
\hfill\break\indent existence; simulations; spectrum shifts}
\begin{abstract}
This work models, analyses and simulates a one-dimensional process of
debonding of a structure made of two viscoelastic bonded slabs that is
described by a rod-beam system. It is motivated, primarily, by the
degradation of adhesively bonded plates in automotive applications and
studies the effects of the humidity, horizontal and vertical vibrations
and temperature on the debonding process.
The existence of a weak solution to the model is established by using
approximate problems, existence theorems for differential inclusions,
and a fixed point theorem. An implicit finite differences algorithm for
the problem is developed and used to simulate the system dynamics.
It is found that the qualitative behavior of the system correlates
well with experimental results. Moreover, it indicates that using the
shifts in the spectrum, as described by the FFT of one component of the
solution, may be used to measure nondestructively the integrity of the
bonds and their deterioration.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{model}[theorem]{Model}
\allowdisplaybreaks
\section{Introduction} \label{sec1}
This work deals with a mathematical model for the process of debonding of
adhesively bonded plates, motivated primarily by its importance
in single lap joints in automotive applications. Whereas
metallic plates are very often joined by welding, nonmetallic or metallic and
nonmetallic plates are often adhesively joined by a thin layer of glue.
It is known, see e.g.\ \cite{MN17, NS15, NM16} and the many references therein,
that the adhesive strength deteriorates as a result of mechanical vibrations,
humidity and temperature,
as well as other mechanical effects. The main interest in
\cite{MN17, NS15, NM16, SN16}
was in the effects of humidity on adhesively bonded parts in vehicles since
they are used in environments of highly varying humidity and possibly temperature.
They constructed a mathematical model using two beams for the processes,
simulated it and compared to their experimental results.
We note that the need to better understand the debonding process can be found in
other industries, in particular, in the Aerospace applications where bonding
of light plates is essential, see, e.g., \cite{SCNM17} and the references therein.
The interest in this work is four-fold: (i) Construction of a model for the process
of debonding of two slabs. (ii) Showing that the model has a weak solution, and
provide for conditions of its uniqueness. (iii) Developing a numerical algorithm
for the discretized problem and its implementation. (iv) Description of
different simulations that show typical debonding processes and
highlight the dependence of the process on the frequency of the driving
traction and the shift in the system's spectrum as debonding progresses.
The model for the debonding of the slabs has been introduced recently in
\cite{sh16}, and its derivation from a full 3D setting is in progress
in \cite{AKS17}. It is considerably simpler than a 2-plate model, nevertheless,
it includes horizontal and vertical tractions in the adhesive region, which
are essential to the process, since it contains four 1D dynamic equations.
The model uses two dynamic rods' equations for the horizontal shear stress
in the (thin) adhesive layer;
two dynamic beams' equations for the vertical shear in the layer; two parabolic
equations for the temperature and humidity in the layer; and a parabolic inclusion
that describes the dynamics of the bonding field. Thus, the system consists six
partial differential equations and a differential inclusion. Moreover, in the beam
equations, to guarantee that the upper beam is above the lower beam, we also
add a set inclusion term. The use of both beams and rods to obtain vertical and
horizontal tractions seems to be new. However, the main novelty is the description
of the debonding process in this setting.
There exists growing recent literature on modeling and analysis
of systems with adhesion where the evolution of the bonding field is
described by parabolic inclusions. These arise from the fact that the bonding field
$\beta=\beta(x, t)$ that measures the fraction density of active bonds
is assumed to be a `damage variable,' which is required to satisfy
$0\leq \beta\leq 1$.
The reader is referred to the monographs \cite{Fr02, MOS13, SST02, SHS06}
for general models using differential inclusions for adhesion of solids,
and the many references therein.
These, naturally, belong to the Mathematical Theory of Contact Mechanics (MTCM),
which has made considerable progress in the last two decades.
Among the many recent papers on the subject of adhesion processes in
mechanical systems, we just mention
\cite{CSS04, HKSS02, JSS01, NAKS05, R11} and the references therein.
We also note that a different approach to humidity related debonding was taken
in \cite{ MN17, NM16}. There, the adhesive layer was considered as an additional
elastic body with humidity dependent elastic-plastic properties, and the
model was static. Their main purpose was to construct a mathematical model for
the prediction of the breaking of
the adhesive. The adhesive layer was assumed to break when the shear or bending
stresses exceeded a prescribed ceiling, the so-called yield limit. The (static)
model was in the form of coupled system of 4th order ordinary differential equations
and the diffusion equation for the humidity. It was solved numerically, and
some of its predictions were compared to experimental results. Here we deal with
the dynamic setting and we introduce the
damage field, which replaces the adhesive layer that is assumed to be very thin.
Since the model constructed in this work is nonlinear and rather complex and
includes two differential inclusions for the bonding field and the motion of the
beams, we first established the existence of its weak solution. To that end,
we used various results from the theory of
differential equations and inclusions for pseudomonotone operators, see, e.g.,
\cite{ berk96, kut99, MOS13}, approximations, a priori estimates and fixed point
arguments. Moreover, under a restriction on the form of the diffusion coefficients
in the heat and humidity equations, we established the uniqueness
of the weak solution. The proof proceeded by first assuming that the bonding,
temperature and the humidity fields were known, established the existence of a
weak solution to the resulting dynamic system with rods and beams, and
then we used fixed point arguments for the full system.
Next, to gain insight into the model solutions, we constructed an algorithm
for the computational approximations of the solutions based on finite differences,
which was fully implicit in time. The scheme was implemented and a numerical study
indicated that it was stable, robust, and seemed to have quadratic convergence.
The proof of convergence is unresolved, yet.
The numerical code was implemented and a number of simulations conducted.
Here, we present results of a few typical simulations. But first, we compare
the exact frequency spectrum of a rod that is held fixed at one end and free
at the other with the computed spectrum. The latter was obtained by using the
Fast Fourier Transform (FFT), and was found to agree very well with
the theoretically calculated frequencies from a simple formula obtained from Fourier
analysis. Then, we computed the FFT spectrum of the system without debonding.
This is used for comparison purposes in the third and fourth
simulations where the FFT spectrum in the cases when the system was
driven by periodic horizontal tractions with frequencies $f=25, 150$ and $ 350\ Hz$.
These are the main computational results in this work. Next, we depict the results
for debonding, when the traction has frequency $ 350\ Hz$, for different values
of the diffusion coefficient in the equation for the bonding field.
Finally, we perform a comparison of the results when the diffusion coefficient
in the equation for humidity is either constant or depends on the bonding function.
We also conduct a numerical study of the convergence of the algorithm. The numerical
solutions of the bonding field at a fixed time and ten decreasing time steps
are presented and compared. It is seen numerically that the algorithm converges
quadratically.
Our main interest in these computer simulations lies in the questions
of how do the traction frequencies affect the debonding process, and how, in turn,
the process affects the shift in the vibration frequencies of the system as
debonding progresses. This may open a way to assess debonding in real systems
by using noninvasive measuring techniques. We return to these points in the
conclusions section.
The rest of the paper is structured as follows. Section \ref{sec2} presents
the derivation of the `classical' model, Model \ref{model1D}, and states clearly
its underlying assumptions.
The weak or abstract formulation of the model is presented in Section \ref{sec4}
within the setting of the appropriate function spaces, Problem \ref{pbmabstract},
and our main theoretical result is stated in Theorem \ref{thm1}.
Using the mathematical preliminaries of Section \ref{sec3}, the existence
of weak solutions is established in Section \ref{sec5}. A numerical algorithm
for the approximate solutions of the model that is based on finite differences,
and which is fully implicit in time, is developed in Section \ref{sec6}. The results
of the simulations are depicted in Section \ref{sec7}. Section \ref{sec8}
presents numerical evidence that the convergence is almost quadratic.
Finally, a summary of the results, some conclusions and a number of
unresolved questions that arise from this work are provided in the last section,
Section \ref{sec9}.
\section{The model} \label{sec2}
We construct a model for the process of debonding of two thin long bonded
slabs, which are assumed to act as rods and beams, that is caused by humidity
and thermal effects. Its derivation from a 3D model, in the limit of two such
long thin slabs, will appear in \cite{AKS17}. General methods for such derivations
can be found in \cite{AV10, TV96}. We use the notion of rods to describe the
horizontal motion of the two slabs and the notion of beams for the vertical motion,
the bending, since these contribute to the debonding process.
The model consists of a nonlinear coupled system of two equations of motion for
the rods, two for the beams, all in terms of the displacements, a parabolic
differential inclusion for the deterioration in the adhesive strength, and two
parabolic equations for the humidity and temperature.
We note that a simpler model has been recently studied in
\cite{MN17, NS15, NM16, SN16}.
There, the authors used a revised Goland-Reissner method for coupled shear
stress-diffusion,
and also a revised Hart-Smith model with diffusion. However, the approach here is
different as we introduce the bonding field explicitly, and we also deal with
the full diffusion equations for humidity and temperature.
The two slabs occupy the intervals $[0, l_2]$ and $[l_1, 1]$,
($00$, and for the sake of simplicity does not depend on $\theta$ or $
\eta$. The temperature is driven from the ends $x=l_1$ and $x=l_2$,
where it is prescribed, hence at the ends $\theta(l_1, t) =\theta _{L}(t)$ and
$\theta(l_2, t) =\theta_{R}(t)$, and initially $\theta=\theta _0(x)$.
Finally, we describe the diffusion of moisture in the adhesive layer, which
is considered one of the main causes of the adhesive deterioration, \cite{NS15, NM16}.
We use the \emph{humidity function} $\eta =\eta (x,t)$, which measures
the water content per unit length, and assume that the diffusion is affected
by the internal strain and the fraction of the active bonds. Therefore, we model the
diffusion as
\begin{equation}
\eta _t-(D\eta _x)_x=0,\quad l_10$, and for the sake of simplicity does not depend on $\theta$ or $
\eta$, too. Humidity diffusion is also driven from the ends where it is that
of the ambient air around the system, so we assume that at the ends $\eta(l_1, t)
=\eta _{L}(t)$ and $\eta(l_2, t) =\eta _{R}(t)$. Although it is usually assumed
that there is no initial water, for the sake of generality, we allow the initial
condition $\eta=\eta _0(x)$.
\vskip4pt
Finally, we must address the constraint that the left beam must be above the
right one, i.e., $w_1\geq w_2$. To that end we introduce
the subdifferential of the indicator function $I_{[0, \infty)}$, which is a
set-valued function given by
\begin{equation*}
\partial I_{[0,\infty)}(r)
=\begin{cases}
0 & \text{if }r>0, \\
(- \infty, 0] & \text{if } r=0, \\
\emptyset & \text{if } r<0.
\end{cases}
\end{equation*}
Then, we add a term with $\partial I_{[0,\infty)}(w_1-w_2)$ to the
equation of motion \eqref{27a} for the beams, thus changing it into a differential
inclusion. These terms guarantee that $w_1\geq w_2$.
Collecting the equations and the conditions above and writing the equations
of motion in terms of the displacements in the rods and beams, leads to the
following \emph{dynamical model for the debonding of two viscoelastic slabs
caused by humidity, heat and vibrations.}
\begin{model}\label{model1D} \rm
Find the functions $u_1,w_1:[0,l_2]\times [0,T]\to \mathbb{R}$,
$u_2,w_2:[l_1,1]\times [0,T]\to \mathbb{R}$, and
$\beta ,\eta ,\theta :[l_1,l_2]\times [ 0,T]\to \mathbb{R}$, such that,
\begin{gather}
u_{1tt}-c_{r1}^2u_{1xx}-\nu _{r1}u_{1txx}
=\beta K_{r1}(u_2-u_1),\label{26} \\
u_{2tt}-c_{r1}^2u_{2xx}-\nu _{r2}u_{2txx} =-\beta K_{r2}(u_2-u_1),
\label{27} \\
\begin{aligned}
&w_{1tt}+c_{b1}^2w_{1xxxx}+\nu _{b1}w_{1txxxx} \\
&\in -\beta K_{b1}(w_1-w_2) &-&\partial I_{[0,\infty )}(w_1-w_2),
\end{aligned}\label{27a} \\
w_{2tt}+c_{b2}^2w_{2xxxx}+\nu _{b2}w_{2txxxx}
= \beta K_{b2}(w_1-w_2), \label{27bb}
\end{gather}
\begin{gather}
\beta _t-k_{\beta }\beta _{xx}+\Phi (|u_2-u_1|, (w_1- w_2),\beta ,\eta ,\theta )
\in -\partial I_{[0,1]}(\beta ), \label{28} \\
\eta _t-(D(|u_2-u_1|, (w_1-w_2),\beta )\eta _x)_x =0, \label{29} \\
\theta _t-(\kappa (|u_2-u_1|, (w_1-w_2),\beta )\theta _x)_x =0;
\label{210}
\end{gather}
\begin{gather}
u_1(0,t)=0, \quad \sigma _{r1}(l_2,t)=0, \label{211} \\
\sigma _{r2}(l_1,t)=0, \quad\sigma _{r2}(1,t)=p(t), \label{212} \\
w_1(0,t)=w_{1x}(0,t)=0, \quad M_1(l_2,t)=\sigma _{b1}(l_2,t)=0,
\label{212a} \\
M_2(l_1,t)=\sigma _{b2}(l_1,t)=0, \quad M_2(1,t)=0,\quad
\sigma_{b2}(1,t)=q(t), \label{212b}
\end{gather}
\begin{gather}
\beta _x(l_1,t)=0, \quad \beta _x(l_2,t)=0, \label{213} \\
\eta (l_1,t)=\eta _{L}(t), \quad \eta (l_2,t)=\eta _{R}(t), \label{214} \\
\theta (l_1,t)=\theta _{L}(t), \quad \theta (l_2,t)=\theta _{R}(t),
\label{215} \\
u_1( \cdot ,0) =u_{10},\quad u_{1t}( \cdot ,0)=v_{10}^{r}, \quad
u_2( \cdot ,0) =u_{20},\quad u_{2t}( \cdot ,0) =v_{20}^{r}, \label{216} \\
w_1( \cdot ,0) =w_{10},\quad w_{1t}( \cdot ,0)=v_{10}^{b}, \quad
w_2( \cdot ,0) =w_{20},\;w_{2t}( \cdot,0) =v_{20}^{b}, \label{216a} \\
\beta ( \cdot ,0) =\beta _0, \quad
\eta ( \cdot ,0) =\eta _0,\quad \theta (\cdot ,0) =\theta _0. \label{217}
\end{gather}
\end{model}
We note that for equations \eqref{26}--\eqref{29} to make sense, we extend
$\beta$ as zero to the rest of the interval $[0, 1]$.
The initial conditions are given in \eqref{216}--\eqref{217}, with $
0<\beta_0(x)\leq 1$ and $\eta_0(x)\geq 0$ on $[l_1, l_2]$. We assume
no flux of $\beta$ at the end points, \eqref{213}. In practice,
$\eta_{L}(t)=\eta_{R}(t)$ and both are given
functions. However, for the sake of generality, we allow them to be
distinct. Similarly, although typically $\theta_{L}(t)=\theta_{R}(t)$, we
allow them to be different.
We note that the system is coupled via $\beta ,K,\Phi ,D$ and $\kappa $, and
we assume that the $K$s are positive constants, although using functions
with appropriate properties leads to similar results.
Finally, we note that $\Phi$ needs to be obtained from experimental data.
\section{Mathematical preliminaries}
\label{sec3}
We now present some mathematical preliminaries used below. First, assume
that $V\subseteq H=H'\subseteq V'$ is a Gelfand triple,
and for $p>1$ we let
\begin{equation*}
\mathcal{V}=L^{p}( [ 0,T] ;V).
\end{equation*}
In the next section we make the spaces concrete.
Next, we have the following theorem, which is a special case of the one in
\cite{berk96,kut99}. Here, $A$ is single valued that is the case in
the last reference, while a generalization of the result can be found in
\cite{kut99}.
\begin{theorem} \label{19novt1s}
Suppose $f\in \mathcal{V}'$ and $\omega \to u_0\in $ $H$.
Also, let $A( u,t) $ be such that $\hat{A}:\mathcal{V}\to \mathcal{V}'$ is monotone,
hemicontinuous and bounded where $\hat{A}$ is the Nemytskii operator
associated with $A( \cdot ,\cdot ) $. That is, $\hat{A}(
u) ( t) $ is defined as $A( u( t) ,t)$.
Then, there exists a solution $u\in \mathcal{V}$ and $u'\in
\mathcal{V}'$ to the abstract initial value problem
\begin{gather*}
u'+z =f\ \text{ in }\mathcal{V}', \\
u( 0) =u_0\ \text{ in }H, \\
z \in \hat{A}( u) .
\end{gather*}
\end{theorem}
In this paper, we are only interested in the case when $p=2$. Then, the
following theorem is a straightforward consequence of Theorem \ref{19novt1s}.
\begin{theorem}\label{6dect1s}
In addition, suppose that $A( \cdot ,t) $ is
strongly monotone, i.e.,
\begin{equation*}
\| \langle A(u, t) -A(v, t)
,u-v\rangle \| \geq \delta \| u-v\|_V^2-\lambda | u-v| _H^2,
\end{equation*}
and hemicontinuous, single valued and bounded. Let $B( \cdot ,t) $
be bounded and satisfy a Lipschitz condition
\begin{equation*}
\| B( u,t) -B( v,t) \| _{V'}\leq K\| u-v\| _V.
\end{equation*}
Then, there exists a unique solution $v\in \mathcal{V}$,
$v'\in \mathcal{V}'$ to the abstract equation
\begin{gather*}
v'+\hat{A}( v) +B( u) =f, \\
v( 0) = _0\in H, \\
u( t) = u_0+\int_0^tv( s)\, ds.
\end{gather*}
\end{theorem}
\begin{proof}
It is a straightforward consequence of a fixed point theorem. Indeed, we fix
$\hat{v}\in \mathcal{V}$ and let $\hat{u}( t) =u_0+\int_0^t
\hat{v}( s) ds$. Then, by Theorem \ref{19novt1s} and standard
monotonicity arguments, there exists a unique solution $v\in \mathcal{V},
v'\in \mathcal{V}'$ such that
\begin{gather*}
v'+\hat{A}( v) +B( \hat{u}) =f, \\
v( 0) =v_0\in H.
\end{gather*}
If $\hat{v},\hat{v}_1$ are two such given functions, then if $v,v_1$ are
the corresponding solutions to the above initial value problem,
\begin{align*}
&\frac{1}{2}| v( t) -v_1( t) |
_H^2+\delta \int_0^t\| v_1( s) -v(
s) \| _V^2ds \\
&\leq \int_0^t\langle B( \hat{u}) -B( \hat{u}
_1) ,v( s) -v_1( s) \rangle ds+\lambda
\int_0^t| u-v| ^2ds \\
&\leq K_{\delta }\int_0^t\| \hat{u}( s) -\hat{u}
_1( s) \| ^2ds+\frac{\delta }{2}\int_0^t\| v( s) -v_1( s) \| ^2ds
+\lambda \int_0^t| u-v| ^2ds.
\end{align*}
Thus,
\begin{align*}
&\frac{1}{2}| v( t) -v_1( t) |
_H^2+\frac{\delta }{2}\int_0^t\| v_1( s)
-v( s) \| _V^2ds \\
&\leq K_{\delta }\int_0^ts\int_0^{s}\| \hat{v}_1( \tau
) -\hat{v}( \tau ) \| _V^2d\tau ds+\lambda
\int_0^t| u-v| ^2ds.
\end{align*}
Using Gronwall's inequality yields
\begin{equation*}
| v( t) -v_1( t) | _H^2\leq
CK_{\delta }e^{\lambda T}\int_0^ts\int_0^{s}\| \hat{v}
_1( \tau ) -\hat{v}( \tau ) \|
_V^2d\tau ds,
\end{equation*}
and then
\begin{equation*}
\int_0^t\| v_1( s) -v( s) \|_V^2ds
\leq C( T,\delta ,\lambda )
\int_0^t \int_0^{s}\| \hat{v}_1( \tau ) -\hat{v}
( \tau ) \| _V^2d\tau ds.
\end{equation*}
Iterating this inequality shows that the map $\theta :\hat{v} \to v$,
just described, is a contraction map on $\mathcal{V}$ for a sufficiently
high power. It follows that there is a unique fixed point that is the unique
solution to the desired initial value problem.
\end{proof}
To these theorems, we add a useful observation. Let $I$ be a time
interval, $H=L^2( I) $ and let $V$ be a closed subspace of $
H^1( I) $ that contains the relevant test functions. Let $\Psi $
be a bounded nonnegative continuous function, and consider the the operator
$A:V\to V'$ given by
\begin{equation*}
\langle Au,v\rangle =( \Psi ( u) u_x,v_x) _H.
\end{equation*}
Then, $A$ is pseudomonotone. This follows from the compactness of the
embedding of $V$ into $C( I) $ so that weak convergence of $u_{n}$
to $u$ yields uniform convergence. Thus
\begin{equation*}
\liminf_{n\to \infty }\int_{I}\Psi ( u_{n})u_{nx}( u_{nx}-u_x)
\geq \liminf_{n\to \infty}\int_{I}\Psi ( u_{n}) u_x( u_{nx}-u_x) =0
\end{equation*}
From this, the $\liminf $ condition for a pseudomonotone operator follows easily.
Finally, we need the following compactness results due to Simon and Lions.
\begin{theorem}[Simon \cite{Sim}] \label{zz23septt3a} %\label{zz23septt3a}%
Let $q>1$ and let $E\subseteq W\subseteq X$, where the injection map is
continuous from $W$ to $X$ and compact from $E$ to $W$. Let $S_{R}$
be defined by
\begin{align*}
S_{R}=\Big\{&
u:\| u( t) \| _{E}\leq R \text{ for all }t\in [ a,b] , \text{ and}\\
&\| u( s) -u( t) \|_{X}\leq R| t-s| ^{1/q}\Big\}.
\end{align*}
Thus, $S_{R}$ is bounded in $L^{\infty }( 0,T,E) $ and the
functions are uniformly H\"{o}lder continuous into $X$. Then,
$S_{R}\subseteq C( [ a,b] ;W) $ and if $\{ u_{n}\}\subseteq S_{R}$,
there exists a subsequence, $\{ u_{n_k}\} $
that converges to a function $u\in C( [ a,b] ;W) $,
\begin{equation*}
\lim_{k\to \infty }| | u_{n_k}-u\| _{C( [ a,b] ;W) }=0.
\end{equation*}
\end{theorem}
\begin{theorem}[Lions \cite{lio69}] \label{zz23septt3aa}
Let $E\subseteq W\subseteq X$, where the injection map
is continuous from $W$ to $X$ and compact from $E$ to $W$.
Let $p\geq 1$, let $q>1$, and define
\begin{align*}
S_{p,q}\equiv \Big\{&u\in L^{p}( [ a,b] ;E) : \|u( t) -u( s) \| _{X}\leq C|
t-s| ^{1/q}\\
&\text{and }\| u\| _{L^{p}([ a,b] ;E) }\leq R\Big\}.
\end{align*}
Thus, $S_{p,q}$ is bounded in $L^{p}( [ a,b] ;E) $ and
Holder continuous into $X$. Then, $S_{p,q}$ is precompact in
$L^{p}([ a,b] ;W) $ and if $\{ u_{n}\} _{n=1}^{\infty}\subseteq S_{p,q}$,
it has a subsequence $\{ u_{n_k}\} $ which converges in $L^{p}( [ a,b] ;W) $.
\end{theorem}
\section{Abstract formulation of the model} \label{sec4}
We now construct a weak and then an abstract formulation of the model
\eqref{26}--\eqref{217}. We use the usual notation for the various Sobolev function
spaces (see, e.g., \cite{Adams}). To that end, we introduce the additional
Sobolev spaces:
\begin{gather*}
V_1 =\{ \phi \in H^1(0,l_2):\phi ( 0) =0\}, \quad
V_2 =H^1( l_1,1) , \\
V_3 =H_0^1( l_1,l_2) , \quad V_4 =H^1(l_1,l_2),
\end{gather*}
and we denote by $V_j'$ the respective dual spaces, where here
and below $j=1,\dots ,4$. We denote by $H_j$ the space $L^2(I_j) $,
where $I_j$ is the interval corresponding to $V_j$.
Additionally, we let
\begin{gather*}
\mathcal{H}_j=L^2( 0,T;H_j) ,\quad \mathcal{V} _j=L^2( 0,T;V_j) ,\\
\mathcal{V}_j'=L^2( 0,T;V_j'), \quad
\mathcal{Z}=L^2( 0,T;H^1( l_1,l_2) ) .
\end{gather*}
We also use
\begin{equation*}
U_1=\big\{ \phi \in H^2(0,l_2):\phi ( 0) =\phi '( 0) =0\},\quad
U_2=H^2( l_1,1),
\end{equation*}
and the spaces $\mathcal{U}_j$ are defined as above.
We expect to have $u_i\in \mathcal{V}_i$, $w_i\in U_i$, for $i=1,2$,
$\beta\in \mathcal{V}_4$, and $\eta, \theta \in \mathcal{V}_3$.
We begin with the rod equations, proceed formally and define the operators
$A_{r1}:V_1\to V_1'$ and $A_{r2}:V_2\to V_2'$ by
\begin{equation*}
\langle A_{r1}\phi, \psi\rangle =\int_0^{l_2}\phi_x\psi_x\,dx,\quad
\langle A_{r2}\phi, \psi\rangle =\int_{l_1}^1\phi_x\psi_x\,dx.
\end{equation*}
Next, we multiply \eqref{26} by a test function $\psi \in V_1$,
and use integration by parts and the boundary conditions \eqref{211}, set
$v_{1r}=u_{1t}$ and recall that $\beta =0$ outside of $[l_1,l_2]$. In
this manner, we obtain the abstract version of equation \eqref{25} in $
\mathcal{V}_1'$, together with the initial and boundary
conditions,
\begin{gather*}
v_{r1}'+c_{r1}^2A_{r1}u_1+\nu _{r1}A_{r1}v_{r1}=\beta K_{r1}( u_2-u_1) , \\
u_1( t) =u_{10}+\int_0^tv_{r1}( s) ds,\quad v_{r1}( 0) =v_{10}^{r}.
\end{gather*}
Similarly, we multiply \eqref{27} by a test function $\psi \in V_2$, use
integration by parts and the boundary conditions \eqref{212} and setting
$v_{r2}=u_{2t}$ we obtain
\begin{gather*}
v_{r2}'+c_{r2}^2A_{r2}u_2+\nu _{r2}A_{r2}v_{r2}=p( \cdot
) \gamma _2^{\ast }-\beta K_{r2}( u_2-u_1) , \\
u_2( t) =u_{20}+\int_0^tv_2( s) ds,\quad v_2( 0) =v_{20}.
\end{gather*}
Here, $\gamma _2$ denotes the trace of a function in $V_2$ at the right
end ($x=1$) and $\gamma _2^{\ast }$ is its adjoint, so that $p( \cdot
) \gamma _2^{\ast }\psi $ is defined as $p( \cdot ) \psi
( \cdot ,1) $. Thus, we can regard $p( \cdot ) \gamma
_2^{\ast }$ as an element of $\mathcal{V}_2'$. We assume for
the sake of simplicity that the traction $p( t) $ is a continuous
and bounded function. \medskip
Next, we consider the beams' inclusions. We multiply the part corresponding
to the spacial derivatives in \eqref{27a} by $\psi \in U_1$, integrate by
parts and using the boundary conditions, we obtain,
\begin{align*}
&\int_0^{l_2}( c_{b1}^2w_{1xxxx}+\nu _{b1}w_{1txxxx}) \psi
\,dx+\int_0^{l_2}r_1^{\ast }\xi\, \psi dx \\
&=c_{b1}^2\langle A_{b1}w_1,\psi \rangle +\nu
_{b1}\langle A_{b1}w_1',\psi \rangle +(\xi,\psi),
\end{align*}
where $r_j:H_j\to L^2( l_1,l_2) $, for $j=1,2$,
are the maps that set each element of $H_j$ as zero off the interval
$( l_1,l_2) $. Thus,
\begin{equation*}
\int_0^{l_2}r_1^{\ast }\xi \psi dx
=( r_1^{\ast }\xi ,\psi ) _{H_1}\equiv ( \xi ,\psi ) _{L^2(l_1,l_2) },
\end{equation*}
where
$\xi \in \partial I_{[0,\infty )}(w_1-w_2)$.
It follows that $\xi _1\in L^2( l_1,l_2) $. The operator
$A_{b1}:U_1\to U_1'$ is given by
\begin{equation*}
\langle A_{b1}w,\psi \rangle =\int_0^{l_2}w_{xx}\psi _{xx}dx.
\end{equation*}
We, next, multiply the part corresponding to the spacial derivatives in
\eqref{27bb} by $\psi \in U_2$, integrate by parts and use the boundary
conditions and find
\begin{align*}
&\int_{l_1}^1( c_{b2}^2w_{2xxxx}+\nu _{b2}w_{2txxxx}) \psi\,dx \\
&=\langle \gamma ^{\ast }q,\psi \rangle +c_{b2}^2\langle
A_{b2}w_2,\psi \rangle +\nu _{b2}\langle A_{b2}w_2',\psi \rangle ,
\end{align*}
where
\begin{equation*}
\langle A_{b2}w,\psi \rangle =\int_{l_1}^1w_{xx}\psi_{xx}dx,\quad
\langle \gamma ^{\ast }q,\psi \rangle =q(t) \psi ( 1),
\end{equation*}
and when $q$ is continuous, or more generally a function in
$L^2(0,T)$, then $\gamma ^{\ast }q\in \mathcal{U}_2'$. Thus, the
two beam equations can be written abstractly in the form
\begin{equation*}
w_1''+c_{b1}^2A_{b1}w_1+\overline{\nu _{b1}}
A_{b1}w_1'+r_1^{\ast }\xi _1=-r_1^{\ast }\beta
K_{b1}(w_1-w_2)_{+}
\end{equation*}
where $\xi _1\in \partial I_{[0,\infty )}(w_1-w_2)$,
\begin{equation*}
w_2''+c_{b2}^2A_{b2}w_2+\nu _{b2}A_{b2}w_2'
+\gamma ^{\ast }q=r_2^{\ast }\beta K_{b2}(w_1-w_2)_{+}\,.
\end{equation*}
As in the equation for $u_i$ it is convenient to write the equations for
$w_i$ as first order equations for a new variable
\begin{equation*}
w_i( t) =w_{i0}+\int_0^ty_i( s) ds.
\end{equation*}
We write the above two equations in terms of this integral and $y_i=w_{it}$
as
\begin{gather*}
y_1'+c_{b1}^2A_{b1}w_1+\nu _{b1}A_{b1}y_1+r_1^{\ast }\xi
_1 =-r_1^{\ast }\beta K_{b1}(w_1-w_2), \\
y_1( 0) =v_{10}^{b}, \quad \xi _i\in \partial I_{[0,\infty
)}(w_1-w_2),\\
y_2'+c_{b2}^2A_{b2}w_2+\nu _{b2}A_{b2}y_2+\gamma ^{\ast }q
=r_1^{\ast }\beta K_{b2}(w_1-w_2), \\
y_2( 0) =v_{20}^{b}.
\end{gather*}
We next note that the set-inclusion for the debonding \eqref{28} is already
written in an abstract form
\begin{equation*}
\beta '-k_{\beta }\beta _{xx}+\partial I_{[ 0,1] }(
\beta ) \ni -\Phi ,\quad \beta ( 0,x) =\beta _0(x) .
\end{equation*}
Here, the debonding source function is $\Phi=\Phi (|u_2-u_1|,
w_1-w_2,\beta ,\eta ,\theta )$, and for the sake of simplicity we assume
that $\beta _0( x) =\beta _0\in ( 0,1] $ is a
constant, although all the results below are valid if $\beta _0$ is a $
C^1$ function with values in $[0,1]$.
Consider next the humidity diffusion equation \eqref{29}. For the sake of
mathematical convenience, we transform the problem so as to have zero
boundary conditions. To that end, we let $\Lambda(x)=(l_2-x)/(l_2-l_1)$,
and let
$\eta _{B}( x,t) =\Lambda ( x) \eta_{L}( t) +( 1-\Lambda( x) ) \eta _{R}(t) $,
which is a known function. Next, we define a new variable $\hat{\eta}$
that vanishes on the boundary of $[l_1, l_2]$ by
\begin{equation*}
\hat{\eta}( x,t) =\eta ( x,t) -\eta _{B}(x,t) .
\end{equation*}
Then, \eqref{29} becomes
\begin{equation*}
\hat{\eta}_t-( D( \hat{\eta}_x+\eta _{Bx}) )_x=-\eta _{Bt},
\end{equation*}
with the modified initial and boundary conditions
\begin{equation*}
\hat{\eta _0}( x,0) =\eta _0( x) -\eta _{B}(x,0) ,\quad
\hat{\eta}( l_1,t) =\hat{\eta}(l_2,t) =0.
\end{equation*}
We recall that $D=D(|u_2-u_1|, w_1-w_2,\beta )$. Now, let
$\psi \in V_3$, and define the operator $N:V_3\to V_3'$ by
\begin{equation*}
\langle N( u_1,u_2, w_1, w_2, \beta ) \hat{\eta}
,\psi \rangle =\int_{l_1}^{l_2} D \hat{\eta}_x\,\psi_x\,dx.
\end{equation*}
We follow the same procedure and multiply the equation for $\hat{\eta}$ with
$\psi $; use integration by parts; use the conditions
$\psi (l_1)=\psi(l_2)=0$; note that
$\eta _{Bx}=(\eta _{R}(t)-\eta _{L}(t))/(l_2-l_1)$
does not depend on $x$, and assume that $\eta _{L}$ and $\eta _{R}$ are
$C^1$ functions, and so $\eta _{Bx}(t)$ is just a known $C^1([0,T])$
function. We now define
$f=f( u_1,u_2, w_1, w_2, \beta ,\eta_{B}) \in \mathcal{V} _3'$, by
\begin{equation}
\langle f( u_1,u_2, w_1, w_2,\beta ,\eta _{B}) ,\psi
\rangle \equiv- \int_{l_1}^{l_2}D\,\eta _{Bx}\psi
_x\,dx-\int_{l_1}^{l_2}\eta _{Bt}\psi \,dx. \label{f}
\end{equation}
To simplify the notation, we use the symbol $\eta $ instead of $\hat{\eta}$
and $N$ and $f$ instead of $N( u_1,u_2, w_1, w_2, \beta )$
and $f( u_1,u_2, w_1, w_2,\beta ,\eta _{B})$,
respectively, and obtain the following problem: Given $\eta _0\in H_3$,
find $\eta \in \mathcal{V}_3$, such that
\begin{gather*}
\eta '+N \eta =f, \\
\eta ( 0) =\eta _0.
\end{gather*}
We now turn to the heat equation \eqref{210}. As in the case of $\eta $, we
transform the problem so as to have zero boundary conditions, so we use
again $\Lambda (x)=(l_2-x)/(l_2-l_1)$, and let
$\theta _{B}(x,t) =\Lambda( x) \theta _{L}( t) +(
1-\Lambda ( x) ) \theta _{R}( t) $, which is a
given function. Next, we define a new variable $\hat{\theta}$ that vanishes
on the boundary by $\hat{\theta}( x,t) =\theta ( x,t)
-\theta _{B}( x,t) $. Then, \eqref{210} becomes
\begin{equation*}
\hat{\theta}_t-\Big( \kappa ( \hat{\theta}_x+\theta _{Bx})
\Big) _x=-\theta _{Bt},
\end{equation*}
with the modified initial and boundary conditions
\begin{equation*}
\hat{\theta _0}( x,0) =\theta _0( x) -\theta
_{B}( x,0) ,\quad \hat{\theta}( l_1,t) =\hat{\theta}
( l_2,t) =0.
\end{equation*}
Next, we define the operator $M:V_3\to V_3'$ by
\begin{equation*}
\langle M( u_1,u_2, w_1, w_2, \beta ) \hat{\theta}
,\psi \rangle =\int_{l_1}^{l_2} \kappa\, \hat{\theta}_x\psi
_x\,dx,
\end{equation*}
where $\kappa=\kappa (|u_2-u_1|, w_1-w_2, \beta)$. We let $h(
|u_2-u_1|, (w_1-w_2), \beta, \theta_{B} ) \in V_3'$
be given by
\begin{equation}
\langle h( |u_2-u_1|, (w_1-w_2), \beta, \theta_{B} )
,\psi \rangle \equiv -\int_{l_1}^{l_2}\kappa \theta _{Bx}\psi
_x\,dx-\int_{l_1}^{l_2}\theta _{Bt}\psi \,dx. \label{g}
\end{equation}
As with $\eta$, we use $\theta $ instead of $\hat{\theta}$, and similarly,
we write $\theta _0( x) $ instead of
$\theta _0(x) -\theta _{B}(x, 0) $.
Thus, the abstract form of the evolution problem for the (scaled)
temperature is to find $\theta \in \mathcal{V}_3$, when
$\theta _0\in H_3$ is given, such that
\begin{gather*}
\theta '+M \theta = h, \\
\theta (0) =\theta _0.
\end{gather*}
Finally, we also need the operator $L:V_4\to V_4'$
defined as
\begin{equation*}
\langle L\beta ,\psi \rangle =\int_{l_1}^{l_2}\beta _x\psi
_x\,dx.
\end{equation*}
Collecting all the equations and conditions above leads to the following
\emph{abstract formulation} of Model \ref{model1D}, \eqref{26} --\eqref{217}.
\begin{problem}[Abstract Formulation] \label{pbmabstract} \rm
Find seven functions $(u_1,u_2, w_1, w_2$,
$\beta ,\eta ,\theta )$ and $v_i=u_i'$, $y_i=w_i'$,
for $i=1,2$, such that
\begin{equation}
\begin{gathered}
v_1'+c_{r1}^2A_{r1}u_1+\nu _{r1}A_{r1}v_1=\beta
K_{r1}( u_2-u_1) \quad \text{in }\mathcal{V}_1',
\\
v_2'+c_{r2}^2A_{r2}u_2+\nu _{r2}A_{r2}v_2=p( \cdot
) \gamma _2^{\ast }-\beta K_2( u_2-u_1) \quad \text{in } \mathcal{V}_2', \\
y_1'+c^2_{b1}A_{b1}w_1+\nu_{b1} A_{b1}y_1+r^{*}_1\xi
=-\beta K_{b1}(w_1-w_2)\quad \text{in }\mathcal{U}_1', \\
y_2'+c^2_{b2}A_{b2}w_2+\nu_{b2} A_{b2}y+\gamma ^{\ast
}q=\beta K_{b2}(w_1-w_2)\\
\text{in }\mathcal{U}_2',\; r^{\ast }_1\xi \in \partial I_{[ 0,\infty ] }(
w_1-w_2), \\
\beta '+\Phi ( |u_2-u_1|, (w_1-w_2),\beta ,\eta
,\theta ) +L\beta \in -\partial I_{K}( \beta ) \quad
\text{in }\mathcal{V} _4', \\
\eta '+N( u_1,u_2, w_1, w_2,\beta ) \eta
=f( u_1,u_2,w_1, w_2,\beta ,\eta _{B}) \quad \text{in }
\mathcal{V}_3', \\
\theta '+M( u_1,u_2, w_1, w_2,\beta ) \theta
=h( u_1,u_2,w_1, w_2,\beta ,\theta _{B}) \quad \text{in }\mathcal{V} _3'.
\end{gathered} \label{43}
\end{equation}
along with the initial conditions
\begin{equation}
\begin{gathered}
u_i( t) =u_{i0}+\int_0^tv_i( s) ds\quad \text{in }V_i, \\
w_i( t) =w_{i0}+\int_0^ty_i( s) ds,\quad \text{in }U_i, \\
v_i( 0) =v_{i0}\in H_i,\quad y_i( 0) =v_{i0}^{b}
\text{in }H_i, \\
\beta ( 0) =\beta _0\in V_4,\quad \beta _0( x)
\in [ 0,1 ] , \\
\eta ( 0) =\eta _0\in H_3,\quad \theta ( 0) =\theta_0\in H_3.
\end{gathered} \label{44}
\end{equation}
where $i=1,2$.
\end{problem}
Thus, we have a nonlinear abstract evolution system with seven coupled
equations and the relevant initial conditions.
The main theoretical result of this work is the following existence and
(partial) uniqueness theorem.
\begin{theorem}\label{thm1}
Assume that the coefficient functions $\kappa $ and $D$ and the
debonding source function $\Phi $ are each bounded and Lipschitz continuous
with respect to all their variables. Then, there exists a solution to
Problem \ref{pbmabstract}. The solution is unique when $D$ and $\kappa $
depend neither on $\eta $ nor on $\theta $.
\end{theorem}
This theorem guarantees the existence and uniqueness (under the above
restrictions) of a weak solution to Model \ref{model1D}. We note here that
the uniqueness of the solution is not known when $D$ and $\kappa $ depend on
$\eta $ or $\theta $.
\section{Existence of weak solutions}\label{sec5}
We establish the existence of a solution to a variational or weak
formulation of Model \ref{model1D}, \eqref{26}--\eqref{217}.
We prove Theorem \ref{thm1}
in steps. For the sake of convenience, we begin with the beam equations in
the case when $\beta$ is a given function. To that end, we introduce an
approximate system in which we replace the subgradient with a regularization
constructed with the projection operator $\mathcal{P}: \mathbb{R} \to
\mathbb{R}_-$, defined by $\mathcal{P}( r) =r$ for $r<0$ and $
\mathcal{P}( r) =0$ if $r\geq 0$.
We have the following result in the case when $\beta$ is given, all the
assumptions above hold true and the subdifferential is approximated with a
sequence of terms $nr_1^{\ast } \mathcal{P}( w_1-w_2)$,
where $n\in \mathbb{N}$. We note that in this case $w_10$,
\begin{equation*}
| y_1-y_2| _{W_1}
\leq \varepsilon |y_1-y_2| _{W_2}+C_{\varepsilon }|y_1-y_2| _{W_0}
\end{equation*}
Using routine manipulations, and Gronwall's inequality and modifying the
constant $C$ leads to an inequality of the form
\begin{align*}
&| y_1( t) -y_2( t) |_{W_0}^2+| ( w_{1xx}-w_{2xx}) ( t)| _{W_0}^2\\
&+\int_0^t| y_1-y_2|_{W_2}^2ds+n\int_{l_1}^{l_2}\Pi ( w_1( t)
-w_2( t) ) dx\leq C.
\end{align*}
This, in turn, implies that $| w_1( t) -w_2(t) | _{W_2}$ is bounded.
Therefore, a more convenient version of the inequality is
\begin{equation}
\begin{aligned}
&| y_1( t) -y_2( t) |_{W_0}^2+| ( w_1-w_2) ( t)| _{W_2}^2 \\
&+\int_0^t| y_1-y_2|_{W_2}^2ds +n\int_{l_1}^{l_2}\Pi ( w_1( t) -w_2( t)) dx
\leq C.
\end{aligned}\label{16maye4}
\end{equation}
Next, we apply \eqref{14may5} to $y_2$, let $H_2=L_2(l_2,1)$, and
obtain
\begin{align*}
&\frac{1}{2}| y_2( t) |
_{H_2}^2+c_{b}^2| w_{2xx}( t) |
_{H_2}^2-c_{b}^2| w_{20xx}| _{H_2}^2
+\nu_{b}\int_0^t| y_{2xx}| _{H_2}^2ds
+\int_0^tq( s) y_2( s,1) ds \\
&=\int_0^t\int_{l_1}^{l_2}\beta K_{b2}( w_1-w_2)_{+}y_2ds.
\end{align*}
From \eqref{16maye1} it follows that
\begin{equation*}
| y_2( t) | _{H_2}^2+|w_{2xx}( t) | _{H_2}^2+\nu_{b}\int_0^t| y_{2xx}| _{H_2}^2ds
\leq C+\int_0^t| y_2| _{H_2}^2ds,
\end{equation*}
where $C$ does not depend on $n$. Then, adding $\nu
_{b}\int_0^t| y_2| _{H_2}^2ds$ to both sides
and using Gronwall's inequality, we obtain the estimate
\begin{equation}
| y_2( t) | _{H_2}^2+\|
w_2\| _{U_2}^2+\int_0^{T}\| y_2\|
_{U_2}^2ds\leq C, \label{16maye2}
\end{equation}
where $C$ does not depend on $n$. Now, we act on $y_1$ with \eqref{14may3}
, let $H_1=L_2(l_2,1)$, and proceeding similarly, we obtain
\begin{align*}
&\frac{1}{2}| y_1( t) |_{H_1}^2+c_{b}^2| w_{1xx}( t) |
_{H_1}^2-c_{b}^2| w_{10xx}| _{H_1}^2
+\nu_{b}\int_0^t| y_{1xx}| _{H_1}^2ds \\
&+n\int_{l_1}^{l_2}\Pi ( w_1( t) -w_2(t) ) dx
+n\int_0^t\int_{l_1}^{l_2}\mathcal{P}(w_1( s) -w_2( s) ) y_2( s) \,dx\,ds\\
&= K_{b1}\int_0^t\int_{l_1}^{l_2}( -\beta ) (
w_1( s) -w_2( s) ) _{+}( y_2(s) ) \,dx\,ds+\frac{1}{2}| v_{10}^{b}|_{H_1}^2.
\end{align*}
This simplifies to an inequality of the form
\begin{equation}
\begin{aligned}
&| y_1( t) | _{H_1}^2+c_{b}^2|
w_{1xx}| _{H_{11}}^2+\nu _{b}\int_0^t|
y_{1xx}| _{H_1}^2ds+n\int_{l_1}^{l_2}\Pi (w_1( t) -w_2( t) ) dx\\
&+\int_0^t\int_{l_1}^{l_2}( K_{b1}\beta ) (w_1( t) -w_2( t) ) _{+}y_2( s)\,dx\\
&+n\int_0^t\int_{l_1}^{l_2}\mathcal{P}( w_1( s)
-w_2( s) ) y_2( s) \,dx\,ds\\
&\leq c_{bs}^2| w_{10xx}| _{H_1}^2+|
v_{10}^{b}| _{H_1}^2.
\end{aligned} \label{16maye3}
\end{equation}
Now, using the estimate \eqref{16maye2} for $y_2$, we find
\begin{align*}
&\int_0^t\int_{l_1}^{l_2}\mathcal{P}( w_1( s)
-w_2( s) ) y_2( s) \,dx\,ds \\
&\geq -\int_0^t\int_{l_1}^{l_2}\frac{1}{2}( w_1-w_2)
_{-}^2\,dx\,ds-\frac{1}{2}\int_0^t\int_{l_1}^{l_2}|
y_2( s) | ^2\,dx\,ds \\
&\geq -\frac{1}{2}\int_0^t\int_{l_1}^{l_2}( w_1-w_2)
_{-}^2\,dx\,ds-CT \\
&\geq -\int_0^t\int_{l_1}^{l_2}\Pi ( w_1-w_2)
\,dx\,ds-CT.
\end{align*}
Next, \eqref{16maye3} implies
\begin{equation}
\begin{aligned}
&| y_1( t) | _{H_1}^2+c_{b}^2|
w_{1xx}| _{H_{11}}^2+\nu _{b}\int_0^t|
y_{1xx}| _{H_1}^2ds+n\int_{l_1}^{l_2}\Pi (
w_1( t) -w_2( t) ) dx \\
&+\int_0^t\int_{l_1}^{l_2}( K_{b1}\beta ) (
w_1( t) -w_2( t) ) _{+}y_2( s) dx\\
& \leq c_{b}^2| w_{10xx}| _{H_1}^2+|
v_{10}^{b}| _{H_1}^2+n\int_0^t\int_{l_1}^{l_2}\Pi
( w_1-w_2) \,dx\,ds+CT.
\end{aligned}
\end{equation}
Estimate \eqref{16maye4} and the continuity of the embedding of
$H^2(l_1,l_2) $ into $C( [ l_1,l_2] ) $ imply
that the last term on the left-hand side of the inequality above is bounded
below by some constant $-C$, which does not depend on $n$. Therefore, an
estimate of the following form is obtained after adjusting the constant $C$.
\begin{equation}
\begin{aligned}
&| y_1( t) | _{H_1}^2+c_{b}^2|
w_{1xx}| _{H_1}^2+\nu _{b}\int_0^t|
y_{1xx}| _{H_1}^2ds+n\int_{l_1}^{l_2}\Pi (
w_1( t) -w_2( t) ) dx \\
&\leq c_{b}^2| w_{10xx}| _{H_1}^2+|
v_{10}^{b}| _{H_1}^2+n\int_0^t\int_{l_1}^{l_2}\Pi
( w_1-w_2) \,dx\,ds+CT.
\end{aligned}
\end{equation}
Now, applying Gronwall's inequality and adjusting the constant $C$, we
obtain
\begin{equation}
\begin{aligned}
&| y_1( t) | _{H_1}^2+c_{b}^2|
w_{1xx}| _{H_1}^2+\nu _{b}\int_0^t|
y_{1xx}| _{H_1}^2ds+n\int_{l_1}^{l_2}\Pi (
w_1( t) -w_2( t) ) dx \\
&\leq c_{b}^2| w_{10xx}| _{H_1}^2+|
v_{10}^{b}| _{H_1}^2+CT.
\end{aligned}
\end{equation}
Then, in view of the boundary condition in $U_1$, if we define a new norm
$| | | w| | | $ by
$| w_{xx}| _{H_1}$ then it is equivalent to the norm
$\| w\| _{U_1}$. Then, the above implies
\begin{equation}
\begin{aligned}
&| y_1( t) | _{H_1}^2+c_{b}^2\|
w_1\| _{U_1}^2+\nu _{b}\int_0^t\| y_1(
s) \| _{U_1}^2ds+n\int_{l_1}^{l_2}\Pi (w_1( t) -w_2( t) ) dx \\
&\leq c_{b}^2\| w_{10}\| _{U_1}^2+| v_{10}^{b}| _{H_1}^2+CT.
\end{aligned}\label{16maye7}
\end{equation}
We have now the necessary a priori estimates, \eqref{16maye7} and (\eqref
{16maye2}. We restore the index $n$ to the solutions $y_{1n},w_{1n}$ of
equations \eqref{14may3} and $y_{1n},w_{1n}$ of equations \eqref{14may5},
and it follows that there exists a subsequence, still denoted with the
subscript $n$, such that all of the following convergences hold true. We
note that the strong convergence comes from Theorem \ref{zz23septt3a}.
\begin{gather*}
y_{in}\to y_i\quad \text{weak*} \text{ in }L^{\infty }( [ 0,T] ,H_i), \\
w_{in}\to w_i\quad \text{weak* in }L^{\infty }( [ 0,T] ,H_i),
y_{in}\to y_i\quad \text{weakly in }\mathcal{U}_i, \\
w_{in}\to w_i\quad \text{strongly in }C( [ 0,T],C( I_i) ), \\
A_{bi}w_{in}\to A_{bi}w_i\quad \text{weakly in }\mathcal{U}_i',\\
A_{bi}y_{in}\to A_{bi}y_i\quad \text{weakly in }\mathcal{U}'.
\end{gather*}
Furthermore, it follows from the inequality \eqref{16maye1} and the above
strong convergence, that
\begin{equation*}
\int_{l_1}^{l_2}\Pi ( w_1( t) -w_2( t)) dx=0,
\end{equation*}
which shows that $w_1( t) \geq w_2( t) $ for all
$x\in [ l_1,l_2] $. Also, \eqref{14may3} implies that
\begin{equation*}
nr_1^{\ast }\mathcal{P}( w_{1n}-w_{2n}) \to \xi \in
\mathcal{U}_1'\text{ weakly}.
\end{equation*}
Now, for $z\in U_1$, we have
\begin{equation}
\begin{aligned}
&\langle nr_1^{\ast }\mathcal{P}( w_{1n}-w_{2n})
,z-( w_{1n}-w_{2n}) \rangle \\
&= n\int_{l_1}^{l_2}( w_{1n}( x) -w_{2n}( x)
,z( x) -( w_{1n}( x) -w_{2n}( x)) ) dx \\
&\leq n\int_{l_1}^{l_2}\Pi ( z) -n\int_{l_1}^{l_2}\Pi
( w_{1n}( x) -w_{2n}( x) ) .
\end{aligned}\label{6dece1s}
\end{equation}
Since we have strong convergence in $C( I_i) $ uniformly in $t$,
for each $x$,
\begin{equation*}
\liminf_{n\to \infty }n\Pi ( w_{1n}( x)
-w_{2n}( x) ) \geq 0=I_{[ 0,\infty ] }(
w_1( x) -w_2( x) )
\end{equation*}
Then, passing to the limit by taking $\limsup $ of \eqref{6dece1s} and
using the strong convergence, if $z\geq 0$, then
\begin{align*}
\langle \xi ,z-( w_1-w_2) \rangle
&\leq 0-\int_{l_1}^{l_2}\liminf_{n\to \infty }n\Pi (
w_{1n}( x) -w_{2n}( x) ) dx \\
&\leq 0-\int_{l_1}^{l_2}I_{[ 0,\infty ] }( w_1(x) -w_2( x) ) dx.
\end{align*}
We note that if $z<0$ on a set of measure zero, the right-hand side above
would be replaced with $\infty $, so the inequality is preserved. It follows
that
\begin{equation}
r_1^{\ast }\xi \in \partial I_{[ 0,\infty ] }(
w_1-w_2) \quad \text{a.e. for each }t. \label{16maye10}
\end{equation}
Thus, passing to the limit in \eqref{14may3} and \eqref{14may5} and using
the strong convergence to deal with the nonlinear terms, we obtain
\begin{gather}
\begin{gathered}
y_1'+c_{b}^2A_{b1}w_1+\nu _{b}A_{b1}y_1+\xi
= -\beta K_{b1}(w_1-w_2)_{+}, \\
y_1( 0) = v_{10}^{b}\text{ in }H_1,
\end{gathered} \label{16maye8} \\
\begin{gathered}
y_2'+c_{b}^2A_{b2}w_2+\nu _{b}A_{b2}y_2+\gamma ^{\ast }q
=\beta K_{b2}(w_1-w_2)_{+}, \\
y_2( 0) =v_{20}^{b},
\end{gathered} \label{16maye9}
\end{gather}
where for all $z\geq 0,z\in U_1$, \eqref{16maye10} holds. Moreover,
\begin{equation}
w_1( t) \geq w_2( t) \quad \text{on }[ l_1,l_2] . \label{16maye11}
\end{equation}
This establishes the following lemma.
\begin{lemma} \label{lem5.2}
Let $\beta $ be a given function in $L^2( [ 0,T]
,V_4) $ that has values in $[ 0,1] $ and extended by zero
off $[l_1l_2]$ Then there exists a unique solution to the system (\eqref
{16maye8}, \eqref{16maye9} such that \eqref{16maye10} holds and also the
inequality \eqref{16maye11} is satisfied.
\end{lemma}
The uniqueness follows from standard monotonicity arguments, since $\beta$
is given.
We turn to the full model. First, to deal with the set-inclusion in the
equation for $\beta $ in \eqref{43}, we replace it with a sequence of
approximate problems involving penalization of the subdifferential
$\partial I_{[ 0,1] }( \beta ) $. To that end, let
$P(\beta ) $ be equal to 0 on $[ 0,1] $ and be piecewise
linear function on $\mathbb{R}$ with slope equal to 1 for
$\beta \notin [ 0,1]$, and let
$\Psi ( \beta ) =\int_0^{\beta }P(r) dr$. Thus, $\Psi ( \beta ) =0$ on
$[ 0,1] $ and is positive off this interval. Then, we substitute for
$\partial I_{K}( \beta ) $ the penalization operator
$nP( \beta ) $, where $n\in \mathbb{N}$. Also, we let
\begin{equation*}
\tau ( r) =\begin{cases}
1 & \text{if }r>1, \\
r & \text{if }r\in [ 0,1] , \\
0 & \text{if }r<0.
\end{cases}
\end{equation*}
The construction of the approximate problems follows. The existence of
solutions to these problems is straightforward to show and then we obtain
the necessary estimates and pass to the limit $n\to \infty $, which
yields a solution to the model with the differential inclusion.
Consider for each $n\in \mathbb{N}$, the following penalized problem, for
$v_{1n},u_{1n},v_{2n}$,
$u_{2n},y_{1n},w_{1n},y_{2n},w_{2n},\beta _{n},\eta _{n}$ and $\theta _{n}$.
However, to simplify slightly the notation, we
omitted the subscript $n$ from the dependent variables.
\begin{subequations}
\label{eqn36}
\begin{gather}
v_1'+c_{r1}^2A_{r1}u_1+\nu _{r1}A_{r1}v=\tau ( \beta
) K_{r1}( u_2-u_1) , \quad \text{in }\mathcal{V}_1', \\
v_2'+c_{r2}^2A_{r2}u_2+\nu _{r2}A_{r2}v_2
=p( \cdot) \gamma _2^{\ast }-\tau ( \beta ) K_{r2}(u_2-u_1) ,
\quad \text{in }\mathcal{V}_2', \\
\beta '+nP( \beta ) +L\beta =-\Phi , \quad \text{in }
\mathcal{V}_4', \\
\eta '+N\eta =f, \quad \text{ in }\mathcal{V}_3', \\
\theta '+M\theta =h, \quad \text{ in }\mathcal{V}_3', \\
y_1'+c_{b}^2A_{b1}w_1+\nu _{b}A_{b1}y_1+r_1^{\ast }\xi
=-\tau (\beta ) K_{b1}(w_1-w_2)_{+},\quad \text{in }\mathcal{U}_1', \\
y_2'+c_{b}^2A_{b2}w_2+\nu _{b}A_{b2}y+\gamma ^{\ast }q=\tau
( \beta ) K_{b2}(w_1-w_2)_{+}, \quad \text{ in }\mathcal{U}_2', \\
r_1^{\ast }\xi \in \partial I_{[ 0,\infty ] }(
w_1-w_2) , \quad \text{if }z\geq 0,
\end{gather}
\end{subequations}
along with the initial conditions.
\begin{equation}
\begin{gathered}
v_i( 0) =v_{i0},i=1,2,\quad
\beta ( 0) =\beta_0,\;\eta ( 0) =\eta _0,\quad
\theta ( 0) =\theta _0\\
y_i( 0) =v_{i0}^{b},\quad
w_i( t)=w_{i0}+\int_0^ty_i( s) ds,\quad i=1,2.
\end{gathered} \label{eqn37}
\end{equation}
We recall that $\Phi =\Phi ( |u_2-u_1|,(w_1-w_2),\beta ,\eta
,\theta ) $, $N=N( u_1,u_2,w_1,w_2,\beta ) $,
$f=f( u_1,u_2,w_1,w_2,\beta ,\eta _{B}) $,
$M=M(u_1,u_2,w_1,w_2,\beta ) $,
$h=h( u_1,u_2,\beta,\theta _{B}) $.
The continuity of the various functions yields the existence of a solution
to the system \eqref{eqn36}, since it follows from Theorems \ref{19novt1s}
and \ref{6dect1s}. Indeed, it is facilitated by the compactness of the
embedding of $H^1( I) $ into $C( I) $.
For given $u_1,u_2, w_1, w_2, \theta $ and $\eta $ there is at most
one solution to the evolution equation for $\beta$, \eqref{eqn36}c, along
with the initial condition $\beta ( 0) =\beta _0$. Moreover,
the operator $nP( \beta ) +L\beta $ can be considered maximal
monotone since it can be considered the subgradient of a convex proper lower
semicontinuous function and so there exists a unique solution to this
evolution equation such that $\beta ',L\beta \in L^2( [
0,T] ,H_3) $ (where $H_3=L^2( l_1,l_2)$) \cite{bre73}.
Then,
\begin{equation*}
( L\beta ,\beta ') _{H_3}=\int_{l_1}^{l_2}-\beta
_{xx}\beta _t\,dx.
\end{equation*}
Now, we can find a sequence $\{ \beta _k\} $ of functions that
are smooth in $x$, such that $\beta _{kxx}\to \beta _{xx}$ and $
\beta _{kx}\to \beta _x$ in $H_3$. Indeed, we can simply use a
convolution with a suitable mollifier $\{ \phi _k\} $ so that
\begin{equation*}
\beta _k=\int_{\mathbb{R}}\beta ( x-y,t) \phi _k(
y) dy,
\end{equation*}
where here (unlike above) we extend $\beta $ so that
$\beta (x,t) =\beta ( l_1,t) $ if $x\leq l_1$, and
$\beta (x,t)=\beta (l_2,t)$ if $x\geq l_2$ and $\phi _k$ has its support in
$B_k\equiv ( -1/k,1/k) $. Then, if $\psi \in C_0^{\infty}( l_1,l_2) $, we find
\begin{align*}
-\int_{l_1}^{l_2}\beta _k( x,t) \psi _x( x) dx
&=-\int_{l_1}^{l_2}\int_{B_k}\beta ( x-y,t) \phi
_k( y) \psi _x( x) \,dydx \\
&=\int_{l_1}^{l_2}\int_{B_k}\beta _x( x-y,t) \phi
_k( y) \psi ( x) \,dydx,
\end{align*}
for all sufficiently large $k$. Also, $\lim_{k\to \infty }\beta
_{kx}( \cdot ,0) =\beta _{0x}$ in $H_3$. Then, standard
arguments involving convolutions imply the convergence to $\beta _x$ in
$H_3$. The situation is similar for the second derivatives. As to the time
derivative, similar reasoning implies that $\beta _{kt}$ converges to
$\beta_t$. Therefore,
\begin{align*}
\int_0^t( L\beta ,\beta ') _{H_3}\,ds
&=\lim_{k\to \infty }\int_0^t\int_{l_1}^{l_2}-\beta
_{kxx}\beta _{ks}\,dx\,ds \\
&=\lim_{k\to \infty }\int_0^t\frac{d}{ds}\int_{l_1}^{l_2}
\Big( \frac{1}{2}\beta _{kx}^2( x,s) \Big) \,dx\,ds \\
&=\frac{1}{2}\| \beta _x( \cdot ,t) \| _{H_3}^2-
\frac{1}{2}\| \beta _{0x}\| _{H_3}^2.
\end{align*}
Then, it follows that
\begin{align*}
&\int_0^t\| \beta '\|
_{H_3}^2ds+\int_{l_1}^{l_2}n\Psi ( \beta ( t)
) dx+\frac{1}{2}\| \beta _x( \cdot ,t) \|
_{H_3}^2-\frac{1}{2}\| \beta _{0x}\| _{H_3}^2 \\
& \leq C\int_0^t\Phi ^2ds+C\int_0^t\| \beta \|
_{H_3}^2ds+\frac{1}{2}\int_0^t\| \beta '\|
_{H_3}^2ds
\end{align*}
where $\Phi $ depends on various variables as above and is Lipschitz in
$\beta $ and also $C$ is a suitable constant.
Then, it follows from elementary arguments, including Gronwall's
inequality and the Lipschitz condition of $\Phi $ in $\beta $, that
\begin{equation}
\| \beta ( t) \| _{H_3}+\| \beta '\| _{\mathcal{H}_3}
+\int_{l_1}^{l_2}n\Psi ( \beta ( x,t) ) dx
+\| \beta ( t) \|_{V_3}\leq C( \beta _0) . \label{eqn39}
\end{equation}
The other evolution equations are somewhat easier to handle. Thus,
straightforward manipulations yield
\begin{equation}
\begin{aligned}
&\| v_1\| _{\mathcal{V}_1}+\| v_2\|_{ \mathcal{V}_2}
+\| v_1'\| _{\mathcal{V}_1'}+\| v_2'\| _{\mathcal{V}_2'}
+\| y_1'\| _{\mathcal{U}_1'}+\| y_2'\| _{\mathcal{U}_2'} \\
&+\| u_1( t) \| _{V_1}+\| u_2(t) \| _{V_2}
+\| \eta '\| _{ \mathcal{V}_3'}+\| \theta '\| _{ \mathcal{V}_4'}
+\| \eta ( t) \| _{H_3}+\| \theta( t) \| _{H_4}\\
&+\| \eta \| _{\mathcal{V}_3}+\| \theta
\| _{\mathcal{V}_4}+\| y_1\| _{\mathcal{ U}_1}
+\| y_2\| _{\mathcal{U}_2}\leq C,
\end{aligned} \label{12nove3}
\end{equation}
where the constant $C$ is independent of $n$. Now, add to the notation of
each of these functions the subscript $n$.
Then, the estimates and Theorems \ref{zz23septt3a},
\ref{zz23septt3aa} imply the existence of a subsequence
that converges as follows:
\begin{gather*}
\beta _{n}\to \beta \quad \text{weak* in }L^{\infty }([ 0,T] ,H_3) , \\
\beta _{n}'\to \beta '\quad \text{weakly in }\mathcal{H}_4, \\
\beta _{n}\to \beta \quad \text{weak* in }\mathcal{V}_4, \\
\beta _{n}\to \beta \quad \text{strongly in }C( [ 0,T],H_3)
\text{ and pointwise}, \\
v_{ni}\to v_i\quad \text{weakly in }\mathcal{V}_i, \\
u_{ni}\to u_i\quad \text{strongly in }C( [ 0,T],H_i) \text{ and pointwise},\\
v_{ni}'\to v_i'\quad \text{weakly in }\mathcal{V}_i', \\
\eta _{n}\to \eta \quad \text{weakly in }\mathcal{V}_3, \\
\eta _{n}\to \eta\quad \text{weak* in }L^{\infty }( [0,T] ,H_3) , \\
\theta _{n}\to \theta \quad\text{weak* } \text{ in }L^{\infty }([ 0,T] ,H_3) , \\
\theta _{n}\to \theta \quad \text{weakly in }\mathcal{V}_3, \\
\eta _{n}'\to \eta '\quad \text{weakly in }\mathcal{V}_4', \\
\theta _{n}'\to \theta '\quad \text{weakly in }\mathcal{V}_4', \\
\theta _{n}\to \theta\quad \text{strongly in }\mathcal{H}_3\text{ and pointwise}, \\
\eta _{n}\to \eta \quad \text{strongly in }\mathcal{H}_3\text{ and pointwise}, \\
A_iu_{ni}\to A_iu_i\quad \text{weakly in }\mathcal{V}_i'\\
A_iv_{ni}\to Av_i\quad \text{weakly in }\mathcal{V}_i' \\
\Phi (u_{1n},\dots,\theta _{n}) \to \Phi (|u_2-u_1|, (w_1-w_2),\beta ,\eta ,\theta )
\quad \text{strongly in } \mathcal{H}_3, \\
N( u_{1n},\dots,\beta _{n}) \eta _{n}\to N(u_1,u_2,w_1, w_2, \beta) \eta \quad
\text{weakly in } \mathcal{V} _3' \\
M( u_{1n},\dots,\beta _{n}) \theta _{n}\to M(
u_1,u_2,w_1, w_2,\beta ) \theta \quad \text{ weakly in }\mathcal{V} _4', \\
nP( \beta _{n}) \to \xi \quad \text{weakly in }\mathcal{V}_4',\\
y_{ni}'\to y_i'\quad \text{weakly in }\mathcal{U}_i', \\
y_{ni}\to y_i\quad \text{weakly in }\mathcal{U}_i, \\
w_{ni}\to w_i\quad \text{ weak* in }L^{\infty }(0,T;U_i'), \\
w_{ni}\to w_i\quad \text{strongly in }C( [ 0,T];C( I_i) ),
\end{gather*}
where $I_i$ is either $[ 0,l_2] $ or $[ l_1,1] $. Furthermore,
\begin{gather*}
A_{bi}y_{ni}\to A_{bi}y_i\quad\text{weakly in } \mathcal{U}_i', \\
A_{bi}w_{ni}\to A_{bi}w_i\quad\text{weakly in } \mathcal{U}_i',\\
\tau ( \beta _{n}) K_{b2}(w_{n1}-w_{n2})_{+}\to \tau
( \beta ) K_{b2}(w_1-w_2)_{+}\quad \text{strongly in }L^2(0,T;H_2), \\
\tau ( \beta _{n}) K_{b1}(w_{n1}-w_{n2})_{+}\to \tau
( \beta ) K_{b1}(w_1-w_2)_{+}\quad \text{strongly in }L^2(0,T;H_1),
\end{gather*}
Also, for each $t$, the above convergences show that there is a constant $C$,
independent of $n$, such that
\begin{equation*}
n\int_{l_1}^{l_2}\Psi ( \beta _{n}( t) ) dx\leq C.
\end{equation*}
Using Fatou's lemma yields
\begin{equation*}
\int_{l_1}^{l_2}\Psi ( \beta ( x,t) ) dx=0.
\end{equation*}
Hence, $\beta (x,t) \in [ 0,1] $ for a.e. $x$.
We turn to the limit in the equation for $\beta_{n}$,
\[
\beta _{n}'+nP( \beta _{n}) +L\beta _{n}=-\Phi (
|u_{2n}-u_{1n}|,(w_{1n}-w_{2n}), \beta _{n},\eta _{n},\theta _{n}) \quad
\text{in }\mathcal{V}_3'.
\]
Using monotonicity considerations and the fact that
$P( \beta ) =0 $ lead to
\begin{align*}
&\langle \beta ',\beta _{n}-\beta \rangle +\langle
nP( \beta _{n}) -nP( \beta ) ,\beta _{n}-\beta
\rangle +\langle L\beta _{n},\beta _{n}-\beta \rangle \\
&\leq ( -\Phi (|u_{n2}-u_{n1}|, (w_{1n}-w_{2n}), \beta _{n},\eta
_{n},\theta _{n}) ,\beta _{n}-\beta ) _{\mathcal{H}_3},
\end{align*}
and so
\begin{equation*}
\langle \beta ',\beta _{n}-\beta \rangle _{\mathcal{V}_4}
+\langle L\beta _{n},\beta _{n}-\beta \rangle _{\mathcal{V}_4}
\leq ( -\Phi (\dots) ,\beta _{n}-\beta ) _{\mathcal{H}_3}.
\end{equation*}
Thus, $\limsup_{n\to \infty }\langle L\beta _{n},\beta_{n}-\beta \rangle \leq 0$
and since $L$ is monotone, hemicontinuous
and bounded as a map from $V_4$ to $V_4'$ we obtain that for
all $\delta \in \mathcal{V}_4$,
\begin{equation*}
\liminf_{n\to \infty }\langle L\beta _{n},\beta _{n}-\delta
\rangle \geq \langle L\beta ,\beta -\delta \rangle .
\end{equation*}
Now, let $\delta ( x,t) \in [ 0,1] $, with
$\delta \in \mathcal{V}_4$, then,
\begin{align*}
&\langle \beta _{n}',\beta _{n}-\delta \rangle
+\langle nP( \beta _{n}) ,\beta _{n}-\delta \rangle
+\langle L\beta _{n},\beta _{n}-\delta \rangle \\
& =( -\Phi ( \dots) ,\beta _{n}-\delta ) _{
\mathcal{H}_3}.
\end{align*}
Now, monotonicity considerations lead to
\begin{align*}
&\langle \beta ',\beta _{n}-\beta \rangle +\langle
\beta _{n}',\beta -\delta \rangle +\langle nP(
\beta _{n}) -nP( \delta ) ,\beta _{n}-\delta \rangle
+\langle L\beta _{n},\beta _{n}-\delta \rangle \\
&\leq ( -\Phi ( \dots) ,\beta _{n}-\delta ) _{\mathcal{H}_3},
\end{align*}
and
\begin{equation*}
\langle \beta ',\beta _{n}-\beta \rangle +\langle
\beta _{n}',\beta -\delta \rangle +\langle L\beta
_{n},\beta _{n}-\delta \rangle \leq ( -\Phi ( \dots)
,\beta _{n}-\delta ) _{\mathcal{H}_3}.
\end{equation*}
Passing to the $\liminf_{n\to \infty} $ we obtain
\begin{equation*}
\langle \beta ',\beta -\delta \rangle +\langle
L\beta ,\beta -\delta \rangle \leq ( -\Phi ( |u_2-u_1|,
(w_1-w_2), \beta ,\eta ,\theta ) , \beta -\delta ) _{
\mathcal{H}_3}.
\end{equation*}
Hence, for the choice of $\delta $ as above,
\begin{equation*}
\langle \beta '+L\beta +\Phi
(|u_2-u_1|,(w_1-w_2), \beta ,\eta,\theta ) ,\delta -\beta \rangle \geq 0,
\end{equation*}
and so, for arbitrary $\delta \in \mathcal{V}_4$,
\begin{equation*}
\langle \beta '+L\beta +\Phi
(|u_2-u_1|,(w_1-w_2), \beta ,\eta ,\theta ) ,\delta
-\beta \rangle _{\mathcal{V}_4}\geq I_{K}( \beta )
-I_{K}( \delta ) .
\end{equation*}
Therefore,
\begin{equation*}
\beta '+L\beta +\Phi ( u_2-u_1,\beta ,\eta ,\theta )
\in -\partial I_{K}( \beta ) .
\end{equation*}
The above convergences are also sufficient to pass to a limit in all the
other equations. Since $\beta ( x,t) \in [ 0,1] $, it
follows that $\tau $ is irrelevant. This yields the existence of a solution
to \eqref{43} and \eqref{44}. This concludes the existence part of
Theorem \ref{thm1}.
\begin{proof} (Uniqueness). It remains to verify the uniqueness part of
Theorem \ref{thm1}. By assumption, the functions $D$ and $\kappa$ are Lipschitz
continuous and independent of $\eta$ and $\theta$. We consider the case of $
\eta$ and $D$, since the case with $\theta$ and $\kappa$ is very similar.
Suppose $\eta ,\hat{\eta}$ are two such solutions with the corresponding
dependent variables without and with hats. To simplify the notations, we let
\begin{gather*}
N=N( u_1,u_2, w_1, w_2, \beta),\quad \hat{N}
= N(\hat{u}_1,\hat{u}_2,\hat{w}_1, \hat{w}_2, \hat{\beta}), \\
D=D( u_1,u_2, w_1, w_2, \beta),\quad \hat{D}= D(
\hat{u}_1,\hat{u}_2,\hat{w}_1, \hat{w}_2, \hat{\beta}).
\end{gather*}
Then
\begin{align*}
&\langle N \eta - \hat{N} \hat{\eta},\eta -\hat{\eta}\rangle = \langle N
\eta -N\hat{\eta}, \eta -\hat{\eta}\rangle +\langle ( N -\hat{N}))
\hat{\eta},\eta -\hat{\eta}\rangle \\
&\geq \delta \| \eta -\hat{\eta}\| _{V_3}^2-| \langle ( N -\hat{
N}) \hat{\eta},\eta -\hat{\eta}\rangle | .
\end{align*}
It follows from the equation for $\eta $ that there is a constant $C$ such
that $\| \eta \| _{\mathcal{V}}$ $\leq C$ for $\eta $.
Thus, the above term on the right-hand side is dominated by
\begin{align*}
&\int_{l_1}^{l_2}| D -\hat{D} | | \hat{\eta}_x|
| ( \eta -\hat{\eta})_x | dx \\
&\leq C_{\delta }\int_{l_1}^{l_2}| D-\hat{D} |
^2| \hat{\eta} _x| ^2dx+\frac{\delta }{2}
\int_{l_1}^{l_2}| ( \eta -\hat{\eta})_x
| ^2dx.
\end{align*}
Hence,
\begin{align*}
&\langle N \eta -\hat{N} \hat{\eta},\eta -\hat{\eta}\rangle \\
&\geq \frac{\delta }{2}\| \eta -\hat{\eta}\|
_{V_3}^2-C_{\delta }\| D -\hat{D} \| _{L^{\infty
}( l_2,l_1) }^2\| \hat{\eta}\| _{V_4}^2 \\
&\geq \frac{\delta }{2}\| \eta -\hat{\eta}\|
_{V_3}^2-C_{\delta }Lip( D) ^2( \| u_1-\hat{u }_1\| _{E}^2+\| u_2-\hat{u}_2\|
_{E}^2) \| \hat{\eta}\| _{V_4}^2,
\end{align*}
where used the facts that $V_3$ embeds continuously into $E$, which embeds
continuously into $C( [ l_1,l_2] ) $. This and a
similar inequality for $\theta $ allow us to use standard arguments and show
that the solution is unique in this special case of Lipschitz continuity and
lack of dependence on $\theta$ and $\eta $.
\end{proof}
When either $D$ or $\kappa $ depend on either $\theta$ or $\eta$, the uniqueness
of the solutions is not known.
\section{Numerical algorithm}\label{sec6}
We present in this section a fully implicit finite difference algorithm for the
approximate solutions of the model. In the next section, we describe its
implementation and depict three representative examples and related numerical results.
For the sake of simplicity, we do not take into account thermal effects,
since in our setting they are similar to those of humidity.
This simplifies slightly the presentation of the algorithm and the
simulations' results.
For the sake of generality, we used the humidity diffusion coefficient
\[
D=d+ d_\beta(1-\beta),
\]
where $d$ and $d_\beta$ are two positive constants. In this way,
the water diffusion in the glue layer increases as the bonding decreases.
It may describe voids that form in the adhesive as it deteriorates.
Moreover, for the sake of simplicity, we assumed that the effects of viscosity
are negligible, so that the slabs were assumed to be purely elastic.
\subsection{The scheme}
We begin with the time and space discretization of the domain.
We divide the time interval $[0,T]$ into $N+1$ mesh points $t_{n}$, where here
and below,
$n=1,2,\dots, N+1$ with uniform time step $\Delta t =T/N$, and then
$t_{n}=(n-1)\Delta t$.
The spatial domain $[0, 1]$ of each one of the slabs is discretized by the
equidistant mesh points $x_j$, where here and below $j=1,2,\dots S+1$, so that
$\Delta x =1/S$ and $x_k = (k-1) \Delta x$. Moreover, the discretization is
such that the right end
of the first slab is the mesh point $l_2=x_{S+1}$
and the left end of the second slab is the mesh point $l_1=x_{j_2}$. Thus, the
nodes $x_1, \dots ,x_{S+1}$ are in the first slab and the nodes
$x_{j_2}, \dots ,x_{S+ j_2}$ in the second one. The common nodes
\[
x_{j_2},\dots, x_{S+1},
\]
are in the region where the slabs are adhesively bonded.
In our numerical experiment it was
set so the adhesive region occupied approximately $40 \%$ of the first slab,
and both were of the same
length. Therefore, we have the mesh $(x_j, t_n)$ in the domain
$[0,1+l_1]\times[0,T]$.
We discretize a function $\phi(x,t)$ defined on $[0,1+l_1]\times[0,T]$
by using its nodal values
\[
\phi_j^n=f(x_j,t_n).
\]
We use the central difference temporal and spacial discretization of first and
second-order derivatives as
\begin{gather*}
\phi_t(x_j, t_{n+1}) \approx \frac{\phi_j^{n+1}- \phi_j^{n-1}}{(2\Delta t)},\quad
\phi_{tt}(x_j, t_{n+1}) \approx \frac{\phi_j^{n+1}
-2\phi_j^n+\phi_j^{n-1}}{(\Delta t)^2}, \\
\phi_{xx}(x_j, t_{n+1}) \approx \frac{\phi_{j+1}^{n+1}-2\phi_j^{n+1}+
\phi_{j-1}^{n+1}}{(\Delta x)^2},
\end{gather*}
and the fourth order derivatives we approximate as follows:
\[
\phi_{xxxx}(x_j, t_{n}) \approx \frac{w_{1j-2}^{k+1}-4w_{1j-1}^{k+1}
+6w_{1j}^{k+1}-4w_{1j+1}^{k+1}+w_{1j+2}^{k+1}}{(\Delta x)^4}.
\]
Since in most the simulations, at this stage, we assume that $D$ is a positive
constant, we do not present its discretization for the sake of simplifying
slightly the matrices below.
Using these expressions in the equations \eqref{26}--\eqref{210}, and the conditions
\eqref{212a}--\eqref{217}, and letting
$j=1,\dots, S+1$ for $u_1$ and $w_1$, $j=1,\dots, S+1$ for $u_2$
and $w_2$, and $j=j_2, \dots, S-j_2+2$ for $\beta$ and $\eta$, leads to the
following algebraic system of equations:
\begin{gather}
\begin{aligned}
&\frac{u_{1j}^{n+1}-2u_{1j}^n+u_{1j}^{n-1}}{(\Delta t)^2}
- c_{r1}^2 \frac{u_{1j+1}^{n+1}
-2u_{1j}^{n+1}+u_{1j-1}^{n+1}}{(\Delta x)^2} \\
&= \beta_{j-j_2+1}^{n} K_1(u_{2j-j_2+1}^{n}-u_{1j}^{n})\chi_j,
\end{aligned} \label{61}
\\
\begin{aligned}
&\frac{u_{2j}^{n+1}-2u_{2j}^{n}+u_{2j}^{n-1}}{(\Delta t)^2}
- c_{r2}^2 \frac{u_{2j+1}^{n+1} -2u_{2j}^{n+1}+u_{2j-1}^{n+1}}{(\Delta x)^2} \\
&= -\beta_j^{n} K_2(u_{2j}^{n}-u_{1j+j_2-1}^{n}) \chi_j,
\end{aligned} \label{62}
\\
\begin{aligned}
&\frac{w_{1j}^{n+1}-2w_{1j}^n+w_{1j}^{n-1}}{(\Delta t)^2}
+ c_{b1}^2\frac{w_{1j-2}^{n+1}
-4w_{1j-1}^{n+1}+6w_{1j}^{n+1}-4w_{1j+1}^{n+1}+w_{1j+2}^{n+1}}{(\Delta x)^4}\\
&=-\beta_{j-j_2+1}^{n}(w_{1j}^{n}-w_{2j-j_2+1}^{n})\chi_j,
\end{aligned} \label{63}
\\
\begin{aligned}
&\frac{w_{2j}^{n+1}-2w_{1j}^n+w_{1j}^{n-1}}{(\Delta t)^2}
+c_{b2}^2\frac{w_{2j-2}^{n+1}-4w_{2j-1}^{n+1}+6w_{2j}^{n+1}
-4w_{2j+1}^{n+1}+w_{2j+2}^{n+1}}{(\Delta x)^4} \\
&=\beta_j^{n}(w_{1j+j_2-1}^{n}-w_{2j}^{n})\chi_j,
\end{aligned}\label{64}
\\
\begin{aligned}
&\frac{w_{2j}^{n+1}-2w_{1j}^n+w_{1j}^{n-1}}{(\Delta t)^2}
+c_{b2}^2\frac{w_{2j-2}^{n+1}-4w_{2j-1}^{n+1}+6w_{2j}^{n+1}
-4w_{2j+1}^{n+1}+w_{2j+2}^{n+1}}{(\Delta x)^4} \\
&=\beta_j^{n}(w_{1j+j_2-1}^{n}-w_{2j}^{n})\chi_j,
\end{aligned} \label{64b}
\\
\begin{aligned}
& \frac{\beta_j^{n+1}-\beta_j^{n-1}}{2\Delta t}
- k_{\beta}\frac{\beta_{j+1}^{n+1}- 2\beta_j^{n+1}
+\beta_{j-1}^{n+1}}{(\Delta x)^2} \\
&+ \Phi(u_{1j+j_2-1}^{n+1}, u_{2j}^{n+1}, w_{1j+j_2-1}^{n+1}, w_{2j}^{n+1},
\beta_j^{n+1}, \eta_j^{n}) =0,
\end{aligned} \label{65}
\\[4pt]
\frac{\eta_j^{n+1}-\eta_j^{n-1}}{2\Delta t}
-D \frac{\eta_{j+1}^{n+1}-2\eta_j^{n+1}+\eta_{j-1}^{n+1}}
{(\Delta x)^2} =0, \label{66}
\end{gather}
where $\chi_j$ is the characteristic function of the discretized adhesive
region $[x_{j_2}, x_{S+1}]$.
This is equivalent to extending $\beta$ as zero outside this interval and
similarly for $\eta$.
The initial conditions are:
\begin{equation} \label{67}
\begin{gathered}
u_{1j}^1=0, \quad u_{1j}^2-u_{1j}^{0}=0, \\
u_{2j}^1=0, \quad u_{2j}^2-u_{2j}^{0}=0, \\
w_{1j}^1=0, \quad w_{1j}^2-w_{1j}^{0}=0, \\
w_{2j}^1=0, \quad w_{2j}^2-w_{2j}^{0}=0, \\
\beta_j^1=1,\quad \eta_j^1=0.
\end{gathered}
\end{equation}
The boundary conditions are:
\begin{equation} \label{eqn68} %\label{eqn69}
\begin{gathered}
u_{1S+2}^ n-u_{1S}^n=0, \quad
u_{2j_2+1}^{n}-u_{2j_2-1}^{n}=0, \\
\frac{u_{2S+1}^{n}- u_{2S-1}^{n}}{2\Delta x}=p^n, \quad
w_{12}^{n}- w_{10}^{n}=0, \quad
w_{1S+2}^{n}-2w_{1S+1}^{n}+2_{1S}^{n}=0, \\
-w_{1S-1}^{n}+2w_{1S}^{n}-2w_{1S+2}^{n}+w_{1S+3}^{n}=0, \quad
w_{2j_2+1}^{n}-2w_{2j_2}^{n}+w_{2j_2-1}^{n}=0, \\
-w_{2j_2-2}^{n}+2w_{2j_2-1}^{n}-2w_{2j_2+1}^{n}+w_{2j_2+2}^{n}=0, \quad
w_{2S+1}^{n}-2w_{2S}^{n}+w_{2S-1}^{n}=0, \\
\frac{-w_{2j_2-2}^{n}+2w_{2j_2-1}^{n}-2w_{2j_2+1}^{n}+w_{2j_2+2}^{n}}
{\Delta x^3}=q^n, \\
\beta_2^n-\beta_1^n=0, \quad \beta_{S-j_2+2}^n-\beta_{S-j_2+1}^n=0, \\
\eta_1^n=\eta_{L}^n,\quad \eta_{S-j_2+2}^n=\eta_{R}^n, \\
u_{1j}^1=u_{10},\quad \frac{u_{1j}^2- u_{1j}^{0}}{2\Delta t}=v_{10}, \\
u_{2j}^1=u_{20}, \quad \frac{u_{2j}^2- u_{2j}^{0}}{2\Delta t}=v_{20}, \\
w_{1j}^1=w_{10},\quad \frac{w_{1j}^2- w_{1j}^{0}}{2\Delta t}=v_{10}^b, \\
w_{2j}^1=w_{10},\quad \frac{w_{2j}^2- w_{2j}^{0}}{2\Delta t}=v_{20}^b, \\
\beta_j^1=\beta_0,\quad \eta_j^1=\eta_0.
\end{gathered}
\end{equation}
\begin{remark} \rm
We point out that the right-hand sides of \eqref{61}--\eqref{65},
which are nonlinear, are
kept frozen at step $n$, in this way the system is linearized.
\end{remark}
Rearranging the system of equations so that all the variables at time $n+1$
are on the left-hand side, while on right-hand side are all variables at
time $n$ and $n-1$, and introducing the notation
\begin{gather*}
C_{r1}=\Big( c_{r1} \frac{\Delta t}{\Delta x} \Big)^2,\quad
C_{r2}=\Big( c_{r2} \frac{\Delta t}{\Delta x} \Big)^2, \\
C_{b1}=\Big( c_{b1} \frac{\Delta t}{\Delta x^2} \Big)^2, \quad
C_{b2}=\Big( c_{b2} \frac{\Delta t}{\Delta x^2} \Big)^2,
\end{gather*}
and for ease of notation, and since in the simulations $\Phi$ was a linear
function of $\beta$, we let $ \Phi=\beta\cdot \widetilde{\Phi}$, where
\[
\widetilde{\Phi}_j^{n+1}=\widetilde{\Phi}(u_{1j+j_2-1}^{n+1}, u_{2j}^{n+1},
w_{1j+j_2-1}^{n+1}, w_{2j}^{n+1}, \eta_j^{n}).
\]
These lead to the linear system:
\begin{gather}
\begin{aligned}
&- C_{r1}u_{1j-1}^{n+1}+(2C_{r1}+1)u_{1j}^{n+1}- C_{r1}u_{1j+1}^{n+1} \\
& =-u_{1j}^{n-1}+(2- \Delta t^2 \beta_{j-j_2+1}^{n} K_1\chi_j)u_{1j}^{n}
+ \Delta t^2 \beta_{j-j_2+1}^{n} K_1\chi_ju_{2j-j_2+1}^{n},
\end{aligned} \label{eqn610}
\\
\begin{aligned}
& - C_{r2}u_{2j-1}^{n+1}+(2C_{r2}+1)u_{2j}^{n+1}
- C_{r2}u_{2j+1}^{n+1}\\
&=-u_{2j}^{n-1}+(2- \Delta t^2 \beta_j^{n} K_1\chi_j)u_{2j}^{n}
+ \Delta t^2 \beta_j^{n} K_1\chi_ju_{1j+j_2-1}^{n},
\end{aligned} \label{eqn611}
\\
\begin{aligned}
& C_{b1}w_{1j-2}^{n+1}-4C_{b1}w_{1j-1}^{n+1}+(1+6C_{b1})w_{1j}^{n+1}
-4C_{b1}w_{1j+1}^{n+1}+C_{b1}w_{1j+2}^{n+1}\\
& =(2-\beta_{j-j_2+1}^{n} \Delta t^2 \chi_j)w_{1j}^{n}+
\beta_{j-j_2+1}^{n} \Delta t^2 \chi_jw_{2j-j_2+1}^{n},
\end{aligned} \label{eqn612s}
\\
\begin{aligned}
& C_{b2}w_{2j-2}^{n+1}-4C_{b2}w_{2j-1}^{n+1}+(1+6C_{b2})w_{2j}^{n+1}
-4C_{b2}w_{2j+1}^{n+1}+C_{b2}w_{2j+2}^{n+1} \\
&=(2-\beta_j^{n} \Delta t^2 \chi_j)w_{2j}^{n}
+ \beta_j^{n} \Delta t^2 \chi_jw_{1j+j_2-1}^{n},
\end{aligned} \label{eqn613}
\\
-2Ck_{\beta}\beta_{j-1}^{n+1} + (1+4Ck_{\beta})\beta_j^{n+1}
-2Ck_{\beta}\beta_{j+1}^{n+1}+ 2\Delta t \beta_j^{n+1}\widetilde{\Phi}_j^{n+1}
=\beta_j^{n-1}, \label{eqn614}\\
-2CD\eta_{j-1}^{n+1}
+ (1+4CD)\eta_j^{n+1}-2CD\eta_{j+1}^{n+1} =\eta_j^{n-1}, \label{eqn615}
\end{gather}
along with initial and boundary conditions \eqref{67}, \eqref{eqn68}
and \eqref{eqn68}.
Moreover, to take care of the subdifferential conditions in \eqref{27a} and
\eqref{28}, in the calculations we preceded as follows. First, to guarantee
that $w^{n+1}_{1j}\geq w^{n+1}_{2j}$, when
it was found that the computed values were
$\widetilde{w}^{n+1}_{1j}< \widetilde{w}^{n+1}_{2j}$
for $j\in \{j_2, \dots,S+1\}$, we set
\[
w^{n+1}_{1j}=w^{n+1}_{2j}=\frac12 (\widetilde{w}^{n+1}_{1j}
+ \widetilde{w}^{n+1}_{2j}).
\]
Second, when it was found computationally that $\widetilde{\beta}^{n+1}_j<0$,
we set
$\beta^{n+1}_j=0$.
We note that from the structure of \eqref{28} and the assumption that $\Phi>0$,
it follows that $\beta$ is decreasing, computationally, too, so that starting
with $\beta_0\leq 1$
\eqref{eqn614} guaranteed that $\beta^{n+1}_j\leq 1$.
Next, we observe that in the system \eqref{eqn68} and \eqref{eqn68}
we have to deal with the cases when $j=2$
and $j=S$ since the equations involve undefined
quantities that are outside of the spatial domain: $w_{10}^{n+1}$,
$w_{1S+2}^{n+1}$, $w_{20}^{n+1}$ and $w_{2S+2}^{n+1}$. However, it
is straightforward to eliminate these quantities using the initial
and boundary conditions specified by \eqref{67}, \eqref{eqn68} and \eqref{eqn68}.
Finally, we set $\mathbf{v}$ to be the $\Lambda$-dimensional column vector
\begin{align*}
\mathbf{v}&=(u_{11},\dots, u_{1S+1}, u_{21},\dots, u_{2 S+1},w_{11},\dots,
w_{1S+1}, w_{2j},\dots, w_{2 S+1},\\
&\quad \dots, \beta_{j_2},\dots, \beta_{S-j_2+2}, \eta_{j_2},\dots, \eta_{S-j_2+2}),
\end{align*}
where
$\Lambda =6(S+1) - 2(j_2-1)$.
Next, we let $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ be the matrices,
based on the system \eqref{eqn610}--\eqref{eqn615}, so that it may be written as
\begin{equation}\label{main-numerical}
\mathbf{A} \mathbf{v}^{n+1}=\mathbf{B} \mathbf{v}^{n} + \mathbf{C} \mathbf{v}^{n-1}.
\end{equation}
Here,
\[
A =
\begin{bmatrix}
A_{u1} &&&0&&\\
& A_{u2}&&& \\
& & A_{w1} \\
& & & A_{w2}\\
& & & & A_{\beta}\\
&0 & & & & A_{\eta}
\end{bmatrix},
\]
where:
\begin{gather*}
A_{u1}= \begin{bmatrix}
1 & 0 & 0 & \ldots & & &0\\
-C_{r1} & 2C_{r1} + 1 & - C_{r1} & 0 & \ldots && 0\\
0 &- C_{r1} & 2C_{r1} + 1 & - C_{r1} & 0 & \ldots & 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0&\ldots&&\ldots&0 & -C_{r1} & 2C_{r1} + 1 &\\
0 &\ldots & & \ldots & 0& -1 & 1
\end{bmatrix},
\\
A_{u2}= \begin{bmatrix}
-1 & 1 & 0 & \ldots & & &0\\
-C_{r2} & 2C_{r2} + 1 & - C_{r2} & 0 & \ldots && 0\\
0 &- C_{r2} & 2C_{r2} + 1 & - C_{r2} & 0 & \ldots & 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0&\ldots&&\ldots&0 & -C_{r2} & 2C_{r2} + 1 &\\
0 &\ldots & & \ldots & 0& -1 & 1
\end{bmatrix},
\\
{\scriptsize A_{w1}= \begin{bmatrix}
1 & 0 & & \ldots & & & & &0\\
0& 1+7C_{b1} & -4C_{b1} & C_{b1} & 0 & \ldots &&& 0\\
C_{b1} &-4C_{b1} & 1+6C_{b1} & - 4C_{b1} & C_{b1} & 0 & \ldots && 0\\
0 & C_{b1} &-4C_{b1} & 1+6C_{b1} & - 4C_{b1} & C_{b1} & 0 & \ldots & 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0&\ldots&&\ldots&0&C_{b1} & -4C_{b1} & 1+5C_{b1} &-2C_{b1}\\
0 &\ldots &&& \ldots &0& 2C_{b1}& -4C_{b1} & 1+2C_{b1}
\end{bmatrix}},
\\
{\scriptsize A_{w2}= \begin{bmatrix}
1+2C_{b2} & -6C_{b2} & 0 & \ldots & & & & &0\\
-2C_{b2} & 1+6C_{b2} & -5C_{b2} & C_{b2} & 0 & \ldots &&& 0\\
C_{b2} &-4C_{b2} & 1+6C_{b2} & - 4C_{b2} & C_{b2} & 0 & \ldots && 0\\
0 & C_{b2} &-4C_{b2} & 1+6C_{b2} & - 4C_{b2} & C_{b2} & 0 & \ldots & 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0&\ldots&&\ldots&0&C_{b2} & -4C_{b2} & 1+5C_{b2} &-2C_{b2}\\
0 &\ldots &&& \ldots &0& 2C_{b2 }& -4C_{b2} & 1+2C_{b2}
\end{bmatrix},
}\\
A_{\beta}= \begin{bmatrix}
-1 & 1 & 0 & \ldots & & &0\\
-2Ck_{\beta} & 1+4Ck_{\beta}+2\Delta t \widetilde{\Phi}_2^{n+1} &-2Ck_{\beta}& 0 & \ldots &&& 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0 &\ldots & & \ldots & 0& -1 & 1
\end{bmatrix},
\\
A_{\eta}= \begin{bmatrix}
1 & 0 & 0 & 0& \ldots & &0\\
-2DC & 1+4CD &
-2CD&0& \ldots && 0\\
\\
0&-2CD & 1+4CD &
-2CD&0 & \ldots & 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0 &\ldots & & -2CD & 1+4CD& & -2CD\\
0 &\ldots & & \ldots & 0& & 1
\end{bmatrix}.
\end{gather*}
We now describe the algorithm for the system. For the sake of simplicity,
we describe it for the case of zero initial displacements and velocities,
and periodic tractions $p^n=p(t_n)$ and $q^n=q(t_n)$ at the end $x=1$.
The modification of the algorithm for the cases when an impulsive force is
applied at $x=1$ initially are straightforward.
To initialize the process we compute $\mathbf{v}^1$ and $\mathbf{v}^2$ as follows:
\begin{gather}
\begin{gathered}
u_{1j}^1=w_{1j}^1=0, \quad j=1,\dots, S+1; \\
u_{2j}^1=w_{2j}^1=0, \quad j=1, \dots, S+1; \\
\beta_j^1=1,\quad j=1, \dots, S- j_2+2; \\
\eta_j^1=0,\quad j=2, \dots,S- j_2+1, \\
\eta_1^1=\eta_{L}^1,\quad \eta_{S-j_2+2}^1=\eta_{R}^1,
\end{gathered} \label{eqn617}
\\
\begin{gathered}
u_{1j}^2=W_{1jk} w_{1k}^2=0, \quad j=1,\dots, S+1; \\
u_{2j}^2=0,\quad j=1, \dots, S, \\
u_{2S+1}^2=\Delta x \cdot p^2; \\
W_{2jk}w_{2k}^2=R_{2jk} w_{2k}^1,\quad j=1, \dots, S+1, \;
k=1, \dots, S+1, \\
\beta_j^2=1,\quad j=1, \dots,S-j_2+2; \\
\eta_j^2=0 \quad j=2, \dots, S-j_2+1; \\
\eta_1^2=\eta_{L},\quad \eta_{S-j_2+2}^1=\eta_{R}.
\end{gathered} \label{eqn618}
\end{gather}
Here, summation over the index $k$ is implied, and $W_1$, $W_2$ and $R_2$
are given by:
\begin{gather*}
{\scriptsize
W_1= \begin{bmatrix}
1 & 0 & 0 & \ldots & & & & &0\\
0 & 2+7C_{b1} & -4C_{b1} & C_{b1} & 0 & \ldots &&& 0\\
C_{b1} &-4C_{b1} & 2+6C_{b1} & - 4C_{b1} & C_{b1} & 0 & \ldots && 0\\
0 & C_{b1} &-4C_{b1} & 2+6C_{b1} & - 4C_{b1} & C_{b1} & 0 & \ldots & 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0&\ldots&&\ldots&0&C_{b1} & -4C_{b1} & 2+5C_{b1} &-2C_{b1}\\
0 &\ldots &&& \ldots &0& 2C_{b1 }& -4C_{b1} & 2(1+C_{b1})
\end{bmatrix},
}\\
{\scriptsize
W_2= \begin{bmatrix}
2(1+C_{b2}) & -6C_{b2} & 0 & \ldots & & & & &0\\
-2C_{b2} & 2+6C_{b2} & -5C_{b2} & C_{b2} & 0 & \ldots &&& 0\\
C_{b2} &-4C_{b2} & 2+6C_{b2} & - 4C_{b2} & C_{b2} & 0 & \ldots && 0\\
0 & C_{b2} &-4C_{b2} & 2+6C_{b2} & - 4C_{b2} & C_{b2} & 0 & \ldots & 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0&\ldots&&\ldots&0&C_{b2} & -4C_{b2} & 2+5C_{b2} &-2C_{b2}\\
0 &\ldots &&& \ldots &0& 2C_{b2 }& -4C_{b2} & 2(1+C_{b2})
\end{bmatrix},
}\\
R_2= \begin{bmatrix}
0 & 0 & 0 & \ldots & & &0\\
0 & 0& 0 & \ldots &&& 0\\
\\
\vdots& &&\ddots&&&\vdots\\
\\
0 &\ldots & & \ldots & 0& 0 & -\frac{C_{b2}}{k_2I}q^2\Delta x^{3}.
\end{bmatrix}.
\end{gather*}
Once the system was initialized, and at time step $n+1$ the solutions
$\mathbf{v}^1,\dots, \mathbf{v}^{n}$ were found, system
\eqref{main-numerical} was solved for $\mathbf{v}^{n+1}$.
\begin{algorithm}
\caption{Finite Differences for the Model}
\begin{algorithmic}
\STATE Set $\mathbf{v}^1, \mathbf{v}^2$ according to (6.17), (6.18)
\hfill\COMMENT{Initial conditions}
\FOR{$n=1\ldots N$}
\STATE Solve $\bf A \mathbf{v}^n = B \mathbf{v}^{n-1} + C \mathbf{v}^{n-2}$
\hfill\COMMENT{Block diagonal solver}
\FOR{$j=j_2 \ldots S+1$}
\IF{$\tilde w_{1j} \ge \tilde w_{2j}$}
\STATE $w^n_{1j},w^n_{2j} = \tilde w^n_{1j},\tilde w^n_{2j}$
\ELSE
\STATE $w^n_{1j} = w^n_{2j} = \frac12 (\tilde w^n_{1j} + \tilde w^n_{2j})$
\hfill\COMMENT{Smooth out}
\ENDIF
\IF{$\tilde \beta^n_j \ge 0$}
\STATE $\beta^n_j = \tilde \beta^n_j$
\ELSE
\STATE $\beta^n_j = 0$\hfill\COMMENT{Clip}
\ENDIF
\ENDFOR
\ENDFOR
\end{algorithmic}
\end{algorithm}
\section{Simulations}
\label{sec7}
The algorithm of the previous section was implemented and run extensively. Here,
we depict some of the simulation results that we consider of interest as they
allow us to gain understanding of the debonding process dynamics.
Since our main interest was in the debonding process resulting from humidity
and mechanical vibration, we present simulations along those lines.
In particular, we depict the dynamics of the system and the evolution of the
debonding.
\vskip4pt
First, we show typical solutions with oscillating tractions. Our main
interest lies in answering two questions:
\begin{itemize}
\item[(1)] How does the debonding process affect the
vibrations spectrum?
\item[(2)] How do the periodic oscillations of the tractions affect the spatial
distribution of the debonding process?
\end{itemize}
The constants used in the simulations are given in Table \ref{T1}.
\begin{table}[htb]
\caption{Simulations system parameter values}
\label{T1}
\begin{center}
\begin{tabular}{|ccc|}
\hline
Parameter & Parameters value& Units\\
\hline
$l_1; l_2$ & 1;\ 1 & m \\
$E_1; E_2$ & $200 \cdot 10^5$; \ $190 \cdot 10^5$ & Pa \\
$A_1; A_2$ & $10^{-4}$; \ $10^{-4}$ & $m^2$ \\
$B_1; B_2$ & $8 \cdot 10^{-6}$;\ $8 \cdot 10^{-6}$ & m$^{4}$ \\
$\nu_{r1}; \nu_{r2};\ \nu_{b1}; \nu_{b2}$ & 0;0;\ 0;0 & m$^2$/s; m$^{4}$/s \\
$c_{r1}; c_{r2}$ & 50;\ 49.13 & m/s \\
$c_{b1}; c_{b2}$ & 14.14; \ 13.89 & m$^2$/s \\
$K_{r1}; K_{r2}$ & $10^{4}$; \ $9 \cdot 10^{3}$ & kg/s$^2$ \\
$K_{b1}; K_{b2}$ & $5 \cdot 10^{3}$; \ $7 \cdot 10^{3}$ & kg/s$^2$ \\
$k_{\beta}$ & 0.01 & m$^2$/s \\
$d; d_{\beta}$ & 0.01; \ 0.01 & m$^2$/s \\
$\kappa$ & $2.3 \cdot 10^{-5}$ & m$^2$/s \\
$\alpha_{h} ;\ \epsilon_{\eta}$ & $650;\ 0.001$ & 1/(m$\cdot$ s); \ 1 \\
\hline
\end{tabular}
\end{center}
\end{table}
Next, the various functions in the model were chosen as follows.
The tractions at $x=1$ were either zero or
\[
p(t)=0.0053 \cos(2\pi f t) ,\quad q(t)= 5.263 \cos(2\pi f t),
\]
where $f$ is the frequency of the traction which we describe in each case.
The humidity at the ends was assumed, for the sake of simplicity, to be
\[
\eta_{L}(t)=1, \quad \eta_{R}(t)=1.
\]
The water diffusion coefficient in the glue layer was assumed to increase
as the bonding decreases, since the water flow in voids is easier, and was
chosen as
\[
D=d+ d_{\beta}(1-\beta).
\]
Thus, when there is full bonding $D=d$, since $\beta=1$, and as debonding progresses
$D$ increases to the value $D=d+ d_{\beta}$, which is the diffusion coefficient
on the debonded surface.
Finally, the debonding source function was assumed to have the `simple' form
\[
\Phi= \alpha_h \beta (|u_1-u_2| +|w_1-w_2|)(\epsilon_\eta +\eta^2),
\]
with the debonding rate constant $\alpha_h$ and the threshold humidity
constant $\epsilon_\eta$.
In some of the simulations, we used $p=0$ and $q=0$ and $D=d$, and we
indicate this in each case.
\subsection{Simulation 1: Natural frequencies of the first rod}
To gain confidence in the numerical simulations, we first compared the natural
frequencies of the first rod, obtained by a simple Fourier analysis,
with the computed frequencies from the simulations. The left rod, without bonding,
was excited with an initial impulse at the right
end and the natural frequencies were found by using the Fast Fourier Transform
(FFT), a subroutine in MATLAB, that was applied to the vibrations of the
(middle of the) rod.
The theoretical natural frequencies were:
\begin{equation*}
f_{n}=( n-\frac{1}{2} ) \frac{1}{2L}c_{r1},\quad n=1,2,3,\dots
\end{equation*}
where $c_{r1}$ is determined from the Young Modulus $E_1$ and the density $\rho$
of the material, i.e. $c_{r1}=\sqrt{{E_1}/{\rho}}$.
For $L=1\ m$ and $c_{r1}=50\ m/s$. Therefore,
\begin{gather*}
f_1=12.5,\quad f_2=37.5, \quad f_3=62.5,\quad f_4=87.5,\\
f_5=112.5, \quad f_6=137.5\quad f_{7}=162.5 \; Hz
\end{gather*}
\begin{figure}[htb]
\centering
\includegraphics[width=3.8in, height=2in]{fig2} % fourierRu1.png
\caption{Natural modes of a single rod (without bonding).}
\label{fnatu1}
\end{figure}
The results of the simulations are depicted in Figure \ref{fnatu1}.
It was found, using
the FFT with moderately small time steps, that the computed frequencies were the
same up to the second place after the decimal point,
\begin{gather*}
\widetilde{f}_1=12.50,\quad \widetilde{ f}_2=37.50, \quad \widetilde{f}_3=62.50,\quad
\widetilde{f}_4=87.50,\\
\widetilde{f}_5=112.50, \quad
\widetilde{f}_6=137.50,\quad \widetilde{f}_{7}=162.50 \; Hz.
\end{gather*}
It is seen that the correspondence between the theoretical and computed
frequencies is excellent.
\subsection{Simulation 2: Natural frequencies of the rigid system}
We studied the natural frequencies of the whole system when
adhesion was full and not changing, i.e., without deterioration, so we considered a
fixed bonding field with $\beta(x,t)\equiv 1$. An impulse was applied to the right
end of the second rod, the tractions were $p=q=0$, $D=d$, and
the natural frequencies of the (rigid) system $f_1,
\dots, f_{7} (Hz)$, depicted in Figure \ref{fourieru1}, were found to be:
\[
f_1=6,\quad f_2=18, \quad f_3=28,\quad f_4=39,\quad f_5=75,\quad
f_6=88,\quad f_{7}=100 \; Hz.
\]
\begin{figure}[htb]
\centering
\includegraphics[width=3.8in, height=2in]{fig3} % system_modes.png}
\caption{The natural frequencies of the fully bonded system, detected in the
middle of the adhesive region $x=l_*$}
\label{fourieru1}
\end{figure}
The system (with $\beta\equiv 1$) was solved and the FFT was used to
obtain the spectrum above from the vibrations at the center of the adhesive region
$x=l_*=(l_1+l_2)/2$. It is seen that the amplitudes of the first two frequencies,
$f_1$ and $f_2$, dominate, and the amplitudes of the fourth and higher resonance
frequencies are orders of magnitude smaller.
\subsection{Simulation 3: Debonding with $25Hz$ traction}
In this simulations a horizontal periodic traction $p(t)$, with period $25Hz$,
was applied to the right rod (with $q=0$), and the deterioration of the bonding was
computed. The choice of $\alpha_h $ was such that almost complete debonding
happened in less than 3 sec, which was not realistic in most applications,
but allowed us to run many simulations. In particular, we addressed the first
question raised at the beginning of the section about the spectrum shift.
First, the displacements of the middle of the adhesive layer vs. time,
$u_1(l_*, t)$, are depicted in Figure \ref{u1f25T}. It is seen clearly that
the oscillations were quite
complicated and changed as the debonding progressed. As we show below, the bonding
was insignificant at $2.4$ sec, and complete debonding happened before $3$ sec.
Then, at about 2 sec, when debonding was weak, the first rod settled into periodic
oscillations in its first free resonance frequency. The second rod became almost
free a bit before $3$ sec, actually by then the bonding field was negligibly small.
\begin{figure}[ht]
\centering
\includegraphics[width=3.8in, height=1.8in]{fig4} % u1_25hz.png
\caption{The displacement $u_1(l_*, t)$ vs. time in the first rod. Debonding
is complete at about 2 s.}
\label{u1f25T}
\end{figure}
Figure \ref{difference} illustrates the
amplitude difference $|u_1-u_2|$ at different times (indicated in color)
in the adhesive region of the slabs.
We recall that this term is a factor in $\Phi$ and, therefore, directly affects the
debonding process. It is seen that the differences were larger at the left end of the
adhesive region and tapered off toward the right end, where they
were an order of magnitude smaller.
\begin{figure}[htb]
\centering
\includegraphics[width=3.8in, height=1.8in]{fig5} % u1u2c.png
\caption{$|u_1-u_2|$ at different times}
\label{difference}
\end{figure}
The evolution of the bonding field $\beta$,
and of the humidity are depicted in Figure \ref{beta_eta_times} at different times
(indicated in color). It is seen that the bonding function decreases monotonically
in time while the humidity increases monotonically. The bonding function almost
vanishes at $2.4$ sec, only a small region in the middle is not negligible, but very
close to being zero.
\begin{figure}[htb]
\centering
\includegraphics[width=2.4in, height=1.3in]{fig6a} % betaTc.png
\includegraphics[width=2.4in, height=1.3in]{fig6b} % etaTc.png
\caption{The bonding field $\beta$ (L) and humidity $\eta$ (R) at different
times (color). The bonding decreases monotonically while the humidity
increases monotonically.}
\label{beta_eta_times}
\end{figure}
\begin{figure}[htb]
\centering
\includegraphics[width=2.4in, height=1.3in]{fig7a} % fu105sN.png
\includegraphics[width=2.4in, height=1.3in]{fig7b} \\ % fu115sN.png
\includegraphics[width=2.45in, height=1.3in]{fig7c} % fu13sN.png
\caption{Spectrum in the first rod detected after 0.5 sec (UL), after
1.5 sec (UR) and after 3 sec (B).}
\label{fu105s}
\end{figure}
We turn to one of the two main findings in our computer simulations, namely,
the shift of the spectral frequencies of the system as debonding progresses.
It was found, as one would expect, that as the debonding
process advances the vibrations of the bonded system change and once
there is full debonding, the vibrations frequencies are those of the the
two free rods.
Indeed, this is depicted in Figure \ref{fu105s}. The peak at the
driving frequency $25\ Hz$ when bonding was essentially full can be seen
clearly in the upper left (UL) figure; when bonding became weak (after 1.5 sec),
the spectrum widened, became less pronounced and shifted to the left
toward the first natural frequency; then the $25\ Hz$
peak completely disappeared after 3 sec, since the bonds were
broken and the first natural frequency $f_1=12.5\ Hz$ (see section 7.1)
was all that remained. Indeed, as can be seen in
Figure \ref{fu13s}, where the three results were combined in one figure,
the frequencies shifted down and converged to the fundamental frequency
of a single rod.
\begin{figure}[htb]
\centering
\includegraphics[width=3.75in, height=1.8in]{fig8} % fu_timesN.png
\caption{The spectrum in the first rod at the three times -- detecting the shift}
\label{fu13s}
\end{figure}
We return to this point in the conclusions section, Section \ref{sec9},
since it may be of considerable interest to use such frequency shifts to
detect the level of deterioration of the adhesive in the system using
nondestructive testing by employing induced vibrations
of controlled frequencies.
\subsection{Simulation 4: Debonding with periodic tractions: $150Hz$ and $350Hz$}
We next describe two simulations of the debonding process induced by
tractions with periods of $f=150Hz$ and $f=350Hz$ in $p$ and $q$. We are
interested in the second question raised above on how does the frequency of
the applied tractions affect the spatial distribution of the adhesive.
It was found in the simulations that the evolution of the bonding field $\beta$
depended strongly on the frequency of the applied tractions. Here, we present two
such simulation results with two different frequencies of the
tractions $p$ and $q$. Again, since the rate of change of $\beta$ is affected by
the term $| u_1-u_2 |$, the dependence is guaranteed, but the spatial
shape of $\beta$ at different times is interesting, and very important in
applications.
Indeed, as can be seen in Figs.\,\ref{absf150} and \ref{absf350}, at the two
frequencies standing waves were formed that strongly affected the shapes of
the adhesive spatial distribution. The frequencies of these standing waves
were those of the applied traction, and the corresponding wavelengths were
\begin{gather*}
\lambda = \frac{c_1}{26 \cdot \Delta x}=0.325 \text{m}, \quad f=153.8 \; Hz, \\
\\
\lambda = \frac{c_2}{11 \cdot \Delta x}=0.1375 \text{m}, \quad f=363 \; Hz.
\end{gather*}
It is noted that when the frequency of the applied traction increased, while keeping
the amplitude the same, the bonding field's deterioration rate became slower.
This seems to be related to the fact that at higher frequencies the wavelengths
were shorter and there were more nodes where there was no deterioration
since $| u_1-u_2 |=0$ at those nodes.
\begin{figure}[htb]
\centering
\includegraphics[width=4in, height=1.5in]{fig9a} % beta_150hz_t.png
\includegraphics[width=4in, height=1.5in]{fig9b} \\ % abs_u1_u2_f150.png
\caption{Evolution of the bonding field $\beta(x,t)$ (U), and the difference
$| u_1-u_2 |$ (B), at different times. The debonding was
much faster at the positions where the difference $| u_1-u_2 |$
was the largest, and was very slow where it was small. The traction force frequency
was $f=150\; Hz$.}
\label{absf150}
\end{figure}
\begin{figure}[htb]
\centering
\includegraphics[width=4in, height=1.4in]{fig10a} % betas_f350.png
\includegraphics[width=4in, height=1.4in]{fig10b} \\ % abs_u1_u2_f350.png
\caption{Evolution of bonding field $\beta(x,t)$ (U), and the difference
$| u_1-u_2 |$ (B), at different times. The debonding was
much faster at the positions where the difference $| u_1-u_2 |$
was the largest, and was very slow where it was small. The traction force frequency
was $f=350\; Hz$.}
\label{absf350}
\end{figure}
Moreover, it was found that as the frequency increased, to obtain smoother curves
the discretization of the spatial domain had to be considerably refined, which in
turn led to much smaller time steps and longer runs. So the same discretization
was kept throughout all of the presented results, except those in Section 8.
\subsection{Simulation 5: Debonding with periodic traction $350Hz$ and different
values of the diffusion coefficient $k_{\beta}$}
To gain insight into the dependence of the debonding process on the diffusion
coefficient $k_{\beta}$, which as was noted above is hard to measure and so it must
be estimated, we run simulations of the debonding process induced by a
horizontal traction with frequency $f=350Hz$ with high and low values of the
coefficient.
We depict in Figure \ref{figkb} the results with two values:
$k_{\beta}=10^{-3}$ (U), and
$k_{\beta}=10^{-5}$ (B). The figure illustrates how a higher value of the
diffusion coefficient causes the smoothing of the debonding process.
Moreover, it also slows the process, as can be seen from the curves in the
figure. It is seen clearly if
one compares the fourth (pink) curve in each case.
\begin{figure}[htb]
\centering
\includegraphics[width=4in, height=1.5in]{fig11a} \\% betas_kb_H.png
\includegraphics[width=4in, height=1.5in]{fig11b} % betas_kb_L.png
\caption{Evolution of the bonding field $\beta$ with $k_{\beta}=10^{-3}$ (U),
and $k_{\beta}=10^{-5}$
(B), at different times.}
\label{figkb}
\end{figure}
It may be of interest to estimate the coefficient $k_{\beta}$ from the form of
the debonding
field found experimentally. Indeed, in some of the results obtained in \cite{NS15}
the debonding was found to be wavy, and it may be possible to correlate these
results to an estimation of the coefficient, using a parameter identification
and optimization methods.
\subsection{Simulation 6: The cases $D=const$. and $D=D(\beta)$}
The final simulations address the issue of the effects when the humidity diffusion
coefficient $D$ is not a constant but depends on $\beta$, thus making the whole
system fully coupled (since when $D=const$. the humidity equation is not coupled
to the rest of the model and can be solved separately). Therefore, we compare
the results in both cases of $D=const$. and $D=D(\beta)$.
We assumed, for the sake of simplicity, that the diffusion coefficient had the
following form,
\[
D=d+d_\beta(1-\beta).
\]
That is, a linear function that was increasing with the advance of debonding.
A comparison of the evolution of humidity in case when the diffusion coefficient
is constant and when it is given above is depicted in Figure \ref{compD}.
The coefficients were chosen as $d=0.01$ and $d_{\beta}=0.01$.
The case with constant $D$ is shown in red
and when $D$ depends on $\beta$ in blue, at five different times.
It is seen that the predictions with such coefficients were very similar
both qualitatively and quantitatively. Moreover, it is
very likely that the diffusion coefficients are much smaller in practice,
which would make the difference even smaller.
However, we note that in the second case the humidity diffusion is slightly faster.
\begin{figure}[htb]
\centering
\includegraphics[width=4in, height=2.0in]{fig12} compD.png
\caption{Humidity $\eta$ with diffusion coefficients $D=d=0.01$
(red), and $D=0.01+0.01(1-\beta)$
(blue), at five different times.}
\label{compD}
\end{figure}
We conclude that using $D$ that depends on $\beta$ may not be very helpful, since
this complicates the model without any clear benefits.
\section{Numerical convergence} \label{sec8}
The theoretical study of the rate of convergence of the algorithm was deemed too
complicated, at this stage, in view of the complex nonlinear structure of the model,
and so was left open. However, to gain additional confidence in the computer
simulations from a different perspective than in Section 7.1, we performed
a numerical study of the scheme's convergence. To that end, we run the simulations
with ten different time steps, each one was half of its predecessor.
Choosing $T=1$, the number of the steps was
\[
N_1=5\cdot10^3, \quad N_2=2\cdot N_1=10^{4},\dots, N_{10}=2^{10}
\cdot N_1=2.56\cdot 10^{6},
\]
with corresponding time steps
\[
\Delta t_1=\frac{T}{N_1}=2\cdot 10^{-4}, \dots ,
\Delta t_{10}=3.9\cdot 10^{-7}.
\]
We assumed, as is customary, that the numerical solution for the smallest time step
($\Delta t_{10}=T/N_{10}$) represents the `true solution.' Then, as a measure of
convergence, we calculated the $l^2$ norm of consecutive differences between
the solution $\beta=\beta^{N_i}$ with $\Delta t_i$ and the solution
$\beta=\beta^{N_{10}}$ with $\Delta t_{10}$, i.e.,
\[
\|\beta^{N_i}-\beta^{N_{10}}\|_2^2=\sum_{l=1}^{S+1}|\beta^{N_i}_{j_2+l}-
\beta^{N_{10}}_{j_2+l}|^2,
\]
for $ i=1,\dots,9$, and also the ratios
\[
R_i=\frac{\|\beta^{N_i}-\beta^{N_{10}}\|_2}{\|\beta^{N_{i+1}}-\beta^{N_{10}}\|_2},
\]
for $ i=1,\dots,8$,
The norm and the ratio are summarized in Table \ref{T2}.
It is seen that as the number of time steps increases the
ratio is converging to $4$, which implies second order convergence of the
algorithm.
We conclude that the algorithm was very stable and robust.
\begin{table}[htb]
\caption{Numerical Errors with Respect to $l^2$ Norm}
\label{T2}
\renewcommand{\arraystretch}{1.2}
\label{eqn618b}
\begin{tabular}{|l|l|l|l|}
\hline
Times Step $*10^3$ & Difference in $l^2$ Norm & Ratio R \\ \hline \hline
5 & 0.00669 & 1.37 \\
10 & 0.004865 & 2.00\\
20 & 0.0024266 & 2.73 \\
40 & 0.0008865 & 3.28 \\
80 & 0.0002695 & 3.61 \\
160 & 0.00007451& 3.80 \\
320 & 0.00001960 & 3.90 \\
640 & 0.000005019& 3.95 \\
1280 & 0.000001270 & \\
2560 & `true solution' & \\ \hline
\end{tabular}
\end{table}
Finally, Figure \ref{betaConv} depicts the behavior of the bonding field at time
$t=1$ for the ten different time steps. It visually confirms the conclusion based
on the table.
\begin{figure}[htb]
\centering
\includegraphics[width=3.75in, height=2.0in]{fig13} % betaCon_f25.png
\caption{$\beta$ for ten different time steps, at time $t=1$. The traction
was $f=25\ Hz$}
\label{betaConv}
\end{figure}
\section{Concluding remarks} \label{sec9}
We presented a `simple' model for the process of debonding of two bonded slabs,
in the so-called single lap joint, as a result of mechanical vibrations,
temperature fluctuations and spread of humidity in the adhesive layer.
The model was based on two 1D beam-rod
systems taking into account the horizontal shear
(rods) and the vertical forces (beams) in the adhesive layer.
Our main interest was in the model,
its analysis and in studying the dependence of the debonding process humidity and
on the frequencies of the mechanical vibrations and the related shift is the
spectrum as debonding progresses.
The debonding process was described by the evolution of the bonding function
$\beta$, the temperature $\theta$ and the diffusion of the humidity function
$\eta$, all three defined on the adhesion region. This gave the system an unusual
form, that of a system of dynamic equations for the slabs coupled
over the adhesive region, together with the evolution inclusion for the
bonding field and the diffusion equations for temperature and humidity.
Since the system was nonlinear and of unusual form, the first step was
to analyze it. The existence of weak solutions to the system was
established by using approximate problems, passing to the limits and
a fixed point argument. Indeed, use was made of a number of recent
results and tools from the theory of differential inclusions with
pseudomonotone operators.
To gain insight into the possible evolution of the system, slightly simplified by
omitting the temperature field, and possibly
creating more accurate models for the prediction of the debonding process,
we turned to computer simulations. A computer algorithm was developed
for the system based on fully implicit time discretization
and a standard spatial discretization. The algorithm was implemented and many
simulations were conduced to gain insight into the model behavior. As
was shown in Section 7.1 and Section 8, the numerical solutions seemed
to be accurate and the algorithm robust and efficient, and it was found
(numerically) to have almost quadratic convergence.
Then, we presented a few simulations' results. Many other simulations of the
system behavior under various conditions and various assumptions on the
problem data, especially the form of the debonding source function $\Phi$,
can be found in \cite{PawelDist}.
The first simulation (Section 7.1) dealt with the spectrum of an excited single
rod and comparison with known natural frequencies. The comparison was found
to be very good. Then, in Section 7.2 the spectrum of the whole system, when
fully bonded, was studied numerically using the FFT. That was the baseline
simulation to which the evolution of the bonding field, actually the
debonding (the decrease in $\beta$) was compared to.
In the third simulation, Section 7.3, we studied the debonding process
caused by the application of a vibrating traction with frequency of $25\ Hz$
and given humidity at the ends. The main interest here was the discovery that the
vibration spectrum of the system changed, moved to lower frequencies and
broadened for some time, as debonding progressed from full bonding to
almost full debonding. This resulted in free rod vibrations when
the debonding was essentially complete, Figs.\,\ref{fu105s} and \ref{fu13s}.
The fourth simulations, Section 7.3, studied a similar setting but with tractions
that had frequencies of $150$ and $350\ Hz$. There, it was found that because
of standing waves in the system, which depended on the applied frequency, the
debonding was slower and exhibited a wavy form, which may be interest in
applications, since it may indicate that regions with high debonding were
separated from regions with low debonding, Figs.\,\ref{absf150} and \ref{absf350}.
Finally, we investigated the dependence of the model on either a constant
humidity diffusion coefficient or a coefficient that depends on the bonding,
i.e., $D=D(\beta)$. It was found that for the values chosen above, there
were some quantitative differences, but qualitatively the solutions looked
very similar in form. However, this waviness, which was also found experimentally
in \cite{NS15, SN16}, may be used to estimate the bonding diffusion coefficient
$k_{\beta}$.
The value of the model, and of this work, lies in the fact that the model
was found simple enough to analyze, but sufficiently complex to
allow for considerable insight into its predictions. Indeed, it allowed us to
run many simulations, since the run times were in minutes.
The main findings are the dependence of the debonding spatial form
on the applied frequency and the shift in the spectrum as debonding progresses.
This clearly indicated that there may be ways to find the bonding/debonding
state of a bonded system by nonintrusive measurements of the spectrum
by externally exciting it. Indeed, it may be possible to excite the system
externally and by measuring the resulting spectrum to correctly estimate
the extent of the debonding.
To continue the line of research begun here, there are four steps
that need to be completed to make the model useful in real applications.
The first step is to find from experimental data and general engineering
approaches an appropriate form of the debonding source function $\Phi$.
Without reasonably accurate $\Phi$, the predictions are likely to be only
qualitative. It may also be of interest to use different functions with
different forms and compare their effects on the debonding process.
The second need is to derive the system, or some related form, from a 3D model
in the limit of slabs that are long and thin. Some progress in this direction has been
made in \cite{AKS17}. A general way of obtaining such models can be found in
\cite{AV10, TV96}.
Thirdly, there is a need to use the model to study the debonding process when
tractions with many frequencies are applied.
Finally, there is a need to introduce randomness to some of the model parameters,
especially those that are difficult to obtain experimentally, and study its effects
on the system predictions. More specifically, to find out to which of the model parameters
it is sensitive and need to be found precisely and to which it is not sensitive so that
their approximate values should be sufficient for reliable predictions.
\begin{thebibliography}{00}
\bibitem{Adams} R. Adams;
\emph{Sobolev Spaces}, Academic Press, New York, 1975.
\bibitem{berk96} J. Berkovits, V. Mustonen;
\emph{Monotonicity methods for nonlinear evolution equations},
Nonl. Anal., \textbf{27} (1996), 1397--1405.
\bibitem{bre73} H. Brezis;
\emph{Op\'{e}rateurs maximaux monotones et
semigroupes de contractions dans les espaces de Hilbert},
Math Studies \textbf{5}, North Holland, 1973.
\bibitem {CSS04} O. Chau, M. Shillor, M. Sofonea;
\emph{Dynamic frictionless contact with adhesion,}
Z. angew. Math. Phys. (ZAMP) \textbf{55} (2004), 32--47.
\bibitem{Fr02} M. Fr\'emond;
\emph{Non-Smooth Thermomechanics}, Springer, Berlin, 2002.
\bibitem{HKSS02} W. Han, K. L. Kuttler, M. Shillor, M. Sofonea,
{\it Elastic beam in adhesive contact,}
Intl. J. Solids and Structures \textbf{39}(5) (2002), 1145--1164.
\bibitem{JSS01} L. Jianu, M. Shillor, M. Sofonea;
\emph{A viscoelastic frictionless contact problem with adhesion,}
Applicable Anal., \textbf{80}(1-2) (2001), 233--255.
\bibitem{kut99} K. L. Kuttler, M. Shillor;
\emph{Set valued pseudomonotone mappings and degenerate evolution inclusions,}
Comm. in Contemporary Mathematics \textbf{1}(1) (1999), 87--123.
\bibitem{lio69} J. L. Lions;
\emph{Quelques Methodes de R\'{e}solution des
Probl\`{e}mes aux Limites Non Lin\'{e}aires}, Dunod, Paris, 1969.
\bibitem{PawelDist} P. Marcinek;
\emph{Computational Topics in Nonlinear
Models for Debonding, Approximation of Noisy Signals and Curve Fitting,}
PhD Dissertation, Oakland University, in preparation.
\bibitem{MN17} E. Mazhari, S. A .Nassar;
\emph{A coupled peel and shear stress-diffusion
model for adhesively bonded single lap joints},
ASME J. Manuf. Sci. Eng., \textbf{139}(9) (2017), 091007-1-9.
\bibitem{MOS13} S. Migorski, A. Ochal, M. Sofonea;
\emph{Nonlinear Inclusions and Hemivariational Inequalities:
Models and Analysis of Contact problems},
Advances in Mechanics and Mathematics \textbf{16}, Springer 2013.
\bibitem{NAKS05} S. A. Nassar, K. T. Andrews, S. Kruk, M. Shillor;
\emph{Modelling and simulations of a bonded rod},
Math. Comput. Modelling, \textbf{42} (2005), 553--572.
\bibitem{NS15} S. A. Nassar, K. Sakai;
\emph{Effect of cyclic heat, humidity,
and joining method on the static and dynamic performance of lightweight
multimaterial single-lap joints,} ASME J. Manuf. Sci. Eng.,
\textbf{137}(5) (2015),
051026. doi:10.1115/1.4030080.
\bibitem{NM16} S. A. Nassar, E. Mazhari;
\emph{A coupled shear stress-diffusion model for
dhesively bonded single lap joints,}
ASME J. Appl. Mech., \textbf{83}(10) (2016), 101006.
\bibitem{R11} M. Raous;
{Surface mechanics: facts and numerical models
Interface models coupling adhesion and friction},
Comptes Rendus M\'ecanique \textbf{339}(7D8) (2011), 491--501.
https://doi.org/10.1016/j.crme.2011.05.007
\bibitem{AV10} A. Rodriguez-Aros. J. M. Via\~{n}o;
\emph{Mathematical justification of viscoelastic beam models by asymptotic methods,}
J. Mathematical Anal. Applic. \textbf{370}(2) (2010), 607--634.
https://doi.org/10.1016/j.jmaa.2010.04.067
\bibitem{AKS17} A. Rodriguez-Aros, K. L. Kuttler, M. Shillor;
\emph{Derivationof a 1D model for humidity caused debonding of rods and beams,}
in preparation.
\bibitem{SN16} K. Sakai, S. A. Nassar;
\emph{Failure analysis of composite-based lightweight
multimaterial joints in tensile-shear tests after cyclic heat at
high-relative humidity,}
ASME J. Manuf. Sci. Eng. \textbf{139}(4) (2016), 041007. doi:10.1115/1.4034888.
\bibitem{SCNM17} G. Scarselli, C. Corcione, F. Nicassio, A. Maffezzoli;
\emph{Adhesive joints with improved mechanical properties for aerospace
applications,} Intl. J. adhesion and adhesives \textbf{75}(6) (2017), 174--180.
\bibitem{sh16} M. Shillor;
\emph{Models of debonding caused by vibrations, heat
and humidity,} Chapter 15 in Mathematical Modelling in Mechanics,
Advanced Structured Materials, Vol 69, Ed. F. dell'Isola, M. Sofonea and D.
Steigmann, Springer, Singapore, 2017, p. 233--250.
\bibitem{SST02} M. Shillor, M. Sofonea, J. J. Telega;
\emph{Models and Analysis of Quasistatic Contact Variational Approach},
Lecture Notes in Physics \textbf{655}, Springer, Berlin, 2004.
\bibitem{Sim} J. Simon;
\emph{Compact sets in the space $L^{p}( 0,T;B) $,} Ann. Mat. Pura. Appl.
\textbf{146} (1987), 65--96.
\bibitem{SHS06} M. Sofonea, W. Han, M. Shillor;
\emph{Analysis and Approximations of Contact Problems with Adhesion or Damage},
Pure and Applied Mathematics \textbf{276}, Chapman \& Hall/CRC
Press, Boca Raton, Florida, 2006.
\bibitem{TV96} L. Trabucho, J. M. Via\~{n}o;
\emph{Mathematical Modelling of Rods,}
in Handbook of Numerical Analysis, FEM (Part 2) PG Ciarlet and JL Lions (Eds.)
Elsevier, 1996.
\end{thebibliography}
\end{document}