Mircea Sofonea

Mikael Barboteu

Weimin Han
email: weimin-han@uiowa.edu

Numerical solution of a contact problem with unilateral constraint and history-dependent penetration

Keywords: Augmented Lagrangian method, Finite element scheme, History-dependent stiffness coefficient, Normal compliance, Numerical simulations, Unilateral constraint

diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Contact phenomena involving deformable bodies lead to nonsmooth, nonlinear mathematical problems. Analysis of the problems, including existence and uniqueness results, is carried out in a large number of works, see for instance [START_REF] Eck | Unilateral Contact problems: variational methods and existence theorems[END_REF][START_REF] Han | Quasistatic contact problems in viscoelasticity and viscoplasticity[END_REF][START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF][START_REF] Shillor | Models and variational analysis of quasistatic contact[END_REF] and the references therein. Numerical analysis of the problems, including error estimation for discrete schemes and numerical simulations, can be found in [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Hlaváček | Solution of variational inequalities in mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][8][START_REF] Wriggers | Computational contact mechanics[END_REF]. The state-of-the-art in the field, including applications in engineering, could be found in the recent special issue [START_REF]Special issue on contact mechanics[END_REF].

In this paper, we provide the numerical solution of a frictionless contact problem for viscoplastic materials. It models the quasistatic contact of a body with an obstacle made of a hard material covered with a thin layer of soft material, in the framework of the small deformation theory. The challenging feature of the problem is mainly caused by the forms of the constitutive law and the boundary condition on the contact surface. We describe the constitutive law of the material by a rate-type relation of the form

σ = Eε( u) + G(σ , ε(u)), (1.1)
The boundary Ŵ is partitioned into three disjoint measurable parts Ŵ 1 , Ŵ 2 , and Ŵ 3 , such that meas (Ŵ 1 )>0, and we denote that [0, T ] be the time interval of interest with some T > 0. The body is clamped on Ŵ 1 × (0, T ),i s submitted to the action of body forces, and surface tractions on × (0, T ) and Ŵ 2 × (0, T ), respectively, and is in contact with a deformable foundation on Ŵ 3 × (0, T ). As a result, its mechanical state evolves, and the evolution is governed by the following mathematical model.

Problem P Find a displacement field u : ×[0, T ]→R d and a stress field σ : ×[0, T ]→S d such that

σ = Eε( u) + G(σ , ε(u)) in × (0, T ), (2.1)
Div σ + f 0 = 0 in × (0, T ), (2.2)

u = 0 on Ŵ 1 × (0, T ), (2.3
)

σν = f 2 on Ŵ 2 × (0, T ), (2.4
)

u ν ≤ g,σ ν + K (ξ(u ν )) p(u ν ) ≤ 0, (u ν -g)(σ ν + K (ξ(u ν )) p(u ν )) = 0 on Ŵ 3 × (0, T ), (2.5
)

σ τ = 0 on Ŵ 3 × (0, T ), (2.6
)

u(0) = u 0 , σ (0) = σ 0 in . (2.7)
Let us comment on the equations and boundary conditions in (2.1)-(2.7) in which, for simplicity, we do not indicate explicitly the dependence of various functions on the variables x or t. Also, here and below we use the index ν and τ to denote the normal and tangential components of vectors and tensors, respectively. Equation (2.1) represents the viscoplastic constitutive law of the material, already introduced in Sect. 1. Equation (2.2) is the equilibrium equation in which f 0 represents the density of body forces; we use it here since the process is assumed to be quasistatic. In (2.1)-(2.2) and below, ε and Div are the deformation and the divergence operators, respectively, defined by

ε(u) = (ε ij (u)), ε ij (u) = 1 2 (u i, j + u j,i ), Div σ = (σ ij, j ),
where an index that follows a comma represents the partial derivative with respect to the corresponding component of the spatial variable, e.g., u i, j = ∂u i /∂ x j . Conditions (2.3) and (2.4) are the displacement and traction boundary conditions, respectively, with f 2 being the density of surface tractions. Condition (2.6) shows that the tangential stress on the contact surface vanishes, corresponding to a frictionless contact. In (2.7), u 0 and σ 0 are the initial displacement and the initial stress field, respectively. We turn to the contact condition (2.5). First, the penetration is limited by a maximal depth g > 0. In other words, at any time t, the normal displacement on the contact surface u ν (t) satisfies the inequality:

u ν (t) ≤ g.
(2.8)

Next, we assume that the normal stress has an additive decomposition of the form:

σ ν (t) = σ D ν (t) + σ R ν (t).
(2.9)

The part σ D ν (t) satisfies a normal compliance contact condition:

-σ D ν (t) = Kp(u ν (t)). (2.10)
Here p is a given nonnegative function which vanishes for negative arguments, and K represents the stiffness coefficient of the foundation. The normal compliance contact condition was used in many publications, see, e.g., [START_REF] Klarbring | Frictional contact problems with normal compliance[END_REF][START_REF] Klarbring | On friction problems with normal compliance[END_REF][START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF][START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] and references therein. The part σ R ν (t) of the normal stress satisfies the Signorini unilateral condition:

σ R ν (t) ≤ 0,σ R ν (t)(u ν (t) -g) = 0. (2.11)
We combine (2.9) and (2.10) to see that

σ R ν (t) = σ ν (t) + Kp(u ν (t)).
(2.12)

Then, we substitute equality (2.12)in(2.11) and use (2.8), to obtain the condition:

u ν (t) ≤ g,σ ν (t) + Kp(u ν (t)) ≤ 0,( u ν (t) -g) σ ν (t) + Kp(u ν (t)) = 0. (2.13)
This condition can be naturally interpreted as follows: the foundation is made of a hard material covered by a thin layer of a softer material with thickness g; the soft material has an elastic behavior, i.e., it is deformable and allows penetration; the contact with this layer is modeled with normal compliance; the hard material is perfectly rigid, and, therefore, it does not allow penetration; the contact with this material is modeled using the Signorini contact condition. To conclude, condition (2.13) models a foundation which has a rigid-elastic behavior: its elastic behavior is given by the layer of the soft material, while its rigid behavior is given by the hard material. Now, consider the more complicated situation where the stiffness coefficient can depend on the accumulation of penetration, denoted as ξ(u ν ):

K = K (ξ(u ν )(t)), (2.14)
where

ξ(u ν )(t) = t 0 u + ν (s) ds, t ∈[0, T ], (2.15)
where u + ν is the positive part of u ν . Combining (2.13) and (2.14) leads to the contact condition (2.5). Note that the dependence K = K (ξ(u ν )(t)) allows us to take into account the variation of the elastic modulus of the thin layer, as the cycles of contact and separation proceed and, therefore, it models the hardening or softening phenomenon. Practical examples of surface hardening or softening abound in industry and various engineering settings, see, e.g., [START_REF] Piotrowski | Smoothing dry friction damping by dither generated in rolling contact of wheel and rail and its influence on ride dynamics of freight wagons[END_REF].

To proceed with a variational formulation of the contact problem P, we need some function spaces. We use the standard notation for Sobolev and Lebesgue spaces associated to and Ŵ. In addition, let

V = u = (u i ) : u i ∈ H 1 ( ), u i = 0a .e. on Ŵ 1 , Q = σ = (σ ij ) : σ ij = σ ji ∈ L 2 ( ) , Q 1 = σ ∈ Q : Div σ ∈ L 2 ( ) d .
These are real Hilbert spaces endowed with the canonical inner products:

(u, v) V = ε(u) • ε(v) dx, (σ , τ ) Q = σ • τ dx, (σ , τ ) Q 1 = σ • τ dx + Divσ • Div τ dx.
The associated norms on these spaces are denoted by • V , • Q , and • Q 1 , respectively. Recall that completeness of the space (V, • V ) follows from the assumption meas(Ŵ 1 )>0, which allows the use of Korn's inequality. For an element v ∈ V , we still write v for the trace of v, and we denote by v ν and v τ the normal and tangential components of v, respectively, on Ŵ given by v ν = v • ν, v τ = vv ν ν. Also, for a regular stress function σ ,weuse the notations σ ν and σ τ , respectively, for the normal and the tangential traces, i.e., σ ν = (σν)•ν and σ τ = σν-σ ν ν.

For every Banach space (X, • X ), we use the notation C([0, T ]; X ) for the space of continuous functions defined on [0, T ] with values on X . It is known that C([0, T ]; X ) is a real Banach space with the norm:

v C([0,T ];X ) = max t∈[0,T ] v(t) X .
Moreover, for a subset K ⊂ X , we still use the symbol C([0, T ]; K ) for the set of continuous functions defined on [0, T ] with values in K .

In the study of the mechanical problem (2.1)-(2.7), we assume that the elasticity tensor E, the constitutive function G, the stiffness coefficient function K , and the normal compliance function p satisfy the following conditions:

