Bruno Zanuttini

Stanislav Živný
email: stanislav.zivny@comlab.ox.ac.uk

A note on some collapse results of valued constraints

Keywords: Valued constraint satisfaction problems, Expressive power, Max-closed constraints, Theory of computation

Valued constraint satisfaction problem (VCSP) is an optimisation framework originally coming from Artificial Intelligence and generalising the classical constraint satisfaction problem (CSP). The VCSP is powerful enough to describe many important classes of problems. In order to investigate the complexity and expressive power of valued constraints, a number of algebraic tools have been developed in the literature. In this note we present alternative proofs of some known results without using the algebraic approach, but by representing valued constraints explicitly by combinations of other valued constraints.

Introduction

The Valued Constraint Satisfaction Problem (VCSP) [START_REF] Schiex | Valued constraint satisfaction problems: hard and easy problems[END_REF] is a general framework to model various optimisation problems [START_REF] Bistarelli | Semiring-based CSPs and valued CSPs: Frameworks, properties, and comparison[END_REF][START_REF]The Handbook of Constraint Programming[END_REF], which generalises the classical constraint satisfaction problem (CSP) [START_REF] Montanari | Networks of constraints: Fundamental properties and applications to picture processing[END_REF][START_REF] Creignou | Complexity Classification of Boolean Constraint Satisfaction Problems[END_REF][START_REF] Dechter | Constraint Processing[END_REF][START_REF]The Handbook of Constraint Programming[END_REF].

Informally, in the VCSP framework, an instance consists of a set of variables, a set of possible values, and a set of (soft) constraints. Each constraint has an associated cost function which assigns a cost (or a degree of violation) to every possible tuple of values for the variables in the scope of the constraint. The goal is to find an assignment of values to all of the variables which has the minimum total cost. The set of cost functions used in the description of the problem is called the valued constraint language. Infinite costs can be used to indicate infeasible solutions, and if the range of all cost functions is {0, ∞} (so-called hard constraints), we get the standard CSP framework as a special case.

An important problem is determining which additional constraints can be expressed by a given valued constraint language. The notion of expressibility has been a key com-ponent in the analysis of complexity for the classical CSP model [START_REF] Jeavons | On the algebraic structure of combinatorial problems[END_REF][START_REF] Bulatov | Classifying the complexity of constraints using finite algebras[END_REF]. It was also a major tool in the complexity analysis of a wide variety of Boolean constraint problems carried out by Creignou et al. [START_REF] Creignou | Complexity Classification of Boolean Constraint Satisfaction Problems[END_REF], where it was referred to as implementation. Expressibility is a particular form of problem reduction: if a constraint can be polynomially expressed in a given constraint language, then it can be added to the language without changing the computational complexity of the associated class of problems. Hence determining what can be expressed in a given valued constraint language is a fundamental step in the complexity analysis of valued constraint problems.

A number of algebraic tools have been developed to understand the complexity and expressive power of constraint languages [START_REF] Jeavons | On the algebraic structure of combinatorial problems[END_REF][START_REF] Cohen | The complexity of soft constraint satisfaction[END_REF][START_REF] Cohen | An algebraic characterisation of complexity for valued constraints[END_REF]. Using the algebraic approach, Cohen et al. showed [START_REF] Cohen | The expressive power of valued constraints: Hierarchies and collapses[END_REF] that some valued constraint languages, including the set of all monotonic valued constraints, can be expressed using only valued constraints of a fixed finite arity. They also showed some other classes of valued constraints, including the set of all monotonic valued constraints with finite cost values, which cannot be expressed by a subset of any fixed finite arity, and hence form an infinite hierarchy [START_REF] Cohen | The expressive power of valued constraints: Hierarchies and collapses[END_REF].

In this paper we consider the class of monotonic, so-called max-closed valued constraints. Our main contributions are new proofs of certain recent results on the expressive power of max-closed valued constraints [START_REF] Cohen | The expressive power of valued constraints: Hierarchies and collapses[END_REF]. The proofs presented here use a very different approach: they do not rely on the algebraic approach, but are based on an explicit representation of valued constraints. Although some of the proofs presented here are simple and considered folklore, we also present new and more involved proofs, thus providing a logical approach to the expressive power of VCSP.

VCSP instances with max-closed valued constraints are well-known to be solvable in polynomial time (generalising the tractable language for the standard Satisfiability problem consisting of Horn clauses) [START_REF] Jeavons | Tractable constraints on ordered domains[END_REF][START_REF] Cohen | The complexity of soft constraint satisfaction[END_REF], but much less is known about their expressive power. The new proofs presented here have a number of interesting new features:

First, the alternative proofs are simpler, not requiring the machinery of the algebraic approach. Moreover, they are constructive and thus give new techniques for studying valued constraint satisfaction problems from the logical point of view, as has proven useful in the propositional case [START_REF] Creignou | Complexity Classification of Boolean Constraint Satisfaction Problems[END_REF].

Second, recent progress on the expressive power of submodular valued constraint has shown that both the algebraic approach and the non-algebraic approach, presented in this paper, can play important roles in investigating the expressive power of valued constraints. In particular, a combination of these two approaches has proved particularly effective.

Submodular constraints are a key concept in operational research and combinatorial optimisation, see for example [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF][START_REF] Topkis | Supermodularity and Complementarity[END_REF][START_REF] Schrijver | Combinatorial Optimization: Polyhedra and Efficiency[END_REF]. As submodular valued constraints can be characterised as min-max-closed valued constraints [START_REF] Cohen | The complexity of soft constraint satisfaction[END_REF], understanding the expressive power of max-closed valued constraints is a natural first step towards complete understanding of the expressive power of submodular valued constraints.

Moreover, the techniques presented in this paper have been recently used successfully. First, using explicit representations of submodular valued constraints, Živný and Jeavons showed [START_REF] Živný | Classes of Submodular Constraints Expressible by Graph Cuts[END_REF] a new class of submodular valued constraints of arbitrary arities expressible by binary submodular valued constraints, thus describing a new class of VCSP instances which can be solved efficiently using Min-Cut. This class was independently discovered in [START_REF] Zalesky | Efficient determination of Gibbs estimators with submodular energy functions[END_REF]. Next, using a combination of the explicit representation of submodular valued constraints and certain algebraic results, Živný et al. showed [START_REF] Živný | The Expressive Power of Binary Submodular Functions[END_REF] a novel class of submodular valued constraints which can be solved efficiently using Min-Cut. Furthermore, they also showed that there are submodular valued constraints which cannot be expressed by binary submodular valued constraints [START_REF] Živný | The Expressive Power of Binary Submodular Functions[END_REF], thus answering an open problem which has been considered within several different contexts in computer science, including computer vision, artificial intelligence, and pseudo-Boolean optimisation. This provides an example how important it is to understand various aspects of the expressive power of valued constraints. We believe that understanding both the algebraic and non-algebraic aspects of the expressive power of valued constraints can contribute to the development of better algorithms for discrete optimisation problems.

Background

In the original definition of the valued constraint satisfaction problem (VCSP) given in [START_REF] Schiex | Valued constraint satisfaction problems: hard and easy problems[END_REF], costs were allowed to lie in any positive totally ordered monoid called a valuation structure. For our purposes it is sufficient to consider costs which lie in the set N consisting of all natural numbers (including zero) together with infinity3 .

Given some fixed set D, a function from D k to N will be called a cost function. If the range of φ lies within {0, ∞}, then φ is called a crisp cost function. Note that crisp cost functions correspond precisely to relations. If the range of φ lies entirely within N, then φ is called a finite-valued cost function. Finally, if φ takes both finite and infinite values, then φ is called a general cost function. Definition 1 An instance I of VCSP is a triple V, D, C , where V is a finite set of variables, which are to be assigned values from the set D, and C is a set of valued constraints. Each c ∈ C is a pair c = v, φ , where v is a tuple of variables of length m, called the scope of c, and φ : D m → N is a cost function. An assignment for the instance I is a mapping s from V to D. We extend s to a mapping from V k to D k on tuples of variables by applying s component-wise. We denote by A the set of all assignments. The cost of an assignment s is defined as follows:

Cost I (s) = v,φ ∈C φ(s(v)).
A solution to I is an assignment with minimum cost.

The VCSP is a very general framework which allows us to describe many optimisation problems, including many NP-hard problems [START_REF] Creignou | Complexity Classification of Boolean Constraint Satisfaction Problems[END_REF][START_REF] Khanna | The approximability of constraint satisfaction problems[END_REF][START_REF] Dechter | Constraint Processing[END_REF][START_REF]The Handbook of Constraint Programming[END_REF]. Various restrictions which give rise to classes of problems solvable in polynomial time have been studied in the literature [START_REF] Cohen | The complexity of soft constraint satisfaction[END_REF][START_REF] Grohe | The complexity of homomorphism and constraint satisfaction problems seen from the other side[END_REF].

In any VCSP instance, the variables listed in the scope of each valued constraint are explicitly constrained, in the sense that each possible combination of values for those variables is associated with a given cost. Moreover, if we choose any subset of the variables, then their values are constrained implicitly in the same way, due to the combined effect of the valued constraints. This motivates the concept of expressibility for cost functions, which is defined as follows: Definition 2 For any VCSP instance I = V, D, C , and any tuple v of m variables of I, the projection of I onto v, denoted π v (I), is the m-ary cost function defined as follows:

π v (I)(x) = min{Cost I (s) | s ∈ A, s(v) = x}.
We say that a cost function φ is expressible over a valued constraint language Γ if there exists an instance I ∈ VCSP such that all cost functions in I are from Γ and a tuple v of variables of I such that π v (I) = φ. We call the pair I, v a representation of φ by Γ. 4Example 3 Let φ be the ternary cost function which returns the sum of its arguments. For instance, the assignment defined by s(v

1 ) = 1, s(v 2 ) = 2, s(v 3 ) = 2 has cost φ( 1, 2, 2 ) = 1 + 2 + 2 = 5. Now let I = V, D, C with V = {v 1 , v 2 , v 3 }, D = {1, 2, 3} and C contains the con- straint v 1 , v 2 , v 3 , φ and the (crisp) constraint v 1 = v 3 . The projection of I onto v = v 1 , v 2 is the binary cost func- tion such that π v (I)(x 1 , x 2 ) = x 1 + x 2 + 1 if x 1 = 1, and π v (I)(1, x 2 ) = x 2 + 3.
The intuition is that if x 1 = 1, then the assignment to v 1 , v 2 can be completed with v 3 = 1, and otherwise it can be completed with v 3 = 2 (since 1 violates v 1 = v 3 and thus yields an infinite cost, which is bigger than the finite cost obtained by assigning v 3 = 2).

Note that the notion of expressibility for crisp cost functions (=relations) corresponds to the standard notion of expressibility using conjunction and existential quantification (primitive positive formulas) [START_REF] Bulatov | Classifying the complexity of constraints using finite algebras[END_REF]. We denote by Γ the expressive power of Γ which is the set of all cost functions expressible over Γ up to additive and multiplicative constants. 5From now on, we assume D is totally ordered with some relation <, and we use < to define a partial order on tuples component-wise. We also use the function max componentwise on tuples.

