Analysis of a debonding model of two elastic 2D-bars

This work establishes the existence of a weak solution to a new model for the process of debonding of two elastic 2D-bars caused by humidity and vibrations. A version of the model was first presented in the PCM-CMM-2019 conference in Krakow, Poland, and was published in (Shillor in J. Theor. Appl. Mech. 58(2): 295–305 2020). The existence of a weak solution is proved by regularizing the problem and then setting it in an abstract form that allows the use of tools for pseudo-differential operators and a fixed point theorem. Questions of further analysis of the solutions, effective numerical methods and simulations, as well as possible controls, are unresolved, yet.


Introduction
This work establishes the existence of a weak solution to a system of partial differentia equations that models the process of debonding, caused by humidity and vibrations, of two elastic 2D-bars.
In many industrial settings, in transportation and everyday life, metallic plates are usually joined by welding. However, when the plates are made of different or nonmetallic materials, they have to be joined with screws or adhesively, or both. This is the case in planes, cars, and other vehicles, where increasingly more parts are made of plastics. One of the issues with car plates or other parts that are glued together is the process of debonding, that is, the process by which the adhesive is loosing its strength over time. This is noticeable in places where there is high humidity and considerable vibrations (driving off the road) and also large variations in temperature. The debonding of plates that are adhesively bonded too often has negative or even catastrophic consequences.
There is a very large mechanical engineering literature on debonding, for obvious reasons, see, e.g., [12,15,17] and the references therein. There is, also, a growing mathematical literature dealing with models of adhesive contact that include the processes of debonding, because of their relations to contact processes. There exist many publications in the Mathematical Theory of Contact Mechanics (MTCM) [6-8, 11, 16, 20, 22] that deal with contact processes with adhesion, and the host of references therein, and note