⎧ ⎨ ⎩ (a) E = (E ijkl ) : × S d → S d . (b) E ijkl = E kli j = E jikl ∈ L ∞ ( ), 1 ≤ i, j, k, l ≤ d. (c) There exists m E > 0 such that Eτ • τ ≥ m E τ 2 for all τ ∈ S d , a.e. in .
(2.16)

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ (a) G : × S d × S d → S d . (b) There exists L G > 0 such that G(x, σ 1 ,ε 1 ) -G(x, σ 2 ,ε 2 ) ≤L G ( σ 1 -σ 2 + ε 1 -ε 2 ) for all σ 1 , σ 2 ,ε 1 ,ε 2 ∈ S d , a.e. x ∈ . (c) The mapping x → G(x, σ ,ε)is measurable on , for all σ ,ε ∈ S d . (d)
The mapping x → G(x, 0, 0) belongs to Q.

(2.17)

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ (a) K : Ŵ 3 × R + → R + . (b) There exists L K > 0 such that |K (x, r 1 ) -K (x, r 2 )|≤L K |r 1 -r 2 | for all r 1 , r 2 ∈ R + , a.e. x ∈ Ŵ 3 . (c) The mapping x → K (x, r ) is measurable on Ŵ 3 , for all r ∈ R + . (d) The mapping x → K (x, 0) belongs to L 2 (Ŵ 3 ). (2.18) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (a) p : Ŵ 3 × R → R + . (b) There exists L p > 0 such that | p(x, r 1 ) -p(x, r 2 )|≤L p |r 1 -r 2 | for all r 1 , r 2 ∈ R, a.e. x ∈ Ŵ 3 . (c) The mapping x → p(x, r ) is measurable on Ŵ 3 , for all r ∈ R. (d) ( p(x, r 1 ) -p(x, r 2 ))(r 1 -r 2 ) ≥ 0 ∀ r 1 , r 2 ∈ R, a.e. x ∈ Ŵ 3 .
(e) p(x, r ) = 0 for all r ≤ 0, a.e. x ∈ Ŵ 3 .