Definition 4 A cost function φ : D k → N is called max- closed if for every u, v ∈ D k , 2φ(max(u, v)) ≤ φ(u) + φ(v).
The class of max-closed relations was first introduced in [START_REF] Jeavons | Tractable constraints on ordered domains[END_REF] and shown to be tractable. In other words, any VCSP instance with max-closed relations over a finite set was shown to be polynomial-time solvable. The class of general and finite-valued max-closed cost functions was introduced and shown to be tractable in [START_REF] Cohen | The complexity of soft constraint satisfaction[END_REF]. Lemma 5 ([4, Lemma 6.14]) If φ is max-closed, then φ is finitely antitone, that is, for all tuples u, v with φ(u), φ(v) < ∞, u ≤ v implies φ(u) ≥ φ(v).

Results

Theorem 7 ([5], alternative proof ) For every d ≥ 3,

(i) R 2 = R 2,3 ; (ii) R max 2 = R max 2,3 ; (iii) R d = R d,2 ; (iv) R max d = R max d,2 ; (v) G max d ⊆ R max d,2 ∪ F max d,1 ; (vi) G max 2 ⊆ R max 2,3 ∪ F max 2,
1 . Theorem 7 was originally proved in [START_REF] Cohen | The expressive power of valued constraints: Hierarchies and collapses[END_REF] using algebraic properties of cost functions: in the case of relations using so-called polymorphisms [START_REF] Bulatov | Classifying the complexity of constraints using finite algebras[END_REF] and in the case of general cost functions using so-called fractional polymorphisms [START_REF] Cohen | An algebraic characterisation of complexity for valued constraints[END_REF]. Here we give a different proof using explicit representations of cost functions which does not rely on these algebraic properties.

PROOF.

(i) It is well-known that any relation R ∈ R 2 can be expressed as a propositional formula ψ in conjunctive normal form: simply a conjunction of clauses which disallow tuples not in R . By the standard SAT to 3-SAT reduction [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF], which replaces a clause (l

1 ∨ l 2 ∨ . . . ∨ l k ) with (l 1 ∨ l 2 ∨ ¬y 1 ) ∧ (y 1 ∨ l 3 ∨ ¬y 2 ) ∧ . . . ∧ (y k-3 ∨ l k-1 ∨ l k )
using new variables y 1 , . . . , y k-3 , ψ is equivalent to ∃y 1 , . . . , y k-3 ψ ′ where ψ ′ consists of clauses of length at most three. Hence ψ ′ , l 1 , . . . , l k is a representation of ψ by R 2,3 . (ii) It is known [START_REF] Jeavons | Tractable constraints on ordered domains[END_REF][START_REF] Gil | Efficient algorithms for description problems over finite totally ordered domains[END_REF] that a relation is max-closed if and only if it can be represented by a conjunction of anti-Horn clauses, that is, clauses of the form (

x 1 ≥ a 1 ∨• • •∨x k ≥ a k ) or (x 1 ≥ a 1 ∨• • •∨x k ≥ a k ∨y ≤ b),
where a 1 , . . . , a k , b are domain values. In the Boolean domain this corresponds to the condition that every clause has at most one negated literal. It is easy to verify that the standard SAT to 3-SAT reduction in the proof of (i) preserves the anti-Horn form of clauses (consider the possibly one negated variable as the last literal). Therefore, the proof of (i) proves (ii) as well. (iii) Because every relation is logically equivalent to some conjunction of clauses [START_REF] Jeavons | Tractable constraints on ordered domains[END_REF][START_REF] Gil | Efficient algorithms for description problems over finite totally ordered domains[END_REF], the proof of (i) gives a weaker result that R d ⊆ R d,3 . We need to show how to express a clause C of length 3 (over the domain D, where d = |D| ≥ 3) as a conjunction of clauses of length 2. Let D = {1, . . . , d} and C = (U 1 (x 1 ) ∨ U 2 (x 2 ) ∨ U 3 (x 3 )) for some literals (unary relations)

U i , 1 ≤ i ≤ 3. We claim that C is equivalent to ∃yC ′ = (U 1 (x 1 )∨N 1 (y))∧(U 2 (x 2 )∨N 2 (y))∧(U 3 (x 3 )∨ N 3 (y))
where y is a new variable and N 1 (y) = D \{1} ("not 1"), N 2 (y) = D \ {2} and N 3 (y) = {1, 2}. It is not difficult to see that a satisfying assignment of C can be extended to a satisfying assignment of C ′ and conversely, a satisfying assignment of C ′ gives a satisfying assignment of C. (iv) It is enough to show that any clause of the form (

x 1 ≥ a 1 ∨ • • • ∨ x k ≥ a k ) or (x 1 ≥ a 1 ∨ • • • ∨ x k ≥ a k ∨ y ≤ b)
, where a 1 , . . . , a k , b are domain values, can be expressed by a conjunction of anti-Horn clauses over at most two variables.

Let C be a clause in R max d :

C = (x 1 ≥ a 1 ∨ x 2 ≥ a 2 ∨ • • • ∨ x k ≥ a k ∨ y ≤ b)
(the case without the y literal is even easier and can be handled similarly). For all i = 1, . . . , k-1, let y i be a fresh variable with d values, so at least three values, say 1, 2, 3 with the natural order. We define the following conjunction of clauses ψ, where y i ∈ {1, 3} is used as a shorthand for y i ≥ 3 ∨ y i ≤ 1 (possible values less than 1 or greater than 3 do not matter) as follows:

ψ = (x 1 ≥ a 1 ∨ y 1 ∈ {1, 3}) ∧ (y 1 ≤ 2 ∨ x 2 ≥ a 2 ) ∧ ∧ k-1 i=2   (y i-1 ≥ 2 ∨ y i ∈ {1, 3}) ∧ (y i ≤ 2 ∨ x i+1 ≥ a i+1 )   ∧ (y k-1 ≥ 2 ∨ y ≤ b).
The intuition is given by reading the second clause as (y 1 ≥ 3 → x 2 ≥ a 2 ) and the third one as (y 1 ≤ 1 → y 2 ∈ {1, 3}). Since the first clause reads "either x 1 ≥ a 1 or y 1 ≥ 3 or y 1 ≤ 1", together with the above implications this gives "either x 1 ≥ a 1 or x 2 ≥ a 2 or y 2 ∈ {1, 3}". Iterating this reasoning, one can see intuitively why the construction works.

More formally, we show that C is logically equivalent to ∃y 1 . . . y k-1 ψ. First, let t be a tuple satisfying ψ. Then if t satisfies x 1 ≥ a 1 , we are done. Otherwise, because of the first clause in ψ, t must satisfy (1) y 1 ≥ 3 or (2) y 1 ≤ 1. In case [START_REF] Bistarelli | Semiring-based CSPs and valued CSPs: Frameworks, properties, and comparison[END_REF], because of the second clause in ψ, t must satisfy x 2 ≥ a 2 and we are done. In case [START_REF] Bulatov | Classifying the complexity of constraints using finite algebras[END_REF], because of the third clause in ψ, t must satisfy y 2 ≥ 3 ∨ y 2 ≤ 1, and we proceed by induction.

Conversely, let t be a tuple satisfying C. We show that t can be completed into a model of ψ by assignments to the y i 's.

Assume first that t satisfies x 1 ≥ a 1 . Then completing t with t(y i ) = 2 for all i = 1, . . . , k -1 yields a model of ψ whatever the values assigned by t to x 2 , . . . , x k , y. This can be seen by examining each clause in ψ. Now assume that t satisfies x i0 ≥ a i0 for some i 0 ∈ {2, . . . , k}. Then completing t with t(y i ) = 1 for all i = 1, . . . , i 0 -2, t(y i0-1 ) = 3 and t(y i ) = 2 for all i = i 0 , . . . , k -1 again yields a model of ψ.

Finally, assume that t satisfies y ≤ b. Then completing t with t(y i ) = 1 for all i = 1, . . . , k -1 yields a model of ψ, which finishes the proof.

Note that this proof makes clear why the same argument does not work for d = 2. Indeed, the "hole" in literal y i ∈ {1, 3} is necessary, since otherwise this literal would be tautologous and thus, so would every second clause be in ψ. We showed in (ii) that in the Boolean case, ternary relations are sufficient and it is well-known (see [START_REF] Cohen | The expressive power of valued constraints: Hierarchies and collapses[END_REF]) that they are indeed necessary. (v) Let φ be an m-ary general max-closed cost function, and write x 1 , . . . , x m for the variables. Let y 1 , . . . , y K be variables (with d values, say 1, . . . , d), and let K = max{φ(x) | φ(x) < ∞} be the biggest finite cost in the range of φ. Intuitively, a cost of k for a tuple will be encoded by y 1 , . . . , y k assigned a cost of 1 (and the others 0). We first encode infinite costs. Let φ R be the relation {u | φ(u) < ∞}, that is, the crisp cost function assigning a cost of 0 to exactly these tuples assigned a finite cost by φ (and ∞ to the others). It turns out that this relation is max-closed. Indeed, for all u, v ∈ φ R we have φ

(u), φ(v) < ∞ by definition of φ R . Since φ is max-closed we have 2φ(max(u, v)) ≤ φ(u)+φ(v) < ∞, so φ(max(u, v)) < ∞ and thus, max(u, v) ∈ φ R . So φ R is max-closed, that is, φ R ∈ R max d .
We now encode finite costs. For an m-tuple t with φ(t) < ∞, write k t for φ(t). We let ψ t be the anti-Horn formula kt j=1 ((

m i=1 x i > t[i])∨y j ≤ 1). Observe that this formula reads x ≤ t → y 1 ≤ 1 ∧ • • • ∧ y kt ≤ 1.
Finally, we define the anti-Horn formula ψ to be ψ R ∧ t∈D m ψ t , where ψ R is an anti-Horn formula equivalent to φ R . By (iv), this formula can be expressed over R max d,2 . The formula ψ encodes the cost of every tuple as a number of y j 's assigned 1. We thus add, to every variable y j , j = 1, . . . , K, the cost function µ defined by µ(1) = 1 and µ(2) = • • • = µ(d) = 0. Clearly, this function is max-closed and therefore in F max d,1 . We now show that ψ, x 1 , . . . , x m is a representation of φ. Let t be an m-ary tuple. Assume first that k t = φ(t) is finite. Then ψ contains the subformula x i > t[i]) ∨ y j ≤ 1). Since obviously t[i] > t[i] holds for no i, every assignment which satisfies ψ sets variables y 1 , . . . , y kt to 1 and thus, has a cost of at least k t . Now let s t be the assignment which is equal to t over x 1 , . . . , x m and which assigns 1 to y 1 , . . . , y kt and 2 to y k+1 , . . . , y K . We show that s t satisfies ψ, which gives an assignment of cost at most k t .

