A Rothe method for a viscoelastic contact problem involving time-fractional derivatives in locking materials

Guangwang Su; Mustapha Bouallala; Boling Chen; Van Thien Nguyen

doi:10.1515/dema-2025-0194

Article Open Access

A Rothe method for a viscoelastic contact problem involving time-fractional derivatives in locking materials

, , and

Published/Copyright: December 4, 2025

Published by

Become an author with De Gruyter Brill

Author Information Explore this Subject

From the journal Demonstratio Mathematica Volume 58 Issue 1

Abstract

This study examines a contact problem involving viscoelastic materials interacting with a rigid foundation. The constitutive relationship is derived from a time-fractional Kelvin–Voigt model. The contact is represented through the normal compliance condition, and friction is modeled using a non-local Coulomb friction law. The variational formulation is established as a hemivariational inequality that incorporates fractional-order derivatives. The existence of a solution to the variational problem is proven using the Rothe method. The second part of the study examines the discretized problem through Euler schemes and finite element methods, offering error estimates for the approximate solutions.

Keywords: fractional viscoelastic constitutive law of Caputo type; contact problem; hemivariational inequality; weak solution; finite element method; error estimate

MSC 2020: 74M15; 74S40; 47J20; 35D30; 74S05

1 Introduction

In deformation theory, locking materials (also known as incompressible locking materials) refer to materials that, under mechanical loading, reach a limit configuration beyond which no further deformation can occur. This phenomenon typically arises due to high internal stresses or specific physical conditions that prevent any additional deformation. The foundational studies on variational problems in interlocking materials were first introduced in [1], [2], [3], [4]. Further developments and extensions of these works can be found in [5]. Reference [6] introduces a mathematical model describing a shape memory elastic locking material. In [7], the authors analyze a static problem involving unilateral contact between a blocking material and a rigid foundation, incorporating a non-local Coulomb friction law. More recently [8], presents a mathematical and numerical study of an electro-elastic locking material model, considering both contact and friction effects with a conductive foundation.

Variational and hemivariational inequalities serve as a fundamental theoretical foundation for examining frictional contact problems. These approaches are especially effective in describing materials exhibiting ideal locking behavior, such as elastic and viscoelastic substances. By using these formulations, it is possible to represent the intricate interactions between a solid and a contact interface, incorporating both frictional effects and the non-linearity of deformations. The theory of hemivariant and variational-hemivariant inequalities was introduced thanks to the work of Panagiotopoulos in [9], [10], [11]. On the other hand [12], [13], [14], focuses on the study of variational and hemivariational elliptic inequalities associated with the displacement field. These references also deal with the variational and hemivariational elliptic equations associated with the displacement field. In [15], the integration of variational-hemivariational inequalities with contact problems for locked materials is discussed. This approach is employed to model Signorini friction, where friction is characterized by a non-monotonic multivalued subdifferential condition.

Fractional differential equations are widely used in the modeling of various physical and mathematical problems, particularly those involving anomalous diffusion, viscoelasticity, fractal media and memory effects, where classical integer-order differential equations may not accurately describe the behavior of the system. These equations offer a more generalized approach, enabling a deeper understanding of complex phenomena in diverse fields such as physics, engineering and finance.

Fractional calculus applications span physical systems [16], [17], [18], [19], [20], viscoelastic mechanics [21], 22], image processing [16], electrical engineering [23], and biological modeling [24].

Recent investigations by Bouallala and collaborators [25], [26], [27], [28] have examined solution existence for innovative contact problem formulations in viscoelasticity, considering both frictional and frictionless scenarios. These studies analyze thermo-viscoelastic materials interacting with heat-conducting foundations, employing the Kelvin–Voigt constitutive relation combined with a fractional time derivative extension of Fourier’s law.

In separate developments [17], 29], researchers have investigated time-fractional hemivariational inequalities and their application to contact problems involving friction. The mathematical treatment of these problems combines Rothe’s temporal discretization technique with monotone operator theory.

In [30], the authors study a viscoelastic contact problem with a rigid foundation, using a fractional Kelvin–Voigt model. The contact condition is a variant of Signorini’s, with non-local Coulomb friction. They establish existence results and propose finite element schemes with error estimates.

The main focus of this study is the analysis of a new contact problem with friction for viscoelastic materials in contact with rigid foundations. This is accomplished by modeling their behavior through the Kelvin–Voigt constitutive law, which includes time fractional derivatives. By analyzing this combination, it becomes clear that the model can be re-expressed as an inclusion problem involving a fractional integral operator and a nonlinear differential equation.

In this model, the contact interaction is characterized by a penalized normal compliance condition, where the penalty coefficient is represented by ϵ. Frictional effects are incorporated using Coulomb’s law. By analyzing the problem under the given boundary conditions, we derive a nonlinear variational formulation in fractional time. The existence of a weak solution is established through monotone operator theory and Rothe’s method.

The present study investigates an alternative computational approach combining:

Spatial discretization via finite element methods.
Complete temporal discretization scheme.

The Caputo fractional derivative formulation serves as the basis for approximation. Furthermore, we establish rigorous error estimates to validate the numerical scheme’s precision and computational reliability.

This paper is organized as follows: Section 2 presents our viscoelastic frictional contact model with locking effects, featuring a nonlinear constitutive law incorporating time-fractional derivatives. We establish the mathematical framework, including notation, data assumptions, and variational formulation. Section 3 proves the existence of weak solutions through rigorous mathematical analysis. The numerical implementation is developed in Section 4, where we:

Design a fully discrete approximation scheme.
Derive optimal error estimates.
Establish convergence properties.

The Appendix provides essential background on fractional calculus and discrete analysis techniques.

2 Mathematical formulation of the fractional-order problem

We study an elastic solid initially occupying a bounded open domain Ω ⊂ R 2 with Lipschitz boundary Γ = ∂Ω. The boundary is decomposed into three pairwise disjoint measurable parts: the Dirichlet boundary Γ_D, Neumann boundary Γ_N, and contact boundary Γ_C, where meas(Γ_D) > 0.

Furthermore, we analyze the system over a finite time horizon [0, T], where T > 0. The displacement field vanishes on Γ_D × (0, T), corresponding to clamped boundary conditions. The body is subjected to:

A distributed body force f ₀ acting throughout Ω × (0, T).
Surface tractions f ₂ applied on Γ_N × (0, T).

The boundary portion Γ_C represents a potential contact surface where the deformable body may interact with a rigid foundation.

Following the notations introduced in [28], we consider the definitions of S 2 and R 2 , along with their respective inner products ⟨⋅, ⋅⟩ and norms ‖⋅‖.

Furthermore, we employ the following notation conventions:

ɛ denotes the linearized strain tensor,
Div represents the divergence operator for tensor fields,

(2.1) ε ( w ) = ε i j ( w ) , ε i j ( w ) = 1 2 w i , j + w j , i , D i v σ = σ i j , j .

where w _i,j: = ∂w _i/∂x _j.

Furthermore, let ν = (ν _i) represent the outward unit normal. We introduce the following functions:

(2.2) w ν = w ⋅ ν , w τ = w − w ν ν , σ ν = σ ν ⋅ ν , σ τ = σ ν − σ ν ν .

The complete mathematical description of the contact mechanics problem can be stated as:

Problem (P): Find a displacement field w : Ω × ( 0 , T ) → R d and a stress tensor σ : Ω × ( 0 , T ) → S d , a.e. t ∈ [0, T] and α ∈ (0, 1] such that