(2. [START_REF] Klarbring | Frictional contact problems with normal compliance[END_REF] For the body forces and tractions densities, we assume

f 0 ∈ C([0, T ]; L 2 ( ) d ), f 2 ∈ C([0, T ]; L 2 (Ŵ 2 ) d ).
(2.20)

For the initial data, we assume

u 0 ∈ V, σ 0 ∈ Q. (2.21)
Next, we introduce the set of admissible displacements defined by

U = v ∈ V : v ν ≤ g on Ŵ 3 . (2.22)
Moreover, we define the functions q : R → R and f :[0, T ]→V by the equalities:

q(x, s) = s 0 p(x, r ) dr ∀ r ∈ R, a.e. x ∈ Ŵ 3 , (2.23)
(f (t), v) V = f 0 (t) • v dx + Ŵ 2 f 2 (t) • v da ∀ v ∈ V, t ∈[0, T ]. (2.24)
Assume now that (u, σ ) are sufficiently regular functions which satisfy (2.1)-(2.7), and let t > 0 be given. We integrate Eq. (2.1) and make use of the initial conditions (2.7) to obtain

σ (t) = Eε(u(t)) + t 0 G(σ (s), ε(u(s))) ds + σ 0 -Eε(u 0 ).
(2.25)

Then, we use integration by parts, the equilibrium equation (2.2), the boundary conditions (2.3)-(2.6), and notation (2.23), (2.24) to see that

(σ (t), ε(v) -ε(u(t))) Q + Ŵ 3 K (ξ(u ν )(t)) (q(v ν (t)) -q(u ν (t))) da ≥ (f (t), v -u(t)) V ∀ v ∈ V. (2.26)
Finally, we note that the first inequality in (2.5) combined with the definition (2.22)impliesthat

u(t) ∈ U. (2.27)
We gather (2.25)-(2.27) together to obtain the following variational formulation of Problem P.

Problem P V Find a displacement field u :[0, T ]→U and a stress field σ :[0, T ]→Q such that

σ (t) = Eε(u(t)) + t 0 G(σ (s), ε(u(s))) ds + σ 0 -Eε(u 0 ), (2.28)
(σ (t), ε(v) -ε(u(t))) Q + Ŵ 3 K (ξ(u ν )(t)) (q(v ν (t)) -q(u ν (t))) da ≥ (f (t), v -u(t)) V ∀ v ∈ U, (2.29)
for all t ∈[0, T ].
The unique solvability of Problem P V is given by the following result, proved in [START_REF] Sofonea | A viscoplastic contact model with normal compliance, unilateral constraint and history-dependent stiffness coefficient[END_REF].

Theorem 2.1 Assume (2.16)-(2.21). Then Problem P V has a unique solution. Moreover, the solution satisfies

u ∈ C([0, T ]; U ), σ ∈ C([0, T ]; Q 1 ). (2.30)
We assume in what follows that (2.16)-(2.21) hold, and we denote by (u, σ ) the solution of Problem P V stated in Theorem 2.1. For each ρ>0, let K ρ be a perturbation of the function K which satisfies the condition (2.18) with a positive constant L K ρ . Consider the following perturbation of Problem P V .

Problem P ρ V Find a displacement field u ρ :[0, T ]→U and a stress field σ ρ :[0, T ]→Q such that σ ρ (t) = Eε(u ρ (t)) + t 0 G(σ ρ (s), ε(u ρ (s))) ds + σ 0 -Eε(u 0 ), (2.31)
(σ ρ (t), ε(v) -ε(u ρ (t))) Q + Ŵ 3 K ρ (ξ(u ρν )(t)) (q(v ρν (t)) -q(u ρν (t))) da ≥ (f (t), v -u ρ (t)) V ∀ v ∈ U, ( 2
∈ C([0, T ]; U ), σ ρ ∈ C([0, T ]; Q 1 ).
In addition, we have the following result, as a consequence of Theorem 5.1 in [START_REF] Sofonea | A viscoplastic contact model with normal compliance, unilateral constraint and history-dependent stiffness coefficient[END_REF].

Corollary 2.2 Assume that there exists a function G

: R + → R + such that (a) |K ρ (r ) -K (r )|≤G(ρ) ∀ r ∈[0, gT], for each ρ>0. (b) G(ρ) → 0 as ρ → 0.
(2.33)

Then the solution (u ρ , σ ρ ) of Problem P ρ V converges to the solution (u, σ ) of Problem P V :

u ρ → u in C([0, T ]; V ), σ ρ → σ in C([0, T ]; Q 1 ) as ρ → 0. (2.34)
Note that Corollary 2.2 states that the weak solution of the viscoplastic contact problem P depends continuously on the stiffness coefficient.

Numerical solution

In this section, we discuss the numerical solution of Problem P. To this end, we use arguments similar to those used in [START_REF] Barboteu | Analysis of quasistatic viscoplastic contact problems with normal compliance[END_REF][START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF][START_REF] Barboteu | Analysis of a contact problem with normal compliance, finite penetration and nonmonotone slip dependent friction[END_REF], based on an adapted combination of the penalty method and the augmented Lagrangian method for the numerical treatment of the specific contact condition (2.5). The starting point of our method is an alternative variational formulation of Problem P by considering a Lagrange multiplier associated with the normal contact stress. Then, an augmented Lagrangian method is naturally used to solve the resulting alternative hybrid variational problem.

A hybrid variational formulation

In order to derive the hybrid variational problem, we introduce a trace space

X ν ={v ν | Ŵ 3 : v ∈ V }. We denote by X ′ ν the dual of X ν and let •, • X ′ ν ,X ν be the duality pairing. Introduce an operator K : X ν → X ′ ν and a function I : X ν → (-∞, +∞] by K u ν ,v ν X ′ ν ,X ν = Ŵ 3 K (ξ(u ν )) p(u ν )v ν da ∀ u ν ,v ν ∈ X ν , I(u ν ) = Ŵ 3 I R -(u ν -g) da ∀ u ν ∈ X ν ,
where I R -represents the indicator function of the set R -= (-∞, 0]. We note that, for all t ∈ (0, T ), the contact condition (2.5) is equivalent to the subdifferential inclusion:

-σ ν (t) ∈ ∂I(u ν | Ŵ 3 (t)) + Ku ν | Ŵ 3 (t) in X ′ ν , (3.1)
where ∂I denotes the subdifferential of I in the sense of the convex analysis. Inclusion (3.1) allows us to supplement the unknowns of the problem with an additional one, the so-called Lagrange multiplier. This leads to the following hybrid variational formulation of the contact problem P.

Problem P V Find a displacement field u :[ 0, T ]→V , a stress field σ :[ 0, T ]→Q and a Lagrange multiplier

λ ν :[0, T ]→X ′ ν such that σ (t) = Eε(u(t)) + t 0 G(σ (s), ε(u(s))) ds + σ 0 -Eε(u 0 ), (3.2)
(σ (t), ε(v)) Q -λ ν (t), v ν| Ŵ 3 X ′ ν ,X ν = (f (t), v) V ∀ v ∈ V, (3.3) -λ ν (t) ∈ ∂I(u ν | Ŵ 3 (t)) + Ku ν | Ŵ 3 (t) (3.4) for all t ∈[0, T ].
Note that Problem P V represents an alternative to the variational Problem P V , and is formulated in terms of three unknown fields.

Numerical approximation

To describe the numerical method for the hybrid variational Problem P V , we first introduce some preparatory material. We use a finite element method for the spatial discretization of the domain and a uniform discretization of the time interval [0, T ]. Assume that is a polyhedral domain, and we denote by {T h } a regular family of triangular finite element partitions of that are compatible with the boundary decomposition:

Ŵ = Ŵ 1 ∪ Ŵ 2 ∪ Ŵ 3 .
Here and below, h represents the spatial discretization parameter. In the numerical examples presented in the next section, we approach the space V by the finite dimensional space of continuous piecewise affine functions, denoted V h . The space Q is approximated by the finite element space of piecewise constants, denoted Q h .F o r the discretization of the Lagrange multiplier λ, we consider a discrete space Y h ν ⊂ X ′ ν ∩ L 2 (Ŵ 3 ). For the time discretization, we use evenly spaced nodes t n = nk,0≤ n ≤ N , where N > 0 is an integer and k = T N is the time step size. For a continuous function v(t) with values in a function space, we write v j = v(t j ) for 0 ≤ j ≤ N . Details about the discretization step can be found in [START_REF] Wriggers | Computational contact mechanics[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers[END_REF][START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF].

We now describe the numerical solution of the hybrid variational Problem P V . The numerical treatment of the condition (3.4) is based on a combination of the penalty method for the normal compliance contact with an adapted augmented Lagrangian argument for the unilateral condition. To develop the Lagrangian approach, we need to introduce additional fictitious nodes for the Lagrange multiplier in the initial mesh. The introduction of these nodes leads to specific contact elements in relation to the geometrical discretization of the interface Ŵ 3 . For the numerical examples in the next section, we consider a "node-to-rigid" contact element, which is composed of one node of Ŵ 3 and one Lagrange multiplier node. Details on this construction can be found in [START_REF] Barboteu | Analysis of quasistatic viscoplastic contact problems with normal compliance[END_REF][START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF][START_REF] Barboteu | Analysis of a contact problem with normal compliance, finite penetration and nonmonotone slip dependent friction[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]. Then, the numerical approximation of Problem P V leads at each time step n to the solution of a system of nonlinear equations of the form:

R(u, λ) = G(u) + F(u, λ) = 0, (3.5)
where the operator G comes from the discretization of the elastic terms of the alternative variational Problem P V and the operator F comes from the discretization of contact terms. In the nonlinear system (3.5), the unknowns are

the generalized discrete displacement field u ={ u i } N h tot i=1 ∈ R d•N h
tot , and the Lagrange multiplier generalized vector

λ ={λ i } N h Ŵ 3 i=1 ∈ R d•N h
Ŵ 3 in which N h tot and N h Ŵ 3 are the total number of nodes issued from the discretization of and Ŵ 3 , respectively. Here u i represents the value of the function u h at the i-th node of T h , and λ i denotes the value of the function λ h at the i-th node of the discretized contact interface. In addition,

G(u) ∈ R d•N h tot × R d•N h Ŵ 3 represents
the generalized elastic term and is defined by

G(u) = (G(u), 0 d•N h Ŵ 3
), where

0 d•N h Ŵ 3 is the zero element of R d•N h Ŵ 3
and the term G(u) ∈ R d•N h tot is given by the equality:

(G(u) • v) R d×N h tot = (σ h ,ε(v h )) Q -(f , v h ) V ∀ v ={v i } N h tot i=1 , ∀ v h ∈ V h ,
where σ h is related to u h by the discrete constitutive law (3.2). The contact operator F(u, λ) in equation (3.5)is related to the contact condition, (3.4). It is defined by

(F (u, λ) • (v, γ )) R d•N h tot ×R d•N h Ŵ 3 = Ŵ 3 (λ h ν + r (u h ν -g)) + ν • v h da + Ŵ 3 -1 r λ h -(λ h ν + r (u h ν -g)) + ν • γ h da + Ŵ 3 K (ξ(χ [0,g] (u h ν )u h ν ) p(χ [0,g] (u h ν )u h ν )ν • v h da for u, v ∈ R d•N h tot , λ, γ ∈ R d•N h Ŵ 3 , u h , v h ∈ V h , and λ h , γ h ∈ Y h ν ,
where r is a positive penalty coefficient, and χ E represents the characteristic function of the set E. In contrast to the problem studied in [START_REF] Barboteu | Analysis of quasistatic viscoplastic contact problems with normal compliance[END_REF], here we have to consider the discretization of the history-dependent term ξ(χ [0,g] (u h ν )u h ν ) and, to this end, we used an implicit rectangular method by approximating ξ(r )(t n ) with n i=1 k (r + ) i . The nonlinear system (3.5) is solved by a generalized Newton method. The pair ũ = (u, λ) is computed simultaneously by using the following linearized iterative scheme:

For i = 0,..., Solve (K i + T i ) ũi+1 =-R( ũi ) ũi+1 = ũi + ũi+1 , (3.6)
until convergence.

At each Newton iteration of index i, the linear symmetric system (3.7) is solved by using a conjugate gradient method with efficient preconditioners, in order to overcome the poor conditioning of the system matrix due to the contact terms. For instance, we use an ILU preconditioner based on an Element-By-Element procedure, see [START_REF] Alart | Solutions of frictional contact problems using an EBE preconditioner[END_REF] for details. The coefficient matrix of the linear system (3.7) is the sum of an elastic stiffness matrix K i and a contact tangent matrix T i . Note that K i = D ũG(u) represents the differential of the functions G with respect to the variable ũ, and T i ∈ ∂ ũF ( ũi ) represents the generalized Jacobian of F at ũi . Due to the fact that the measure of the set of the points in which the function F is nondifferentiable is zero, T i is reduced to a single-valued classical Jacobian matrix. In the case of the "node-to-rigid" contact element composed of one elastic node of Ŵ 3 and one contact Lagrange multiplier node, the elementary contact tangent matrix T ec i takes the form:

T ec i = r (M + N) M M 1 r (M -I) .
Here I is the 2 × 2 identity matrix, N and M are 2 × 2 block matrices which are defined according to the contact state:

N = 0 if u ν < 0, u + ν ν ⊗ D u K (ξ ) + K (ξ )ν ⊗ ν if 0 ≤ u ν < g, and M = 0 if u ν < g, ν ⊗ ν if u ν ≥ g,
where 0 is the 2 × 2 zero matrix, and u ⊗ v represents the tensor product of two vectors u and v, defined by (u ⊗ v) l,m = u l v m in which l, m = 1,...,d. The global contact tangent matrix T i is obtained by a finite element assembly procedure from all the contact elements

T j i ,1≤ j ≤ N h Ŵ 3 .

Numerical simulations

The aim of this section is to present some results of numerical simulations. We pay particular attention to the mechanical interpretation of the contact condition (2.5). We complete them with two parametric studies which illustrate the dependence of the numerical solution with respect to the stiffness coefficient and the numerical discretization parameters.

The physical setting

A representative academic example of compression of a ball against a foundation is considered. Due to the symmetry of the problem, we only need to consider the physical setting depicted in Fig. 1. There,

= (x 1 , x 2 ) ∈ R 2 : x 1 > 0, x 2 1 + (x 2 -1) 2 < 1 ,
Ŵ 1 ={0}×[0, 2], Ŵ 2 = (x 1 , x 2 ) ∈ R 2 : x 1 ≥ 0, x 2 ≥ 1, x 2 1 + (x 2 -1) 2 = 1 , Ŵ 3 = (x 1 , x 2 ) ∈ R 2 : x 1 ≥ 0, x 2 < 1, x 2 1 + (x 2 -1) 2 = 1 .
The domain represents the cross section of a three-dimensional deformable body subject to the action of tractions in such a way that a plane stress hypothesis is valid. The horizontal component of the displacement field vanishes on Ŵ 1 and vertical tractions act on Ŵ 2 . No body forces are assumed to act on the body during the process. The body is in frictionless contact with an obstacle on the part Ŵ 3 of its boundary. For the discretization in Fig. 1,w eu s e 30,808 elements with 128 elements containing sides on the contact boundary. We assume that the mechanical response of the ball is purely elastic, i.e., we model the material's behavior with a constitutive law of the form (1.1) in which the function G vanishes. In addition, we assume that the material is homogeneous and isotropic. We denote by E and κ the Young modulus and the Poisson ratio of the material, respectively, and use δ ij for the Kronecker symbol. Then, the elasticity tensor E is given by

(Eτ ) ij = Eκ (1 + κ)(1 -2κ) (τ ii )δ ij + E 1 + κ τ ij , 1 ≤ i, j ≤ 2, (4.1)
where the summation convention is used. For the numerical experiments, we use the following data:

E = 1000 N/m 2 ,κ = 0.3, f 0 = (0, 0) N/m 2 , f 2 = (-1000 × t, 0) N/mo n Ŵ 2 , p(r ) = r + , K (ξ ) = c ν (1 + αξ) β where ξ(r )(t n ) = n i=0 k (r + ) i , c ν = 100 N/m 2 ,α = 20,β =±1, g = 0.2m, T = 1.1s, N = 11, k = 0.1.
We note that for β = 1 the function ξ → K (ξ ) is increasing and, therefore, the stiffness coefficient increases with the history of the penetration. This behavior models the hardening process of the foundation. In contrast, for β =-1, the function ξ → K (ξ ) is decreasing, and, therefore, the stiffness coefficient decreases with the history of the penetration. This behavior models the softening process of the foundation.

t = 0.8s

Fig. 3 Evolution of deformed meshes and contact interfaces forces in the case β =-1 each time increment of the process, for both cases of history dependence. We note that the coefficient grows from 100 N/m 2 to 323 N/m 2 for the case β = 1, while the coefficient decreases from 100 N/m 2 to 21 N/m 2 for the case β =-1.

We now compare the contact condition (2.5) with two other contact conditions, found in the literature, which can be obtained as limit cases of (2.5). The first one is a contact condition with normal compliance and unilateral constraint, introduced in [START_REF] Jarušek | On the solvability of dynamic elastic-visco-plastic contact problems[END_REF] and studied in [START_REF] Barboteu | Analysis of quasistatic viscoplastic contact problems with normal compliance[END_REF][START_REF] Barboteu | On the behaviour of the solution of a viscoplastic contact problem[END_REF]; it can be obtained from the condition (2.5) in the case where K is a constant, and it models the contact with an elastic-rigid foundation, without hardening or softening. The second one is the classical normal compliance condition introduced in [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] and studied, for instance, in [START_REF] Klarbring | Frictional contact problems with normal compliance[END_REF][START_REF] Klarbring | On friction problems with normal compliance[END_REF][START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF]; this condition models the contact with an elastic foundation; it is characterized by the fact that the penetrations are not limited, and it can be obtained from the contact condition (2.5) in the limit g →∞.

To perform our comparison, in the upper two graphs in Fig. 5, we plot the deformed configurations as well as the contact interface forces corresponding to Problem P V at the final time t = 1.1 s, in the case β = 1 and β =-1.

In the lower two graphs, we plot the deformed configurations as well as the contact interface forces corresponding to Problem P V with α = 0, g = 0.2 m and α = 0, g = 2000 m, respectively, which correspond to the contact condition with normal compliance and unilateral constraint and the classical normal compliance contact condition, respectively. Observe that the penetrations obtained in the case with increasing stiffness coefficient (the case β = 1) are smaller than the penetrations obtained in the case where the stiffness coefficient is constant (the case of normal compliance with unilateral constraint and the case of normal compliance). The penetrations obtained in the case with decreasing stiffness coefficient (the case β =-1) are larger than that obtained in the case where the stiffness coefficient is constant (the case of normal compliance with unilateral constraint and the case of normal compliance). This behavior is, again, a consequence of the hardness and softness of the foundation. It is confirmed in Fig. 6, where the penetrations and the normal contact stresses are plotted at each time increment of the process. The magnitude of the corresponding normal stresses related these four cases of contact conditions, at each time increment of the process, is plotted in Fig. 7.

Parametric studies

We now proceed with two parametric studies which highlight the behavior of the numerical solution of the contact problem with respect to the stiffness coefficient, on one hand, and illustrate the numerical convergence of the discrete scheme, on the other hand. First, we present numerical evidence of the convergence result given in Corollary 2.2. To this end, we consider the discrete solution (u hk ρ , σ hk ρ ) of the perturbated Problem P ρ V constructed with the perturbed stiffness coefficient K ρ (ξ ) = ρc ν (1 + αξ) β , with ρ being a small positive parameter which converges to zero. We denote by (u hk * , σ hk * ) the discrete solution of the contact Problem P V with K ≡ 0. Obviously, the functions K ρ and K ≡ 0 satisfy condition (2.33); therefore, Corollary 2.2 applies. Also, note that for K ≡ 0, the contact conditions (2.5) reduce to the classical Signorini condition in a form with the gap g, which models the contact with a perfectly rigid foundation. We denote by P * V this particular form of Problem P and stress that it describes a limit situation when the layer of soft materials disintegrates by itself and does not offer resistance to penetration.

In Fig. 8, we plot the numerical estimations of the following difference:

u hk ρ -u hk * V + σ hk ρ -σ hk * Q
at the time t = 1.1 s, for various values of the parameter ρ. We change ρ from 1 to 10 -6 . We note that these numerical estimations converge to zero when ρ tends to zero. This provides numerical evidence of the convergence (2.34). This convergence is also illustrated in Fig. 9, in which we plot the deformed configurations as well as the contact interface forces for four values of the parameter ρ. We note that for ρ = 1, all the contact nodes are in the state of normal compliance, while for ρ = 0.01, almost all the contact nodes are in the state of unilateral contact with a gap g = 0.2m.

we provide a parametric study with respect to the discretization parameters h and k in order to investigate the numerical convergence order of the fully discrete scheme. To this end, we compute a sequence of numerical solutions obtained by using uniform partitions of the time interval [0, T ] and quasi-uniform triangulations of the domain . We start with the values h = 1/4 and k = 1/4 which are successively halved. The parameter h represents the average size of a finite element of the mesh, and k is the time step. The numerical results reported in Figs. 1, 2, 3, 4, 5, 6, and 7 correspond to h = 1/128 and k = 1/128. In addition, the numerical solution for h = 1/512 and k = 1/512 is taken as the "true" solution u, and we presume that this discretization corresponds to a problem with 472,074 degrees of freedom and 470,252 elements. Then, we compute the numerical errors uu hk V at t = 1.1s for several values of the discretization parameters h and k. The numerical results are presented in Fig. 10, where the dependence of the error estimates on the sum h + k is plotted. We conclude from here that our numerical method has a linear asymptotical convergence behavior with respect to the discretization parameters h and k. A rigorous proof of this statement represents an open question which, clearly, deserves to be investigated in the future.

Discussion

Studies of contact problems related to modeling, mathematical analysis, numerical analysis, and numerical simulations have attracted steady attention from researchers in engineering and mathematics. This paper provides a numerical study of a newly formulated model in contact mechanics. The model concerns frictionless contact for viscoplastic materials. The process is quasistatic, and the contact is modeled with normal compliance and unilateral constraint; the novelty of the model is due to the fact that the stiffness coefficient is allowed to depend on the history of the penetration. The problem is highly nonlinear, and it is represented by a system of nonlinear, multivalued equations. The main contribution of the paper is the introduction of an efficient numerical method to solve the problem. Numerical results on some examples are reported to show the good performance of the method.

Further work along the lines of the subject of the paper are (a) a rigorous error analysis of the numerical method, including convergence and error estimates; and (b) adaptive solution algorithms. These tasks are by no means trivial, due to the form of the boundary integral term in the variational inequality of the system. To develop adaptive solution algorithms, it is necessary to first derive a posteriori error estimate for the numerical method (cf. [START_REF] Bartels | Averaging techniques yield reliable a posteriori finite element error control for obstacle problems[END_REF][START_REF] Bostan | Recovery-based error estimation and adaptive solution of elliptic variational inequalities of the second kind[END_REF][START_REF] Bostan | A posteriori error analysis for a contact problem with friction[END_REF][START_REF] Bostan | A posteriori error estimation and adaptive solution of elliptic variational inequalities of the second kind[END_REF][START_REF] Han | A posteriori error analysis via duality theory, with applications in modeling and numerical approximations[END_REF][START_REF] Moon | A posteriori error analysis for parabolic variational inequalities, M2AN Math[END_REF][START_REF] Nochetto | Pointwise a posteriori error control for elliptic obstacle problems[END_REF] on a posteriori error analysis and adaptive solution algorithms for simpler models of elliptic or parabolic variational inequalities). Success of the work on these two aspects will not only lead to more efficient and effective numerical methods for solving the contact problem of this paper, but also be useful in numerically solving other history-dependent contact problems.

Fig. 1

 1 Fig. 1 Physical setting: compression of a half-ball against a foundation

Fig. 4 Fig. 5 Fig. 6 Fig. 7

 4567 Fig.[START_REF] Shillor | Models and variational analysis of quasistatic contact[END_REF] Evolution of the history-dependent stiffness coefficient during the time process

Fig. 8

 8 Fig. 8 Numerical estimates of at t = 1.1 s for various values of the parameter ρ

Fig. 9 Fig. 10

 910 Fig. 9 Deformed mesh and contact interface forces for Problem P ρ V at t = 1.1 s, for various parameters of ρ

Acknowledgments This work was supported by a the Marie Curie International Research Staff Exchange Scheme Fellowship within the 7th European Community Framework Programme under Grant Agreement No. 295118. The work of the second author was also supported by grants from the Simons Foundation.

 -7. Below we provide some comments on these results, with emphasis on their mechanical interpretations.

First, in Fig. 2, we show the deformed meshes and the contact interface forces with β = 1, at times t = 0.2s, t = 0.5s, t = 0.8 s and t = 1.1s. In Fig. 3, we show the deformed meshes and the contact interface forces at the same moments of time, in the case where β =-1. A comparison of these figures reveals that, at each time moment, the penetrations obtained with β = 1 are smaller than that with β =-1. This agrees with the fact that the case β = 1 corresponds to a hardening of the layer of thickness g whereas β =-1 corresponds to a softening of the layer.

Note that in Fig. 2, all the contact nodes are in the state of normal compliance, since the penetration does not reach the limit g. In contrast, in Fig. 3, a large proportion of the contact nodes is in the state of unilateral contact; there, the complete flattening of the asperities of size g = 0.2 m is reached, due to the decreasing of the stiffness coefficient. This behavior is confirmed in Fig. 4 in which the evolution of the stiffness coefficient is plotted at


Santiago Alvarez-Guerrero

Jose Ordonez-Miranda

Romeo De Coss

Juan Jose Alvarado-Gil

Effect of the metal-insulator transition on the thermoelectric properties of composites based on 𝐁𝐢 𝟎.𝟓 𝐒𝐛 𝟏.𝟓 𝐓𝐞 𝟑 with 𝐕𝐎 𝟐 nanoparticles

Keywords: Vanadium dioxide nanoparticles, Vanadium dioxide composites, Thermoelectric figure of merit

) is analyzed across the metal-insulator transition of VO ) and for temperatures up to 550 K. For temperatures higher than 350 K, it is shown that VO ) nanoparticles with a concentration of 34% enhance the electrical conductivity and ZT of the matrix by about 16% and 10%, respectively, while the Seebeck coefficient remains pretty much constant. This indicates that VO ) nanoparticles provide an effective way to enhance the thermoelectric efficiency of Bi *.+ Sb ,.+ Te -materials. The calculated ZT values for VO ) are in good agreement with the experimental data reported in the literature for temperatures higher than 350 K. The thermal conductivity values obtained for VO ) in the insulating phase are in good agreement with the experimental data reported in the literature, which are used to calculate the interface thermal resistance between Bi *.+ Sb ,.+ Te -and VO ) . Furthermore, the ratio 𝑘 '(( 𝑇𝜎 '(( ⁄ is found to be higher than the Lorenz number obtained for pure metals, such that its values increase with temperature and the VO ) concentration, for temperatures higher than the transition temperature (342.5 K) of VO ) .