First let ψ t ′ ∈ ψ, and recall that ψ t ′ reads x ≤ t ′ → y 1 ≤ 1 ∧ • • • ∧ y φ(t ′ ) ≤ 1. If t ≤ t ′ , then since φ is max-closed and both costs are finite (by definition of ψ t ′ ), we have φ(t) ≥ φ(t ′ ) by Lemma 5. It follows {y 1 , . . . , y φ(t ′ ) } ⊆ {y 1 , . . . , y kt }, so s t assigns 1 to y 1 , . . . , y φ(t ′ ) and thus satisfies ψ t ′ . Otherwise, if t ≤ t ′ , then t[i] > t ′ [i] for some i and thus s t satisfies ψ t ′ (it does not satisfy its premises). Finally, s t satisfies ψ t ′ for all t ′ . Now s t satisfies ψ R by definition of φ R , since s t equals t over x 1 , . . . , x m and φ(t) < ∞ by assumption. We finally have that for all t with finite cost under φ, s t satisfies ψ and thus, the projection of ψ assigns a cost of at most k t to t. Since the cost is at least k t as shown above, we have the result. Now if t has infinite cost under φ, then by definition of φ R we have that t does not satisfy ψ and thus, has infinite cost under the projection of ψ as well, as desired.

(vi) The same as the proof of (v). In this case ψ is a relation over the Boolean domain, and therefore ψ can be expressed, by (ii), over R max 2,3 .

Definition 6

 6 For every d ≥ 2, we define the following: -R d,m denotes the set of all relations of arity at most m over a domain of size d, and R d = ∪ m≥0 R d,m ; -R max d,m denotes the set of all max-closed relations of arity at most m over an ordered domain of size d, and R max d = ∪ m≥0 R max d,m ; -F d,m denotes the set of all finite-valued cost functions of arity at most m over a domain of size d, and F d = ∪ m≥0 F d,m ; -F max d,m denotes the set of all finite-valued max-closed cost functions of arity at most m over an ordered domain of size d, and F max d = ∪ m≥0 F max d,m . -G d,m denotes the set of all general cost functions of arity at most m over a domain of size d, and G d = ∪ m≥0 G d,m ; -G max d,m denotes the set of all general max-closed cost functions of arity at most m over an ordered domain of size d, and G max d = ∪ m≥0 G max d,m .

See[START_REF] Cooper | High-order consistency in valued constraint satisfaction[END_REF] for a discussion of why limiting ourselves to the N valuation structure is not a severe restriction.

The pair I, v is also known as a gadget for expressing φ over Γ[START_REF] Cohen | The expressive power of valued constraints: Hierarchies and collapses[END_REF], or a constraint network[START_REF] Dechter | Constraint Processing[END_REF], or an implementation[START_REF] Creignou | Complexity Classification of Boolean Constraint Satisfaction Problems[END_REF].

This is only due to technical reasons, but it does not change the complexity of the language, that is, L and L remain polynomialtime equivalent to each other may additive and multiplicative constants be allowed or not[START_REF] Cohen | An algebraic characterisation of complexity for valued constraints[END_REF].

Acknowledgements The authors would like to thank the anonymous reviewers for useful comments on an earlier draft of this paper. Stanislav Živný would like to thank Peter Jeavons for many helpful discussions and ongoing support and encouragement.

