A general thermodynamical model for adhesive frictional contacts between viscoelastic or poro-viscoelastic bodies at small strains
Tom\'a\v{s} Roub\'i\v{c}ek

TL;DR
This paper introduces a comprehensive thermodynamical model for adhesive and frictional contacts between viscoelastic or poro-viscoelastic bodies, including a stable numerical scheme and an extension to porous media with diffusant effects.
Contribution
It presents a unified model for adhesive and frictional contacts in viscoelastic bodies, with a stable, energy-conserving numerical discretization and an extension to porous media with diffusant influence.
Findings
Developed a semi-implicit energy-conserving numerical scheme
Extended the model to include porous media with diffusant effects
Demonstrated stability and convergence of the proposed method
Abstract
A general model covering a large variety of the so-called adhesive or cohesive, possibly also frictional, contact interfaces between visco-elastic bodies with inertia considered in a thermodynamical context is presented. A semi-implicit time discretisation which conserves energy, is numerically stable and convergent, and which advantageously decouples the system is devised. An extension to porous media with adhesive contact influenced by a diffusant concentration is devised, too.
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A general thermodynamical model for adhesive
frictional contacts between viscoelastic or
poro-viscoelastic bodies at small strains
Tomáš Roubíček
Institute of Thermomechanics, Czech Academy of Sciences
Dolejškova 5, CZ-182 08 Praha 8, Czech Republic
Abstract
A general model covering a large variety of the so-called adhesive or cohesive, possibly also frictional, contact interfaces between visco-elastic bodies with inertia considered in a thermodynamical context is presented. A semi-implicit time discretisation which conserves energy, is numerically stable and convergent, and which advantageously decouples the system is devised. An extension to porous media with adhesive contact influenced by a diffusant concentration is devised, too.
2010 Mathematics Subject Classification: 35R45, 65K15, 74A55, 74F10, 80A17.
Keywords: Contact mechanics; tribology; thermomechanics; poromechanics; dry friction; adhesion; evolution variational inequalities; numerical approximation.
1 Introduction
Contact mechanics is an important and widely developed part of continuum mechanics of solids, intensively studied during past decades from the modelling, computational, and analytical aspects, too [58, 60, 59, 66, 68]. We focus on a class of models for contacts with adhesion (cohesion) and friction in the thermodynamical context. Our goal is to formulate a rather general multi-purpose model unifying a relatively wide menagerie of such particular adhesion-friction models and to discuss its thermodynamics, outline its conceptually efficient discretisation amenable for calculation of vibrations/waves.
It is certainly impossible to present in a comprehensive way the relevant literature to the part of contact mechanics covered by the general model developed in the present paper. Anyhow, from the engineering point of view, a detailed analysis of related variational formulations of interface cohesive models was developed in [21]. Another related and thermodynamically based approach for general imperfect interfaces was proposed and numerically studied in [24, 30, 29], considering e.g. an interface temperature different from the temperatures of adjacent bulks, cf. also [13, 14, 15] for analytical considerations. Nevertheless, except [15], friction dissipation was not taken into account in these works. Coupling of cohesive models and friction contact was developed after the pioneering work [62], which assumed the friction is activated after complete cohesive failure, in a series of relevant contributions [2, 9, 19, 22, 27, 36, 37, 41, 42, 53, 57, 61, 64] among others, providing, for isothermal configurations, a continuous transition from the cohesive state to the frictional contact state. However, no rigorous proof of the solution existence or the convergence of a numerical scheme to the exact solution was provided in these works. From the mathematical point of view, overviews of existing partial models in literature can be found, e.g., without Coulomb friction () in [20, 49, 52] (isothermal), without inertia/friction and interfacial plasticity but with fracture-mode-mixity sensitive debonding (sometimes called also delamination) in [31] (isothermal) or [43, 44] (anisothermal); see also a review chapter [50].
Models of adhesive contacts are to a great extent inspired by bulk damage possibly combined with plasticity. The role of a scalar damage variable is then replaced by an interface damage (a concept invented by M. Frémond [25, 26], also called debonding or delamination, the term adopted hereinafter) variable and plasticity can be transferred to a plastic slip, cf. Remark 2.4 for more details. Like in the bulk models, there are many options:
(i)
delamination can influence the elasticity of the adhesive (through decaying elastic moduli) or the plasticity (through decaying plastic yield stress for the plastic slip), or both;
(ii)
vice versa, plasticity can influence delamination indirectly through influencing the stress and displacement jump, or directly through influencing activation threshold for delamination;
(iii)
delamination evolution can be either unidirectional or reversible (i.e. admitting a so-called healing);
(iv)
delamination and interfacial plasticity can be considered rate-independent (and then choosing an appropriate concept of solution becomes a vital part of the model itself and should carefully be considered, cf. [38]) or rate-dependent (visco-plasticity, viscous delamination), i.e. 4 options altogether, or more if also healing is involved either as rate-dependent or rate-independent;
(v)
interfacial plasticity can be either with hardening or without hardening;
(vi)
length scale (gradients) can be considered in interfacial damage or/and plasticity.
In contrast to the mentioned bulk models, the interfacial analog can additionally be combined with friction. We can then identify other options (beside the trivial option just ignoring friction):
(vii)
friction acts on the displacement velocity, or can be realized through the interfacial plastic slip;
(viii)
the contact is unilateral or bilateral, in the former case the adhesion being mode sensitive or insensitive, possibly in the normal-compliance variant counting with a small penetration (microscopically related to various asperities, cf. e.g. [63] for a detailed discussion).
Of course, as a general feature in (contact) mechanics, there is still an option as far as temperature variation and related heat transfer concerns:
(ix)
either to ignore it by considering the model isothermal, or to consider heat transfer in the bulk and through the contact interface but ignoring the heat capacity and the heat transfer in the adhesive itself, or even to involve this latter heat problem in the adhesive, too.
In addition, various non-linearities can be considered at various parts of the models, leading e.g. to various cohesive-contact models etc. Eventually, there are standard options as far as the bulk models: purely elastic or visco-elastic with various rheologies (possibly involving creep) or some additional inelastic processes like plasticity, purely quasistatic models versus inertial effects involved, simple or nonsimple materials, cf. [44], linear or nonlinear materials, purely isothermal versus anisothermal with or without thermal expansion, etc.
Contact mechanics has many applications in engineering or in (for example geo-) physics, and various combinations of the above identified options (i)–(ix) may serve for various purposes. Of course, not all combinations give sensible models or are amenable for a rigorous analysis supporting also executable numerical schemes with guaranteed stability and convergence.
The above outlined rich menagerie of inelastic-contact models however should not create a false impression that it covers everything. In particular, in this paper we ignore phenomena like fatigue or wear, and we confine ourselves to small strains and, in addition, we also assume small displacements (and thus a prescribed contact zone, avoiding e.g. contact of rolling bodies or a varying Hertz’ contact zone).
Moreover, there is also a lot of contact-mechanical models based on hemivariational inequalities involving nonsmooth nonconvex functions, cf. e.g. [8, 4, 39]. Often, this mathematically involved and numerically more difficult setting results rather from lacking of suitable internal variables. In contrast, the definite advantage of our general model is that it relies on classical variational inequalities involving convex functions only.
The plan of this paper is the following: In Section 2 we formulate the general model coupling the adhesive layer with the viscoelastic bulk. Then, in Section 3 we devise its “monolithic” formulation treating the bulk and the adhesive in a unified way, and discuss its energetics and thermodynamical consistency. Then, in Section 4, a semi-implicit energy-conserving time discretisation allowing for efficient computer implementation is devised and its numerical stability and convergence is outlined. Eventually, in Section 5, we outline the extension towards poro-viscoelastic bulk materials and adhesives, devising again a monolithic formulation. For readers’ convenience, let us summarize the basic notation used in what follows in Table 1.
2 The general model
We consider a bounded Lipschitz domain split to two (or more) bodies by some internal interface(s) , cf. Figure 1. We use a concept of an “interface” temperature in an infinitesimally thin but still heat conductive adhesive on the -dimensional contact interface which glues two -dimensional visco-elastic bodies, cf. [13, 14, 16] or [38, Sect.5.3.3.3] for such a model, or in combination with a friction also [15]. For a relation to models which avoid this interface-temperature concept see Remark 2.2. Rather for notational simplicity, we will consider no more complicated rheology than the linear Kelvin-Voigt one (so that, in particular, no additional internal variables in the bulk will be considered).
The unit normal vector to is denoted by . Hence on , while we denote just by on , Fig. 1. A convention for definition of jumps across for traces of fields defined in and is used as well as their normal and tangential components, namely
[TABLE]
Thus, e.g. \mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}>0 means that the boundaries of and at separate (we refer to such situation as opening) whereas for \mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}<0 overlapping of these subdomains takes place.
We further consider the boundary of split to two parts, and , where the Dirichlet and the Neumann boundary conditions will be prescribed. Considering a fixed time horizon , we will use the notation
[TABLE]
One of the main ingredient of the model will be the force balance on the contact interface :
[TABLE]
where denotes the pair of the traces of a possibly discontinuous temperature field from both sides of the interface , and , and . Thus, the expression means with and denoting temperature fields in the two domains adjacent to the interface and and the bulk-to-adhesive heat coefficients from the corresponding bulk domains.
In (2.3f-h), “” denotes the gradient over -dimensional surface (interface) , i.e. the tangential derivative defined as for defined around with (or, equally, ); note that takes into account the profile of in a neighbourhood of although eventually a.e. depends only on the trace on . Further the corresponding surface divergence “” is the trace of ; thus is the so-called Laplace-Beltrami operator.
In the bulk, we consider the standard thermo-visco-elasticity with thermal expansion consisting of the force balance coupled with the heat-transfer equation:
[TABLE]
with being the applied bulk force (typically gravitational); let us say that we consider constant over each subdomain so that the reference temperature occurs explicitly only in (2.3) but not in (2.4a).
Further, we supplement (2.4) with standard boundary conditions counting also with the heat transfer from the adhesive:
[TABLE]
with denoting the normal to the -dimensional boundary of the interface , cf. Figure 1.
In addition, we consider an initial-value problem by prescribing the following initial conditions
[TABLE]
The conceptual rheological diagram corresponding to the part of the model, namely (2.3)–(2.4a), is shown in Fig. 2. The modelling assumption is that, under compression, the normal compliance dominates the response of the adhesive, i.e. \gamma_{\scriptscriptstyle\textrm{C}}({\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}})\gg\gamma_{{}_{\rm N}}(\alpha,{\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}}) for {\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}}<0, otherwise rather should enter (2.3b-d), cf. also Remark 4.6. The standard symbols are used for elastic and viscous elements (springs and dashpots), respectively, and for the dry-friction-type elements.
Remark 2.1** (More general and .)**
Some applications may need rather and . The -dependent interface heat capacity needs generalization of (3.10b) below so that and yields some adiabatic-like heat source/sink contribution to the heat equation, namely .
Remark 2.2** (Models without interfacial temperature:
and .)**
A particular case arises when the adhesive has zero heat conductivity and capacity . Then, it is reasonable to consider and the interfacial temperature can be eliminated. This conceptually reflects the modelling assumption about zero thickness of the adhesive. Then the heat-transfer problem (2.3h) degenerates to the transient conditions:
[TABLE]
on , which reflects that the heat generated by delamination is distributed with equal proportions and into the two subdomains adjacent to . Cf. also [38, Rem. 5.3.23]. Such a model has been considered e.g. in [43]. This enhances the concept of a so-called weakly conducting (Kapitza) interface, cf. [30], which would be obtained for when the normal heat flux would be continuous across .
Remark 2.3** (Models with only interfacial temperature
.)**
In contrast to Remark 2.2, sometimes an idea of ‘an external constant temperature reservoir’ from [10] is a relevant approximation – then is considered as a fixed constant and the heat-transfer problem (2.4b)–(2.5c,d) is omitted, cf. also [45].
Remark 2.4** (Role of the plastic slip .)**
The plastic slip may play various roles:
To distinguish delamination in Mode I and Mode II under gradually increasing load as in [40, 47] and, after the (irreversible) delamination process is completed, it will stop evolving, or
To facilitate healing (=reversible delamination) and, allowing for free tangential shift in the delaminated state, to make sure that the healing will aim the shifted configuration instead of the original one, cf. [49].
If \mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm T}}}}{[[u]]_{{{}_{\rm T}}}}{[\![u]\!]_{{{}_{\rm T}}}}{[\![u]\!]_{{{}_{\rm T}}}}\sim\pi can be expected during the whole process, an approximation of friction can be realized through the activation threshold for evolving .
3 A monolithic form and thermodynamics of the model
Let us briefly present the thermodynamics of the boundary-value problem (2.3)–(2.5). This justifies the problem physically. Moreover, adapting [38, Sect. 5.3] for our situation, a “monolithic” (in the sense of simultaneous treatment of the -dimensional interface and -dimensional bulk) description will show a (relatively) simple structure of the problem, cf. (3.9) below.
First, the underlying overall Helmholtz free energy \varPsi=\varPsi(e,{\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]}},\pi,\alpha,\theta_{\scriptscriptstyle\textrm{A}},\theta_{\scriptscriptstyle\textrm{B}}) has a bulk (=volume) and an interface parts, i.e.
[TABLE]
The terms represent a (temperature-dependent) energy related to debonded adhesive, based on the phenomenon that every damage microscopically means more microcracks/voids which always bears some additional energy.
The typical choice for and describing a linearly responding (for a fixed delamination variable ) elastic adhesive
[TABLE]
while the typical choice for is a unilateral normal-compliance:
[TABLE]
Yet, it should be mention that the normal compliance is more natural is our a setting than the Signorini condition which is a bit artificial because, when the adhesive is rigid, the issue is of fracture and not debonding.
The other underlying ingredient of the model is the overall pseudo-potential of the dissipative forces which also has a bulk and a contact-interface contributions and , namely:
[TABLE]
where
[TABLE]
the -term representing the specific dissipation rate of the delamination process on the contact boundary while the quadratic terms describes viscosity of the adhesive. Furthermore, the overall kinetic energy is the quadratic functional of bulk velocity , i.e.
[TABLE]
with the mass density. Eventually, we still need the (linear) functional of external mechanical loading which, after the standard transformation to time-constant (homogeneous) Dirichlet boundary condition by the shift with being an extension of from (2.5a) on , is given by
[TABLE]
For simplicity and without excluding interesting applications, we assume that and are far from each other so that one can assume and thus this shift transformation does not influence the flow rules on .
In terms of these functionals, using also the linear operator
[TABLE]
the mechanical part can be written as:
[TABLE]
where we have used the notation for an abstract “gradient/difference” operator and a heat-conductivity operator defined respectively as
[TABLE]
where denotes the surface Lebesgue measure on .
The inclusions in (3.9a,b,c) are related with nonsmoothness of if and if , and if a_{1}({\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]}},\alpha,\theta_{\scriptscriptstyle\textrm{A}},\cdot) is nonsmooth, respectively. The abstract heat-transfer equation (3.9d) in the notation (3.10) is to be understood in the weak form as:
[TABLE]
for any with , where the bilinear form is defined by
[TABLE]
with and , and with and .
Introducing the specific entropy , we can write the abstract entropy equation in the form:
[TABLE]
This means two coupled equations for the interfacial entropy and the bulk entropy :
[TABLE]
here and are the heat-production rates in the adhesive and in the bulk due to the dissipative process, i.e. . Considering also the boundary conditions (2.5c) reflected in and (2.5d), from (3.14b) we obtain the overall entropy balance in the bulk as
[TABLE]
where the last integral means with again meaning the traces of temperature from the two regions adjacent to at a spot in question, while from (3.14a) we obtain the overall entropy balance in the adhesive
[TABLE]
Altogether, summing (3.15) and (3.16) and using the elementary algebra , we obtain the overall entropy balance
[TABLE]
Considering thermally isolated system (i.e. ), we can see that the overall entropy is nondecreasing in time, which represents the 2nd law of thermodynamics. Of course, we relied on positive-definiteness of and and non-negativity of , cf. (4.5) below. Furthermore, we relied on non-negativity (or rather positivity) of temperatures and , which can be ensured by suitable initial/boundary conditions in combination with the adiabatic terms and in (2.3f) and (2.4b), which may alternate sign (and hence cause both heating and cooling) but such possible cooling is “switched off” if the temperature approaches zero. In other words, our system is consistent also with the 3rd law of thermodynamics.
The mere mechanical energy balance can be obtained by testing (3.9a,b,c) respectively by , , and . We first realize the special form of
[TABLE]
where is purely mechanical part of the free energy with defined in (3.8), the linear functional {\boldsymbol{b}}$$(\boldsymbol{E}u,\pi,\alpha):\boldsymbol{\theta}=(\theta_{\scriptscriptstyle\textrm{A}},\theta_{\scriptscriptstyle\textrm{B}})\mapsto-\int_{\mathchoice{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}}b_{0}(\alpha)\theta_{\scriptscriptstyle\textrm{A}}\,\mathrm{d}S-\int_{\varOmega{\setminus}\mathchoice{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}}\theta_{\scriptscriptstyle\textrm{B}}\mathbb{B}{:}e\,\mathrm{d}x facilitates the thermo-mechanical coupling, and where the purely thermal part is . The mentioned test then gives:
[TABLE]
where R(\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]},\alpha,\boldsymbol{\theta};\boldsymbol{E}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}},\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\pi}}},\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\alpha}}})=\langle{\boldsymbol{r}},(1,1)\rangle=\langle\partial_{(\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{e}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{e}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{e}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{e}}},{\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}\big{]}\big{]}}{[[\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}]]}{[\![\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}]\!]}{[\![\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}]\!]}},\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{\pi}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{\pi}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{\pi}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{\pi}}},\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{\alpha}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{\alpha}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{\alpha}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{\alpha}}})}\varXi(\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]},\alpha,\boldsymbol{\theta};\boldsymbol{E}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}},\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\pi}}},\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\alpha}}}),(\boldsymbol{E}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}},\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\pi}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\pi}}},\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\alpha}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\alpha}}})\rangle denotes the overall dissipation rate.
Denoting by the specific internal energy and by the specific free energy, we have the Gibbs’ relation . We can write the abstract specific-internal-energy balance . We denote the overall internal energy by
[TABLE]
where is the total heat energy. The total-energy balance can be obtained by testing the heat equation by 1, i.e. substituting into (3.11) and summing it with (3.19). Realizing that , it gives
[TABLE]
4 A semi-implicit fractional-step time discretisation
An efficient way for time-discretisation guaranteeing numerical stability and convergence and, after a further spatial discretisation, allowing for a constructive solution is to decouple the system by a suitable fractional-step method, in engineering literature usually referred to as a staggered scheme. Beside being computationally amenable and efficient conceptual algorithm, it can simultaneously served analytically to a proof of a mere existence of a weak solution. For it, one can benefit from the monolithic formulation and use abstract results that are essentially already available in literature (in particular in [32, 38]) with possibly a little modification, combining the fully-implicit time-discretisation usually used for purely analytical purposes with some other [48], more suitable to capture dynamical effect, as used also below. In the isothermal case, also other abstract-operator-based techniques available for dynamic frictional normal-compliance contact between visco-elastic bodies exist, cf. [5, 34, 35] and likely they can be expanded if the tangential slip is involved, too.
The general strategy for splitting the variables is to respect the assumed separate convexity (or at least semi-convexity) of the free energy and additive splitting of the dissipation potential . Here, we can apply splitting of the state to and and , assuming that
[TABLE]
In addition, we introduce a new variable in the position of a velocity, and in the limit indeed satisfy . This transforms the 2nd-order system into the 1st-order one. When combined with a mid-point time discretisation of the Crank-Nicolson type, it allows for energy-conserving scheme and for varying time step, which just for notational simplicity is considered fixed, denoted by , so that equidistant partitions of the time interval are used. Let us point out that the fully implicit discretisation of the inertial term, i.e. , would serve for analytical purposes similarly good but computationally it would lead to a scheme with spurious numerical attenuation which practically does not allow for computation of vibrations or waves.
This mid-point formula applied to the dynamical equations transformed into the form of the 1st-order system can be understand as a particular case of the celebrated Hilber-Hughes-Taylor formula [28], sometimes called a Simo’s scheme [56, Sect. 1.6]. Its combination of fractional-step split and 1st-order-system transformation with mid-point formula, called sometimes the Yanenko formula [67], leads to a three-step decoupled scheme, written in terms of the notation from (3.9) as follows. Having, , we first solve:
[TABLE]
In (4.2a) and (4.2d), we used the notation for the time-difference defined as
[TABLE]
assuming that is continuously differentiable. If is a quadratic form with symmetric, then gives exactly the midpoint difference. Cf. also [48] for such a general discretisation scheme in the isothermal variant. Here, if (3.2) is adopted, it is important that the dependence of on is quadratic up the normal-compliance boundary term (3.3) which however depends on the scalar-valued normal displacement, so that the fraction in (4.3) has a good sense, and also the dependence on is allowed generally nonlinear because it is scalar-valued like the mentioned normal displacement. For example, if the dependence on would be quadratic, then \varPsi^{\circ}_{\alpha}(\boldsymbol{E}u_{\tau}^{k},\pi_{\tau}^{k},\,\cdot\,,\theta_{\scriptscriptstyle\textrm{A},\tau}^{k-1})\big{]}(\alpha_{\tau}^{k},\alpha_{\tau}^{k-1}) in (4.2d) would simplify to .
Note also that we used a regularization making from (4.2f) bounded in terms of the rates (=the time differences) so that is bounded in an -space (in contrast to which altogether is bounded only in an -space). Similar remark concerns also . Of course, this regularization is devised to disappear in the limit for so that (4.2) has a capacity to approximate really the original continuous system (3.9), which is indeed ensured by the regularizing positive coefficients , , , , and , which are of different physical dimensions.
The decoupled system (4.2) is to be solved recursively for , starting for with
[TABLE]
Speaking generally, desired attributes of any efficient and reliable discretisation are the following:
(A1)
existence and algorithmic amenability of solutions for the discrete problems,
(A2)
energetic consistency (here even energy conservation) for any discretisation (i.e. here for any the time step ),
(A3)
numerical stability (= a-priori estimates) for any discretisation (i.e. any time step ),
(A4)
convergence (at least in terms of subsequences) to a suitably defined solutions of the original continuous problem for .
Let us emphasize that not all these attributes are always sufficiently justified in engineering literature.
The full model is very general and satisfaction of the above attributes (A1)–(A4) needs relatively strong assumptions which can be weakened in particular special cases. Considering the general dimension , a rather general variant of the basic data qualification can be cast as follows:
[TABLE]
where is the extension of used in (3.7) and is the set of non-negative reals. Here the symmetry (4.5a) means that , and the same for . We have used the standard notation for the function spaces on and their norms. Namely, denotes the Lebesgue space of measurable functions whose -power in integrable while denotes the Sobolev space of functions which are together with all their -order derivatives in . Moreover, we abbreviate , as usual. Below in (4.14g,h), refers to the dual space. For vector or matrix valued cases, we will write e.g. or etc. Moreover, we use the Bochner spaces of Banach-space-valued functions on the time interval , denoted by .
In fact, (4.14g) requires controlled growth of the boundary conditions, in particular (4.5e), and excludes the mentioned Signorini contact in the dynamical case (where and thus ) while the bilateral contact (3.4) works if the underlying space is replaced by the corresponding closed linear subspace respecting (3.4).
The following assertion concerns the attributes (A1) and (A2). Considering a fixed time step , we define the piecewise-constant and the piecewise affine interpolants respectively by
[TABLE]
Similar meaning has also , , , , etc.
Proposition 4.1** (Existence and energetics of approximate solutions.)**
Let (4.5) holds. Then, for any , there exists the solution to the scheme (4.2)–(4.4). Moreover, the discrete mechanical energy balance, i.e. an analog of (3.19) integrated over time, holds
[TABLE]
and, at least if there is no adiabatic coupling , the analog of the total energy balance(3.20) integrated over time holds at least as an inequality
[TABLE]
for any with . Moreover, the non-negativity of on holds, too.
Sketch of the proof. As to the attribute (A1), note that (4.5a,b,e,g) implies (4.1). Further, it is important that the scheme (4.2) has a variational character for all three sub-problems for (u_{\tau}^{k},{\color[rgb]{0,0,0}v_{\tau}^{k}},\pi_{\tau}^{k}), for , and for . These potentials are coercive and weakly lower semicontinuous, which ensures by the standard direct method that a solution of this system always does exist. Noteworthy, even though the system for the first step (4.2a-c) is nonsymmetric, it indeed has a potential after elimination of , namely
[TABLE]
if would be quadratic (i.e. ) while for the nonquadratic case the primitive of the difference quotients (4.3) are to be involved, cf. [48, Formula (4.4)]. Also, note that, due to the regularization of the heat-sources in (4.2e-g), we can use the conventional -theory for the elliptic boundary-value problems.
Testing the particular equations (4.2a,b,c) respectively by , , and , we obtain
[TABLE]
For this we used, among others, the binomial formulas
[TABLE]
a.e. on and on , respectively. Then, testing (4.2d) by , gives
[TABLE]
Summing (4.10) and (4.11) up, we can enjoy cancellation of the terms and, when still sum for , we obtain
[TABLE]
Thus (4.7) has been obtained.
Next, we test the heat-transfer problem (4.2e-h) by . When added to (4.12), this “nearly” cancels the dissipation -term up to the regularizing -coefficients; i.e. the cancellation would really occur if these coefficients would vanish, while otherwise this regularized dissipation heat is always dominated by the -term in (4.12) and thus it gives the inequality in (4.8).
Eventually, we have to show the non-negativity of a.e. on . This can be performed standardly by the test of the negative part of , still belonging to , hence being a legitimate test function.
Moreover, if (4.1) holds and in addition also
[TABLE]
then all the underlying potentials are convex, which eliminates the usual algorithmic difficulties with finding critical points of nonconvex functionals and thus makes the numerical solution relatively simple. Note that (4.1c) ensures that and , determining in (4.2e) and defined by (3.10b), are monotone so that their primitive functions occurring in the potential underlying for (4.2d-f) are convex, and similarly the contribution coming from (4.2g) is convex (and even quadratic) if (4.13) holds. Often, these minimization problems have (possibly after a Mosco-type transformation [38, Sect. 3.6.3]) a structure of linear-quadratic programming or second-order cone programming problems for which efficient computational algorithms are widely available. If healing is not allowed (thus always ), one can weaken (4.13) by allowing depending on in a way so that \theta_{\scriptscriptstyle\textrm{A}}\mapsto\theta_{\scriptscriptstyle\textrm{A}}\partial_{\alpha\theta_{\scriptscriptstyle\textrm{A}}}^{2}\psi_{\scriptscriptstyle\textrm{\hskip 0.25002ptA}}\big{(}\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]},\pi,\alpha,\theta_{\scriptscriptstyle\textrm{A}}\big{)} is nondecreasing, which still makes the contribution to the underlying potential coming from (4.2g) convex (although not quadratic). In more general cases when (4.13) is not satisfied, the potential of (4.2d-f) may be nonconvex and finding its critical point (not necessarily a global minimum) by some iterative routine is necessary.
Another important attribute of the scheme (4.2) is . This is facilitated by the adiabatic heat sources involving instead of occurring in (4.2a,b) which, on the other hand, needs a bit stronger qualification of the data as far as growth conditions than physically necessary, namely the heat capacities are required to depend on temperature with an at least linear growth, see (4.5c).
As to the attribute (A3), one can imitate the procedure which led to (3.19) and (3.20) when replacing time derivatives by time differences. It allows for optimal data qualification when rather technical interpolation estimates for the adiabatic term is used, cf. e.g. [43, 44]. Alternatively, we use here a simpler based on (4.15) below, although it needs a data qualification (4.5c,d) which is a bit stronger in comparison what would be needed when used the mentioned interpolation.
The general physically motivated strategy for (A3) is to get the mechanical and the heat energies bounded uniformly in time and the overall energy dissipated during the process bounded, both uniformly in , too. Based on estimation of , we have:
Proposition 4.2** (Numerical stability)**
Let (4.5) hold. With and independent of , the following a-priori estimates hold:
[TABLE]
where denotes the continuous piece-wise affine interpolant from the values and similarly also , , etc.
Sketch of the proof. To show numerical stability of the scheme is to test the particular equations/inclusions in (4.2) successively by , , , , and (meaning a constant 1/2 on both and ) instead of used for (3.20); cf. [38, Sect.5.3]. Under the assumption (4.1a,b), one can see a telescopic cancellation effect and the resulted estimate:
[TABLE]
where , , and is from (4.2h). Here we have used the special ansatz leading to (3.18). By using the discrete Gronwall inequality for (4.15) and by the coercivity of the mentioned energies, it then yields the a-priori estimates (4.14a-d).
Moreover, special nonlinear tests by with allow for important estimates of the temperature gradients, namely (4.14e-f), cf. [12] or also e.g. [43, 38]. There are several variants [32, Sect. 8.3] how to perform this test. The simplest is to substitute for and to eliminate by this way. Then we can rely on (4.14d) to derive (4.14e-f).
Eventually, by comparison from the equations, we obtain still the “dual” estimates (4.14h).
Furthermore, the desired attribute (A4), i.e. convergence, is based on several carefully assembled arguments to be applied in a proper order. The definition of a weak solution is to take into account that (3.9a,b,c) are inclusions involving convex subdifferentals if and if and if is nonsmooth, which results to variational inequalities rather than equalities. More specifically, (3.9a) says that F(t)-M^{\prime}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}-E^{*}\partial_{(\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{e}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{e}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{e}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{e}}},{\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}\big{]}\big{]}}{[[\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}]]}{[\![\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}]\!]}{[\![\mathchoice{{\buildrel\hskip 0.70004pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\Large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}{{\buildrel\hskip 0.70004pt\text{\large.}\over{u}}}]\!]}})}\varXi\big{(}\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]},\alpha,\boldsymbol{\theta};\boldsymbol{E}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}\big{)}-E^{*}\partial_{(e,{\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]}})}\varPsi(\boldsymbol{E}u,\pi,\alpha,\boldsymbol{\theta}) belongs to the subdifferential of the convex functional \mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}\mapsto\varXi(\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]},\alpha,\boldsymbol{\theta};\boldsymbol{E}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}), which straightforwardly (after integrating the term by part in time) leads to the weak formulation of (3.9a) as
[TABLE]
for all \mathchoice{\text{\small\widetilde{\text{\normalsize}}\hskip 0.26999pt}}{\text{\small\widetilde{\text{\normalsize}}}}{\widetilde{\pi}\hskip 0.29999pt}{\widetilde{\pi}}\in L^{2}(I;H^{1}(\mathchoice{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}};{\mathbb{R}}^{d})) and \mathchoice{\text{\small\widetilde{\text{\normalsize}}\hskip 0.26999pt}}{\text{\small\widetilde{\text{\normalsize}}}}{\widetilde{\alpha}\hskip 0.29999pt}{\widetilde{\alpha}}\in L^{2}(I;H^{1}(\mathchoice{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}}{\varGamma_{\mbox{\tiny\rm C}}})), while the heat equation (3.9d)–(3.10) with the initial condition from (2.6) has been weakly formulated already as (3.11), and eventually the resting initial conditions are satisfied, i.e. , , and .
Proposition 4.3** (Convergence)**
Let again (4.5) hold and let be some approximate solution from Proposition 4.2. Then there is a subsequence with converging weakly in the topologies from (4.14a-c,e-g) to some limit . Moreover, each of these limits is a weak solution in the sense (4.16) with (3.11) to the original problem (3.9) with the initial conditions (2.6).*
Sketch of the proof. The general strategy is, after formulating the discrete scheme (4.2) in the form like (4.16) and (3.11) by using the by-part summation in time and after selecting weakly* convergent subsequences by Banach’s selection principle, to make the limit passage in the mechanical part by using also the compactness arguments through the Aubin-Lions theorem, and further to show the mechanical-energy conservation and to make limit passage in the total dissipated energy as an equality, which further shows the strong convergence of the dissipative heat. This eventually allows for a limit passage in the heat-transfer equation. We refer e.g. to a rather abstract scheme [38, Sect.5.3.2] which covers our general model under the assumptions (4.5).
Remark 4.4** (Numerical implementation.)**
Coming back to (A1), the numerical implementation of (4.2) needs still a spatial discretization. To this goal, a finite-element method is well applicable. One can consider P1 or Q1 finite elements in bulk for discretisation of , , and , and their interface variants for discretisation of , , , and ; cf. e.g. [23] for coupled bulk-surface FEM. In any case, in coupled anisothermal situations, the adiabatic-like effects caused by (4.2g) makes the non-negativity of temperature in spatially discretised problem questionable even on acute triangulations and only a successive limit passage seems to work, i.e. first the space discretisation and afterwards the time discretisation; cf. also the discussion in [6, 7].
Remark 4.5** (The role of and .)**
The dissipation and stored interfacial energies and should be distinguished in particular for anisothermal or reversible-delamination situations. Typically, involves rate-dependent effects reflecting the idea that fast processes leads to some extra microscopical atomic vibrations, i.e. some heat production, while reflects the idea that degradation of the adhesive creates microscopical new interface or microcracks in the adhesive, which deposits some energy without causing a heat production and which can be get back if healing is allowed. A specific example might be
[TABLE]
with some small and causing small heat production during fast delamination and very slow healing, respectively.
Remark 4.6** (Some variants of the model.)**
Actually, there are several variations of the presented model. E.g., considering the Signorini contact (resulting by sending \gamma_{\scriptscriptstyle\textrm{C}}({\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}})\to+\infty if {\color[rgb]{0,0,0}\mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}}<0 in (3.3)) does not comply with (4.5e) and needs either ignoring heat transfer or inertial effects, cf. e.g. [43, 44, 52]. Engineering models typically do not allow for healing, i.e. in (4.17b), which needs weakening of the continuity assumption (4.5h) but still allows for a limit passage in the conventional weak formulation of the delamination flow rule, cf. again [43, 44]. If the adiabatic effects on would be suppressed (by setting ), then the delamination might be allowed to evolve rate-independently by considering both and in (4.17b). The viscosities in the adhesive and the bulk (cf. Remark 4.7 below) then avoid the question of a suitable concept of weak solutions. Let us point out that, otherwise, a combination of two rate-independent processes governed not by a jointly convex potential leads to various quite distinct solution concepts, cf. [38] for a thorough discussion; actually our semi-implicit formula (4.2) would lead to a certain stress-driven-like weak solution, cf. [65].
Remark 4.7** (Hidden rate-dependency of plastic slip.)**
The reason for considering viscosity in the adhesive (i.e. the dumpers and in Fig. 2) is to control the slip rate in even though there is no explicit rate dependence in the flow rule (2.3g); note that \mathchoice{\big{[}\big{[}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}\big{]}\big{]}}{[[\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}]]}{[\![\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}]\!]}{[\![\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}]\!]} is well controlled in due to the viscosity in the bulk. Otherwise, would be controlled only in or, more precisely, in measures on the closure of , which would make analytical troubles unless would not be constant or, at least, unless and would depend only on , which is then in , but not on .
Remark 4.8** (Thermally expansive adhesive.)**
Like the term in (3.1b), the term in (3.1d) could be replaced by \theta_{\scriptscriptstyle\textrm{A}}b_{0}(\alpha)\cdot\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]} with some . This modification could model a thermally expansive adhesive while not substantially affecting the analysis.
5 Extension towards poro-viscoelastic media and adhesives
The set of internal variable (so far composed from the interfacial damage and plastic-like interfacial slip) can be enhanced in particular applications. One possible way is to consider damage and plastic strain also in the bulk . Another enhancement may consider also other physical processes both in the bulk and in the adhesive on the interface. E.g. delamination of piezoelectric materials has been extensively studied in literature.
Here we want to illustrate the applicability of the structure (3.9) on processes with nonlocal dissipation potentials and which are undergoing not only in the adhesive surface but possibly also in the bulk. More specifically, both humidity/moisture (or water pressure) and temperature are know to affect the debonding processes in specific applications, cf. e.g. [1, 3, 17, 33, 54] or also [55, Chap. 6]. Then, beside heat transfer, also a water transport in the porous adhesive layer and/or in the poro-viscoelastic bulk is to be considered. In fact, the diffusant need not be just water (as typically in civil engineering constructions or in rock or soil mechanics) but some solvents in polymers or hydrogen in metals undergoing the so-called metal/hydride transition, cf. e.g. [51]. In some situations, only solvent transport over the adhesive interface is relevant, while in other situations the transport inside the bulk accompanied with pronounced swelling effects may lead to delamination, cf. Remark 5.2 below. Of course, both mentioned situations can occur simultaneously, which will be covered by the model presented below.
Models used in engineering are often very phenomenological and without rational-mechanical structure. In contrast, consistently with the previous part, here we want to cast a thermodynamically consistent model where the diffusion is governed by the chemical potential (containing also a pressure) so that we cover both Fick and Darcy law. Other concepts used are the Biot model for saturated flow [11] and capillarity due to the gradient of concentration involved in the free energy, leading to the Cahn-Hilliard model [18].
We denoting by and the content of diffusant in the interfacial adhesive layer and in the bulk, respectively. Then the specific bulk free energy (3.1b) expands as:
[TABLE]
with as in (3.1b) and with the so-called Biot modulus, the Biot coefficient, a coefficient (related but not identical with the bulk modulus) and the equilibrium content.
Analogously we also cast the interface model. In contrast to (3.1d), the analog of the Biot term couples the normal and the tangential components of \mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]} because the analog of the strain trace, i.e. with denoting the identity matrix, is the sum of components of \mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]}, i.e. \mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm N}}}}{[[u]]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}{[\![u]\!]_{{{}_{\rm N}}}}+{\rm I}_{d-1}\cdot\mathchoice{\big{[}\big{[}u\big{]}\big{]}_{{{}_{\rm T}}}}{[[u]]_{{{}_{\rm T}}}}{[\![u]\!]_{{{}_{\rm T}}}}{[\![u]\!]_{{{}_{\rm T}}}}={\rm I}_{d}\cdot\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]} with the vector . Then the specific contact interface energy (3.1d) is enhanced as
[TABLE]
Like in (3.1a), the overall free energy is
[TABLE]
The driving force for diffusion is the gradient of the chemical potential which is a derivative of the free energy according concentration. Here, we have both the bulk and the adhesive chemical energies, denoted respectively by and . Because of the gradients of and of , we must consider the functional derivative, which gives
[TABLE]
The diffusion equations in the bulk and in the adhesive layer then represent conservation of the mass of the diffusant:
[TABLE]
The right-hand side of (5.5a) is the source/sink of the diffusant in the adhesive layer due to the transfer from the adjacent bulk region, and means the trace of in on and analogously for . This flux occurs naturally in the boundary condition (5.5c) on . In principle, the coefficient may depend on , , and on the traces of on , or may even be different on the boundary of and . Of course, this system is to be completed by the initial conditions, namely
[TABLE]
This diffusion problem (5.5) is coupled with the original model thermomechanical model (2.4) and its boundary condition (2.5) as well as the contact interfacial model (2.3) by considering the stress taking into account (5.1), which now gives
[TABLE]
with denoting the spherical (sometimes called volumetric) part of the strain tensor . This is now to be used in (2.3a), and thus in (2.4a) and in (2.5b), too.
Consistently with the previous notation, we denote and and formulate the diffusion problem (5.5) in a monolithic way. Again we use the gradient/difference operator from (3.10a) and, like , we define the monolithic mobility operator
[TABLE]
Then the diffusion problem (5.4)–(5.5) can be written simply as
[TABLE]
where we used already the special choice (5.2) which causes, in particular, that does not depend on and .
The energetics of this model (5.9) coupled with (3.9) can be revealed by testing the equations in (5.9) respectively by and , and the equations in (3.9a,b,c,d) again respectively by , , , and . Instead of (3.1), the overall free energy now is from (5.3), while the overall dissipation rate (3.5) is enhanced by an extra-dissipation due to diffusion processes so that altogether it reads as:
[TABLE]
Written monolithically, it means
[TABLE]
where the bilinear form is again from (3.12). The specific dissipation rates and in (3.14a,b) are then expanded by the interfacial heat source and by the bulk heat source , and also part (say 1/2) of the heat production due to the last term in (5.10), so that
[TABLE]
while the resting heat production due to the last term in (5.10) should go into the adjacent bulk domains through enhancement of the right-hand side in the boundary condition (2.5d) by .
Actually, one expects the dissipation rate to depend on rates, i.e. here on and rather than on the gradients and or the difference . Indeed, the chemo-mechanical part of the problem essentially enjoys the general structure (3.9) with the rate-dependent (pseudo)potential of dissipative forces which is now nonlocal, however. To see this quite well-known fact, we adapt a bit technical calculation [46] to our monolithic bulk/interfacial model. We realize the structure of (5.9)
[TABLE]
with from (5.8), and where denotes a quadratic potential
[TABLE]
Thus, the Gâteaux derivative of in (5.13) “monolithically” expresses the combination of the Beltrami-Laplace operator with the mobility matrix on and the Laplace operator with the mobility matrix on , and the transition conditions between and the adjacent bulk. More specifically, is the linear positive-definite operator on denoting the Hilbert space modulo constant functions, or essentially equivalently . Here and are assumed symmetric to ensure existence of a potential of this linear operator.
The convex conjugate of is also quadratic and defines by means of its Gâteaux differential the inverse operator which is nonlocal. Realizing the definition of the chemical potentials, (5.13) then rewrites as
[TABLE]
which already fits with the form of (3.9). The dissipation rate (5.10) can be then equivalently expressed in terms of rates as
[TABLE]
where denotes the mentioned convex conjugate functional.
The entropy balance (3.17) now uses the enhanced heat production (5.12) together with the entropy production on the boundaries due to the mentioned additional heat production in the boundary conditions (2.5d).
As the free energy is now convex (even quadratic) in (u,\mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]},\pi,\zeta_{\scriptscriptstyle\textrm{A}},\zeta_{\scriptscriptstyle\textrm{B}}), i.e. (4.1a,b) now reads as
[TABLE]
the system (3.9) enhanced by the diffusion equation (5.13) bears the energy-conserving discretisation [46] similarly as (4.2) as a three-step decoupled scheme. The first sub-problem now handles . More specifically, having in mind the monolithic form (5.15), thanks to (5.17a), the first step (4.2a-c) will be now made jointly with
[TABLE]
where that the difference formula (4.3) simplifies to (5.18b) because is quadratic in the variable involved in . Of course, the recursive scheme (5.18) should start for by using the initial conditions . Like (4.9), the resulted recursive boundary-value problems have potentials. More specifically, in (4.9) now is to contain also the argument and the whole functional (4.9) is to be expanded by the term , cf. also [46, Formula (17a)] for the bulk part.
The energetics and the a-priori estimates related to this diffusion enhancement can be obtained by testing (5.18a) by and (5.18b) by , which yields
[TABLE]
which is to be added to (4.10) and then to (4.12), too. In addition to the a-priori estimates (4.14), assuming and positive definite, , , and , and also and , we now obtain also
[TABLE]
This guarantees numerical stability of the staggered discretisation enhanced by diffusion, and also allows for convergence towards the weak solution now enhanced by suitable weak formulation of the initial-boundary-value problem (5.4)–(5.5)–(5.6), cf. [46] for details as far as the bulk model concerns while the -interfacial variant is analogous.
Remark 5.1** (More coupling.)**
The dissipation mechanism in the thermo-mechanical part can easily be made dependent. For example, the friction coefficient often depend on the diffusant content . Also, the coefficient can depend on and on \mathchoice{\big{[}\big{[}u\big{]}\big{]}}{[[u]]}{[\![u]\!]}{[\![u]\!]}. In some applications it is also reasonable to make elastic moduli in the bulk or in the adhesive and dependent on the diffusant content or , respectively. Yet, this latter coupling would affect the chemical potentials and and might make the analysis and computer algorithmic implementation more complicated.
Remark 5.2** (Models with .)**
Some engineering models consider moisture propagation only along in the adhesive, while the bulk is pores-free. This situation is covered by the above model simply by putting the initial bulk content and the transient adhesive-bulk coefficient .
Acknowledgment
The author is deeply thankful for many inspiring discussions with Vladislav Mantič and Roman Vodička. Also, many comments and conceptual suggestions of two anonymous referees are truly acknowledged. This research has been partly covered by the Czech Science Foundation through the grants 16-34894L “Variational structures in continuum thermomechanics of solids” and 17-04301S “Advanced mathematical methods for dissipative evolutionary systems”. Eventually, the hospitality of the University of Seville within several stays during 2015–2018 as well as the institutional support RVO: 61388998 (ČR) is acknowledged, too.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] G. Alfano, S. Marfia, and E. Sacco. A cohesive damage-friction interface model accounting for water pressure on crack propagation . Comput. Methods Appl. Mech. Engrg. , 196:192–209, 2006.
- 2[2] G. Alfano and E. Sacco. Combining interface damage and friction in a cohesive-zone model. Int. J. Numer. Method. Eng. , 68:542–582, 2006.
- 3[3] L.E. Asp. The effects of moisture and temperature on the interlaminar delamination toughness of a carbon/epoxy composite . Composites Sci. Technology , 58:967–977, 1998.
- 4[4] C.C. Baniotopoulos, J. Haslinger, and Z. Morávková. Mathematical modeling of delamination and nonmonotone friction problems by hemivariational inequalities . Applic. Math. , 50:1–25, 2005.
- 5[5] M. Barboteu, K. Bartosz, and P. Kalita. A dynamic viscoelastic contact problem with normal compliance, finite penetration and nonmonotone slip rate dependent friction . Nonlinear Anal., Real World Appl. , 22:452–472, 2015.
- 6[6] S. Bartels and T. Roubíček. Thermo-visco-elasticity with rate-independent plasticity in isotropic materials undergoing thermal expansion. Math. Modelling Numer. Anal. , 45:477–504, 2011.
- 7[7] S. Bartels and T. Roubíček. Numerical approaches to thermally coupled perfect plasticity. Numer. Meth. for Partial Diff. Equations , 29:1837–1863, 2013.
- 8[8] K. Bartosz. Hemivariational inequalities modeling dynamic contact problems with adhesion . Nonlinear Analysis T.M.A. , 71:1747–1762, 2009.