(2.3) σ ( t ) ∈ A ε ( w ( t ) ) + B ε D t α 0 C w ( t ) + ∂ I A ( ε ( w ( t ) ) ) in Ω × ( 0 , T ) ,

(2.4) D i v σ ( w ( t ) ) + f 0 ( t ) = 0 in Ω × ( 0 , T ) ,

(2.5) w ( t ) = 0 on Γ D × ( 0 , T ) ,

(2.6) σ ( t ) ⋅ ν = f 2 ( t ) on Γ N × ( 0 , T ) ,

(2.7) σ ν w ν ( t ) = − 1 ϵ w ν ( t ) + , ϵ > 0 on Γ C × ( 0 , T ) ,

(2.8) σ τ ( t ) ≤ π w τ ( t ) B σ ν ( w ( t ) ) , σ τ ( t ) < π w τ ( t ) B σ ν ( w ( t ) ) ⇒ w τ ( t ) = 0 , σ τ ( t ) = π w τ ( t ) B σ ν ( w ( t ) ) ⇒ ∃ κ ≠ 0 , such that σ τ ( t ) = − κ w τ ( t ) on Γ C × ( 0 , T ) .

(2.9) w ( 0 ) = w 0 in Ω .

Where I _A is defined by

(2.10) I A ( r ) = 0 , if r ∈ A , + ∞ , otherwise ,

of closed and convex set A given by

(2.11) A = { β ∈ S d , ‖ β ‖ ≤ M } .

We recall that equation (2.3) describes the constitutive law governing viscoelastic materials with locking constraints, incorporating a time-fractional derivative in the sense of Caputo. Here, A = a ijkl and B = b ijkl correspond to the elasticity and viscosity tensors, respectively. Equation (2.4) expresses the equilibrium condition for the stress-displacement field. The boundary conditions for displacement and traction are given in equations (2.5) and (2.6), respectively. The normal compliance contact condition on Γ_C is described by equation (2.7), while the initial condition is formulated in equation (2.8). Lastly, Coulomb’s friction law is presented in equation (2.9), where π denotes the friction coefficient and B represents a regularization operator. For further details, refer to [13], 14].

We work within the following Hilbert space framework:

W = L 2 ( Ω ) 2 = w = ( w i ) ∣ w i ∈ L 2 ( Ω ) , W 1 = H 1 ( Ω ) 2 , W = σ = ( σ i j ) ∣ σ i j = σ j i ∈ L 2 ( Ω ) .

These spaces are equipped with the respective inner products:

( w , φ ) W = ∫ Ω w i φ i d x , ( σ , τ ) W = ∫ Ω σ i j τ i j d x , ( w , φ ) W 1 = ( w , φ ) W + ε ( w ) , ε ( φ ) W .

Accounting for the boundary condition (2.5), we define the subspace:

X = w ∈ W 1 ∣ w = 0 on Γ D ,

and the corresponding closed convex set:

X a = w ∈ X ∣ | ε ( w ) | ≤ M a.e. in Ω .

The set X _a represents the set of displacement fields w such that the strain ɛ(w) does not exceed a critical value M in norm, at each point (or almost everywhere) in the material. This models an elastic limit or a non-damage constraint on the material.

In the case of locking materials (such as rigid-plastic materials or some polymers), this condition expresses a rigid behavior below a strain threshold (no significant strain as long as ‖ɛ(w)‖ ≤ M), followed by a plastic or damage response beyond this threshold. The value M acts as a locking threshold.

We equip the space V with the energy inner product defined by

( w , φ ) X : = ε ( w ) , ε ( φ ) W ,

which generates the energy norm  ‖ w ‖ X : = ( w , w ) X = ‖ ε ( w ) ‖ W .

Remark 2.1.

The norm ‖⋅‖_X is equivalent to the standard Sobolev norm  ‖ ⋅ ‖ W 1 on the subspace X, making (X, ‖⋅‖_X) a complete Hilbert space.

The condition meas(Γ_D) > 0 ensures the validity of Korn’s inequality [31]:

(2.12) ‖ ε ( w ) ‖ W ≥ c 0 ‖ w ‖ W 1 ∀ w ∈ X ,

where c ₀ > 0 is a constant depending Ω and Γ_D.

Applying the Sobolev trace theorem yields a constant c ₁ = c ₁(Ω, Γ_C, Γ_D) > 0 such that:

(2.13) ‖ w ‖ L 2 ( Γ C ) 2 ≤ c 1 ‖ w ‖ X ∀ w ∈ X .

To establish the mathematical foundation for analyzing Problem (P), we introduce the following working hypotheses:

(HP1) The elasticity tensor A : Ω × S d → S d satisfies the following standard properties:
1. Symmetry and regularity:
  a ijkl = a jikl = a klij ∈ L ∞ ( Ω ) .
2. Uniform ellipticity:
  ∃ m a > 0 such that a ijkl ( x ) ξ i j ξ k l ≥ m a ‖ ξ ‖ 2 ∀ ξ ∈ S d , a.e. x ∈ Ω .
3. Boundedness:
  ∃ M a > 0 such that | ( A ε ( w ) , ε ( φ ) ) W | ≤ M a ‖ w ‖ X ‖ φ ‖ X ∀ w , φ ∈ X .
(HP2) The viscosity operator B : Ω × S d → S d satisfies the following properties:
1. Symmetry and regularity:
  b ijkl = b jikl = b klij ∈ L ∞ ( Ω ) .
2. Uniform ellipticity (coercivity):
  ∃ m b > 0 such that b ijkl ( x ) ξ i j ξ k l ≥ m b ‖ ξ ‖ 2 ∀ ξ ∈ S d , a.e. x ∈ Ω .
3. Boundedness:
  ∃ M b > 0 such that | ( B ε ( w ) , ε ( φ ) ) W | ≤ M b ‖ w ‖ X ‖ φ ‖ X ∀ w , φ ∈ X .
(HP3) The forces, traction, and initial conditions adhere to the following regularities:
f 0 ∈ C 0 , T ; L 2 ( Ω ) 2 , f 2 ∈ C 0 , T ; L 2 ( Γ N ) 2 , w 0 ∈ X a .
(HP4) The friction coefficient π : Γ C × R + → R + satisfies the following conditions:
1. Lipschitz continuity: There exists L _π > 0 such that for all x , y ∈ R + ,
  | π ( ⋅ , x ) − π ( ⋅ , y ) | ≤ L π | x − y | a.e. on Γ C .
2. Measurability: For each y ∈ R + , the mapping
  x ↦ π ( x , y ) is measurable on Γ C .
3. Uniform boundedness: For each y ∈ R + , the mapping
  x ↦ π ( x , y ) is π * − bounded a.e. on Γ C ,
  where the bound is given by
  π * = sup y ∈ R + ‖ π ( ⋅ , y ) ‖ L ∞ ( Γ C ) .
(HP5) The mapping B : W Γ C * → L ∞ ( Γ C ) is linear, compact, and continuous, and it satisfies
c B = ‖ B ‖ L ∞ ( Γ C ) .

By applying Riesz’s representation theorem, we can uniquely identify elements λ ∈ V through the following characterization:

(2.14) λ ( t ) , φ X * × X = ∫ Ω f 0 ( t ) ⋅ v d x + ∫ Γ N f 2 ( t ) ⋅ v d Γ .

We introduce the bilinear forms a : X × X → R and b : X × X → R , defined respectively by:

