\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}