Introduction

Thermoelectric (TE) materials enable direct and reversible heat-into-electricity conversion, for constituting the basis for green and environment-friendly energy technology. [1] Bi ) Te --based materials are known as narrow-gap semiconductors, with a bandgap of about 0.13 eV. [START_REF] Austin | Proc. Phys. Soc[END_REF] These compounds give rise to bipolar thermal conduction (𝑘 . ) at high temperatures which contributes to the total thermal conductivity according to the equation: 𝑘 = 𝑘 ' + 𝑘 / + 𝑘 . , (1) where 𝑘 ' and 𝑘 / are the electronic and lattice thermal conductivity, respectively. The increase in thermal conductivity leads to a poor thermoelectric efficiency. The efficiency of TE materials is determined by the dimensionless figure of merit 𝑍𝑇 = 𝑆 ) 𝜎𝑇/𝜅, where 𝑆 is the Seebeck coefficient, 𝜎 is the electrical conductivity, 𝑇 is the average temperature, and 𝜅 is the thermal conductivity. Unfortunately, these parameters are interdependent, improving one of them generally unbalances or diminishes the others, resulting in a limited energy-conversion efficiency. [START_REF] Ge | [END_REF] Considerable effort has been devoted to improving the thermoelectric properties and overcome this challenge.