(2.15) a w ( t ) , φ = : A ε ( w ( t ) ) , ε ( φ ) W , b w ( t ) , φ = : B ε ( w ( t ) ) , ε ( φ ) W .

We introduce the functional j c : X × X → R and j f : X × X → R as follows:

(2.16) j c w ( t ) , φ = 1 ϵ ∫ Γ C w ν ( t ) + φ ν d Γ = 1 ϵ w ν ( t ) + , φ ν Γ C , ∀ w , φ ∈ X ,

(2.17) j f w ( t ) , φ = ∫ Γ C π w τ ( t ) B σ ν ( w ( t ) ) ‖ φ τ ‖ d Γ , ∀ w , φ ∈ X .

Now, from (2.3), we have

(2.18) σ ( t ) = A ε ( w ( t ) ) + B ε D t α 0 C w ( t ) + Y ε ( w ( t ) ) ,

where

(2.19) Y ( ε ( w ( t ) ) ) ∈ ∂ I X a ( ε ( w ( t ) ) ) in Ω × ( 0 , T ) .

Due to the convexity of I X a , we obtain

(2.20) Y ( ε ( w ( t ) ) , ε ( φ − w ( t ) ) ) ≤ I X a ( ε ( φ ) ) − I X a ( ε ( w ( t ) ) ) = 0 , ∀ φ ∈ X a .

We now consider sufficiently regular functions w and σ that satisfy the following variational formulation of Green’s formula:

(2.21) σ ( t ) , ε ( φ ) W + D i v ( σ ( t ) ) , φ W = ∫ Γ σ ( t ) ν ⋅ φ d Γ , for all φ ∈ W 1 .

Combining (2.4), (2.15), (2.18) and (2.19), we have

(2.22) a w ( t ) , φ − w ( t ) + b D t α 0 C w ( t ) , φ − w ( t ) ≥ f 0 ( t ) , φ − w ( t ) X * × X + ∫ Γ σ ( t ) ν ⋅ ( φ − w ( t ) ) d Γ .

By incorporating the contact conditions (2.7) and (2.8) with the constitutive relations (2.16)– (2.18), and employing the Einstein summation convention, we derive the following variational formulation:

Problem (WF): Find a displacement field w : Ω × ] 0 , T [ → R 2 a.e. t ∈ ]0, T[, for all φ ∈ X _a and α ∈ (0, 1] such that:

(2.23) a w ( t ) , φ − w ( t ) + b D t α 0 C w ( t ) , φ − w ( t ) + j c w ( t ) , φ − w ( t ) + j f w ( t ) , φ − j f w ( t ) , w ( t ) ≥ λ ( t ) , φ − w ( t ) X * × X ,

(2.24) w ( 0 , x ) = w 0 .

3 Existence result

The existence of solutions to Problem (WF) is established through the following result:

Theorem 3.1.

Assume that (HP1)–(HP5) hold. There exists at least one solution of Problem (WF) satisfies the following regularity

w ∈ W 1,2 0 , T ; V .

The proof relies on the utilization of the nonlinear operator, specifically the Caputo fractional derivative, in conjunction with the Rothe method. It unfolds through several distinct steps.

Step 1: Let J be the functional defined by J : X × X → R as follows:

(3.1) J ( w ( t ) , φ ) : = j c ( w ( t ) , φ ) + j f ( w ( t ) , φ ) , for all φ ∈ X .

We now present the corresponding result:

Lemma 3.2.

The functional J : X × X → R satisfies the following properties:

For each fixed φ ∈ X, the mapping J(⋅, φ) is Lipschitz continuous on X.
The Clarke generalized gradient ∂J exhibits sublinear growth:

(3.2) ‖ ∂ J ( w , φ ) ‖ X * ≤ c J ( 1 + ‖ φ ‖ X ) ∀ w , φ ∈ X ,

where c _J > 0 is a constant.

Here, ∂J(w, φ) means Clarke’s subdifferential of J with respect to w, while φ is considered as a fixed parameter.

Proof.

1) Let w ₁, w ₂ be the elements in X for all t ∈ [0, T] and φ ∈ X.

By (3.2), we have that

(3.3) J ( w 1 ( t ) , φ ) − J ( w 2 ( t ) , φ ) = j c ( w 1 ( t ) , φ ) + j f ( w 1 ( t ) , φ ) − j c ( w 2 ( t ) , φ ) − j f ( w 2 ( t ) , φ ) .

Using (2.16), and this inequality [ w 1 ] + − [ w 2 ] + ≤ w 1 − w 2 , we conclude that

(3.4) j c ( w 1 ( t ) , φ ) − j c ( w 2 ( t ) , φ ) ≤ c 1 2 ϵ w 1 ( t ) − w 2 ( t ) X ‖ φ ‖ X .

Another word by (HP4)–(HP5), we can see that

(3.5) j f ( w 1 ( t ) , φ ) − j f ( w 2 ( t ) , φ ) ≤ c 1 2 c B L π + π * c B * w 1 ( t ) − w 2 ( t ) X ‖ φ ‖ X .

Such that c B * is derived from the continuity of the operator B.

This means that for each φ ∈ X fixed, J(⋅, φ) is Lipschitz continuous on X.
From [32], Corollary 4.15(iv) and Theorem 7.3], we find that

(3.6) ∂ J ( w ) ⊂ ∂ j f ( w ) + ∂ j c ( w ) .

It follows from [32], Corollary 4.15(v)] that

(3.7) ∂ J ( w , φ ) X * ≤ m J 1 + ‖ φ ‖ X ,

with m _J = m _c + m _f, where m c = 1 ϵ and m f = c B L π + π * c B * .□

Step 2: Now give, y ( t ) = D t α 0 C w ( t ) and w ( t ) = I t α 0 y ( t ) + w 0 . By using the subdifferential ∂J, the variational inequality (2.23) can be written as finding w(t) such that:

λ ( t ) − b ( y ( t ) ) − a ( w ( t ) ) ∈ ∂ J ( w ( t ) ) .

By substituting w ( t ) = I t α 0 y ( t ) + w 0 , we obtain:

λ ( t ) − b ( y ( t ) ) − a I t α 0 y ( t ) + w 0 ∈ ∂ J I t α 0 y ( t ) + w 0 ,

which gives the following inequality

(3.8) b ( y ( t ) ) + a I t α 0 y ( t ) + w 0 + ∂ J I t α 0 y ( t ) + w 0 ∋ λ ( t ) .

Consider a fixed integer N ∈ N and let δ = T/N be the temporal step size. We investigate the following discrete approximation of the fractional integral operator I t n α 0 y ( t ) :

(3.9) I ̃ t n α 0 y ( t ) = 1 Γ ( α ) ∑ i = 1 n ∫ t j − 1 t j ( t n − s ) α − 1 y ( t i ) d s = δ α Γ ( 1 + α ) ∑ i = 1 n y ( t i ) ( n − i + 1 ) α − ( n − i ) ,

where t _k = kδ. Also, we define λ k δ as follows:

(3.10) λ i δ = δ − 1 ∫ t i − 1 t i λ ( s ) d s , for i = 1 , … , N .

By employing Rothe’s discretization method to equation (3.8), we obtain the following time-discrete fractional problem:

Problem (FR): Find y k δ ⊂ X a for k = 1, …, N such that

(3.11) b y k δ + a w k δ + ∂ J w k δ ∋ λ k δ ,

where

(3.12) w k δ = w 0 + δ Γ ( 1 + α ) ∑ i = 1 k y i δ ( k − i + 1 ) α − ( k − i ) α .