(Stanislav Živný 1 This work was supported by ANR grant CANAR (ANR-06-BLAN-0383-02). 2 This work was supported by EPSRC grant EP/F01161X/1.


Krzysztof Nobis

Ewaryst Wierzbicki

Czeslaw Wozniak

On the physical interpretation of nonstandard methods in mechanics


Consuelo Granata
email: granata@isir.upmc.fr

Joseph Salini
email: salini@isir.upmc.fr

Ragou Ady
email: ady@isir.upmc.fr

Philippe Bidaud
email: philippe.bidaud@isir.upmc.fr

Human whole body motion characterization from embedded Kinect

Non-verbal communications such as kinesthetics, or body language and posture are important codes used to establish and maintain interpersonal relationships. They can also be utilized for safe and efficient human robot interactions. A correct interpretation of the human activity through the analysis of certain spatio-temporal and dynamic parameters represent an outstanding benefit for the quality of human machine communication in general. This paper presents an effective markerless motion capture system provided by a mobile robot for sensing human activity, in non-invasive fashion. We present a physical model based method exploiting the embedded Kinect. Its performances are evaluated first comparing the results to those obtained with a precise 3D motion capture marker based system and to data obtained from a dynamic posturography platform. Then an experiment in real life conditions is performed to assess the system sensitivity to some gait disturbances.

I. INTRODUCTION

Social robots must be able to interact effectively with humans, to understand their needs and to interpret theirs orders. These challenging tasks can obviously benefit from considering the human body language (posture, gestures, etc.). Two thirds of communication between humans is non-verbal [START_REF] Hogan | Can't Get Through: 8 Barriers to Communication[END_REF]. The analysis of non-verbal cues provide access to a deeper reading of emotions (e.g. boredom, interest , sadness, etc.) and even intentions of the person. Body language creates a very important part of the communication. The body then gives a certain "form" in mind. Beyond signs, being able to understand and anticipate human activity directed towards a goal is an important element to optimize, through interaction, the robot perception and motor activity. To reach this understanding, we can try to rely on how the brain generates predictions. A hypothesis on the "predictive coding" suggests that the brain continuously generates expectations on sensory input from motor activities. This hypothesis derives from the ideomotor theory [START_REF] Elsner | Effect anticipation and action control[END_REF]. It considers that actions are cognitively represented in terms of their perceived effects. This principle of "ideomotor" argues that during the execution of a particular action, a motor model (reduced order model) is automatically associated with the input of perception representing the perceptual effects of the action [START_REF] Jeannerod | The cognitive neuroscience of action[END_REF]. Perception and action planning are considered similar processes based on the characteristics of activation codes that represent external events. Having such a model must allow to design the robot functions (communication, perception, action) in a form directed by the interaction as any other form of human-machine interface.

We focus here on the analysis of human locomotor activity, that is generally based on the spatio-temporal parameters (Center of Mass (CoM) trajectory, step length, step frequency, step width, walking speed, trunk acceleration etc.) considered to be relevant for gait evaluation [START_REF] Webster | Validity of the GAITRite R walkway system for the measurement of averaged and individual step parameters of gait[END_REF], [START_REF] Zijlstra | Assessment of spatio-temporal gait parameters from trunk accelerations during human walking[END_REF], but also on the parameters characterizing the state of dynamic balance of the postural system. Typically we use the Centre of Pressure (CoP) which is by definition identical to the to the Zero Moment Point (ZMP) [START_REF] Sardain | Forces acting on a biped robot. center of pressure-zero moment point[END_REF]. The CoP/ZMP is the point where the Ground Reaction Force (GRF) vector intersects the plane of Base Of Support (BOS) of the feet on the ground. This reaction force represents the resultant of the contact force distribution on the ground. In the last two decades, many technical solutions were developed to access these different parameters of human movement, e.g. treadmills, instrumented surfaces such as force platforms or particular shoes [START_REF] Brenière | Why we walk the way we do[END_REF], [START_REF] Schepers | Ambulatory estimation of center of mass displacement during walking[END_REF]. To capture the kinematics movement, the most reliable results were achieved by systems such as Motion Capture (MoCap) which exploits high speed digital cameras to capture the 3D motion performed by a subject. Several systems can be used but most commonly the subject is fitted with either passive or active markers [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF]. The passive markers are often covered with an infrared reflective material and then attached to the subject on predefined anatomical landmarks. The estimation of 3D human motions has been recently considered by using single or multi video cameras [START_REF] Sminchisescu | Estimating articulated human motion with covariance scaled sampling[END_REF], [START_REF] Ren | Learning silhouette features for control of human motion[END_REF] (for more details about the motion capture methods we refer to the Poppe's review [START_REF] Poppe | Vision-based human motion analysis: An overview[END_REF]). The collected data are used to animate a model reproducing the anthropometric characteristics of the subject. The kinematic solver is formulated as a frame-by-frame weighted least-squares problem that minimizes the differences between the measured marker locations and the model's virtual marker locations. By introducing thus obtained joint motions, the dynamic results from the multi-body motion with the help of a forward dynamic model [START_REF] Sandholm | MotionLab: a Matlab toolbox for extracting and processing experimental motion capture data for neuromuscular simulations[END_REF]. These technologies are designed to be used in controlled laboratory conditions (with specific environmental set-up or instrumented clothes worn by the person) and cannot be employed for activity assessment in real life environments. Gonzáles et al. [START_REF] González | Estimation of the Center of Mass with Kinect and Wii balance board[END_REF] proposed a non-invasive (marker-less) portable system based on the combination of 2 sensors, a Wii balance board and a Kinect to compute the CoM. This system is unsuitable for the analysis of locomotion because of the use of the balance board. During locomotion, the most challenging phases for postural balancing system are beginning and end the walk, turning, avoiding obstacles (changing the length of the pitch, changing direction, stepping over objects, etc.) and bumping into people and objects. A system able to quantify human responses to perturbations induced by the daily life environment can have a use-value. To our knowledge there is no system that can be used for an in-depth analysis of human activity (especially of posture balance and stability) in non controlled environmental conditions, in a non-invasive fashion (marker-less method) and with an accuracy comparable with marker based technologies. We approach this need with a robot which can perceive and track humans while moving and with the embedded Kinect sensor to capture in real-time multi-segmental human motion by using the Microsoft Kinect SDK. The paper describes a system developed to analyze a person's activity by using a non-invasive (marker-less) motion capture system. A method for determining the individual's anatomical parameters and the algorithms developed for the digital animation model based on the measurement of a number of characteristic points obtained by a Kinect are described. These algorithms are validated experimentally on different exercises during which the positions of the CoP/ZMP, measured and calculated using a posturography platform, are compared. The proposed system is also assessed analyzing the gait of a subject followed by the robot. It shows strong sensitivity to some disturbances obtained manipulating artificially the lower limb joints.

II. EXPERIMENTAL SYSTEM

The robot used for the implementation is a commercial robot, Kompaï (Figure 1) specifically designed for services in private or public spaces. In order to raise the potential of this platform, we added a Microsoft Kinect and an extra laptop for computation. We developed and integrated a framework for the synthesis of smart person following behavior composed of a multimodal human detector and a fuzzy logic based decisional engine (Figure 1). For more details we refer to our previous work [START_REF] Granata | A framework for the design of person following behaviors for social mobile robot[END_REF].

III. FILTERING METHOD

The Kinect sensor belongs to the class of devices known as depth cameras. These cameras interpret 3D scene information based on projected infrared light system called light coding that generates a 3D point cloud. The Kinect sensor has an approximate resolution of 0.013 m per pixel and works at 30 Hz. The 3D point cloud is converted into a depth map, and from a single depth map, body parts are inferred using a randomized decision forest-based approach, learned from over 1 million training examples. A mean shift algorithm is used to robustly compute 3D positions of body joints from the modes of discrete probability distributions. More details can be found in Shotton et. al ( [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF]). The accuracy of the joints position provided by the Kinect skeletal detection algorithm has been investigated in [START_REF] Obdrzálek | Accuracy and Robustness of Kinect Pose Estimation in the Context of Coaching of Elderly Population[END_REF] where the inter-joint distance (the limb length) variations have been estimated as up to more than 0.1 m. For facing this issue an improved skeletonization solution to a space-time constraint problem is required. In previous works, the pose correction was approached in different ways: according to physical-model [START_REF] Tak | Physically based motion retargeting filter[END_REF], by using Kalman-like filters [START_REF] Lee | Motion fairing[END_REF], by adopting different kinds of regressors (random forest regressor [START_REF] Shen | Exemplarbased human action pose correction and tagging[END_REF], nearest neighborhood regressor [START_REF] Duda | Pattern Classification and Scene Analysis[END_REF], Gaussian process regression [START_REF] Schölkopf | New support vector algorithms[END_REF] etc.). We propose a method for pose correction that takes into account both the coherence of the anthropomorphic model of the skeleton and the spatial-temporal motion consistency. This method is composed of 3 steps (Figure 2):

Step 1: Initialization When the person is in the Kinect field of view, his skeleton is detected and the 3D positions of all his joints are provided (Figure 6 left). During the initialization phase, we use a first set of Kinect measurements to define a reference model of the involved body. In particular, the system computes the anthropometric kinematic model of the body according to the mean distances between the consecutive joints.

Step 2: Physical model calibration This step is required to keep the human physical skeleton consistent with the reference model (in anthropomorphic point of view). Our approach relies on a standard mathematical solution, constrained optimization. It consist in minimizing the scalar value quadratic objective function f (X) while respecting a certain number of constraints g(X) [START_REF] Hogan | Can't Get Through: 8 Barriers to Communication[END_REF]. minimize f (X) subject to g(X) = 0

(1)

In the objective function f (X) (Equation 2), X is a vector residing in (3 * nJ)-dimensional space (3D coordinates for nJ joints), a i , b i , c i are the weights associated with the coordinates of the joint i (in our case we consider equal weights on the 3 dimensions).

f (X) = nJ i=1 a i (x i -x i ) 2 + b i (ȳ i -y i ) 2 + c i (z i -z i ) 2 (2)
(x i , ȳi , zi ) and (x i , y i , z i ) are respectively the coordinates of the joint i provided by the skeleton detector algorithm and the coordinates that we are searching for to minimize the objective function while respecting the constraints in Equation 3.

g j (X) = d 2 (J j ; J j+1 ) -d 2 (M j ; M j+1 ) = = (x j -x j+1 ) 2 + (y j -y j+1 ) 2 + + (z j -z j+1 ) 2 -d 2 j,j+1 = 0 (3)
with j = 1, ..., nJ -1, (J j ; J j+1 ) consecutive detected joints, (M j ; M j+1 ) consecutive joints in the model, and d j,j+1 the distance between M j and M j+1 .

All the constraints are equality constraints and are fixed according to the reference model (Initialization step). The optimization of the objective function is made by using the quadratic interior point method (QIPM), which is based on the improvement of initial conditions (measurements) for solving quadratic programming problems.

Step 3: Model Simulation and Parameters Analysis Finally, the 3D trajectories of the kinect skeleton joints virtual markers are used to animate a physical model of the subject on a dynamic simulator Arboris-Python [START_REF] Barthélemy | Arboris-Python[END_REF]. Arboris-Python is an open-source constrained multibody dynamics simulation software written in Python language. It includes a generic and easily extensible set of joints (singularity-free multi-degree of freedom joints, non-honolomic joints, etc.) for the design and modeling of tree structure mechanisms with a minimal set of state variables. It gives access to the completed mechanical properties of the system as well as to the constraints and to the controllers implemented to get the desired behavior of the virtual human. Various control algorithms have been implemented in the Arboris-Python software, from the proportional-integral-derivative controllers to the predictive model based controllers which are used for the control of the locomotion and postural balance task or the interaction tasks with adaptive impedance. The equations of motion of these multibody systems are obtained with the Boltzmann-Hamel formalism [START_REF] Duindam | Lagrangian dynamics of open multibody systems with generalized holonomic and nonholonomic joints[END_REF] from which the first-order approximation of the model is computed. The resulting equations are then integrated using a time-stepping method and a semi-implicit Euler integration scheme. In this way, it is possible to introduce and solve additional constraints, i.e. the kinematic loops, which can be either unilateral (contact) or bilateral (joint), with the help of a Gauss-Seidel algorithm [START_REF] Liu | Computation of three dimensional rigid body dynamics of multiple contacts using time-stepping and Gauss-Seidel[END_REF]. For a complete description of this software you may refer to [START_REF] Salini | Dynamic control for the task/posture coordination of humanoids: toward synthesis of complex activities[END_REF]. Once the model of the subject generated from a generic virtual human model (with 36 degrees of freedom) is instantiated with the kinematic data retrieved from the calibration phase and by inferring the anthropometric table Leva [START_REF] De Leva | Adjustments to Zatsiorsky-Seluyanov's segment inertia parameters[END_REF], the simulation is ran using the 3D cartesian points as target points for the selected joint axis through PID controllers. The input torque vector producing the virtual human motion in accordance with the one tracked by the Kinect are computed by solving a Linear Quadratic Program (LQP) that optimizes a set of weighted tasks (virtual joint marker trajectories and postural control) subject to equality and inequality constraints translating the physical limitations to implicitly satisfy the human motion. For the LQP problem formulation we refer to [START_REF] Salini | Synthesis of complex humanoid whole-body behavior: a focus on sequencing and tasks transitions[END_REF]. A number of parameters characterizing walking or steady state of the person can be obtained from the simulation of the virtual human (spatio-temporal gait parameters, energy, postural balance etc.). In the experiments presented hereafter we consider more particularly the CoP/ZMP. The CoP/ZMP is computed as follows: c = J com (q)q + Jcom (q, q) q where J com is the Jacobian of the CoM, h is its height, z and c are respectively the horizontal position vectors of the CoP/ZMP and of the CoM, and g is the gravity value. These equations directly refer to the Linear Inverted Pendulum Model (LIPM), an unstable system whose control must consider the prediction of the future states of the system. z and c would respectively represent the foot of the pendulum and the projection of the point mass onto the ground. Others parameters such as ankles and knees trajectories are proven be distinguishing factors useful to characterize the postural balance and to detect deviations during walk.

z = c - h g c (4)

IV. SYSTEM VALIDATION

A. Material and procedure

In order to assess the reliability of the system, a preliminary set of experiments were made in a laboratory setting (Fig. 6 right). Five healthy subjects (see Table I for physical body data) were asked to execute 3 different movements (arm movement, rocking movement and side steps) on a posturography platform while wearing 13 CodaMotion markers, to validate the consistency of our system (see Fig. 6 center for marker placement). The Kinect embedded in the robot was simultaneously used for skeleton detection. The goal of this set of experiments was to assess the CoP/ZMP trajectory provided by Arboris by comparing it with the CoP/ZMP trajectory measured with a posturography platform (ground truth). The CoP/ZMP was chosen because it's the only parameter that can be directly measured by a posturography platform. Both CodaMotion and Kinect data are collected and replayed by dynamic simulation and the resulted CoP/ZMP trajectory was recorded and compared with the posturography platform measurements. All data was filtered using a second order lowpass Butterworth filter with cutoff frequency of 10Hz. For Kinect data, the pose correction method described in Section III was applied.

B. Results

To evaluate the consistency of our system with respect to the ground truth, the CoP/ZMP trajectories computed by Arboris replaying both the Kinect and the CodaMotion data, and the CoP/ZMP trajectories measured by the posturography platform were compared. In Figure 3, the 3 CoP/ZMP trajectories related to arm movement executed by one subject are shown. Figures 4 and5 report the CoP/ZMP trajectories corresponding to rocking movement and to lateral steps for the same subject. The max mean error between the trajectories provided by Arboris by using only the Kinect data (with the correction pose algorithm) and the ground truth was lower than 0.08 m. This result proves the accuracy of the system estimating some important parameters related to human activity. Not only the overall motion was well preserved during the simulation, but also the posture assessment showed good results. . Fig. 6. The anatomical joints detected by the Microsoft Kinect SDK algorithm (left), the CodaMotion markers placement on the subjects' bodies (center) and the experimental setting (right).

V. WALKING ANALYSIS DURING PERSON FOLLOWING

A. Material and procedure

System capabilities to extract some walking spatiotemporal and dynamics parameters of a subject have been evaluated in a real life situation. For this, the robot followed a walking person tracked him by using the on board sensors under optimal conditions. The experiments were performed with a subject walking normally, then the mobility of the lower limb joint was artificially manipulated in a controlled way. The objective of these experiments was to evaluate the sensitivity of the system to the disturbances, aiming the analysis of pathological walking activities later. Three experiments were conducted to produce permanent changes by mechanical effects on lower limb:

1) Subject walked in his comfortable walking speed (CWS) wearing his usual shoes; 2) Subject walked in his CWS wearing a brace on the right leg; 3) Subject walked in his CWS wearing a ski boot on the right foot.

Experiments 1 was the reference for the normal gait to which the results of the other experiments are compared. Experiments 2 and 3 aimed to evaluate the system accuracy in detecting anomalies in the cases of loss of mobility of the knee and ankle joints, respectively. The experiments were performed in the subject's CWS to avoid the influence of velocity constraints in the deviations effects. The mechanical disturbances on the subject are shown in Figure 7. The subject's movements were sensed by using only the Kinect embedded on the robot and running the skeleton detector algorithm provided by the Microsoft Kinect SDK. The collected data were then processed following the filtering method described in Section III.

B. Results

The CoP/ZMP trajectory during walking is an important parameter for postural balance assessment. In Figure 8 the CoP/ZMP trajectory and the feet position during the normal walk (Experiment 1) are shown. The CoP/ZMP trajectory is computing by using the LIPM.

Overall, the results presented below show that the system is accurate enough to detect both mechanically induced gait disturbances. The 2 parameters that proved key factors to label walking deviations were: the knees flexions and the ankle lateral movement. As expected, when the subject walked wearing his usual shoes, the average flexion of right and left knees were almost equal (Figure 9). Wearing the brace in the right leg, the mean flexion of the right knee during the support phase was lesser than the mean flexion of the left knee (-0.17 rad of right knee versus -0.24 rad of the left). The use of the ski boot entailed a slight flexion difference as well (-0.23 rad of right knee versus -0.21 rad of the left).

The 2 types of lower limb disturbances tested are distinguishable observing the side movement of the ankles, e.g. the movement orthogonal to the CoM/ZMP during the swing phase. As shown in Figure 10, in the case of normal gait, both movements of right and left ankles were equal. The use of the brace entailed an asymmetry of side movement. The mean amplitude of the lateral movement observed for the right ankle was bigger than the mean value observed for the left ankle (0.07 m versus 0.12 m). The same phenomenon was observed when the subject walked wearing the ski boot (0.07 m versus 0.13 m).

VI. CONCLUSION