There are currently several methods for improving the performance of thermoelectric materials which are related to quantum approaches known as band engineering and phonon engineering. The pillars of these theories were established by Dresselhaus and Hicks in the early 1990s. They proposed the idea of size quantization with the possibility of improving the thermoelectric properties of materials. [4,5] Dresselhaus suggested the use of low dimensionality for decreasing the phonon thermal conductivity and increasing the product 𝑆 ) 𝑛, where 𝑛 is the charge carrier concentration. Dresselhaus and Hicks [4] continued this idea and calculated ZT values for a twodimensional thin film. Their results showed that the 2D structure had an enhanced ZT coefficient compared with the 3D system of the same material. Following the same line of investigation, they also calculated ZT values for a quantum wire (1D) and confirmed that the corresponding ZT values for the one-dimensional system were enhanced with respect to the 2D quantum well. [5] In brief, they demonstrated that the nanostructures confine the electron wavefunction in one or two dimensions improving the thermoelectric efficiency. Therefore, size quantization acts on electrons and enables the so-called band structure engineering that is ascribed to the change in the density of states. The phonon engineering is related to obtain a reduced phonon mean free path in the material by means of interface scattering and hence a reduction in the lattice thermal conductivity. [5,6] For the case of a material in which 𝑘 / of Eq. 1 may be negligible, it could reach the upper limit of thermoelectric efficiency. [7] The fundamentals of the enhancement of Seebeck coefficient by means of an increase in the density of states (DOS) establishes that the valence or conduction band of the host semiconductor resonates with one energy level of a localized atom in a semiconductor matrix. [START_REF] Mahan | Proc. Natl. Acad. Sci[END_REF] Additionality, we can find alternative methods to boost the Seebeck coefficient that may include, enhancement of Seebeck coefficient by carrier energy effects, [START_REF] Faleev | [END_REF]10] inhibition of the intrinsic excitation effect, [11] or by embedding superparamagnetic nanoparticles into p-type matrices to improve the thermoelectric performance, [12] among other approaches.

There are interesting reports on the thermoelectric properties of different composite materials based on the Bi *.+ Sb ,.+ Te -matrix, [10][11][12] for instance, Kim et al. [10] reported an enhancement of Seebeck coefficient in Bi *:+ Sb ,:+ Te -with high-density tellurium nanoinclusions. They reported experimental evidence that the addition of Te nanoparticles enhances the electron scattering by the carrier energy filtering effect due to the barrier potential formed by the Te/Bi *:+ Sb ,:+ Te -interfaces. It was found that the addition of 15% in volume of tellurium into the Bi *:+ Sb ,:+ Te -films increases the Seebeck coefficient significantly by more than 45% at 320K. However, the electrical conductivity decreases with increasing density of Te nanoparticles throughout the temperature range of 320 to 520 K.

More recently, Zhu et al. [11] investigated the addition of Cu impurity into a BST matrix. According to these authors, Cu impurities induce an enhancement in the thermoelectric behavior with an ZT average value of 0.85 across the entire temperature range due to synergistic optimization of the electronic and thermal transport properties. This ZT value was larger than that of the BST matrix of 0.52.

The effect of introducing magnetic nanoparticles in BST has also been explored. Li et al. [12] demonstrated that introducing superparamagnetic nanoparticles is an effective approach to simultaneously enhance the thermoelectric and cooling performance of p-type BiSbTe-based alloys. The authors incorporated Fe -O 1 nanoparticles into a commercial Bi *.+ Sb ,.+ Te -matrix and found that this induces the enhancement of the Seebeck coefficient due to the carrier multiple scattering. It also induces a significant decrease in thermal conductivity due to the phonon scattering caused by magnetic moment fluctuations. The composite reached a maximum ZT corresponding to 32% higher than that of the matrix.

It is well known that thermoelectric properties of strongly correlated electron systems could contribute to the semiconductor to metal transition (SMT). [13] Vanadium dioxide (VO ) ) is a strongly correlated oxide, [14] which shows a SMT close to room temperature, which induces remarkable changes in its structural, optical, and electrical properties. These characteristics have made vanadium dioxide the basis in the development of a great diversity of technological applications, such as temperature sensors, [15] energy storage media, [START_REF] Paone | Opt. Model. Meas. Sol. Energy Syst. III[END_REF] thermally-driven radiative diodes, [START_REF] Forero-Sandoval | [END_REF] among others. However, VO ) has received little attention in the field of thermoelectricity. [14] Few applications for VO ) in thermoelectricity has recently emerged. [18,19] Kosta et al. [18] studied VO ) nanoparticles embedded into thermoelectric n-type Mg ) Si *.222 Sn *., Sb *.*,) alloy. They found that the addition of VO ) nano inclusions provide an extra scattering for the low frequency phonons, in addition to that caused by the point defects and dislocations that appear during the ball milling process, which can scatter high/medium-frequency phonons. As consequence, the effective thermal conductivity diminishes and the ZT increases by 38%.

Back et al. [19] investigated the thermoelectric properties of MgO/VO ) BST composites by extrinsic phase mixing of MgO and VO ) nanoparticles in the BST matrix. It is important to note that BST is a narrow gap semiconductor and presents the bipolar diffusion effect, which represents a critical point to increase the performance of the semiconductor. As is known, when the energy gap of a material is small, some electrons can be thermally excited to the next higher energy band leaving the corresponding empty states (holes). [START_REF] Chen | Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons[END_REF] Consequently, the total thermal conductivity is enhanced according to Eq. 1. Thus, the authors found that MgO and VO ) distribution in the BST effectively scatters phonons, inducing a significant reduction of thermal conductivity by decreasing lattice and bipolar thermal conductivity. These phenomena allowed the enhancement of the ZT values for the MgO/VO ) BST composites over a wide range of temperatures.