We now establish the following existence result for the discretized problem:

Lemma 3.3

(Existence for Discrete Problem). There exists a critical step size δ ̄ > 0 such that for all δ ∈ ( 0 , δ ̄ ) , the discrete problem (3.11) admits at least one solution.

Proof.

We assume that y k δ k = 0 n − 1 are provided, and we will select y n δ ∈ X satisfying (3.11) and (3.12). To accomplish this, we introduce the following multivalued operators: M: X → X and M ₀: X → X* such that

(3.13) M ( y ) : = M 0 ( y ) + ∂ J ( y ) ,

and

(3.14) M 0 ( y ) : = b ( y ) + a δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ ( n − i + 1 ) α − ( n − i ) α + δ α y Γ ( 1 + α ) + w 0 .

Our objective is now to prove that the operator M is surjective.

To begin, we show that M is coercive.

By applying Lemma 3.2 together with equation (3.14), we derive the following:

(3.15) ∂ J δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ ( n − i + 1 ) α − ( n − i ) α + δ α Γ ( 1 + α ) + w 0 X * ≤ m J δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ X ( n − i + 1 ) α − ( n − i ) α + δ α Γ ( 1 + α ) ‖ y ‖ X + 1 + ‖ w 0 ‖ X ≤ m J 1 + m δ + δ α Γ ( 1 + α ) ‖ y ‖ X ,

where

(3.16) m δ = δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ ( n − i + 1 ) α − ( n − i ) α .

Based on hypothesis (HP2) and equation (2.13)

(3.17) M ( y , y ) ≥ m b ‖ y ‖ X 2 − M a m J + δ α Γ ( 1 + α ) ‖ y ‖ X ‖ y ‖ X − m J 1 + m δ + δ α Γ ( 1 + α ) ‖ y ‖ X ≥ m b − δ α ( M a + m J ) Γ ( 1 + α ) ‖ y ‖ X 2 − M a c δ + m J δ α Γ ( 1 + α ) ‖ y ‖ X − m J 1 + m J .

We define δ ̄ = Γ ( 1 + α ) m b δ α ( M a + m J ) to establish the coervcivity of the operator M.

Taking into account the conditions imposed on a and b, we deduce that

(3.18) M 0 ( y , y ) ≥ m b − δ α M a Γ ( 1 + α ) ‖ y ‖ X 2 .

As a result, M ₀ can be identified as a pseudo-monotone operator.

We now proceed to show that the operator M ₁: X → X*, defined by

(3.19) M 1 ( y ) : = ∂ J δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ ( n − i + 1 ) α − ( n − i ) α + δ α Γ ( 1 + α ) y + w 0 ,

is pseudo-monotone.

Thanks to the properties of J and the reflexivity of the space X, it follows that M(y) is nonempty, convex, and weakly compact for every y ∈ X.

Moreover, according to Lemma 3.2, the operator M is bounded.

Let us now consider a sequence {y _m} ⊂ X such that y _m ⇀ y weakly in X as m → ∞, and

(3.20) ξ m ∈ ∂ J δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ ( n − i + 1 ) α − ( n − i ) α + δ α Γ ( 1 + α ) y m + w 0 .

Since the operator ∂J is bounded, the sequence {ξ _m} remains bounded in X*.

Possibly after extracting a subsequence, we assume that y _m ⇀ y weakly in X* as m → ∞.

Moreover, taking into account that the graph of the multivalued operator

(3.21) y ↦ ∂ J δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ ( n − i + 1 ) α − ( n − i ) α + δ α y Γ ( 1 + α ) + w 0 ,

is closed with respect to the X × X* topology (see [32], Proposition 3.23(v)]), we conclude that

(3.22) ξ ∈ ∂ J δ α Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ ( n − i + 1 ) α − ( n − i ) α + δ α y Γ ( 1 + α ) + w 0 .

Subsequently, it becomes evident ξ ∈ M(y), and we define

(3.23) ξ m , y m → ξ , y X * × X , as m → ∞ .

By applying Definition 4.5, we conclude that the operator M possesses the property of pseudo-monotonicity.

As a result, M is pseudo-monotone, which ensures the existence of at least one solution to Problem (FR).□

Step 3: At this stage, we introduce the sequence of solutions associated with the fractional Rothe problem (3.11).

The following result can then be established:

Lemma 3.4

(Uniform Bounds). Under assumptions (HP1)–(HP5) and equation (3.16), there exist positive constants δ ₀ > 0 and c > 0, independent of the discretization parameter δ, such that for all 0 < δ < δ ₀, the solution triple y k δ , w k δ , ξ k δ of Problem (FR) satisfies:

Uniform boundedness:

(3.24) max 1 ≤ k ≤ N ‖ y k δ ‖ X + ‖ w k δ ‖ X + ‖ ξ k δ ‖ X * ≤ c .

Discrete variational inclusion:

(3.25) b y k δ + a w k δ + ξ k δ = λ k δ , ξ k δ ∈ ∂ J y k δ ,

for all time steps k = 1, …, N.

Proof.

Let 1 ≤ n ≤ N, multiplying equation (3.25) by y n δ , we have

(3.26) b y n δ , y n δ + a w n δ , y n δ + ξ n δ , y n δ X * × X = λ n δ , y n δ X * × X .

Using to (3.12), (HP1)–(HP2) and Lemma 3.2, we obtain

(3.27) λ n δ , y n δ ≥ m b ‖ y n δ ‖ X 2 − M a ‖ w n δ ‖ X ‖ y n δ ‖ X − m J 1 + ‖ w n δ ‖ X ‖ y n δ ‖ X ≥ m b ‖ y n δ ‖ X − M a ‖ w 0 ‖ X + δ α Γ ( 1 + α ) ∑ i = 1 n ‖ y i δ ‖ X ( n − i + 1 ) α − ( n − i ) α ‖ y n δ ‖ X − 1 + ‖ w 0 ‖ X + δ α Γ ( 1 + α ) ∑ i = 1 n ‖ y i δ ‖ X ( n − i + 1 ) α − ( n − i ) α ‖ y n δ ‖ X ≥ m b ‖ y n δ ‖ X 2 − M a δ α Γ ( 1 + α ) ‖ y n δ ‖ X 2 − M a ‖ w 0 ‖ X ‖ y n δ ‖ X − m J ‖ y n δ ‖ X − m J ‖ w 0 ‖ X ‖ y n δ ‖ X − m J δ α Γ ( 1 + α ) ‖ y n δ ‖ X 2 − M a δ α Γ ( 1 + α ) ∑ i = 1 n − 1 ( n − i + 1 ) α − ( n − i ) α ‖ y i δ ‖ X ‖ y n δ ‖ X − m J δ α Γ ( 1 + α ) ∑ i = 1 n − 1 ( n − i + 1 ) α − ( n − i ) α ‖ y i δ ‖ X ‖ y n δ ‖ X .

Therefore, based on the preceding analysis, we deduce that

(3.28) λ n δ X + δ α ( M a + m J ) Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ X ( n − i + 1 ) α − ( n − i ) α + m J + m J + M a ‖ w 0 ‖ X ≥ m b − δ α ( M a + m J ) Γ ( 1 + α ) y n δ X .

Selecting δ ̄ = m b Γ ( 1 + α ) 2 ( M a + m J ) 1 / α , we deduce that m b − δ α ( M a + m J ) Γ ( 1 + α ) ≥ m b 2 for all 0 < δ < δ ̄ . Thus