The dynamics of two bonded 2D-bars
This section presents a highly nonlinear and nonsmooth system that models the dynamics of two elastic 2D-bars that are bonded. However, here we take into account the effects of mechanical vibrations and humidity on the debonding process. The setting is depicted in Fig. 1.
We let u(x, y, t) denote the horizontal displacement of the 2D-bar; and using symmetry and a simple assumption, the vertical displacement of the central horizontal axis w(x, t) depends only on x. Thus, this is actually a "1.5D" system, since w depends only on x. Below, we use the subscripts x, y, t to denote partial derivatives. The model we study describes the debonding process of two such 2D-bars, caused by mechanical vibrations and humidity diffusion. Our main interest lies in the dynamics of the bonding field β = β(x, t) that is defined on the contact surface, B in Fig. 1(red).
The setting is depicted in Fig. 1. The top 2D-bar occupies the domain 1 = {0 < x < l 2 , 0 < y < 2h} and the bottom 2D-bar occupies 2 = {l 1 < x < L, -2h < y < 0}, where 0 < l 1 < l 2 < 1, and, for the sake of simplicity, we scaled all the distances so that the length is L = 1. The quantities with index 1 refer to the top bar and those with 2 to the bottom bar, thus for i = 1, 2, we refer to bar i . The horizontal displacements u i = u i (x, y, t) are defined on i . The vertical displacement of the central line w 1 = w 1 (x, t) is defined on 1 ∩ {y = h}, and w 2 = w 2 (x, t) on 2 ∩ {y = -h}. We denote by ρ i the 2D material density, E i the Young modulus, and G i the shear modulus of bar i = 1, 2. These are positive mechanical, found from measurements.
The debonding process takes place on the surface B = {l 1 < x < l 2 , y = 0} that is occupied by the glue (red segment in Fig. 1), where the bonding field β = β(x, t) and the humidity field η = η(x, t) are defined. We assume that there are no body forces acting in the system; the vertical tractions p + = p + (x, t) act on the top boundary of 1 while it is clamped on D , and free at x = l 2 and on the bottom (0 < x < l 1 ). Similarly, the vertical tractions p -= p -(x, t) act on the bottom boundary of 2 , it is free on x = l 1 , and on the top for l 2 ≤ x < 1, and, finally a traction p * (t) acts at x = 1.
The bonding of the two bars is modeled by the bonding field β = β(x, t), which is defined on B . It represents the fraction of active surface bonds and has the character of a damage variable (see, e.g., [2,3,8,20]), namely, 0 ≤ β ≤ 1; indeed, when β = 1, all the adhesive bonds at the point are active; when β = 0, all the bonds are severed; and when 0 < β < 1, the fraction β are active bonds, which describes the decrease in the load carrying capacity of the adhesive.
The debonding process on B is controlled by the debonding source function , which depends on the bonding field, humidity, and the mechanical tractions. Mathematically, we write it as (x, t) = |u 2u 1 |, (w 1w 2 ), β, η , on B . This function contains the information about the processes and must be determined in conjunction with experimental data, as was noted in [13], but in a simpler setting. We assume that the debonding process also depends on diffusion as the field at a point is affected by what happens in neighboring points. To take into account the constraints 0 ≤ β ≤ 1 and w 2 (x, t) + h ≤ w 1 (x, t)h on B (so there is no interpenetration there), we proceed as follows. We let I [0,1] (β) and I [0,∞) (r) be the indicator functions of the intervals [0, 1] and [0, ∞), with subdifferential ∂I [0,1] (β) and ∂I [0,∞) (r), respectively.
To prevent the top bar from penetrating the lower bar, we add the term to the right-hand side of the equation of motion for w i , for i = 1, 2. This guarantees that when w 1 (x, t)w 2 (x, t) -2h > 0, so there is a gap between the bars, this term vanishes. When the top bar tries to interpenetrate the bottom bar, w 1 (x, t)-w 2 (x, t)-2h = 0, the term generates enough resistance force to prevent the penetration. More precisely, there exists a fictitious force ξ ∈ ∂I [0,∞) (w 1 (x, t)w 2 (x, t) -2h) such that ξ = 0 when w 1 (x, t)w 2 (x, t) -2h > 0, and when w 1 (x, t)w 2 (x, t) -2h = 0 then -ξ > 0 is pushing the top bar up, and ξ < 0 pushing the bottom bar down, preventing penetration. We note that ξ need not be a function, only a distribution, which complicates the mathematical analysis considerably. The evolution equation of the bonding field, see, e.g., [8,20], is assumed to have the form of a constrained diffusion equation, or set-inclusion, where k β is the bonding diffusion coefficient, and the subdifferential on the right-hand side of (2.1) guarantees that 0 ≤ β ≤ 1.
Next, we assume that the humidity function evolves by nonlinear diffusion from the edges (x = l 1 , y = 0) and (x = l 2 , y = 0), where the adhesive layer is exposed to the environment. The coefficient of diffusion D is assumed to depend on the bonding field, stress, and humidity density in the adhesive layer. Thus, The boundary conditions are were η L (t) and η R (t) are given humidity densities at x = l 1 and x = l 2 , respectively. In most applications it is reasonable to assume that η L = η R = η amb , where η amb is the ambient humidity density. However, these may be different and both varying with time. Also, 0 ≤ η 0 (x) is the initial humidity density in the adhesive, which in practice is likely to be unknown. The model for the debonding of two adhesively joined 2D-bars caused by vibrations and humidity is the following.
Problem 1 Find the functions u i , defined on i , i = 1, 2; the functions w 1 defined on 1 ∩ {y = h} and w 2 defined on 2 ∩ {y = -h}; and the functions β, η defined on B , such that the following equations and conditions hold for 0 < t ≤ T: The boundary conditions for β and η are 8) and the initial conditions The balance of vertical and horizontal tractions on B , due to the adhesive, σ βv and σ βh , is assumed to be Here, K βv and K βh are the adhesive stiffnesses, which depend on the humidity. Next, the tractions acting on the system are and zero tractions everywhere else, except for D , where Finally, the initial conditions for u i , w i are We note that the interesting aspects of the physics of the process enter via the debonding source function and the humidity diffusion coefficient D, both of which must be found from experimental data in conjunction with numerical optimization. This is a major issue in the applications of the model. However, even qualitatively simple forms of and D may provide very useful information in concrete applications.