This paper proposed a system for the analysis of human activity in everyday life contexts using an autonomous robotic system. We introduced a method for determining the individual's anatomical parameters and the algorithms developed for the digital animation model based on the measurement of a number of characteristic points obtained by the Kinect embedded in the robot. The system was assessed experimentally on different exercises during which the positions of the CoP/ZMP, calculated and measured using a posturography platform, were compared. The potential of the system was also tested in real life conditions exploiting the mobility of the robot and its ability to track persons. These experiments were performed with a subject walking normally, then the mobility of the lower limb joint was artificially manipulated in a controlled way. The system showed sensitivity to the disturbances meaning that it can be exploited to detect pathological walking activities. This research is partially supported by the project ANIPEV funded by European Regional Development Fund (ERDF) and the e-CareBot project funded the DREAM-IT Foundation. .

Fig. 1 .

 1 Fig. 1. Overall system architecture.

Fig. 2 .

 2 Fig. 2. Example of pose correction for a person in standing position. From the noisy data collected by the Kinect to the animation of the dynamic mannequin.

Fig. 3 .

 3 Fig. 3. CoP/ZMP trajectory results of arms movement for a subject. The trajectories issued from Kinect and CodaMotion data are calculated by Arboris.

.

Fig. 4 .

 4 Fig. 4. CoP/ZMP trajectory results of rocking movement for a subject. The trajectories issued from Kinect and CodaMotion data are calculated by Arboris.

.

Fig. 5 .

 5 Fig. 5. CoP/ZMP trajectory results of lateral steps movement for a subject. The trajectories issued from Kinect and CodaMotion data are calculated by Arboris.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Mechanical disturbances on the subject: normal (left), right knee immobilization (center) and ski boot on right foot (right).

Fig. 9 .Fig. 10 .

 910 Fig. 9. Knees angular positions: normal walking (left), right knee immobilization (center) and ski boot on right foot (right).

TABLE I

 I

	.	MAIN PHYSICAL BODY CHARACTERISTICS OF THE SYSTEM
		VALIDATION EXPERIMENT PARTICIPANTS.
	.
		Subject	Height (m)	Weight (kg)	Gender
		P1	1.66	60	F
		P2	1.78	78	M
		P3	1.8	71	M
		P4	1.75	88	M
		P5	1.86	79.5	M


Dr Z Naniewicz

Dr Cz Wolniak

Keywords: In diesem Aufsatz sehlagen wir ein mathematisches Modell yon Verbundwerkstoffen, die nnr begrenzte Belastungen fibertragen, vor. Es wird bewiesen, dab die Randwertprobleme zur gewissen Semi-Variations-Ungleichungen ffihren

On a quasi-stationary model of debonding processes in layered composites Z. 5~aniewicz and Cz. Wo~niak, Warszawa

Introduction

Throughout the paper we deal with composites made of a finite number of sheets (laminae) consisting of parallel fibres embedded in a matrix. It is assumed that every sheet is modelled as a transversely isotropic linear elastic material [START_REF] Jones | Mechanics of composite materials[END_REF]. The sheets are bonded together by very thin layers of the bonding material. These bonding layers are modelled as surfaces, i.e., their thickness is neglected. The proposed approach is based on the assumptions that the bonding material can sustain only restricted values of the interlaminar tractions, and that it behaves elastic prior to possible delamination. Hence during the deformation discontinuities in displacement between the layers can occur. The objective of the contribution is to propose and discuss a mathematical model of debonding processes under the assumption of small displacement and velocity gradients. We show that the displacement jumps across the interlaminar surfaces can be related to the interlaminar stresses by means of some multivalued nonmonotone operators derived from the non-smooth and non-convex strain energy function via the concept of Clarke's generalized gradient [START_REF] Clarke | Generalized gradients of Lipschitz functionals[END_REF].

At the same time the irreversibility of the debonding process is taken into account. It is also proved that the pertinent boundary value problems lead to certain hemi-variational inequalities. In the general case the governing relations are non-local in time. It is shown that after time discretization ~he debonding process is described by the sequence of hemi-variational inequalities and the succesive problems can be treated as time independent. The existence of solutions is established by making use of the theory of generalized pseudo-monotone set valued mappings introduced by Browder and Hess [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF]. The obtained results can be used as the basis for the formulation of homogenized models for debonding problems in layered composites. The proposed approach is an alternative to that given in [START_REF] Wo~niak | Mathematical modelling of delamination effects[END_REF].

Mathematical tools

Throughout the paper we use some concepts of the functional analysis which are summarized below; for particulars the reader is referred to [START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF]. We define ~

indD(w) ----[ 0 if wE J), t +~ if wE W\D
is an example of a convex lowersemicontinnous function defined on W. In the paper we deal also with Lipschitz continuous functions j: W-~ X which may be neither differentiable nor convex and are defined on the Hilbert space W. For such functions we use the concept of the generalized Clarke's gradient [START_REF] Clarke | Generalized gradients of Lipschitz functionals[END_REF]. To this aid we define first the generalized Clarke's directional differential ]~ v) of ]: W --> R at w E W and in the direction v E W by means of

i(w + h + ;~v) -j(w + h) j~ ; v) ~ lira sup h-->0 /~ ~'~0
Then the generalized gradient ~(w) of j(.) at w E W is defined as ~i(w) ~ {~ ~ W: j0(w; v) >= (v, ~), V v E WI where (., .) denotes the scalar product in the Hilbert space W. As it is known [START_REF] Clarke | Generalized gradients of Lipschitz functionals[END_REF] ~](w) coincides with the subdifferentia] ~](w) of ](.) at w E W, if ](-) is convex. The function ](.) is said to be regular in the sense of Clarke, if j0(w; v) = i'(w; v) holds for every w E W, v E W with

i'(w; v) : lim ](w + ,~v) --i(w)

as the one-sided directional differential of j(.) at w E W in the direction v E W.

Basic assumptions and the general statement of the problem

Let ~2 be the regular region in ]R 3 occupied by the undeformed composite which is assumed to be made of S disjointed sheets (laminae). Hence D = u AK (g = 1,..., •), A z being the region occupied by the K-th undeformed sheet, see Fig. 1. Thus the interlaminar surfaces separating the adjacent undeformed sheets are given by HK = (OA~: n ~AK+I) \ ~D where K = 1 ..... S --1. We use the notations The deformations of the composite are analysed in the time interval [0, v I] with the initial time instant related to the undeformed (and unstressed) body. We assume that for every 3 E [0, 3z] the body forces (per mass density) b(x, 3) at x C A and the surface tractions p(x, v) at x ~ F are known where F is the corresponding part of the boundary ~D. We assume also that on F0 -------0D \ F the traces u0(x, r), x E Fo, of the displacement fields u(., ~) : A ~ IR a are know~ for every r ~ [0, Tz] and mes(Fo) > 0. By T(x, T) we denote the Cauchy stress tensor at x ~ A, 3 E [0, @, and by t(z, T) the interlaminar stress vectors at z ~ H, 3 E [0, D] which for z E HK c:: H are acting on the K-th sheet AK across the interface HK. Every sheet is assumed to be modelled as linear-elastic homogeneous material. By e(x), C(x), we denote the known mass density and the tensor of elastic moduli at x E A, respectively, which are constartt in every Ax (K = 1 ..... S). We introduce the strain tensor E(u) ~. (Vu + Vur)/2 with V(.) as the gradient operator and (.)~ for transposition. Under the assumptions that in all problems the inertia terms can be neglected and that no external forces are applied to the interfaces between the sheets we have (

) 1
The Eqs. (1) have to be considered together with the interrelation between the interlaminar displacement jumps [u] (z, v) and the interlaminae stresses t(z, v) where z E H, v C [0, 3f], Here this interrelation will be based on the following physical assumptions:

1. The interlaminar bonding material (which in the undeformed state is modelled by the surface H = u ]/g (K ~-1 ..... S-1) posseses linear elastic properties before possibly debonding, 2. The local debonding conditions (i.e. the debonding conditions for the fixed z E H) are based on the fact that the interlaminar bonding materiM can sustain only restricted values of interlaminar stress.

3. After debonding only unilateral contact between the adjacent sheets can be realized; for the sake of simplicity we neglect the effect of friction.

provided that the inequality n0(z, ~-u]l (z, 0)) < 0 holds for every o C [0, v). After the debonding only the unilateral contact between the layers can be realized and the strain energy of the bonding is assumed to be equal to zero. Thus in the modelling of the interlaminar relation the irreversibility of the debonding process has to be taken into account. To this aid we introduce the function n:

here and in the sequel the operator ~ is taken with respect to ~-u]] (z, ~) for fixed z E/7 and 0(z, z) E {0, 1}. The parameter 0(z, v) is determined by (4) which describes the irreversibility of the debonding process. It has to be emphasized that for the local debonding condition (given by (3) with 0(z, ~) = 1) the values of s(z, z) are not uniquely determined by [START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF].

The formulae (1)--( 5) represent the general governing relations of the debonding process for the composites under consideration. The basic unknown is u : A X [0, vl] -+ ]pa. The problem of finding u(x, ~) for x E A, ~ E [0, v~], such that (1)--( 5) hold, can treated as quasi-stationary since no time derivatives enter into the governing relation. However, the general statement of the debonding problems given above is non-local in time and hence rather complicated to be the basis for more detailed discussion. That is why we propose in the next section the time-local approximation of the foregoing statement of the debonding problem.

Time-local approximation

Now we are showing that after time discretization of the governing relations (1)--( 5) we arrive at some formulation of the time-nonlocal debonding problem which, from the formal point of view, can be treated as a sequence of time-local problems. Setting Zo ~ 0, v~ ~ T I where M is a sufficiently large positive integer, and introducing the finite sequence ~o, ~1 .... ,3M of time instants such that ~o < ~ < ... < ZM, we assume that the general relations (1)--( 5) presented in Sect. and (4) has to be replaced by [START_REF] Wo~niak | Constraints in constitutive relations of mechanics[END_REF]. All aforementioned relations have to hold for m = 1,2 ..... M. Now observe that for m = 1 we have 0(z, 3,,-i) = 0(z, 0) = 1 for z E/7, and the basic unknown in ( 7) is the displacement field u(., zl). Assume that u(., ~1) has been obtained as a specific solution of (7) (in Sect. 6 we shall prove that such solution exists) ; then from ( 6) we can calculate the function 0(., vl). Hence for m = 2 the basic unknown in ( 7) is the displacement field u(-, v2) since now 0(-, ~) is the known function. Using this procedure we can calculate u(., Tin) for m = 1, 2 ..... M by finding the successive solutions to [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF] (provided that such solutions exist) and using successively [START_REF] Wo~niak | Constraints in constitutive relations of mechanics[END_REF]. All these problems can be treated, from the formal point of view, as time-local problems. Moreover, in every such problem the time instant ~ plays the role of ~he parameter.

Variational formulation

Taking into account the results of Sect. 4 we restrict ourselves to an arbitrary but fixed time instant z = vm and we assume that 0(z, ~_,) for z ~ H, has been calculated previously. Then, from the formal point of view, [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF] for v = z~ can be treated as the governing relations of some time independent problem for the function u(., vm) : A --~ ]R a. For the sake of simplicity we omit the argument ~,~ in [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF] and write fi(z, [u~ (z), 0(z)) instead of fi(z, ~-u~ (z, zm), 0(z, ~m-1)).

Let V be the space of the sufficiently regular test functions v: A --~ ]R 3 such that every vlA~ (K = 1, ..., S) has well defined traces on OA K and v(x) = 0 for x E/~0. We assume also that for every v E V the divergence theorem can be applied to v(x) 9 Div T(x), x E A. Taking into account the possible discontinuity of v(.) across H = u HE (K = 1 ..... S --1) we obtain