(3.29) 2 λ n δ X m b + m J + m J ‖ w 0 ‖ X + M a ‖ w 0 ‖ X m b + 2 δ α ( M a + m J ) m a Γ ( 1 + α ) ∑ i = 1 n − 1 y i δ X ( n − i + 1 ) α − ( n − i ) α ≥ y n δ X .

Utilizing hypothesis (HP3), for every δ > 0 and n ∈ N . There exists a positive constant c _λ > 0 such that λ n δ X ≤ c λ .

Let’s denote m ̃ = 2 m b c λ + M a + ‖ w 0 ‖ X + m J ‖ w ‖ X + m J .

By employing the generalized discrete Gronwall inequality, we find that there exists a positive constant m ₀ such that

(3.30) λ n δ X ≤ m 0 ⁡ exp 2 δ α ( M a + m J ) m b Γ ( 1 + α ) ∑ i = 1 n − 1 ( n − i + 1 ) α − ( n − i ) α ≤ m 0 ⁡ exp 2 T α ( M a + m J ) m b Γ ( 1 + α ) ≤ c .

Combing (3.12)– (3.30), we derive the following result

(3.31) w n δ X = w 0 + δ α Γ ( 1 + α ) ∑ i = 1 n y i δ ( n − i + 1 ) α − ( n − i ) α X + ‖ w 0 ‖ X + ∑ i = 1 n t n − i + 1 α − t n − i α ≤ ‖ w 0 ‖ X + c T α Γ ( 1 + α ) ≤ c .

Finally, using Lemma 3.2, we derive the following estimate for λ n δ by

(3.32) λ n δ X ≤ m J 1 + ‖ w n δ ‖ X ≤ m J ( 1 + c ) .

Thus, Lemma 3.4 is established.□

Step 4: In this phase, we aim to demonstrate the existence of a solution to Problem (PV).

Proof of Theorem

(3.1). We consider the sequence δ n such that δ _n → 0, as n → ∞.

Based on the estimate (3.24), the sequence y ̄ δ , w ̄ δ and λ ̄ δ , which interpolate to y δ , w δ and ξ δ respectively, are bounded for k = 1, …, N.

Therefore there exist y ∈ X, w ∈ X and ξ ∈ X* such that

(3.33) y ̄ δ → y weakly in X , as δ → 0 ,

(3.34) w ̄ δ → w weakly in X , as δ → 0 ,

(3.35) ξ ̄ δ → x weakly in X * , as δ → 0 .

Using [33], Lemma 4(a)], we derive that

(3.36) I t α 0 y ̄ δ → I t α 0 y δ weakly in X , as δ → 0 ,

and by (3.33), it follows that

(3.37) w ̄ δ ( t ) − w 0 − I t α 0 y ̄ δ ( t ) X = δ α Γ ( 1 + α ) ∑ i = 1 n y i δ ( n − i + 1 ) α − ( n − i ) α − 1 Γ ( α ) ∫ 0 t ( t − s ) α − 1 y ̄ δ ( s ) d s X ≤ c Γ ( α ) ∫ t n − 1 t n ( t n − s ) α − 1 d s + ∫ 0 t ( t − s ) α − 1 − ( t n − s ) α − 1 d s ≤ c Γ ( α ) ( t n − t ) α + t α + ( t n − t ) α − t n α .

for all t ∈ t n − 1 , t n . Then

(3.38) w ̄ δ − w 0 − I t α 0 w ̄ δ ( t ) → 0 strongly in X , as δ → 0 ,

Combining this result with equation (3.37), we obtain the weak convergence:

(3.39) w ̄ δ ( t ) ⇀ w 0 + 0 I t α y ( t ) weakly in X as δ → 0 .

Since the subdifferential mapping w ↦ ∂J(w) is upper semicontinuous from X to X* (endowed with weak topologies), and using (3.35) together with [32], Theorem 3.13], we deduce:

(3.40) ξ ( t ) ∈ ∂ J w 0 + 0 I t α y ( t ) for a.e. t ∈ ( 0 , T ) .

We now introduce the corresponding Nemytskii operators:

(3.41) b ̄ y ( t ) = b y ( t ) , and a ̄ y ( t ) = a w 0 + I t α 0 y ( t ) .

Given assumptions (HP1)–(HP2), as well as equations (3.32) and (3.34), we have

(3.42) b ̄ y ̄ δ → b ̄ y weakly in X , as δ → 0 ,

(3.43) a w 0 + I t α 0 y ̄ δ ( t ) → a w 0 + I t α 0 y ( t ) weakly in X , as δ → 0 ,

Derive from (3.24) and (HP1), it follows that

(3.44) ∫ 0 T a I t α 0 y ̄ δ ( t ) + w 0 X ≤ M a c Γ ( 1 + α ) ∫ 0 T t α d t + T M a ‖ w 0 ‖ X ≤ M a c T 1 + α Γ ( 2 + α ) + T M a ‖ w 0 ‖ X .

Applying the Lebesgue Dominated Convergence Theorem yields the following representation:

(3.45) lim δ → 0 a y ̄ δ , φ = lim δ → 0 ∫ 0 T a I t α 0 + y ̄ δ ( t ) + w 0 , φ ( t ) d t = ∫ 0 T lim δ → 0 a I t α 0 y ̄ δ ( t ) + w 0 , φ ( t ) d t = ∫ 0 T lim δ → 0 a I t α 0 y δ ( t ) + w 0 , φ ( t ) d t = a ( y , φ ) .

Conversely, invoking [34], Lemma 3.3], we obtain the strong convergence:

(3.46) λ δ → λ strongly in X as δ → 0 .

By combining the convergence results (3.42)– (3.44) and passing to the limit in inequality (3.11), we conclude that y ∈ L ²(0, T; X) solves problem (3.8). Consequently, the function w ∈ W ^1,2(0, T; X) defined by

(3.47) w ( t ) = w 0 + 0 I t α y ( t ) for a.e. t ∈ ( 0 , T ) ,

provides a solution to Problem (WF).□

4 Complete spatio-temporal discretization

In this section, we develop a fully discrete approximation method for Problem (P) and establish its optimal convergence rate. The scheme is constructed using the finite element space V ^h and a partition of the time domain [0, T] into discrete points 0 = t ₀ < t ₁ < … < t _N = T.

For each n = 1, …, N, the local time step is defined as k _n = t _n − t _n−1. The method accommodates variable time steps, with the global step size given by k = max 1 ≤ n ≤ N k n .

Given a continuous function w(t) belonging to an appropriate function space, we adopt the notation w _i = w(t _i) for each discrete time point i = 0, 1, …, N. The forward difference operator is defined as Δw _n = w _n − w _n−1, with the associated difference quotient given by δw _n = Δw _n/k _n.

Consider T h , a shape-regular family of triangular finite element subdivisions of the closure Ω ̄ , which respects the boundary partition Γ = Γ ̄ C ∪ Γ ̄ D ∪ Γ ̄ N . The corresponding finite element spaces for approximating the displacement field w are constructed as:

X h = φ h ∈ [ C ( Ω ̄ ) ] d φ h | Tr ∈ [ P 1 ( Tr ) ] d ∀ Tr ∈ T h , φ h = 0 on Γ ̄ D

with X a h = X a ∩ X h .

The fully discrete numerical scheme employs an implicit Euler time discretization and takes the following form:

Problem (PVHK): Find a displacement field w n h k ⊂ X a h for all φ ^h ∈ X ^h and n = 1, …, N