Existence
The system is highly nonlinear and nonsmooth, so it is unlikely to have any classical solutions, so, we reformulate it and look for weak solutions. Moreover, generally, the solutions are unlikely be unique, therefore, we deal only with proving existence. However, we replace the the nonpenetration condition in (2.5), represented by the subdifferential for w i , which physically is an idealized condition, with a more physically reasonable one. For the sake of simplicity, we rescale the variables so that ρ 1 = ρ 2 = 1. The system is rather complex, involving six dependent variables, and so we do not provide all the details in the proof of existence. We just outline the approach, pointing out that existence can be proved in this manner.
We first introduce the spaces and subspaces where we seek the various functions, others being chosen by analogy. Let, which are also the spaces of test functions. Here, H 1 denotes the usual Sobolev space over the relevant domain. We let H 1 = L 2 ( 1 ) and H 2 = L 2 ( 2 ).
We start with the equation for u 1 , and let ψ ∈ V 1 . It follows from (2.4), using ψ as a test function and integrations by part, that To write this equality as an abstract equation, we define the operator L 1 by and note that it is linear, coercive, and bounded as a map from V 1 to V 1 . Next, let B 1 be a Lipschitz function, given by Therefore, we write the equation for u 1 , see (2.4) (i = 1), in abstract form as Here, the prime denotes a weak time derivative.
The same reasoning applies to u 2 , see (2.4) (i = 2), with ψ ∈ V 2 , so we let and note again that it is linear, coercive (up to a constant), and bounded as a map from V 2 to V 2 . Next, let B 2 be a Lipschitz function, given by Then, the abstract form of (2.4) (i = 2) is We note that, thanks to the compact embedding theorems in Sobolev spaces, if w in → w i in L 2 (0, T; W i ) and w in is bounded in L 2 (0, T; W i ), one may conclude that We turn to the equations for w 1 and w 2 . Since the Signorini condition of perfect nonpenetration enters as a body force into the equations, rather than as a boundary condition, the standard methods used previously, do not apply. This may reflect the fact that such a fictitious perfect force does not make sense physically and, indeed, there is no such thing as a perfectly rigid body. Therefore, we replace it with a more physically realistic condition. To that end, we introduce the function R such that R(s) = 0 if s ≥ 0 and R(s) = s if s < 0. Then, we replace ∂I [0,∞) (s) with 1 ε R(s), assuming that 0 < ε 1. Thus, we allow some interpenetration, but penalize it heavily. Then we replace the term (w 1w 2 ), which originally is nonnegative, with |w 1w 2 |, which is nonnegative.
To proceed, we multiply (2.5) (i = 1) with ψ ∈ W 1 such that ψ(l 2 ) = 0, integrate over 1 , use integration by parts, note that ψ and w 1 do not depend on y, and use algebraic Here, we used the fact that ψ(l 2 ) = 0, hence, 1 u 1xy ψ dx dy = - This is an important observation because in the weak solutions u xy is not well defined. We define the operators L 1 and B 1 by Thus, the abstract formulation of (2.5) (i = 1) is Similarly for (2.5) (i = 2), using ψ ∈ W 2 , such that ψ(l 1 ) = ψ(1) = 0, we define L 2 and B 2 by Thus, the abstract formulation of (2.5) (i = 2) is We note that the L i are coercive, bounded, and linear.
To complete the mechanical model, the initial conditions need to be added. We do it below, when we discuss the smoothed mechanical problem.
Next, we consider the set-inclusion for the bonding field β, with the initial condition β(x, 0) = β 0 (x), for x ∈ [l 1 , l 2 ]. The boundary conditions are β x = 0. Let ψ ∈ H 1 (l 1 , l 2 ), we define the operator L β : H 1 (l 1 , l 2 ) → H 1 (l 1 , l 2 ) by and note that the operator L β , is the subgradient of the function Next, multiplying the equation with ψ as a test function, and using integration by parts, we obtain the weak form We assume that is bounded and Lipschitz in all its variables. Also, we let H = L 2 (l 1 , l 2 ). We have the following preliminary result.
It follows from the theory of maximal monotone operators (see Brezis [5]) that for each n, there exists a solution β n to the problem. This solution satisfies β n ∈ L 2 ((0, T) : L 2 (l 1 , l 2 )) and β nxx ∈ L 2 ((0, T) : L 2 (l 1 , l 2 )). Multiplying both sides by β n and integrating yields (H ≡ L 2 (l 1 , l 2 )) where C is a positive constant that depends on . Hence, we obtain the bound Similarly, multiplying by β n and integrating yields where C (β 0 ) depends also on β 0 . Therefore, we find Thus, there is a subsequence, still denoted by β n , such that β n → β weakly in L 2 (0, T); H , Passing to the limit as n → ∞, we obtain the desired solution. However, we lose regularity of β because of the subgradient term. The weak solution is unique because of monotonicity considerations.
We turn the equation for η (see (2.7)) with the boundary conditions (2.8), η(l 1 , t) = η L , η(l 2 , t) = η R . To obtain zero boundary conditions, we let η = η -( x-l 1 l 2 -l 1 η R + l 2 -x l 2 -l 1 η L ), and the equation in terms of η has the form Thus, without loss of generality, we assume that the initial-value problem has zero boundary conditions, and the equation is replaced with where η ∈ H 1 0 (l 1 , l 2 ), and the hat is omitted. Again, we assume that D, are Lipschitz continuous and bounded. Also, we assume that D(·) ≥ δ > 0, so the diffusion coefficient does not vanish and thus the problem is nondegenerate. Assume, in addition, that the functions a, b ∈ L ∞ ([l 1 , l 2 ] × [0, T]), and a(t, x) ≥ δ > 0, and define the operator A as Proof This is a standard result since A is monotone, bounded, hemicontinuous, and coercive as a map from L 2 ((0, T) : H 1 0 (l 1 , l 2 )) to its dual space. To include the dependence of D on η, one can use the Schauder fixed point argument or to retard in time the functions on which the diffusion coefficients depend and then passing to a limit using compactness.
We conclude that the initial value problem for η has a weak solution. We now revisit the mechanical system and modify it by adding viscosity, which acts as smoothing operator of the solutions, and is needed for the a priori estimates below. The modified system is along with initial conditions. Here, γ is a small positive number, representing material viscosity. Also, K i ,K i are suitable linear maps from an appropriate Sobolev space S to its dual space S , which ensure that u i and w i have two square integrable derivatives. This Sobolev space is a closed subspace of H 2 ( i ) such that u ix = 0 on the boundary of I i , and is dense in V i . Then, the K i can be considered as Riesz maps from these new spaces to their duals. We assume sufficient regularity on the initial conditions to justify all of the following argument and, in particular, the initial data u i0 lie in S. The regularized abstract form of the full model follows.