f v(x). Div T(x) dV = --f tr [Vv(x)T(x)] dV + f v(x). T(x) n(x)

dA A A F S--1 + E f [v-(z) 9 T-(z) --v+(z) 9 T+(z)] N(z) dA K=I HK
where we have used the notations introduced in See. 3. Then it can be shown that [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF] for ~ = v~, in the simplified notation introduced above, is equivalent to the equations s(z) (z),

f tr [T(x) Vv(x)] dV + f t(z). Iv] (z) dA = f Q(x) b(x). v(x) dV A H A + f p(x) 9 v(x) dA for V v E V, (8)
In this section we are to show that under proper regularity conditions (8), (9) lead to some hemivariational inequality [START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF]. We bear in mind that ~(x), C(x) are constant in everyAa (K = 1, .. implies that the first condition from (9) holds almost everywhere on H. Now we have to establish the interrelation between the functional s* and the displacement u(.), such that the second condition from (9) holds. To this aim we define the functional +(v)=f~(z,~q(z),~(~))dV, v~ v ]7

bearing in mind that J(v) depends also on the function 0(.). Taking into account the definition of fl(.) we conclude that J(.) is Lipschitz continuous on V and is regular in the sense of Clarke (Sects. 2 and 3). The interrelation between s* and u(.) that we are going to esta.blish will be determined by the Lemma. Let s* E ~J(u). Then there exists s C (L2(//)) a such that (i) <v, ~*> = f s(z). [q (z) dA for V v C V, H (ii) s(z) E @(z, [u] (z), ~(z)) for almost every z ~ 1t.

Pro@ Let us assume that for the fixed 0(.) the function fi: H XIRS-+ IR given by fl(z, w) fl(z, w, O(z)) has the properties 1. fl(., w) is measurable for any w E IR 3, 2. fl(., 0) E El(//), 3. Ifl(z, wl) --fl(z, w2)l < k(z) []wl --w~ll, k(z) > 0, for every wl, w~ E IR 3, z E/I with k(-)

L2(H).

  ~-]P, u { t e o} u {--~} and ~,+ ~ ll(+ v {0} t~ (~-~}. Let ~: W --> ]l( be a convex lowersemicontinuous function defined on a linear topological space W and let W* be a dual of W, i.e., the space of all linear continuous functionals defined on W. The value of w* E W* on the element w E W is denoted by (w, w*) and the bilinear form (-, .): W• W* --~ ~ determines a dual pairing between W and W*. The set of all w* E W*, such that (v -w, w*) <= ~(v) -~(w), V v c W, and ~(w) ~ +~, is a subdifferential of r at w E W and is denoted by ~(w). If D is a closed convex set in W, then the indicator function of D given by

7 Fig. 1 ,

 71 Fig. 1, Scheme of the undeformed composite

  DivT(x, 3)+0(x)b(x,r)~0 at xEA, T(x, 3) = C(x) E(u(x, 3)) at x ~ A, T(x, r) n(x) = p(x, 3) st x ~ F, 9 ~ [0, TI] n(x, 3)=Uo(X, 3) at XCPo, T+(z, 3) I~(z) --T (z, r) :N(z) = t(z, 3) at z C H.

  In order to specify the aforementioned assumptions and to establish the pertinent interlaminar relation we introduce the decompositiont(z, v) = s(z, ~) + r(z, v), r(z, z) : N(z) r~(z, v) for z C H(2)where s(z, v) represents the interlaminar stress (before the local debonding of sheets) which is due to the elastic properties of the bonding material and where r(z, T) : N(z) ru(z, v) is the reaction due to the impenetrability of the adjacent sheets. The impenetrability condition has the form of the uniteral constraints [u~ (z, ~) ~ 0 (z E//, v E [0, vl]), and [u] (z, v) : 0 implies rx(z, ~) ~ 0 as well as ~u]]~ (z, v) > 0 implies r~(z, v) : 0. It can be easily veryfied that the aforementioned conditions can be written in the simple form r~(z, 3) E e ind~+([U]N (z, v)) for z E//, v E [0, vii which also follows immediately from the general constraint approach [6]. Setting we assume that before possible debonding the strain energy function of the bonding materialfor every z E/7 is given by the function n o : H • R a ~ R defined as YN, YT, ~ are positive constants. Here y~, YT can be treated as moduli (longitudinal and shear, respectively) which characterize the properties of the bonding material. The constant characterizes the local debonding condition which is postulated in the form ~z0(z, [u]] (z, z)) :0 for zE/7, v~[0,~i]

/ 7 •

 7 If( ~ --> ]l( given by n(z,w)~~ if no(Z,w)<O, v(o))~ {~ if ~(z,v(o))<0 for every o E [0, v), 0=<a<~ if ~(z,v(o)) ----0 for some o E [0, v) defined for every z E//, 9 E (0, vZ] on the set of (sufficiently regular) functions v : [0, v) -+ ]l( for an arbitrary time instant v E [0, vii and in an arbitrary deformation process the strain energy of the interlaminar bonding is determined by the function fl: H X~ a X {0, 1} -+ R, given by (z, 0(z, 3)) n(z, in-1 (z, v)) 0(z, with n(.) and 0(.) as introduced above. The functions fl(z,., 0(z, T)) defined on R 3 are non-smooth and non-convex but for every z s H, 0(z, z) E {0, 1} are regular in the sense of Clark (Sect. 2). The local debonding takes place, if 0(z, z) = 1 and ~u]] (z, 3) satisfies (3) ; in any other case tile functions fl(z,., 0(z, z)) are smooth and the values of s(z, 3) can be derived from fl(z,., 0(z, v)) as the derivatives with respect to ~'u] (z, 3). However, using the concept of the generalized Clarke's gradient of fi(z,., O(z, ~)) we also obtain the general interrelation between s(z, T) and [u] (z, ~): s(z, ~) ~ @(z, [u~ (z, 3), O(z, ~));

7 )

 7 3 have to hold only for 3 E {v0, ~1 ..... VM}. Let us define 1 if ~r(z, [u~ (z, ~)) < 0 for every ~ ~ {Vo ~1,..., Tm--1} (6) 0(z, Vm_,) ~-0 if z(z, In] (z, a)) = 0 for some a E {To, r,,..., Tm_l}' for every z C H and for m = 1, ..., M. Taking into account the definition of the functionals O(z, [u~ (z, .)) in (4) it is easy to see that after the time diseretization the functions 0(., ~m) have to be replaced by the functions 0(., ~m-~) where m = 1, ..., M. Thus the time discretization implies the approximation O(z, zm)~-~O(z, 3~_l) with m=l,...,M and zE//. Hence (1)--(3) and (5) take the form DivT(x,~m)+O(x)b(x,,~)=O for xCA, T(x, vm) = C(x) E(u(x, ~m)) for x E A, T(x, ~m) n(x) = p(x, vm) for x E F, u(x, ~) = Uo(X, ~,~) for x E 17o, (T+(z, 3~) I~(z) --T-(z, ~m) N(z) --t(z, ~m) for Z E /7, t(z, 3~) = S(Z, Zm) + r(z, z~), r(z, z~o) = ~'(z) ra(z, zm) for z E /7, r~(z, Vm) C Oind~ § (z, vm )) for z C/7, s(z, zm) E ~fl(z, ~'u]] (z, Zm), 0(z, vm-x)) for z C/7

F

  , S), and we assume that b C (L~(sg)) a. p C (Le(I')) ~. Using the known notations for the Sobolev spaces we define V --{v C (L~(~))~: VIA~ C (HI(A~)) ~ (K = 1 ..... S); ~oV = 0} where yov is the trace of v on To, and ZK ~ {T E (L~(AK))ax3: Div T C (L2(AK)) a (K = 1 ..... S)} (for an arbitrary set M the symbol M ~ x 3 stand for a set of all 3 X 3 symmetric matrix elements which belong to M). In order to simplify the analysis, from now on we assume that uo(x, ~) = 0 for x E/'0, T C [0, ~I] (cf. (1)). Then the displacement field u(.) in (8) can be treated as an element of the space V and [u~ E (H1/~(H)) a. Moreover, we assume that T[A~ E Z~ for K= 1 ..... S; then t E (L~(//)) a. Hence all integrals in the variational condition (8) have a well defined sense. Denoting by (., .} the bilinear form defined on V X V* and setting <., v)--f tr [r(~(~))C(~)E(u(x))] dr, A (v,/} ~ f e(x) b(x). v(x) dV § f p(x). v(x) dA, , r*} ~-f r~(z) ~V]]N (Z) dA // we obtain from (8) the variational condition ~(u, v) + (v, s*) + (v, ~*) = (v,/) for V v ~ V 00) in which the functionals s*, r* E V* have to be interrelated with the displacement field u C V by means of the conditions (9). In order to interrelate the functional r* with the displacement field u(.) we introduce the closed convex set ~7 in V by means of IN: ~-{v E V: ~-v]]~ (z) > 0 for aimost every z E H}. Then .r* C ~ind~(u) (nn)

  Under the conditions 1.--3. the functional 2: (L2(H)) 3 -+ ]R given by J(w) =:f fi(z, w(z))dA, w E (L2(T/)) a is Lipschitzian and has the generalized gradient (Sect. 2) ~J(w) ~ {s E (L~(//))8: OT~ v) > (v, w) for every v E (L2(H)) 3} where (., .) is a scalar product in (L2(II)) a and ](w + h + ~v) -Y(w + h

is Clarke's directional differential at w and in the direction v. Iris known [START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF] that, if s E ~J(w), then s(z) C ~fl(z,w(z)) for z E and hence s(z) c (z), 0(z))

for almost every z E//. On the other hand the condition s* E ~J(u) implies the existence of s E (H1/2(II)3) * with (V, 8}B~ : (~'V'~, S)1/2 where (., ")1/2 is the bilinear form determined by the pertinent dual pairing. Since the space H1/2(1I) is dense in L2(II), it is possible to extend the functional s* defined on (HI/2(//)) s onto (L2(H)) a in such a way that s c

By virtue of the meaning of ~J(.), the aforementioned condition leads to the conditions (i)--(ii)

of the Lemma which ends the proof. Now making use of ( 10), (11) and the Lemma formulated above we obtain

where u E V, / E V*. Taking into account that the functional J(.) depends on 0(7) ://--> {0, 1}, J(v) -~ J(v; 0), and that u ----u(7~), / ----/(7m), 0 = 0(~-1) we arrive at the sequence of inequalities

for every v E V, u(Tm) E V, m = 1, 2, ..., M; this represents the variational form of the governing relations [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF] with 0(7m-1)-0(', 7m-1) defined by [START_REF] Wo~niak | Constraints in constitutive relations of mechanics[END_REF] and no(x, 7m)= 0 for every x E Fo and m = 0, 1 .... , M. Let/(Tin) (m = 1, 2 .... , M) be known. Then the successive solutions u(vm) of the inequalities (13) for m = 1, 2 ..... M (provided that they exist) can be obtained from the sequence of the time-local problems for n(7~), since for every v ~ 7m the function 0(V~-l) --~ 6(., 7m-1) has been previously calculated after having solved the problem for v = 7m-1 (cf. the comments at the end of Sect. 4). Hence for every 7 = 7,, (m --~ 1, 2, ..., M) in the problem of finding the displacement field u(7~) = u(-, 7m) from (13) we assume that the functional J(., O(Vm-1)): V -~ ]R is known. Since fl(z, [ui] (z, 7m), 0(z, 7m-,)) ~-a(Z, [nil (Z, 7m)) 0(Z, ~m--1), the aforementioned functional is defined by J(V; 0(7,.-1) ) ~. f Yl(g, ~-V' ~ (Z))0(Z, 7m_1)di H for every v E V and m = 1, 2,..., M. The functional J(.; 0(v~_l) ) will be referred to as the superpotential of the interlaminar bonding tractions. Due to the time irreversibility of the process it depends on the argument 0(V=-l) given by ( 6). The variational formula (13) constitutes the hemivariational inequality [START_REF] Panagiotopoulos | Inequality problems in mechanics and applications[END_REF] for the displacement field u(v,,) E V. This fact makes it possible to establish the existence of solutions to the problems under consideration which will be done in the following section of the paper. Let us remark that after replacing in (13) u(7=) ~nd 0(7~-1) by u [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF] and 0 [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF], respectively, we obtain the variational form of the general governing relations (1)--( 5) (for u0 0) which has to hold for every v E [0, 7]].

Existence

We restrict ourselves to the problem of existence of a solution u(v,~) to the hemi-variational inequality (13) for an arbitrary but fixed m E {1, 2 ..... M} under assumption that the function O(vm_l) has been previously calculated. Hence the superpotential J(.) = J(-, 0(7m_1) ) is assumed to be known and the problem is governed by the hemi-variational inequality of the form (12)

where for the sake of simplicity we have neglected the time parameter vm and the dependence of the superpotential J on the argument 0(., v,~-l). The analysis is carried out under the extra assumption that mes(Fo r, ~AK) > 0 for K --~ 1 ..... S and is based on some facts from the theory of generalized pseudo-monotone set valued mappings introduced in [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF]. The reader who is not familiar with this theory can pass now to the theorem formulated at the end of this section.

Let us introduce the Lamd operator A : V -> V* by means of

It is easy to see that the hemi-variational inequality (12) is equivalent to the multiequation / E Au + ~J(u) + ~indK(u ).

(

Thus in order to prove that for every / E V* there exist the displacement field u E V such that (14) holds, it suffices to show-that Range(A + ~J + ~indK) : V*.

(

First we establish for this purpose the pseudo-monotonicity of the multioperator ~J and then we apply the known theorem on the range of the generalized regular pseudo-monotone mappings [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF]. Using the results of [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF] we show that the following three conditions hold:

(i) ~J(v) is a non-empty, closed, convex and bounded set for every v E V;

(ii) ~J is an upper semi-continuous mapping from the strong topology on V to the weak topology on V* ;

(iii) if (vn) is a sequence in V weakly converging to v, v n E }J(vn) and lira sup (v~ --v, v*) ~ 0, then for any w E V there exists v*(w) E ~J(v) with lim inf (v, --w, v*) >_ (v --w, v*(w)).

It is known that the generalized gradient of a Lipschitz function satisfies (i) and (ii) [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF]. To get (iii) we use the method proposed in [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF]. Let (v,) converge weakly to v and let v* E ~J(v~). Taking into account the semi-continuity of j0(., .) [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF] we obtain, by passing to the limit, for some v*(w) E ~J(v) which ends the proof of (iii).

The generalized gradient ~J being pseudo-monotone satisfies the condition

with a constant c, since J(.) is Lipschitzian on V. Thus, in the terminology of [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF], ~J(.) is a generalized regular pseudo-monotone operator. On the other hand, it is obvious that

E Domain (A -Jr ~indK).

Thus A -+ ~ind K ~-~J is generalized regular pseudo-monotone [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF]. The coercivness of A on V follows from Korn's inequality, from the assumption that mes(F o n ~AK) ~ 0 for K = i ..... S, and from the strong ellipticity of the tensor of elastic moduli C(x), x E A. It implies the coercivness of A + ~ind~ + ~J and by the theorem given in [START_REF] Browder | Nonlinear mappings of monotone types in Banach spaces[END_REF] we arrive at (15). The obtained result can be summarized in the form of the Theorem. Under the regularity conditions introduced in Sect. 5 and the condition mes(F o n AK) > 0 for K ~-1, ..., S, the hemi-variational inequality (12) posseses at least one solution for an arbitrary / E V*.

Final remarks

The quasi-stationary model of debonding processes in layered linear-elastic composites, given by the sequence of hemi-variational inequalities (13) and by the conditions (6), constitutes the starting point for further approximations rather than for the direct analysis of special problems and numerical solutions. As a rule, the layered composites are made of a great number of sheets. Hence the great number of interfaces IIK (K z 1 ..... S --1), across which the displacement field may suffer discontinuity, argues against the introduced model from the point of view of engineering applications. However, the mathematical model of debonding processes proposed in this paper can serve as a foundation for the derivation of some homogenized models, e.g., those based on the microlocal modelling approach [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF]. The derivation of such models will be performed in a forthcoming paper. On the other hand the governing relations obtained in this paper can be generalized by introducing in [START_REF] Wolniak | A nonstandard method of modelling of thermoelastic periodic composites[END_REF] the inertia forces and by interpreting the conditions (5) as the constraints imposed on s(z, ~). On this way, using the method of constraints [START_REF] Wo~niak | Constraints in constitutive relations of mechanics[END_REF], it is possible to formulate the evolutionM nlode]s of the debonding processes [START_REF] Wo~niak | Mathematical modelling of delamination effects[END_REF]. Let us remark that the general line of approach to the debonding problems proposed in this paper can be applied also to layered composites made of inelastic materials.


Stanislaw J Matysiak

Czeslaw Wozniak

On the microlocal modelling of thermoelastic periodic composites


shall assume e << 1. For an arbitrary V-periodic integrable function f = f(x) we introduce the averaging operator

1 6r> = fvf(x) do,
where dv----dXl dx2dx3. By ui, s~j, ti, we denote the displacements, stresses and boundary tractions, respectively. The material properties are described by elastic modulae C~jkt = C,'jkt(X), viscous modulae Dijkt = D;m(x) and mass density p = p(x), which are V-periodic functions of x,

The body forces are denoted by b i and are assumed to be constant.

GOVERNING RELATIONS AND MODELLING HYPOTHESES

The starting point of the proposed method of modelling is the known governing equations of the composite body under consideration, given by the principal of virtual work f sijt~u~jdv=~o ti6ulda+f P(b,-ii,)Su~dv, [START_REF] Bensoussan | Asymptotic Analysis of Periodic Structures[END_REF] which is assumed to hold for every virtual displacement ~U i and the constitutive equation of the Voigt material Sij = CijklU(k,l) "~ Oijk/t(k,t).

(

In micro-periodic composite materials p = p(x), C;jkt = Cijkt(x), Dijkt = Dijkt(X), X e ~, are V-periodic functions, where length dimensions of V are very small as compared with those of fL Hence the direct description of the composite body, in which no restrictions on displacement fields are imposed, leads to the system of partial differential equations for u; with variable highly oscillating (V-periodic) coefficients. The analytic form of these equations is not convenient for a computational analysis of special problems and that is why we look for certain simplified (averaged) mathematical models of the composites under consideration.

The method of modelling proposed below will be based on that introduced by Wo~haiak [4] for the linear elastic composites and generalized by Wierzbicki et al. [5] on the non-linear elastodynamics of composite materials. The model introduces two auxiliary concepts and is based on certain macro-modelling hypotheses. To make the paper self-consistent the foundations of the modelling approach will be recalled below.

The first auxiliary concept is that of the regular macro function. Denoting by AF the accuracy of computations related to the values of a function F(x, t), x ~ f2, t ~ R, we shall refer to this function as the macro function, if for each x, y e f2 and every t, such that x -y ~ V, condition IF(x, t) -F(y, t)l < AF holds. If F is continuous and has continuous derivatives then it will be called the regular macro function provided that similar conditions (with pertinent computation accuracies) also hold for all derivatives of F. Roughly speaking, increments of regular V-macro functions and their derivatives on an arbitrary cell x + V in the domain of their definition can be neglected when compared with values of those functions at x.

The second auxiliary concept is that of the micro-shape function system. This is a system of n, n-> 1, linear independent functions ha(x), a = 1 ..... n which are V-periodic, continuous, have piecewise continuous first order derivatives and satisfy conditions: (ha)= 0, (pha)= O, ha(x) e O(l) and ha.i(x)e O(1), for every x (i.e. values of ha,~ are independent of the microstructure length parameter l).

The avenue leading from equations ( 1) and ( 2) to the proposed macro-description of viscoelastic composite materials is based on the following hypotheses.

The kinematic hypothesis

This hypothesis states that the displacements in the composite body under consideration can be expected in the form ui(x, t) = Ui(x, t) + ha(x)Va(x, t),

x E ~, t -> 0, (3)
where U, Va are certain regular macro functions and h a is the system of micro-shape functions, postulated in every problem under consideration. Functions U~ and Va constitute the new kinematic variables called macro-displacements and correctors, respectively. The choice of the micro-shape functions depends on the character of the material heterogeneity in the representative volume element V and on the expected accuracy of the proposed model. Generally speaking, formulae ui(xo, t)=Ui(xo, t)+ ha(x)V~(xo, t), x E Xo + V(xo) have to represent a certain discretization of an arbitrary but fixed element Xo + V(xo) of the region f£

The virtual work hypothesis

This hypothesis assumes that the prinicple of virtual work (1) holds for

~UI(I ) : ~Ui(x ) -~-ha(x)SVa(x), x ~ ~, (4)
where ~Ui, 6Va a:re arbitrary linear independent regular V-macro functions. Moreover, it is assumed that 6ui(x) = 6U/(x) for x ~ 0f£ This hypothesis is strictly related to the kinematic hypothesis formulated above.

The macro-approximation hypothesis

This hypothesis states that the following will be neglected: (i) Terms O(AF) in formulae involving values of F, where F is an arbitrary regular macro function.

(ii) Terms O(e) in calculations of integrals over ft.

The above hypothesis takes into account the definition of a macro function and the condition e<<l.

REFINED MACRO-MODEL

Substituting the right-hand side of equation (2) into equation (1), using equations ( 3) and ( 4) and applying the macro-approximation hypothesis, after rather lengthy calculations which will be not presented here, we arrive at the system of equations for the new kinematic variables U. Va. This system can be written down in the form of equations of motion 5) it can be seen that we obtain the system of differential equations for ~, V~ with constant coefficients. Hence, the resulting equations represent a certain macro-model of a periodic composite which can be applied to the analysis of special engineering problems. Moreover, this system of equations involves in the explicit form the microstructure length parameter /. The relations obtained describe the microstructure length-scale effect of the body behaviour, which is why equations ( 5)-( 7) will represent what will be called the refined macro-model of the composite under consideration. The characteristic feature of the model is that the correctors V~ are governed by the ordinary differential equations, involving the second and first time derivatives of correctors. It means that the correctors are independent of the boundary conditions and can be interpreted as certain internal variables. For homogeneous bodies, by virtue of (C~jkth~,k)= C~jkt(h~,k)= O, and under homogeneous initial conditions for correctors we obtain VT(x, t) = 0 for every x E ~, t -> 0. Hence the correctors describe the effect of the material inhomogeneity on the behaviour of the composite. For linear-elastic materials D~jkt = 0 and equations ( 5) and (7) reduce to the form proposed by Wofiaiak [4] and representing the linear refined macro-elastodynamics of periodic composites.

Sij,j -(p )Ui "~-(p )bi = O, 12 yab ~T b i dr H a = O,

LOCAL MACRO-MODEL

The local macro-model of the composite under consideration can be derived by neglecting terms O(l 2) in equations (5). In this case the equations of motion reduce to the form

S~j,j -(p)l)i + (p)bi = O, (8)
natural boundary conditions have a form similar to that given by equation ( 6

In this case correctors V~ are the internal variables and equation (10) can be referred to as the evolution equation. It can be observed that the above local model of viscoelastic composites is much simpler and useful to applications than that obtained by the asymptotic homogenization approach, cf. [START_REF] Bensoussan | Asymptotic Analysis of Periodic Structures[END_REF][START_REF] Bakhvalov | Osrednienie Processor v Periodiceskich Sredach[END_REF].

Let us observe that for quasi-stationary problems (where inertia forces are neglected) the refined and local macro-models of composites coincide.

At the end of general considerations it has to be noticed that the solution U~, V7 to the initial-boundary value problems for equation ( 5) and (7) as well as for equations (8)-(10) have a physical meaning only if Ui, V~ are regular V-macro functions. This condition imposes certain restriction on the boundary traction si in equation ( 6) and is also related to the boundary-layer problem for composite materials, cf. [START_REF] Sanchez-Palencia | Homogenization Techniques for Composite Media[END_REF].

LAMINATED MATERIALS

Let us apply the general results of Sections 3 and 4 to the two-constituent periodic laminates. The representative layer of these composites is made of two' homogeneous laminae of thicknesses 1', l", respectively. Hence l = l' + l" is the microstructure length parameter. We shall assume that the xl-axis is normal to the lamina interfaces and introduce one /-periodic and equations (7) yield Hi = (C,,k,(h,,)2)Vk + (Cakth,,)Uk,, + (Dilkl(h,l)2}12k + (Dak,h,1}l[]k,t.

S# = (Ci]kl)Uk,l + (Ci]lth,l)Vl + (Dqkl)(-]k,l + (O,juh,1)~,

(12) Under an extra assumption that every material constitutent is isotropic with the Lam6 modulae A, 1. and the viscosity */we obtain Cijkl : Aaij~kl "~-g(~ik~l + ~il~jk),

Dqk, = 2*/( Sik~]t--~ 8i]tSkl),

where A, 1., 7/take the values A', 1.', */' and A", 1.", */" in the laminae of thicknesses l' and l", respectively. Let c run over A, 1., 7/and define l' l"

e>=Tc' +7 c',

[c] -= c" -c', lc' lc" {c}-~--+--. l' 1"
For isotropic constituents and under the above denotations the constitutive equations ( 12 + ~<~>(0,,i, "}-03,23 ) "]-2V'3 [jl~]V2, , "Jr 2V'-3 [h]V, 2 + 2V"3 ['q](l?2,, 2 • -~ v,.2) <p>02 = o, <2/, + A)U3.33 + (Iz)(U3,n + U3,22) q-(]A -[-/~)(UI.13 q-U2.23 ) -I-<]j~) •3.,i -l-U3,22 q-~ 03,33

1 "a t-~ ('0)(L~/,,13 + 02,23 ) -am 2V3 []./,]V3, , + 2V3 [hlV,.3 + 2V'3 [rl](l?3,, + ~ I71,3) -(O)t)'3 = O, X/3 [2t* + A]U,,, + X/'3 [Al(U2,2 + U3,3) + ~ V~ [r/l(2U,,, -02,2 -/-713,3) + 3{4t, + 2a.}v,8{rl}f', + 112<p)9. ' = O, 2 vr3 []l£](U1, 2 -]-U2,1) q-V~ [n1(~1,2 -]-02,1) q-6{/z}V2 + 6{7/}15"2 + 1 12(p)(22 = O, 2 12 V'3 [/*](U,,3 + U3.,) + V~ [,/](/5"1.3 +/)'3,,) + 6{/*}I/2 + 6{,/}15"3 -~ l (P)?3 = 0. ( 13
)
The above equations describe the laminate under consideration in the framework of the refined model. For the local model we neglect in equations ( 13) the terms involving the length parameter I.

EXAMPLE

Using the equations derived in Section 5 we shall investigate the length-scale effect on the behaviour of the composite body. To this end we shall consider the periodic laminated material occupying the region f~= (0, L,)x (0, L2)× (0, L3) where the microstructure length parameter 1 (the thickness of the repetitive layer) is assumed to be very small compared with L-= min{L,, L2, L3}. Let the body be subjected to the extension 6 in the xraxis direction. Hence, the boundary conditions for U, are given by: U, = 0 for x, = 0, U, = 6 for x, = L,, U~ = 8L-'x, for x2 = 0, L2 and for x3 = 0, L3; at the same time U2 =/-/3 = 0 on f~. The initial conditions will be assumed in the form Ul(x,O)=Sx,/L, U,(x,0)=0, U2(x,0)= U3(x,0)= Uz(x, 0) = U3(x, 0) = 0 and V/(x, 0) = V/(x, 0) = 0 for every x e f~.

Let us define the non-dimensional positive parameter:

It can be shown that the solution to the aforementioned boundary value problem in the framework of the refined model is related to ~: by the following conditions:

Since ~ depends not only on the material properties of constituents but also on the microstructure length parameter then the response of the viscoelastic body in the dynamic processes depends on the size of microstructure. Neglecting in equations ( 13

CONCLUSIONS

In order to show the differences between the proposed length-scale model and the homogenized model let us compare the solution ( 14)-( 16) with that given by equations (17). It can easily be seen that the formulae for stresses S,j obtained in the framework of the local

( 5 )

 5 where jab_ (phahb)l-2 are inertial modulae (the values of which are independent of the microstructure length parameter), natural boundary conditionsSiynj = si(6) and constitutive equations Sit = (Cijkl)Uk, l "~-(Cijlkhak)va + (Oi]kl)~-]k, l + (Oi]klha, k)('ra, Hi = (Cijkth jh ,,)Vk + (Cijkth~,j}Uk,, + (Dqkth~.jhb, jhb,,)(/'b + (Oijklha,])~-]k,I • Fields Sij, Ha are called the averaged stresses and microdynamic forces, respectively. After substituting the right-hand sides of equations (7) into equations (

  ), the constitutive equations are given by SO = (Cokt)Uk,t + (Cq~th~,k)V~ + (Diykl)Uk,t + (Dijklh~,k)?~, (9) and for correctors V~' we obtain the system of ordinary differential equations (Dijkth~hbt) ~ "b + ( Cqkth~.hbt) v b = -(Dol,,h~.)(.lk,t -( C~jkth~)Uk,,.

  micro-shape function h=h~(xl) which in [0,1] satisfies conditions h(O)=h(l)=l/X/3, and is linear in [O,l'], [l', l], where the coordinate planes xl = O, xl = l', h(l') = -l/V~ x2 = l' + l" coincide with laminae interfaces. In this case equations (5) reduce to the form Soj -(p)O~ + (p)b~ = 0 12(p)~"i "Jr H i = 0 (11)

  ) read Sn = 2(1*)Uu + (a)(Ul,l +/-/2 2 +/-/3,3) + 2V~ [2t* + A]V1 + 2 ' 3 (*/)(E,l -/5"2,2 -03,3) + [*/l~, $12 = $21 = (1*)(U,.z + V2a) + 2V3 [1.]V2 + (*/)(01,2 + U2,,) + 2X/3 [,/]122, Sl 3 = S31 : (1*)(U1. 3 -~-U3,1 ) + 2X/3 [/.1I/3 + (*/)(E,3 + U3a) + 2X/3 [*/193, 2 . 4~/3 SEE = 2(1.)U2.2 + (/~)(U1,1 + U2,2 + U3,3) + 2X/3 [/~]V 1 -~-~ (*/)(2U=,= -~/1,1 -03,3) ---[*/]121, U3,3 + (a)(Um + U2,2 +/-/3,3/+ 2X/3 [a]Vl + g (*/)(203,3 -0,,1 -02,2) -4___ [*/]121" In this case the governing equations for U~, V~ obtained from equations (11) and (12*/)(02,21 + 03,13) + ~ [41. + 2A]VI,, + X/3 [2/z](V2,z + V3,3) + X/3 [2*/] Vl,l -t-122,2 -]-~r3,3 --(p)UI = O) ~,)U2,22 -4-(1£)(U2,11 + U2,33 ) + <]d, +/~)(U1,12 + U3,23 ) + ('0) 02,11 +3 02,22 + 02,33

  (i) if 0 < ~ < 1 then V~ [2~ + ,~] 6 ~ / 1 -V~--~exp(-alt) + 1 + V1 -~exp(-~2t)-2 ~ ) Iv', = 12 {4/, + 2hi L, V'I -~¢ ~ ' U, = 6 ( L , )-' x , , U2 = U3 = I/2 = 113 = 0 where 8{n} ~I a, = F-~p), + V1-~:), 8{7} tl "2 = l=-~p>, -Vl -O; 1 [2tz + a] ( 16 £[rl]{n}~ 8 F3 -V~¢ -1 {2/x + a------~ [A] + 3 12(p) ] L,"

  ) the terms involving the length parameter I and introducing similar boundary and initial conditions (in this case initial conditions for V/are not taken into account), we obtain the boundary-value problem formulated in the framework of the local macro-model but describing the same physical situations. It can be verified that the solution to this problem takes the form where [2/, + A] 3 II1 (1 -exp(-3"t)), {4t* + A} L~ G = 8(L,)-'x,, I-]2 = lib = V2 = 1/3 = O, + A} [A] + ~ [171y L--7" (17)

1 .

 1 Diagrams of micro-stresses Sii based on equations (16) and (17) for the length-scale model.

Fig. 2 .

 2 Fig. 2. Diagrams of micro-stresses S;~ for the local model,

Acknowledgement--This research was supported by the Polish Commiteee for Scientific Research (KBN Warsaw) under grant 3 3310 92 03.

(ii) if g = 1 then &l = C + E~ exp(-at) + E]t exp(-at), $22 = $33 = D + F~ exp(-at) +F]t exp(-at),