(4.1) a w n h k , φ h − w n h k + b D t n α 0 C w n h k , φ h − w n h k + j c w n h k , φ h − w n h k + j f w n h k , φ h − j f w n h k , w n h k ≥ λ n , φ h − w n h k X ,

(4.2) w 0 h k = w 0 h ,

where w 0 h ∈ X h is an approximate of w ₀.

To simplify the notation once again, we introduce the velocity

(4.3) u n h k = δ w n h k , n = 1 , … , N , and w n h k = ∑ j = 1 n k j u j h k + w 0 h .

By employing a discrete counterpart of Theorem 3.1, we establish the existence of at least one solution w n h k ∈ X h to the discrete problem (4.1) and (4.2).

We now turn to the analysis of convergence properties for the fully discrete approximation scheme. Specifically, we obtain the following convergence result:

Theorem 4.1.

Under the hypotheses (HP1)–(HP5) and for an initial approximation w 0 h ∈ X h satisfying

(4.4) ‖ w 0 − w 0 h ‖ X → 0 as h → 0 ,

the fully discrete solutions of system (4.1) and (4.2) exhibit the following convergence behavior:

(4.5) max 1 ≤ n ≤ N ‖ w n − w n h k ‖ X → 0 as h , k → 0 .

Proof.

Take φ = w n h k in (2.21) at t = t _n, we find

(4.6) a w n , w n − w n h k ≤ b D t n α 0 C w n , w n h k − w n + j c w n , w n h k − w n + j f w n , w n h k − j f w n , w n − λ n , w n h k − w n X .

Substituting φ ^h by φ n h in (4.1), we obtain

(4.7) a w n h k , w n h k − φ n h ≤ b D t n α 0 C w n h k , φ n h − w n h k + j c w n h k , φ n h − w n h k + j f w n h k , φ n h − j f w n h k , w n h k − λ n , φ n h − w n h k X .

In author word, we have

(4.8) a w n − w n h k , w n − w n h k − a w n h k , φ n h − w n = a w n , w n − w n h k + a w n h k , w n h k − φ n h .

Adding (4.6), (4.7) combined with (4.8), it follows that

(4.9) a w n − w n h k , w n − w n h k ≤ R j f + R 1 + R j c + R 2 ,

where

(4.10) R j f = j f w n , w n h k − j f w n , φ n h + j f w n h k , φ n h − j f w n h k , w n h k ,

(4.11) R 1 = a w n h k − w n , φ n h − w n + b D t n α 0 C w n h k − D t n α 0 C w n , φ n h − w n h k ,

(4.12) R j c = j c w n , w n h k − φ n h − j c w n h k , w n h k − φ n h ,

and

(4.13) R 2 = a w n , φ n h − w n + b D t n α 0 C w n , φ n h − w n + j c w n , φ n h − w n + j f w n , φ n h − j f w n , w n − λ n , φ n h − w n X .

Using the following inequalities

(4.14) x y ≤ β x 2 + 1 4 β y 2 , ∀ β > 0 ,

(4.15) [ x ] + − [ y ] + ≤ | x − y | ,

(4.16) φ n h − w n h k X ≤ φ n h − w n X + w n − w n h k X ,

combined with the hypotheses (HP1)–(HP2), (2.15) and (2.11), there exists a positive constant c such that

(4.17) R 1 ≤ c w n − w n h k X + φ n h − w n X + D t n α 0 C w n − D t n α 0 C w n h k X ,

and

(4.18) R j c ≤ c w n − w n h k X + φ n h − w n X .

Thanks to (2.14) and the hypothesis (HP4)–(HP5), (4.16) and (2.11), we found the following estimate

(4.19) R j f = j f w n , w n h k − j f w n , φ n h + j f w n h k , φ n h − j f w n h k , w n h k ≤ c B L π c 1 2 w n − w n h k X × φ n h − w n X .

Now, using Definition 4.2, we can deduce the discrete approximation for the fractional derivative in the following manner:

(4.20) D t n α 0 C w n h k = 1 Γ ( 1 − α ) ∫ 0 t n δ w n h k t n − s − α d s = 1 Γ ( 1 − α ) ∑ j = 1 n ∫ ( j − 1 ) k j k 1 k w j h k − w j − 1 h k + O ( k ) n k − s − α d s = 1 Γ ( 1 − α ) ( 1 − α ) ∑ j = 1 n 1 k w j h k − w j − 1 h k + O ( k ) ( n − j + 1 ) 1 − α − ( n − j ) 1 − α k 1 − α = 1 k α ( 1 − α ) Γ ( 1 − α ) ∑ j = 1 n 1 k w j h k − w j − 1 h k ( n − j + 1 ) 1 − α − ( n − j ) 1 − α + 1 ( 1 − α ) Γ ( 1 − α ) ∑ j = 1 n ( n − j + 1 ) 1 − α − ( n − j ) 1 − α O k 2 − α .

We now formally introduce the following mathematical definitions that are fundamental to our analysis:

(4.21) m α = ( 1 − α ) Γ ( 1 − α ) k n α − 1 ,

and

(4.22) m j = j 1 − α − ( j − 1 ) 1 − α .

Now that these have been established, we are ready to express:

(4.23) D t n α 0 C w n h k = m α ∑ j = 0 n m j w n − j + 1 h k − w n − 1 h k + n 1 − α ( 1 − α ) Γ ( 1 − α ) O k 2 − α = m α ∑ j = 0 n m j w n − j + 1 h k − w n − 1 h k + O ( k ) .

Subsequently, we infer that

(4.24) D t n α 0 C w n h k − D t n α 0 C w n X = m α ∑ n = 0 n m j w n − j + 1 h k − w n − j − w n − j + 1 h k − w n − j X ≤ m α ∑ j = 0 n m j w n − j − 1 − w n − j − 1 h k X + m α ∑ j = 0 n m j w n − j − w n − j h k X .

For the estimation of the term w n − j − 1 − w n − j − 1 h k X , it can be observed that:

(4.25) w n − j − 1 − w n − j − 1 h k = ∫ 0 t n − j − 1 δ w ( s ) d s + w 0 − ∑ l = 1 n − j − 1 δ w l h k k l − w 0 h = w 0 − w 0 h + ∑ l = 1 n − j − 1 δ w l − δ w l h k k l + ∑ l = 1 n − j − 1 ∫ t l − 1 t l δ w ( s ) d s − δ w l k l .

Also, we find

(4.26) ∑ l = 1 n − j − 1 ∫ t l − 1 t l δ w ( s ) d s − δ w l k l X = ∑ l = 1 n − j − 1 ∫ t l − 1 t l δ w ( s ) − δ w l d s X ≤ ∑ l = 1 n − j − 1 ∫ t l − 1 t l δ w ( s ) − δ w l X d s = I k ( δ w ) .

where

(4.27) I k ( δ w ) ≤ k ‖ w ‖ H 2 ( 0 , T ; X ) .

Then

(4.28) w n − j − 1 − w n − j − 1 h k X ≤ ∑ l = 1 n − j − 1 δ w l − δ w l h k X k l + w 0 − w 0 h X + k ‖ w ‖ H 2 ( 0 , T ; X ) .

We use j ≤ n ≤ N and Nk = T, we have

(4.29) ∑ j = 0 n k w n − j − 1 − w n − j − 1 h k X ≤ c T w 0 − w 0 h X + k 2 ‖ w ‖ H 2 ( 0 , T ; X ) + T ∑ j = 0 n k ∑ l = 1 n − j − 1 δ w l − δ w l h k X .