Problem 4
Find the functions u i , w i , η and β, (for i = 1, 2) such that The initial conditions are: We have the following existence result.

Theorem 5
There exists a solution to the system (3.8), for each 0 ≤ γ . The theorem guarantees the existence of weak solutions to the equations for two bars that are in adhesive contact, provided the subgradient acting on the bars for w i is replaced with a penalization operator approximating the subgradient on B . Thus, we allow for γ = 0, so the materials are purely elastic. We only give a sketch of the proof since there are many dependent variables and the considerations are mostly routine. The proof is based on the following pair of compactness results, due to [14] and [21], and used below.
Theorem 6 Let q > 1 and let E ⊆ W ⊆ X where the injection map is continuous from W to X and compact from E to W . Let S be defined by ∞ (a, b, E)  Proof We provide a sketch of the proof. It is based on time delay. We start with the first two equations in (3.8) for u i when w 1 , w 2 are known, so, we replace the time variable in w i by the delayed time th and for t < h, we use the initial conditions. Therefore, w 1 , w 2 are known on the time interval [t, t + h]. Similarly, in the second pair of equations for w 1 , w 2 , we use w 2 (th) in the equation for w 1 , and w 1 (th) in the equation for w 2 . In the equation for η, we replace t with th in all the variables in D(|u 2u 1 |, |w 1w 2 |, β, η) and use the the initial condition if t ≤ h. Next, we do the same in all the variables in (|u 2u 1 |, |w 1w 2 |, β, η), (t, x, |u 2u 1 |, |w 1w 2 |, β, η). Then, routine considerations, tedious a priori estimates, along with Proposition 2, yield a solution to this time delayed regularized system of equations which has the property that u ix , u iy vanish on the boundaries and w ix , too. Then, we obtain the necessary estimates on this solution by multiplying through by -E i u ixx -G i u yy in the equations for u i and use integration by parts and the first term to argue that there is an estimate on u it (t) V i which is independent of γ . A similar process is applied to the equations for w i .

Then, S R is bounded in L
Next, having the estimates, we pass to the limit as h → 0, using Theorems 6, and 7, and the embedding theorems in Sobolev space, and obtain a solution to the regularized (2022) 2022: 15 Page 13 of 14 system without the time delay, along with estimates on u it (t) V i , which are independent of γ , and similar estimates for w i . Finally, we use these theorems to pass to the limit in a subsequence for which γ → 0. This yields the existence of a weak solution to the original problem in which the idealized subgradient term in the equations for w i are replaced with a penalized version of this subgradient. The passage to the limits is straightforward since the higher-order operators are linear and the boundary operators and other functions which are not linear depend on lowerorder terms and can be estimated using compactness considerations in the two theorems.
We consider the uniqueness of the solution very unlikely, and so we skip it. Now that the existence of weak solutions has been established, some questions of interest for further research are: • Perform additional analysis of the solutions, especially about their regularity.
• Construct an efficient and convergent algorithm for computer simulations.
• Find numerically a typical debonding process behavior.
• Study how the debonding depends on the excitation frequency.
• Get an estimate of the degree of debonding from the shift in the vibrations spectrum.
• Correlate the computer experiments with experimental data to find a possible structure for the debonding source function and the humidity diffusion coefficient D.