Inspired by the above-mentioned works [18,19] and the fact that vanadium dioxide shows a fully reversible first-order metal-to-insulator transition accompanied by a relatively little change in its thermal conductivity and an abrupt jump by nearly two orders of magnitude in its electrical conductivity on heating, [START_REF] Morin | [END_REF] we have investigated the influence of VO ) on the electrical conductivity, thermal conductivity and Seebeck effect for VO ) /BST composites. Recently, Jung et al. [22] obtained an analytical expression for the effective thermoelectric properties and dimensionless figure of merit of a composite considering interfacial resistance. The model was developed via the mean-field homogenization scheme called the Mori-Tanaka method. The results were validated against the Finite Element Methods (FEM) with good agreement for volume fraction of inhomogeneities up to 20%.

Recently, Qiu et al. [23] experimentally studied the thermal conductivity from mesoporous to macroporous SiOC ceramics. Their thermal conductivities were measured in freestanding samples using the 3ω technique, subsequently they employed three typical models (series, Maxwell-Eucken 1, and effective medium theory) assigning adequate weights to each model and combining them to derive an empirical formula for the thermal conductivity. In this work, we combine the Bruggeman scheme, [24] (considered as one of the most accurate for high filler volume fractions) with the accuracy of the finite element methods [25] to evaluate the thermal conductivity, electrical conductivity, and Seebeck coefficient for composites of three different vanadium dioxide concentrations.

Numerical methodology

The methodology for modeling the thermoelectric properties of particulate VO ) /BST composites by means of numerical simulations is presented in this section. This approach is based on the finite element method of a 3D model that consists of a set of spherical particles randomly distributed inside a cubic matrix. The spatial distribution of the inclusions is generated using a MATLAB TM code that ensures the non-overlapping of the inclusions, as detailed in our previous work. [26]

Effective electrical conductivity

A similar procedure to the one described by Nilsson et al. [27] was used to calculate the effective electrical conductivity. The cubic domain (matrix) with VO ) inclusions is illustrated in Fig. 1(a). The bottom side of the domain was electrically grounded at 0 V (𝑉 Г = 0 V), a constant voltage was applied to the upper side (𝑉 Г = 1 V), whereas periodic boundary conditions were established on the remaining four sides (𝑉 Г4 = 𝑉 Г. ). For the particles, the relationship describing the electrical conductivity of VO ) [24] as a function of temperature is as follows:

𝜎 56 ! (𝑇) = (3𝑓 7 -1)𝜎 7 + (3𝑓 8 -1)𝜎 8 + J((3𝑓 7 -1)𝜎 7 + (3𝑓 8 -1)𝜎 8 ) ) + 8𝜎 7 𝜎 8 4 , (2)
where 𝜎 7 = 84175.1 Sm -1 and 𝜎 8 = 4.86 x 10 6 𝑒 9-,-:.,/< Sm -1 are the electrical conductivity in its dielectric (low temperature) and metallic (high temperature) phases, respectively; 𝑇 = = 3626.2 K is the transition temperature. 𝑓 7 and 𝑓 8 = 1 -𝑓 7 are the volume fractions of the metallic and insulating domains during the SMT, both are temperature dependent functions according to the equations:

𝑓 7 = - 𝑇 𝑈 𝑊 * Q- 𝑇 2𝑈 exp T- 𝑇 𝑈 U erfc T 𝑇 = -𝑇 √2∆𝑇 U[, (3)
𝑓 8 = 𝑇 𝑈 𝑊 * Q- 𝑇 2𝑈 exp T- 𝑇 𝑈 U erfc T 𝑇 -𝑇 = √2∆𝑇 U[, (4)
where 𝑈 = 323.7 K is the activation energy, 𝑊 * is the main branch of the Lambert 𝑊 function, and ∆𝑇 = 5.5 K is the standard deviation for the heating process. The main branch of the Lambert function was implemented in MATLAB TM code and called as an external function from COMSOL TM software.

For the Bi *.+ Sb ,.+ Te -material, the electrical conductivity values were taken from the experimental work reported by Li et al. [12] and implemented in COMSOL TM as the interpolation function "sigma(T)". The effective electrical conductivity was calculated by solving the Poisson equation with respect to the electric potential V, using a DC stationary method, [27] ∇. (𝜎∇V) = 0.

(

) 5
The magnitude of the effective electrical conductivity was obtained by integrating the current, 𝐼 = 𝜎∇𝑉 over the whole cubic computational domain and dividing the result by the area, 𝐴 of the grounded face and the applied voltage as follows, Temperature dependence of the electrical conductivity for the VO ) (blue line) and for Bi *.+ Sb ,.+ Te -(black line).

𝜎 '(( = 𝐼 𝐴𝑉 . ( 6

Effective thermal conductivity

The analytical function that describes the thermal conductivity of VO ) , 𝑘 56 ! (𝑇) can be reasonably well described as follows, [28] 𝑘 56 ! (𝑇) = 𝑘 > + 𝑘 7 -𝑘 > 1 + 𝑒 9?(<9< " ) , (7) where 𝑘 > = 3.6 Wm -1 K -1 and 𝑘 7 = 6.0 Wm -1 K -1 are the thermal conductivity in its dielectric (low temperature) and metallic (high temperature) phases, respectively; 𝑇 B = 342.5 K is the transition temperature and 𝛽 = 3.6 K -1 is the phase-transition slope of 𝑘 56 ! (𝑇) at 𝑇 = 𝑇 B . For the Bi *.+ Sb ,.+ Te -material, the thermal conductivity values were taken from the experimental data reported by Back et al. [19] and implemented in COMSOL TM as the interpolation function "kappa(T)". The effective thermal conductivity was calculated, solving the Fourier equation by means of the stationary method, as is described in our previous work. [26] In the 3D model (see Fig. 2(a)), the top surface temperature (red zone) is kept higher than the bottom one (white zone) to induce a heat flux from top to bottom. All the remaining four surfaces are thermally insulated. The average value of conductive heat flux, normal to all three xz surfaces, is used to calculate the 𝑘 '(( by means of the Fourier law of heat conduction,

∇. e𝑘 '(( 𝑇f = 0. ( 8
)
The heat flux 𝑞 in the cubic domain was evaluated as the surface integration of the normal total heat flux in the top-to-bottom direction as follows,

𝑞 = h 𝑞 C (𝑥, 𝑦, 𝑧)𝑑𝑦𝑑𝑧 , (9)
and thus, the effective thermal conductivity, 𝑘 '(( , is calculated as follows,

𝑘 '(( = 𝑞 𝑙(𝑇 DE= -𝑇 BE/> ) , ( 10
)
where 𝑙 is the length of the cubic domain.

Calculation of the interface thermal resistance between 𝐕𝐎 𝟐 and 𝐁𝐢 𝟎.𝟓 𝐒𝐛 𝟏.𝟓 𝐓𝐞 𝟑

Interfacial thermal resistance is among the most relevant parameters determining the thermal conductivity of composites, and it is even present at atomically perfect interfaces. This resistance appears because of the physical differences, between neighboring materials, generating heat carriers scattering. These effects have been extensively explored in the literature and it has been shown that the numerical values appear in a broad range. This property has been specially studied in carbon nanostructures in which it has explored the determinant role of low frequency phonons. [29,30] For our system made up with two different materials, we theoretically estimated the interfacial thermal resistance magnitude between VO ) and Bi *.+ Sb ,.+ Te -based on previous experiments. As we mentioned before, Back et al. [19] investigated the thermoelectric properties of VO ) / Bi *.+ Sb ,.+ Te -composites. They found that the dispersion of VO ) particles effectively scatter the heat-carrying phonons resulting in a reduced thermal conductivity. The idea that the interfacial phonon scattering enhances the ZT figure of merit in thermoelectric materials was pointed out by Hicks and Dresselhaus. [5] According to their calculations, a significant increase in ZT can be achieved by going to lower dimensions, [5,6] particularly, a one-dimensional Bi ) Te -nanowire will undergo an increase in phonon scattering from the surfaces and this leads to a reduction in the lattice thermal conductivity.

The lattice thermal conductivity, 𝑘 / is related to the mean free path, 𝜆 according to the equation:

𝑘 / =
, -𝐶 F 𝑣𝜆, whrere 𝐶 F is the lattice heat capacity and 𝑣 is the sound velocity. We know from Ref. 19 that the mean free path is ~ 70 nm. Thus, if the diameter 𝑑 of the particles is smaller than 𝜆, the scattering at the surface will occur and cause a decrease in 𝑘 / . [5] Therefore, the total thermal conductivity is reduced with respect to the composite material. Based on this fact and assuming that the quantization effects could be observed for nanostructures < 50 nm, we have considered a diameter for the spherical particles of 𝑑 = 6 nm.