When we combine equations (4.9), (4.17), (4.19), (4.24) and (4.29), we obtain the following relation.

(4.30) w n − w n h k X ≤ c φ n h − w n X + w 0 − w 0 h X + R 2 + k 2 ‖ w ‖ H 2 ( 0 , T ; X ) + k ∑ j = 0 n w n − j − w n − j h k X + T ∑ j = 0 n k ∑ l = 1 n − j − 1 δ w l − δ w l h k X .

By using the definition of R ₂, as well as equations (2.2), (2.4) and (2.5), we can then establish the following:

(4.31) R 2 ≤ σ t n , ε φ n h − ε w n + j c w n , φ n h − w n + j f w n , φ n h − j f w n , w n − λ n , φ n h − w n X .

Utilizing (HP3)–(HP5), along with the regularity of u _n and equation (2.11), we can derive the following:

(4.32) R 2 ≤ c φ n h − w n X .

Similarly to the approach used in [35], we find that

(4.33) k ∑ l = 1 n − j − 1 δ w l − δ w l h k X ≤ c k 2 ‖ w ‖ H 2 ( 0 , T ; X ) .

By applying the generalized discrete variant of the Gronwall inequality from Lemma 4.4, we can conclude that:

(4.34) max 1 ≤ n ≤ N e n ≤ c g n + c ∑ j = 0 k e n − j ,

where

(4.35) e n = w n − w n h k X , e n − j = w n − j − w n − j h k X ,

and

(4.36) g n = φ n h − w n X + w 0 − w 0 h X + k 2 ‖ w ‖ H 2 ( 0 , T ; X ) .

Finally, we employ P h ( w n ) as the standard finite element interpolation operator for w. We then derive the interpolation error estimate, as discussed in [35].

(4.37) w n − P h w n X ≤ c h w n X , and w 0 − P h w 0 X ≤ c h w 0 X .

By combining equations (4.34) and (4.37) with condition (4.4), as k → 0, we arrive at the order estimate (4.5).□

Corresponding author: Boling Chen, Center for Applied Mathematics of Guangxi, and Guangxi Colleges and Universities Key Laboratory of Complex System Optimization and Big Data Processing, Yulin Normal University, Yulin 537000, Guangxi, P.R. China, E-mail: bolinchenylnu@163.com

Acknowledgments

The authors would like to thank Chunyu Li, Meiyan Zhou, Haimei Chen and Danlu Feng for their recommendations and remarks aiming at improving the manuscript in terms of clarity. This project has received funding from the Natural Science Foundation of Guangxi Grant Nos. GKAD23026237 and 2025GXNSFGA069001.

Authors contributions: All authors of this manuscript contributed equally to this work.
Data Availability: No data were used to support this study.
Conflicts of Interest: The authors declare that there are no conflicts of interest.
Ethical approval: Not available.

Appendix

We now recall several fundamental concepts from fractional calculus and nonlinear functional analysis, which have been previously developed in [12], 36], 37].

Definition 4.2

(Caputo fractional derivative of order 0 < α ≤ 1). Let X be a Banach space and consider a finite time interval (0, T). For any function f ∈ AC(0, T; X), where AC(0, T; X) denotes the space of absolutely continuous functions mapping (0, T) to X, the Caputo fractional derivative of order α is given by

D t α 0 C f ( t ) = 0 I t 1 − α f ′ ( t ) = 1 Γ ( 1 − α ) ∫ 0 t ( t − s ) − α f ′ ( s ) d s , ∀ t ∈ ( 0 , T ) .

Remark 4.3.

When α = 1, the Caputo derivative coincides with the classical first-order derivative:

D t 1 0 C f ( t ) = f ′ ( t ) for almost every t ∈ ( 0 , T ) .

Lemma 4.4.

Let T > 0 be given. For a positive integer N, define k = T/N. Assume that g n n = 1 N and e n n = 1 N are two sequence of nonnegative numbers satisfying

e n ≤ c g n + c ∑ j = 1 n k e j , n = 1 , … , N ,

for a positive constant c independent of N or k. Then, there exists a positive constant c, independent of N or k, such that

max 1 ≤ n ≤ N e n ≤ c max 1 ≤ n ≤ N g n .

Definition 4.5.

Let X be a reflexive Banach space and consider an operator T : X → 2 X * satisfying:

For each v ∈ X, the image Tv is nonempty, closed, and convex in X*;
T is a bounded operator;
For any weakly convergent sequences v _n ⇀ v in X and v n * ⇀ v * in X* with v n * ∈ T v n , if

lim sup n → ∞ ⟨ v n * , v n − v ⟩ ≤ 0 ,

then v* ∈ Tv and ⟨ v n * , v n ⟩ → ⟨ v * , v ⟩ .

Then the operator A is pseudomonotone.

Definition 4.6

(Clarke generalized directional derivative and generalized gradient). Let J : X → R be a locally Lipschitz function. We denote by J ⁰(u, v) the Clarke generalized directional derivative of J at the point x ∈ X in the direction y ∈ X is defined by

J 0 ( x , y ) = lim sup λ → 0 + , z → x J ( z + λ y ) − J ( z ) λ .

The generalized gradient of J : X → R at x ∈ X is defined by

∂ J ( x ) = ξ ∈ X * , / J 0 ( x , y ) ≥ < ξ , y > X * , X , for all y ∈ X .

References

[1] E.-H. Essoufi and A. Zafrar, Dual methods for frictional contact problem with electroelastic-locking materials. Optimization, Taylor and Francis, 2020, pp. 1–28.10.1080/02331934.2020.1745794Search in Google Scholar

[2] W. Prager, On ideal locking materials, Transactions of the Society of Rheology, J. Rheol. 1 (1957), no. 1, 169–175, https://doi.org/10.1122/1.548818.Search in Google Scholar

[3] S. Migorski, A. Ochal, and M. Sofonea, Nonlinear inclusions and hemivariational inequalities. Models and analysis of contact problems. Advances in Mechanics and Mathematics 26, Springer, New York, 2013.10.1007/978-1-4614-4232-5Search in Google Scholar

[4] S. Bourichi and E. L.-H. Essoufi, Penalty method for an unilateral contact problem with Coulomb’s friction for locking materials, Int. J. Math. Model. Comput. 6 (2016), no. 1, 61–81.Search in Google Scholar

[5] M. Bouallala and E. H. Essoufi, A thermo-viscoelastic fractional contact problem with normal compliance and coulomb’s friction, J. Math. Phys. Anal. Geom. 17 (2021), no. 3, 280–294, DOI https://doi.org/10.15407/mag17.03.280.Search in Google Scholar

[6] W. Prager, Elastic solids of limited compressibility, Proceedings of the 9th International Congress of Applied Mechanics, vol. 5, 1958, pp. 205–211.Search in Google Scholar

[7] W. Prager, On elastic, perfectly locking materials, Appl. Mech. (1966), 538–544.10.1007/978-3-662-29364-5_72Search in Google Scholar

[8] A. G. Radwan, Resonance and quality factor of the RLαCα fractional circuit, IEEE J. Emerg. Sel. Top. Circuits Syst. 3 (2013), no. 3, 377–385, DOI https://doi.org/10.1109/jetcas.2013.2272838.Search in Google Scholar

[9] S. Shen, F. Liu, J. Chen, I. Turner, and V. Anh, Numerical techniques for the variable order time fractional diffusion equation, Appl. Math. Comput. 218 (2012), 10861–10870, DOI https://doi.org/10.1016/j.amc.2012.04.047.Search in Google Scholar