It is well known that the reduction of the effective thermal conductivity caused by the scattering of charge carriers is synonymous with the existence of interfacial thermal resistance. Phonon scattering mainly slants the heat transfer due to a barrier that is formed at the boundary interface yielding thermal boundary resistance (TBR). In addition, the less contact between fillers yields thermal contact resistance (TCR), which results in lower thermal conductivity. Therefore, the interfacial thermal resistance is a contribution of both TBR and TCR. [31] Different studies have experimentally reported the total thermal conductivity values for different vanadium dioxide samples. For thin films of VO ) , Kizuka et al. [32] reported that the total thermal conductivity increases from 3.6 to 5.4 Wm -1 K -1 . Oh et al. [33] reported an increase from 3.5 to 6 Wm -1 K -1 in thin films of VO ) . Meanwhile Chen et al. [34] reported a decrease from 3.5 to 2.2 Wm -1 K -1 for polycrystalline VO ) . The authors explain that this behavior of the thermal conductivity with temperature is caused by the reduction of the lattice thermal conductivity across the phase transition. Therefore, the differences in the variation trends of the thermal conductivity, reported in previous works, might be caused by discrepancies induced by the carrier contribution, which are sensitive to the structure and defects of the samples. Therefore, according to abovementioned experimental data, [32][33][34] it can be inferred that whatever the sample, bulk or thin films, the thermal conductivity of VO ) is ~3.6 Wm -1 K -1 in the insulating phase. In the case under consideration in this work, we followed the results reported by Back et al. [19] These authors synthetized a composite mixing 15 g of Bi *.+ Sb ,.+ Te -(𝑚 GHI ) powder with 5 mol.% of VO ) (𝑛 56 ! ) nanopowders. The volumetric composition, 𝑓 for this composite was obtained using the equation:

𝑓 = 100𝑛 VO 2 𝑀 VO 2 𝜌 BST 𝑛 VO 2 𝑀 VO 2 𝜌 BST + 𝑚 BST 𝜌 VO 2 . ( 11
)
From the density values for BST (𝜌 BST ) and for VO ) (𝜌 56 ! ) and its corresponding molar mass (𝑀 56 ! ) values, we obtain the value 𝑓 = 29.4% of VO ) . Then, the corresponding interfacial thermal resistance is calculated by fitting the experimental data and theoretically predicted effective thermal conductivity for an ITR value that matches in a temperature range in which VO ) undergoes the insulating phase.

To carry out the simulations we considered a thin interface layer between both adjacent different materials (see Fig. 1(a)) where heat cannot flow directly through the layer which results in a temperature jump across the interface. [35] The interfacial thermal resistance, R is defined as

𝑇 EJ= -𝑇 8K = -𝑅𝑞 L . 𝑛, ( 12
)
where T out and T in refer to temperatures at the outer and inner surfaces of the interface, respectively, qs is the heat flux at the interface, and the n is the outward surface normal vector. The results of this approach are shown in section 3.1.

Seebeck coefficient calculation

To know the overall effect of VO ) on the Bi *.+ Sb ,.+ Te -thermoelectric matrix, we have also determined the Seebeck coefficient for the VO ) /Bi *.+ Sb ,.+ Te -composites. Seebeck coefficient (𝑆) is a property of a material depending on its electronic structure near the Fermi level. [36] The sign of 𝑆 indicates which charge carriers dominate the electronic transport in semiconductors as well as in metals. Negative values of 𝑆 indicate that the material is an n-type semiconductor.

We simulated the Seebeck effect by means of a thermoelectric generator system implemented in COMSOL TM . This system allows the Seebeck effect, a phenomenon where the difference in temperature of a material leads to a potential difference. We employed a brick-shaped domain with six sides as shown in Fig. 3(a). The thermoelectric part is made up of Bi *.+ Sb ,.+ Te -and it is capped by two thin copper electrodes while the VO ) spherical particles are embedded inside the matrix.

To estimate the Seebeck coefficient magnitude, we employed the floating potential boundary condition to calculate the potential difference that undergoes the system. Then, we take the ratio of the potential difference or floating potential (that is the thermoelectric voltage seen at the Cu terminals) to the temperature difference as follows,

𝑆 '(( = - ∆𝑉 ∆𝑇 . ( 13
)
A considerably large temperature difference was set on the surface of the electrodes, for example, 100 °C at the left of the boundary of the Cu electrode and 0 °C at the right of the Cu electrode. The relative permittivity values for VO ) were taken from the experimental results reported by Yang et al. [37] and implemented in COMSOL TM as the piecewise function "𝑒𝑝𝑠𝑖𝑙𝑜𝑛(𝑇)". Although the dielectric constant of the Bi *.+ Sb ,.+ Te -matrix is unknown, it is likely lower than that of Sb ) Te - and probably near to 50. [38,39] The Seebeck coefficient values for VO ) were taken from the work reported by Khan et al. [14] and implemented in COMSOL TM as the piecewise function "𝑆𝑒𝑒𝑏𝑒𝑐𝑘_VO ) ". Meanwhile the corresponding Seebeck coefficient values for Bi *.+ Sb ,.+ Te - material were taken from the experimental data reported by Back et al. [19] and implemented in COMSOL TM in the same way with the function "𝑆𝑒𝑒𝑏𝑒𝑐𝑘_BST". Fig. 3(b) shows the Seebeck coefficient for both VO ) and BST as a function of the temperature and in general, the properties employed to the simulations are detailed in Table 1.

Results and discussions

In this section we present the results of the simulations performed for the VO ) /Bi *.+ Sb ,.+ Te - composites. It is important to mention that VO ) are not participating as dopant (which means that vanadium dioxide does not change the crystallinity of the Bi *.+ Sb ,.+ Te -matrix) and our simulations contemplate that VO ) undergoes the "ideal efficiency" in its thermal conductivity (3.6 -6.0 Wm -1 K -1 ) and electrical conductivity (131 -8.4 x 10 4 Sm -1 ). To ensure the physical property calculations accuracy, in each case 3 repetitions were carried out, and the average results are reported as follows.

Electrical conductivity, thermal conductivity and Seebeck coefficient

First, we discuss the results for electrical conductivity. The temperature-dependent electrical conductivity for VO ) nanoparticles embedded in Bi *.+ Sb ,.+ Te -as a function of temperature is shown in Fig. 4. At low temperatures (< 350 K), the effective electrical conductivity tends to decrease, since VO ) is still in the insulating phase and the conductivity of the BST decreases with increasing temperature (Fig. 1(b)). At about 327 K the curves reach a minimum, then the electrical conductivity starts to increase continuously with increasing temperature, ought to the changing of VO ) into the metallic phase until reaching a maximum at ~348 K. It is interesting to note that this maximum value corresponds to the intersection of the conductivity curves of both VO ) and BST as seen in Fig. 1(b). Thereafter, the curve decreases, but the EEC values are significantly higher than those of the BST matrix.

With the increasing of VO ) concentration, the electrical conductivity enhances continuously. As the temperature increases, the electrical conductivity of Bi *.+ Sb ,.+ Te -tends to ~500 Scm -1 meanwhile the effective electrical conductivity for the VO ) /Bi *.+ Sb ,.+ Te -composite with 𝑓 = 34.0% of VO ) reached 625 Scm -1 at 550 K, which was around 16% greater than that of the BST matrix (525 Scm -1 ) at the same temperature. This is relevant considering the range of temperature for thermoelectric applications (T < 1000 °C). The improvement of the electrical conductivity could contribute to the direct enhancement of the thermoelectric efficiency in Bi *.+ Sb ,.+ Te -bulk materials if the thermal conductivity is reduced or kept nearly constant.

Turning our attention to the effective thermal conductivity results. The ETC values as a function of temperature are plotted in Fig. 5. We tested various interface thermal resistance values until we obtained the blue curve depicted in Fig. 5(a), the inset in the figure shows a good agreement with the experimental results reported by Back et al. [19] within the temperature range for T < 350 K. Then, based on the methodology described in section 2.2.1, we estimated the interface thermal resistance between Bi *.+ Sb ,.+ Te -and VO ) as R = 1.91 m 2 KGW -1 . This approximation is reasonable because the VO ) nanoparticles yield large interfacial area and the interfaces among different phases have been shown to play an important role as phonon scattering centers. [39,40] The calculated ETC values increase as the temperature increases, and at ~ 343 K a change in the slope of the curve is observed. In Fig. 5, it can be clearly observed the effect on the thermal conductivity of the VO ) fraction undergoing the phase change from insulator to metal on heating. In other words, our results (blue curve in Fig. 5(a)) show the behavior for the VO ) /Bi *.+ Sb ,.+ Te - composite if the VO ) experiences what we have called the "ideal efficiency" on its thermal conductivity (3.6 -6.0 Wm -1 K -1 ). Although the magnitude of the thermal conductivity of VO ) in the metallic phase is larger than the one for the matrix, the total ETC values for all concentrations of VO ) are lower than that of the pure BST (see Fig. 5(b)), which is favorable for thermoelectric purposes. In the following, the results for the thermoelectric behavior of the VO ) /Bi *.+ Sb ,.+ Te -are presented. We determined the Seebeck coefficient values, in which the simulated thermoelectric generator was tested at various temperature differences (ΔT), while the cold side was maintained at 0 °C. Fig. 6(a) portrays the variation of Seebeck coefficient (S) against temperature (T). Positive Seebeck coefficient values indicate the p-type conduction behavior. For low temperatures, the Seebeck coefficient values for the composite with 𝑓 = 34.0% of VO ) is slightly above the corresponding values of the matrix. This is due to the fact that in the semiconductor phase, VO ) shows large values of S (see Fig. 3(b)) as a consequence of the significant enhance in the average energy of carriers, that takes part in transport phenomenon at low temperatures. [14] When the temperature increases, the total Seebeck coefficient of the composites slightly decreases due to the dominant contribution of the matrix. In other words, the addition of VO ) in the BST matrix does not improve the values of Seebeck coefficient. This can be explained as another consequence of the interface thermal resistance as follows: The Seebeck coefficient refers to the generated voltage difference for a given temperature difference through an electrical path. If a jump of the temperature occurs at the interface, the effective total temperature difference is reduced across the hot to cold side while the electrical current can flow through the inclusion. As a result, the thermoelectric composite produces less voltage as the interface thermal resistance increases. [22] Notice that for T < 350 K, the Seebeck coefficient values for composites are slightly larger than those of the matrix. For T > 350 K, S slightly decreases with respect to the matrix. From this analysis we conclude that this thermoelectric property remains almost unchanged. Fig. 6(b) shows the dimensionless thermoelectric figure of merit (ZT) for pure BST and for the composites as function of temperature. It can be observed that ZT values of all composites increase with increasing temperature to reach a maximum value and then to decrease with further increase in temperature. Moreover, the ZT values significantly increased with the incorporation of VO ) nanoparticles and it is notable that there is a good agreement between the experimentally determined ZT values and our calculations of VO ) /Bi *.+ Sb ,.+ Te -composites for 𝑓 = 29.4% VO ) . The low ZT values of composites at low temperatures is a consequence of changes that VO ) undergoes for both thermal and electrical properties. Nevertheless, the composite with 𝑓 = 34.0% exhibited a maximum ZT value of ~0.7 at 375 K, which is around 10% higher than that of BST matrix. On the other hand, Li-based ternary Heusler compounds reach its maximum TE efficiency of ~0.8 at about 200 K [41] which indicate that they can be adequate for cooling applications. In our case, considering the VO ) /Bi *.+ Sb ,.+ Te -performance makes it suitable for room temperature or intermediate (350 -450 K) applications.

From Fig. 6(b), it is also observed that for higher temperatures, the figure of merit values showed a significant improvement of around 20% compared with pure Bi *.+ Sb ,.+ Te -. Thus, this improvement is primarily due to the enhancement of the electrical conductivity and reduction in thermal conductivity. The trends of the ZT curves for the composites are similar to the Heusler quaternary compound PdZrTiAl, the metallic nature of this alloy in spin down leads to an acceptable extent amount of ZT in the entire temperature range, so it has the ability to use in power generators. [42] In our case, the VO ) /Bi *.+ Sb ,.+ Te -composite material with a ZT range 0.8-0.3 may be suitable for waste heat utilization. [START_REF] Pöttgen | Rare Earth Chemistry[END_REF] The improvement in the ZT values of BST matrix is due to the reduction in the thermal conductivity and the remarkable increase in the electrical conductivity compensating for the slight decrease of Seebeck coefficient caused by the addition of VO ) nanoparticles. Experimentally, it has been also reported that the embedding of nanoparticles in BST matrix significantly enhances the thermoelectric performance. [START_REF] Kim | [END_REF]

Relationship between thermal conductivity, electrical conductivity, and temperature

The relationship between the thermal conductivity and the electrical conductivity is established by the Wiedemann-Franz (WF) law. This law was originally developed for metals based upon the fact that heat and electrical transport both involve the free electrons. The law states that the ratio of the electronic contribution of the thermal conductivity to the electrical conductivity of a metal is proportional to the temperature according to the following relationship,

𝜅 𝜎 = 𝐿𝑇, ( 14
)
where 𝐿 is the constant of proportionality known as Lorenz number. [START_REF] Kittel | Introduction to Solid State Physics[END_REF] For metals, this constant is denoted as 𝐿 * and its value is equal to 2.45 × 10 -8 WS -1 K -2 . For pure materials, the values of 𝐿 are also well known. However, for alloys or composite materials this value is unknown.

Here, it is noteworthy that due to the correlated nature of VO ) it is not allowed to calculate the expected electronic thermal conductivity, using the Wiedemann-Franz approach for the entire range of temperatures. [START_REF] Lee | [END_REF]47] The relationship between the conductivities and temperature is quantified in Fig. 7. It is seen how the values of 𝜅/𝑇𝜎 for composites are lower than those of the matrix. In addition, for all concentration the corresponding values are slightly higher than 𝐿 * , which can be helpful in evaluating the role of electrons as energy carriers in these systems.

Conclusions

The thermoelectric properties of VO ) /Bi *.+ Sb ,.+ Te -composites with volume fractions of 0, 24.7%, 29.4% and 34.0% have been investigated across the metal-insulator transition of VO ) and for temperatures up to 550 K, by means of simulations based on Finite Element Methods. We have shown that VO ) nanoparticles with a concentration of 34.0% enhance the electrical conductivity by approximately 16% for temperatures higher than 350 K, while the Seebeck coefficient remains pretty much constant. The resulting maximum value of ZT = 0.7 shows up at 375 K, which represents an increase of 10% in comparison with the corresponding one of the matrix. The calculated ZT values agree well with experimental data and indicate that VO ) nanoparticles provide an effective way to enhance the thermoelectric efficiency of Bi *.+ Sb ,.+ Te -materials. We have estimated the interface thermal resistance magnitude between Bi *.+ Sb ,.+ Te -and VO ) as R = 1.91 m 2 KGW -1 . Furthermore, it has been found that the ratio 𝑘 '(( 𝑇𝜎 '(( ⁄ is higher than the Lorenz number obtained for pure metals, such that its values increase with temperature and the VO ) concentration, for temperatures higher than the transition temperature (342.5 K) of VO ) . Finally, we want to emphasize that the numerical approach presented here can be useful to predict the effective thermoelectric properties of composites based on VO ) .

)Figure 1 .

 1 Figure 1. (a) 3D domain geometry to simulate the effective electrical conductivity. The colors indicate an electrical potential ranging from V = 1 V (blue zone) to V = 0 V (green zone). (b)

Figure 2 .

 2 Figure 2. (a) Schematic representation of the 3D composite model. (b) Temperature dependence of the thermal conductivity of VO ) (blue line) and for Bi *.+ Sb ,.+ Te -(black line).

Figure 3 .

 3 Figure 3. (a) Schematic representation of the modeled composite made up of a Bi *.+ Sb ,.+ Te - matrix hosting in which particles of VO ) nanoparticles. (b) Seebeck coefficient as a function of temperature of BST (black line) and VO ) (blue line). The blue and red colors in (a) stand for the electrical potential applied to the composite material, blue (1 V) and red (ground potential).

Figure 4 .

 4 Figure 4. Temperature dependence of electrical conductivity for VO ) /Bi *.+ Sb ,.+ Te -composites for different VO ) concentration.

Figure 5 .

 5 Figure 5. (a) The comparison between the theoretical effective thermal conductivity (blue line) predicted by the simulations and the experimental [19] results (pink line). (b) The effective thermal conductivity versus temperature calculated for various concentrations.

Figure 6 .

 6 Figure 6. Thermoelectric properties of VO ) /Bi *.+ Sb ,.+ Te -composite material as a function of temperature (a) Seebeck coefficient, (b) Figure of merit.

Figure 7 .

 7 Figure 7. Temperature dependence of 𝜅/𝑇𝜎 relation for VO ) /Bi *.+ Sb ,.+ Te -composite as temperature function.

Table 1 .

 1 Electrical, thermal, and thermoelectric properties for the components of the VO ) /Bi *.+ Sb ,.+ Te -composite. *We have previously developed a COMSOL TM -readable heat capacity function for VO ) .[26]

	Property	Cu electrodes	Bi *.+ Sb ,.+ Te -	VO )
	Relative permittivity	1	50	𝑒𝑝𝑠𝑖𝑙𝑜𝑛(𝑇)
	Electrical conductivity	5.998 x 10 7 [S/m]	sigma(T)	Eq. 2
	Heat capacity at constant	385[J/(kg*K)]	185.9895[J/(kg*K)]	*
	pressure
	Density	8960[kg/(m*K)]	6840[kg/(m*K)]	4381.161-
				0.1409686*T
	Thermal conductivity	400[W/(m*K)]	kappa(T)	Eq. 7
	Seebeck coefficient	6.5 x 10 -6 [V/K]	Seebeck_BST(T)	Seebeck_VO2(T)

Availability of data:

The data that support the findings of this study are available from the corresponding author upon reasonable request.