[10] C. Zou, L. Zhang, X. Hu, Z. Wang, T. Wik, and M. Pecht, A review of fractional-order techniques applied to lithium-ion batteries, lead-acid batteries, and supercapacitors, J. Power Sources 390 (2018), 286–296, DOI https://doi.org/10.1016/j.jpowsour.2018.04.033.Search in Google Scholar

[11] S. Migorski and A. Ochal, Quasi-static hemivariational inequality via vanishing acceleration approach, SIAM J. Math. Anal. 41 (2009), 1415–1435, DOI https://doi.org/10.1137/080733231.Search in Google Scholar

[12] C. Carstensen and J. Gwinner, A theory of discretization for nonlinear evolution inequalities applied to parabolic Signorini problems, Ann. Mat. Pura Appl. 177 (1999), no. 1, 363–394, DOI https://doi.org/10.1007/bf02505918.Search in Google Scholar

[13] F. Demengel and P. Suquet, On locking materials, Acta Appl. Math. 6 (1986), no. 2, 185–211, DOI https://doi.org/10.1007/bf00046725.Search in Google Scholar

[14] C. Eck and J. I. Jarusek, Existence results for the static contact problem with Coulomb friction, Math. Models Methods Appl. Sci. 8 (1998), no. 03, 445–468, DOI https://doi.org/10.1142/s0218202598000196.Search in Google Scholar

[15] M. Bouallala, Variational and error estimation for a frictionless contact problem in thermo-viscoelasticity with time fractional derivatives, Comput. Anal. Comput. 2 (2024), 1–18, DOI https://doi.org/10.3934/cac.2024001.Search in Google Scholar

[16] J. Cen, M. Bouallala, and A. Bendarag, Modeling and numerical analysis of a viscoelastic contact problem for locking materials with time-fractional derivatives. J. Georgian Math. (2025).10.1515/gmj-2025-2071Search in Google Scholar

[17] I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of their Solution and Some of Their Applications, Elsevier, 1998.Search in Google Scholar

[18] S. Migórski and J. Ogorzały, A variational-hemivariational inequality in contact problem for locking materials and nonmonotone slip dependent friction, Acta Math. Sci. 37 (2017), no. 6, 1639–1652, DOI https://doi.org/10.1016/s0252-9602(17)30097-8.Search in Google Scholar

[19] M. Bouallala, E. H. Essoufi, and S. Liao, A problem involving time-fractional contact with normal compliance and Tresca’s friction within the realm of electro-viscoelasticity, Chaos Solitons Fractals 188 (2024), 115484, DOI https://doi.org/10.1016/j.chaos.2024.115484.Search in Google Scholar

[20] F. Mainardi, Fractional calculus and waves in linear viscoelasticity: An introduction to mathematical models, World Sci. (2022).Search in Google Scholar

[21] M. Sofonea, A nonsmooth static frictionless contact problem with locking materials, Anal. Appl. (World Sci.) 16 (2018), no. 6, 851–874, DOI https://doi.org/10.1142/s0219530518500215.Search in Google Scholar

[22] M. Sofonea and S. Migórski, Variational-Hemivariational Inequalities with Applications, CRC Press, 2017.10.1201/9781315153261Search in Google Scholar

[23] J. A. C. Martins and J. T. Oden, Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws, Nonlinear Anal. Theory Methods Appl. 11 (1987), no. 3, 407–428, DOI https://doi.org/10.1016/0362-546x(87)90055-1.Search in Google Scholar

[24] C. Ionescu, A. Lopes, D. Copot, J. T. Machado, and J. H. Bates, The role of fractional calculus in modeling biological phenomena. A review, Commun. Nonlinear Sci. Numer. Simul. 51 (2017), 141–159, DOI https://doi.org/10.1016/j.cnsns.2017.04.001.Search in Google Scholar

[25] Y. Ayyad, M. Barboteu, and J. R. Fernandez, A frictionless viscoelastodynamic contact problem with energy consistent properties: Numerical analysis and computational aspects, Comput. Methods Appl. Mech. Engrg. 198 (2009), no. 5–8, 669–679, DOI https://doi.org/10.1016/j.cma.2008.10.004.Search in Google Scholar

[26] M. Barboteu, W. Han, and S. Migórski, On numerical approximation of a variational-hemivariational inequality modeling contact problems for locking materials, Comput. Math. Appl. 77 (2019), no. 11, 2894–2905.10.1016/j.camwa.2018.08.004Search in Google Scholar

[27] D. Paola, R. Heuer, and A. Pirrotta, Fractional visco-elastic Euler–Bernoulli beam, Int. J. Solids Struct. 50 (2013), no. 22–23, 3505–3510, DOI https://doi.org/10.1016/j.ijsolstr.2013.06.010.Search in Google Scholar

[28] S. Zeng, Z. Liu, and S. Migorski, A class of fractional differential hemivariational inequalities with application to contact problem, Z. Angew. Math. Phys. 69 (2018), no. 2, 1–23, DOI https://doi.org/10.1007/s00033-018-0929-6.Search in Google Scholar

[29] M. Bouallala, E. L.-H. Essoufi, V. T. Nguyen, and W. Pang, A time-fractional of a viscoelastic frictionless contact problem with normal compliance, Eur. Phys. J. Spec. Top. 232 (2023), 1–10, DOI https://doi.org/10.1140/epjs/s11734-023-00962-x.Search in Google Scholar

[30] P. D. Panagiotopoulos, Hemivariational Inequalities, Applications in Mechanics and Engineering, Springer, Berlin, 1993.10.1007/978-3-642-51677-1Search in Google Scholar

[31] W. Han and M. Sofonea, Evolutionary variational inequalities arising in viscoelastic contact problems, SIAM J. Numer. Anal. 38 (2002), no. 2, 556–579, DOI https://doi.org/10.1137/s0036142998347309.Search in Google Scholar

[32] Z. Naniewicz and P. D. Panagiotopoulos, Mathematical theory of hemivariational inequalities and applications, Dekker, New York, 1995.Search in Google Scholar

[33] S. Zeng and S. Migórski, A class of time-fractional hemivariational inequalities with application to frictional contact problem, Commun. Nonlinear Sci. Numer. Simul. 56 (2018), 34–48, DOI https://doi.org/10.1016/j.cnsns.2017.07.016.Search in Google Scholar

[34] J. T. Oden and E. B. Pires, Nonlocal and nonlinear friction laws and variational principles for contact problems in elasticity, J. Appl. Mech. 50 (1983), no. 1, 7–76, DOI https://doi.org/10.1115/1.3167019.Search in Google Scholar

[35] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo, Theory and Applications of Fractional Differential Equations, vol. 204, Elsevier Science Limited, 2006.Search in Google Scholar

[36] J. Necas and I. Hlavácek, Mathematical Theory of Elastic and Elasto-Plastic Bodies: An Introduction, Elsevier, 2017.Search in Google Scholar

[37] M. Sofonea, History-dependent inequalities for contact problems with locking materials, J. Elasticity 134 (2019), no. 2, 127–148, DOI https://doi.org/10.1007/s10659-018-9684-3.Search in Google Scholar

Received: 2025-04-14

Accepted: 2025-10-20

Published Online: 2025-12-04

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/dema-2025-0194

Keywords for this article

fractional viscoelastic constitutive law of Caputo type; contact problem; hemivariational inequality; weak solution; finite element method; error estimate

Creative Commons

BY 4.0