Josep Rifà
email: josep.rifa@autonoma.edu

Lorena Ronquillo
email: lorena.ronquillo@autonoma.edu

Construction of New Completely Regular Z 2 Z 4 -linear Codes from Old *

Keywords: additive codes, completely regular codes, perfect binary codes, Z2Z4-linear codes

A code C is said to be Z2Z4-additive if its coordinates can be partitioned into two subsets X and Y , in such a way that the punctured code of C obtained by removing the coordinates outside X (or, respectively, Y ) is a binary linear code (respectively, a quaternary linear code). The binary image of a Z2Z4-additive code, through the Gray map, is a Z2Z4-linear code, which is not always linear. Given a perfect Z2Z4-linear code, which is known to be completely regular, some constructions yielding new Z2Z4-linear codes are computed, and the completely regularity of the obtained codes is studied.

Introduction

Let Z α 2 and Z β 4 denote the set of all binary vectors of length α and the set of all quaternary vectors of length β, respectively. Let GF (q) be a Galois field with q elements, q being a power of some prime number, and let GR(q m ) be a Galois ring with cardinality q m , which comes from an extension of the ring Z q . The classical Hamming weight wt(v) of a vector v ∈ GF (q) n is the number of coordinates which are different from zero, and the Hamming distance d(u, v) between two vectors u, v ∈ GF (q) n denotes the weight of their difference.

Any non-empty subgroup C of Z α 2 ×Z β 4 is a Z 2 Z 4 -additive code. Let φ : Z 4 -→ Z 2 2 be the usual Gray map, that is, φ(0) = (0, 0), φ(1) = (0, 1), φ(2) = (1, 1), and φ(3) = (1, 0); and let Φ : Z α 2 × Z β 4 -→ Z n 2 , where n = α + 2β, be the extended Gray map (Id, φ) given by Φ(u 1 , . . . , u α |v 1 , . . . , v β ) = (u 1 , . . . , u α |φ(v 1 ), . . . , φ(v β )).

It is worth noting that the extended Gray map is an isometry which transforms Lee distances defined in a Z 2 Z 4 -additive code C onto Hamming distances defined in the binary code C = Φ(C), where the length of C is n = α + 2β.

A Z 2 Z 4 -additive code C is isomorphic to an abelian structure like Z γ 2 × Z δ 4 . Thus C has |C| = 2 γ 4 δ codewords, where 2 γ+δ of them are of order two. We call such code C a Z 2 Z 4 -additive code of type (α, β; γ, δ), and its binary image C = Φ(C) is a Z 2 Z 4 -linear code of type (α, β; γ, δ), which may not be linear.

Given a Z 2 Z 4 -additive code C, its dual is also a Z 2 Z 4 -additive code, denoted by C ⊥ , and defined as the set of vectors in Z α 2 × Z β 4 which are orthogonal to every codeword in C. We use the following definition (see [START_REF] Borges | Z2Z4-linear codes: generator matrices and duality[END_REF]) of inner product in Z α 2 × Z β 4 :

u, v = 2 α i=1 u i v i + α+β j=α+1 u j v j ∈ Z 4 , (1)
where u, v ∈ Z α 2 × Z β 4 and computations are made considering zeros and ones in the α binary coordinates as quaternary zeros and ones, respectively.

After the Gray map over C ⊥ , the binary code

C ⊥ = Φ(C ⊥ ), of length n = α + 2β, is called the Z 2 Z 4 -dual code of C.
Two Z 2 Z 4 -additive codes C 1 and C 2 of the same length are monomially equivalent if one can be obtained from the other by permuting the coordinates and, if necessary, changing the sign of certain quaternary coordinates.

Although a Z 2 Z 4 -additive code C of type (α, β; γ ′ , δ ′ ) may not have a basis, every codeword is uniquely expressible in the form

γ ′ i=1 λ i u (i) + γ ′ +δ ′ j=γ ′ +1 µ j v (j) , (2)
where

λ i ∈ Z 2 for 1 ≤ i ≤ γ ′ , µ j ∈ Z 4 for γ ′ + 1 ≤ j ≤ γ ′ + δ ′ and u (i) , v (j) are vectors in Z α 2 × Z β 4
of order two and order four, respectively (see [START_REF] Borges | Z2Z4-linear codes: generator matrices and duality[END_REF]). Vectors u (i) , v (j) give us a generator matrix G which can be written as follows:

G = B 2 Q 2 B 1 Q 1 , (3)
where B 2 and B 1 are binary matrices of size (γ ′ ×α) and (δ ′ ×α), respectively; Q 2 is a (γ ′ × β)-matrix whose elements are in {0, 2} ⊂ Z 4 and Q 1 is a quaternary (δ ′ × β)-matrix with row vectors of order four. Refer to [START_REF] Borges | Z2Z4-linear codes: generator matrices and duality[END_REF] for further details.

Depending on the computations we are doing it might be convenient to take a different representation for the above matrix.

Let D be a matrix written as follows:

D = B ′ 2 Q ′ 2 B ′ 1 Q ′ 1 , (4)
where B ′ 2 and B ′ 1 are binary matrices of size (ψ ×α) and (θ ×α), respectively, in which binary zeros and ones have been represented as quaternary zeros and twos; Q ′ 2 is a (ψ × β)-matrix whose elements are in {0, 2} ⊂ Z 4 and Q ′ 1 is a quaternary (θ × β)-matrix with row vectors of order four. We will refer to this matrix as a matrix of type (α, β; ψ, θ). Note that matrix in (4) represents a group homomorphism

D : Z ψ 2 × Z θ 4 -→ Z α 2 × Z β 4
, where for any v = (λ 1 , . . . , λ ψ , µ ψ+1 , . . . , µ ψ+θ ) ∈ Z ψ 2 × Z θ 4 , vector D(v) is computed as in [START_REF] Bassalygo | Remark on uniformly packed codes[END_REF]. From a matrix point of view, this can be done by computing vD, which yields a vector in Z α 2 × Z β 4 whose binary zeros and ones are represented, respectively, by quaternary zeros and twos.

Given a code C ⊥ of type (α, β; γ, δ), corresponding to the Z 2 Z 4 -additive dual code of C, its generator matrix H is a parity check matrix of C. Thus all codewords in C are orthogonal to every row vector in H through the inner product in [START_REF] Bassalygo | Uniformly packed codes[END_REF]. As for the usual case of codes over a finite field, given a vector v ∈ Z α 2 × Z β 4 , we can refer to the syndrome S v as the vector in Z γ 2 × Z δ 4 obtained by computing vH T , where H T is the matrix obtained by writing the rows in H as column vectors, and it is written as in (4).

Let C be a binary code. The distance of any vector

v ∈ Z n 2 to C is d(v, C) = min x∈C {d(v, x)}, and the covering radius of C is ρ = max v∈Z n 2 {d(v, C)}. Let us also define C(i) = {v ∈ Z n 2 : d(v, C) = i}, i = {1, .
. . , ρ}. We say that two vectors u and v are neighbours if d(u, v) = 1. In this paper, the definition of completely regularity given in [START_REF] Neumaier | Completely regular codes[END_REF] will be used. Definition 1. A code C with covering radius ρ is completely regular, if for all l ≥ 0 every vector x ∈ C(l) has the same number a l of neighbours in C(l), the same number b l of neighbours in C(l + 1), and the same number c l of neighbours in C(l-1). Moreover, note that a l +b l +c l = n and c 0 = b ρ = 0. The intersection array of C is then (b 0 , . . . , b ρ-1 ; c 1 , . . . , c ρ ).

It is well known that any linear completely regular code C implies the existence of a coset distance-regular graph. Binary images of Z 2 Z 4 -additive codes are nonlinear binary codes, but we can also construct the quotient graph which is distance-regular when the initial Z 2 Z 4 -linear codes are completely regular.

Given a code C, its external distance is the number of nonzero terms in the MacWilliams transform. From any Z 2 Z 4 -additive code C we can construct the corresponding Z 2 Z 4 -additive dual code C ⊥ and, in this case, their respective weight enumerator polynomials are related by the MacWilliams identity [START_REF] Delsarte | An algebraic approach to the association schemes of coding theory[END_REF]. Thus in the specific cases of Z 4 or Z 2 Z 4 additivity, we can also refer to [START_REF] Rifà | Translation Invariant Propelinear Codes[END_REF] and, hence the external distance of C = Φ(C) coincides with the number of nonzero weights occurring in the distance distribution of

C ⊥ = Φ(C ⊥ ).
Throughout this paper we will use the concept of uniformly packed code (usually called "in the wide sense") defined in [START_REF] Bassalygo | Uniformly packed codes[END_REF].

Definition 2 ([1]

). Let C be a binary code of length n and covering radius ρ. We say that C is uniformly packed "in the wide sense" if there exist rational numbers τ 0 , . . . , τ ρ such that, for any

v ∈ Z n 2 , ρ k=0 τ k A k (v) = 1,
where A k (v) is the number of codewords in C at distance k from v.

In the following Proposition we summarize a few well known results.

Proposition 1. Let C be a code (not necessarily linear) with error correcting capability e ≥ 1, covering radius ρ and external distance s. Then:

(i) [8] ρ ≤ s.
(ii) [START_REF] Bassalygo | Remark on uniformly packed codes[END_REF] ρ = s if and only if code C is uniformly packed (in the wide sense).

(iii) [START_REF] Brouwer | Distance-Regular Graphs[END_REF] If C is completely regular, then it is uniformly packed (in the wide sense). (iv) [START_REF] Delsarte | An algebraic approach to the association schemes of coding theory[END_REF][START_REF] Rifà | Translation Invariant Propelinear Codes[END_REF] If C is a Z 2 Z 4 -linear code, then the external distance s of C coincides with the number of nonzero weights occurring in the

Z 2 Z 4 -dual code of C. (v) [7] If 2e + 1 ≥ 2s -1, then C is a completely regular code. A perfect Z 2 Z 4 -linear code C = Φ(C) of length n = 2 m
-1 is a binary perfect code, that is, a binary code of minimum distance 3, where all vectors in Z n 2 are within distance one from a unique codeword.

It is well known [START_REF] Borges | Rifà: A characterization of 1-perfect additive codes[END_REF] that for any m ≥ 2 and each δ ∈ {0, . . . , ⌊ m 2 ⌋} there exists a perfect Z 2 Z 4 -linear code C of binary length n = 2 m -1, such that its

Z 2 Z 4 -dual code is of type (α, β; γ, δ), where α = 2 m-δ -1, β = 2 m-1 -2 m-δ-1
and m = γ + 2δ (recall that the binary length can be computed as n = α + 2β). This allows us to write, for a given value of δ, the parity check matrix H of any Z 2 Z 4 -additive code C corresponding to a perfect Z 2 Z 4 -linear code C = Φ(C). This parity check matrix can be expressed in form [START_REF] Borges | Z2Z4-linear codes: generator matrices and duality[END_REF] where the first α columns are all possible nonzero vectors in Z γ+δ 2 , and the last β columns are all possible order four vectors in Z γ 2 × Z δ 4 , up to scalar multiples, where binary zeros and ones are respectively represented as quaternary zeros and twos (see [START_REF] Borges | The rank and kernel of extended 1-perfect Z4-linear codes and additive non-Z4-linear codes[END_REF]).

In this paper we construct some completely regular codes by modifying, in various ways, a perfect Z 2 Z 4 -linear code C. In Section 2 we extend, puncture and also shorten code C, while in Section 3 we describe a new method to obtain nonlinear completely regular codes in GF (4) by lifting perfect Z 2 Z 4 -linear codes defined over GF [START_REF] Bassalygo | Remark on uniformly packed codes[END_REF]. The same approach is also taken in Section 4 to obtain uniformly packed codes in GF (4) by lifting extended perfect Z 2 Z 4 -linear codes over GF (2).

2 Extending, puncturing and shortening perfect Z 2 Z 4 -linear codes

Let C be a Z 2 Z 4 -additive code, with parity check matrix H, such that its binary image C = Φ(C) is a perfect Z 2 Z 4 -linear code of length n = 2 m -1, for m ≥ 3, and its Z 2 Z 4 -additive dual code C ⊥ is of type (α, β; γ, δ).

Let Ĉ be the extended code of C, that is, the Z 2 Z 4 -additive code whose parity check matrix Ĥ is obtained from H by first adding a zero column at the beginning of the binary part, and then adding the all-twos row vector, where the twos in the first α + 1 coordinates are representing binary ones. The binary code Ĉ = Φ( Ĉ) has length n = 2 m and every codeword has even weight.

Let Ċ be the punctured code of C, that is, the Z 2 Z 4 -additive code obtained by removing the i-th column of the generator matrix of C or, equivalently, obtained by deleting the i-th coordinate from each codeword in C. Each time a coordinate is deleted in C, the length of the corresponding binary code Ċ = Φ( Ċ) drops by 1 when 1 ≤ i ≤ α, or by 2 when α + 1 ≤ i ≤ α + β.

Let C * be the shortened code of C, that is, the Z 2 Z 4 -additive code obtained by taking all codewords in C having a zero as their i-th component, where 1 ≤ i ≤ α + β, and subsequently delete the i-th component from these codewords. The procedure of shortening a code obviously decreases its length, thus the binary code

C * = Φ(C * ) is of length n = 2 m -2 if 1 ≤ i ≤ α, and of length n = 2 m -3 if α + 1 ≤ i ≤ α + β.
The following propositions can be proved. Using properties (ii) and (v) from Proposition 1, the following proposition is straightforward. Proposition 4. Uniformly packed codes.

(i) Code Ĉ is uniformly packed. (ii) Code Ċ is uniformly packed. (iii) Code C * is uniformly packed only when the i-th shortened coordinate is 1 ≤ i ≤ α.
Finally, we can state the following proposition which is given without proof.

Proposition 5. Completely regular codes.

(i) Code Ĉ is a completely regular code with intersection array (2 m , 2 m -1; 1, 2 m ).

(ii) Code Ċ is a completely regular code with intersection array (2 m -2; 2) when the i-th punctured coordinate is 1 ≤ i ≤ α, and with intersection array

(2 m -4; 4) when α + 1 ≤ i ≤ α + β. (iii) Code C * is completely regular with intersection array (2 m -2, 1; 1, 2 m -2),
when the i-th shortened coordinate is 1 ≤ i ≤ α.

Remark 1. Notice that in those cases where the extended (shortened) codes are completely regular codes, the intersection array coincides, respectively, with that of the extended (shortened) Hamming binary code of the same length. As for the punctured code, its intersection array coincides with that of the punctured Hamming binary code of the same length when we are puncturing a coordinate within the first α binary coordinates, while it coincides with the intersection array of the 2-punctured Hamming binary code when we are puncturing a quaternary coordinate. In other words, we have constructed new completely regular codes, but with the same intersection array as the corresponding linear codes.

3 Lifting perfect Z 2 Z 4 -linear codes

Lifting of binary perfect linear codes was previously studied in [START_REF] Rifà | On lifting perfect codes[END_REF]. In this article we now introduce and study lifted perfect Z 2 Z 4 -linear codes.

As in Section 2, let C be a Z 2 Z 4 -additive code, with parity check matrix H, such that its binary image C = Φ(C) is a perfect Z 2 Z 4 -linear code of length n = 2 m -1, for m ≥ 3, and its Z 2 Z 4 -additive dual code C ⊥ is of type (α, β; γ, δ). We define the code C r over an extension GF (2 r ) of the binary finite field as

C r = C + ξC + ξ 2 C + . . . + ξ r-1 C,
where ξ i C denotes the result of multiplying every codeword in C by ξ i , and ξ is a primitive element in GF (2 r ). Let C r = Φ -1 (C r ) be the corresponding code in GF (2 r ) α × GR(4 r ) β . Code C r has H r as parity check matrix, where H r can be seen as an α × β × r matrix in which each submatrix (h ijk ) for a fixed k ∈ {1, . . . , r} coincides with the parity check matrix H of C. We will say that code C r is obtained by lifting the perfect Z 2 Z 4 -linear code C. Note that code C r may not be linear over GF (2 r ).

Any vector v ∈ C r can be seen as a matrix of size r × (α + β) by representing every coordinate in GF (2 r ) from the first α coordinates and every coordinate in GR(4 r ) from the last β coordinates, as a column vector in GF (2) r and in GR(4) r , respectively. This matrix therefore contains r row vectors in Z α 2 × Z β 4 . Moreover, we will represent the binary zeros and ones in the first α coordinates of every such rows as quaternary zeros and twos, respectively. Let this matrix representation of v be denoted by [v]. Proposition 6. Let C be a Z 2 Z 4 -additive code having a perfect Z 2 Z 4 -linear code binary image C = Φ(C) of length n = 2 m -1, for m ≥ 3; let C r be the code defined over GF (2 r ) α × GR(4 r ) β , obtained by lifting C, and let C r be the corresponding binary image of C r . For r = 2, the covering radius ρ of C r is 2.

Proof. Let v be any vector in GF (4) α × GR(16) β and not in C 2 , and let S v = vH T ∈ GF (4) γ × GR(16) δ be its syndrome. We can then write v = w + e, where w ∈ C 2 and e is a vector in GF (4) α × GR(16) β of minimum weight such that its syndrome is S e = S v . Therefore, the distance from v to C 2 is the weight of vector e, that is d(v, C 2 ) = wt(e).

Let us represent vector e as a matrix [e] of two row vectors in Z α 2 × Z β 4 . Since the covering radius of any perfect Z 2 Z 4 -linear code is 1, we have that the Lee weight of every row in [e] is at most 1. Therefore, we know that wt(e) ≤ 2, thus d(v, C 2 ) ≤ 2 and the covering radius of C 2 is ρ ≤ 2.

However, since C 2 is not a perfect code, we know that ρ cannot be 1. Hence the covering radius ρ of C 2 is 2.

Lemma 1. Let v be a vector in GF (4) α × GR(16) β , and let (x, y) T be a column vector in the matrix representation [v] of v, where x, y ∈ Z 4 . Then,

wt L (v) = α+β i=1 wt L (([v] 1,i , [v] 2,i ) T ),
where wt L ((x, y) T ) = max{wt L (x), wt L (y), wt L (xy)}.

Recall, when computing the Lee weight defined in Lemma 1, that the first α column vectors in [v] contain quaternary zeros and twos which are representing binary zeros and ones, respectively. The lemma above let us go further on the Lee weights of codewords in the dual code of C 2 .

Proposition 7. The nonzero Lee weights of the codewords in C ⊥

2 are 2 m-1 and 2 m-1 + 2 m-2 , where m = γ + 2δ.

From the above proposition it is clear that the external distance of C 2 is s = 2, which leads us to the following proposition. Proposition 8. Code C 2 is uniformly packed.

We include a technical lemma, which will help us later in Theorem 1.

Lemma 2. Let C 2 be the code defined over GF (4) α × GR(16) β , and let C 2 be the corresponding binary image which we know has covering radius 2. Let µ i be the number of cosets in C 2 (i), for 0 ≤ i ≤ ρ. Then, µ 1 = 3n; µ 2 = n(n -1), where n = α + 2β . Now, we can prove the main theorem.

Theorem 1. Let C be a Z 2 Z 4 -additive code, with parity check matrix H, such that its binary image C = Φ(C) is a perfect Z 2 Z 4 -linear code of length n = 2 m -1, for m ≥ 3, and its Z 2 Z 4 -additive dual code C ⊥ is of type (α, β; γ, δ). Let C 2 be the code defined over GF (4) α × GR(16) β , obtained by lifting C, and let C 2 be the corresponding binary image of C 2 . Code C 2 is a completely regular code with intersection array (3n, 2(n -1); 1, 6).

Proposition 2 .Proposition 3 .

 23 External distance. (i) Code Ĉ has external distance s = 2. (ii) Code Ċ has external distance s = 1. (iii) Code C * has external distance s = 2 when the i-th shortened coordinate is 1 ≤ i ≤ α, and s = 3 when the i-th shortened coordinate is α+1 ≤ i ≤ α+β. Covering radius. (i) The covering radius of Ĉ is ρ = 2. (ii) The covering radius of Ċ is ρ = 1. (iii) The covering radius of C * is ρ = 2.

* This work was partially supported by the Spanish MICINN Grants MTM2009-08435, PCI2006-A7-0616, and also by the CUR del DIUE de la Generalitat de Catalunya and the European Social Fund with Grants 2009SGR1224 and FI-DGR.

Proof. Code C 2 is a projective code, hence e ≥ 1. From Proposition 6 and Proposition 8 it is easy to see that code C 2 satisfies item (v) in Proposition 1 and so, code C 2 is completely regular.

For the computation of the intersection array (b 0 , b 1 ; c 1 , c 2 ), as we know the code is completely regular, we easily have that c 1 = 1 and c 2 = 6;

In general, the lifted code C r , when r > 2, is neither completely regular, nor uniformly packed. We can take a specific example for the case r = 3 and show that the corresponding code C 3 has ρ = s. The same happens, in general, for r ≥ 3.

Let C be a Z 2 Z 4 -additive code whose Z 2 Z 4 -additive dual code is of type (3, 6; 0, 2), and having the following parity check matrix:

Code C corresponds to a perfect Z 2 Z 4 -linear code C = Φ(C) of length n = 15. Let C r be the code defined over GF (2 r ) 3 × GR(4 r ) 6 , obtained by lifting C, for r = 3, and let C 3 be its binary image.

Let v be any vector in GF (8) 3 × GR(64) 6 , and not in C 3 . As in the proof of Proposition 6, we have d(v, C 3 ) = wt L (e), where v = w + e, w ∈ C 3 , and e is a vector in GF (8) 3 × GR(64) 6 of minimum weight in the coset of v.

Since the covering radius of any perfect Z 2 Z 4 -linear code is 1, it follows that the Lee weight of every row in the matrix representation [e] of vector e is at most 1. Therefore, the covering radius of C 3 is ρ ≤ 3.

We proceed to compute the external distance of C 3 . Let w be a codeword in C ⊥ 3 , and let [w] be its matrix representation.

Let w 1 = (220|123011), w 2 = (202|011112), w 3 = w 1 + w 2 = (022|130123) and w 4 = 3(w 1 + w 2 ) = (022|310321) be four vectors generated by matrix in [START_REF] Borges | Rifà: A characterization of 1-perfect additive codes[END_REF]. Note that the Lee weight of w is 13 when [w] has vectors w 1 , w 2 and w 3 as rows; it is 15 when [w] has vectors w 1 , w 2 and w 4 ; it is 8 when [w] has vector w 1 and two zero rows and, finally, it is 12 when [w] has vector w 1 , w 2 and one zero row. Thus the external distance s of C 3 is s ≥ 4 and, therefore, C 3 is not a uniformly packed code.

Remark 2.

As already mentioned at the beginning of this article, the quotient graph corresponding to the binary image of a Z 2 Z 4 -additive code is distanceregular when that Z 2 Z 4 -linear code is completely regular.

From the completely regular codes C 2 described in this paper, we obtain distance-regular graphs with classical parameters. The bilinear forms graphs [6, Sec. 9.5] have the same parameters and they were not described, until recently in [START_REF] Rifà | On lifting perfect codes[END_REF], as coset graphs of completely regular linear codes. The interesting point here is that Theorem 1 is giving now the same description but using nonlinear codes.

4 Lifting extended perfect Z 2 Z 4 -linear codes Let C be a Z 2 Z 4 -additive code having a perfect Z 2 Z 4 -linear code binary image C = Φ(C) of length n = 2 m -1, for m ≥ 3, and let Ĉ be the extended code of C, studied in Section 2. Let Ĉr , for r = 2, be the code defined over GF (4) α+1 × GR(16) β obtained by lifting code Ĉ. Code Ĉ2 has Ĥ2 as parity check matrix, where Ĥ2 is an (α + 1) × β × 2 matrix in which each submatrix (h ijk ) for a fixed k ∈ {1, 2} coincides with the parity check matrix of Ĉ, that is Ĥ. Let Ĉ2 be the corresponding binary image of Ĉ2 .

The following proposition can be proved.

Proposition 9. Code Ĉ2 has covering radius ρ = 3 and it is uniformly packed, but not completely regular.


Thomas Blanchet
email: thomas.blanchet@wanadoo.fr

How useful is the Berlin Centre of Competence for Water for the urban management of water supply? A contribution to an inter-organisational analysis in the field of sustainable innovation

Keywords: Knowledge management, public private partnership, sustainable development, urban water, water privatization

come

How useful is the Berlin Centre of Competence for Water for the urban management of water supply? A contribution to an inter-organisational analysis in the field of sustainable innovation

Thomas Blanchet

INTRODUCTION

This paper represents a first outline of my Phd project, which aims to focus on knowledge management in the water sector in a post-privatization context. Over the past two decades, the water sector has been experiencing significant changes. While the majority of water utilities have remained public worldwide, private operators in such services have gained in influence. Considering the private sector participation as the only chance for a great number of countries to obtain access to drinkable water (Camdessus, 2003), the International Monetary Fund and the World Bank encouraged such involvement (Goldman, 2007). In Europe, several actors also supported the liberalization of public services. In order to bring more efficiency to urban services, the European Commission strongly favored this process in various sectors. While the water sector escaped deregulation, the introduction of competition through obligatory invitations to tender was discussed and the cooperation between local authorities and private operators through public-private partnerships encouraged. Taking into account this context, I aim to discuss in this project the following question:

How can cooperation between actors with potential divergent interests and conflicting logic contribute to the development of sustainable urban water management?

The traditional water management in Germany

To answer this question, I will focus on the development of public-private partnerships in Germany. In contrast to many other countries, Germany has a strong institutionalized water sector with a long and solid management tradition and high technical standards. In Germany the task of water supply, as well as several other activities, belongs to the field called "communal economy" (Kommunalwirtschaft). Through the self-administration principle (Selbstverwaltung) municipalities are free to decide how to manage this task. Municipalities generally created utilities in the form of a public company that integrates water distribution, energy and public transportation. More than 6000 German companies are responsible for water supply and sanitation. These utilities, called Stadtwerke, are under the strict control of the municipality. Even when a private operator buys into the capital of such companies, the majority of the shares generally remain in public hands (Drouet, 1988: 45). This organization is the consequence of a weak and decentralized state, which enabled the development of strong cities relying on considerable technical and financial resources. In contrast, French urban services, including water distribution, have followed another path. Instead of having created municipal-owned public companies regulated by legal principles, municipalities have "succeeded in co-operating, delegating and adjusting their actions by enlisting the aid of other partners" (Lorrain, 1992: 80). Private operators have become increasingly powerful, regrouping all competencies necessary to offer turnkey contracts, meaning "with architectural and technical design, realization of secured performance, operational experience and sometimes funding" (Camilleri, 2006: 30). The French context led to the foundation of two "national champions" which have become global players in the water services.

The case of the Berlin Water Utility

Following Reunification, the German model of water management has been subject to various pressures to change. In addition to the financial crisis of the municipalities, the cost of the Reunification and the effect of the "New Steering Model", local authorities had to cope with direct competitive pressures represented by private operators in search of new contracts. In such a context, the Berlin senate sold, in 1999, 49,9% of the shares of the Berlin Water Utility to a consortium composed of Veolia, RWE and Allianz (which in turn sold its shares to the other two partners in 2002).

Among German cities of over 500.000 inhabitants, Berlin is the only one in which a French private operator is involved in water management. The privatization took place under a cooperation model, that is, the private partner and the city cooperate in one organization, where both have a certain amount of shares. One peculiarity of this organizational model is that Berlin Water Utility still remains a company under public law. Through the creation of a holding approved by the Berlin Constitutional Court, the participation of private partners was made possible on October 21, 1999. This organizational change was the result of a long process involving strong competition between private operators, heated debates within the Berlin parliament, and hard bargaining between private operators and the senate of Berlin.

The main goal of this privatization was financial: the finance administration aimed to cope with the structural deficit of the municipality. However, to be politically and socially accepted, the privatization should also lead to social, environmental and technological improvements at the local level. In this context, the senate imposed several conditions for the transaction with the future private partner. One of these was the creation of the Berlin Centre of Competence for Water. This Center was founded in December 2001. Its purpose is to contribute to the development of techniques and solutions in water management through the cooperation between heterogeneous actors: the private multinational Veolia, the Berlin Water Utility, local universities, the Technology Foundation of Berlin (TSB), whose objective is to support innovation. Its official goal is the creation, receipt and dissemination of knowledge on the theme of water. The analysis of this center permits to raise several questions on the function of such an organization in a broader privatization context, and on its activity (scope and limits): Does this organization influence the policy making of the Berlin water management? To what extent does the cooperation between heterogeneous actors with conflicting interests lead to reach the common goal of the organization? What are the power relationships inside this center and how do they influence the performance of the organization?

METHOD

In this PhD project, I aim to adopt a qualitative analysis based on a single-case research design. I chose a qualitative analysis because of the complexity of water management and the multiplicity of actors involved in it. I opted for a case-study design for three reasons. Firstly, "case study analysis focuses on a small number of cases that are expected to provide insight into a causal relationship across a larger population of cases" (Gerring, 2007: 86). Secondly, I aim to analyze a phenomenon in its real life context where the relationship between the context and the phenomenon under scrutiny is not so clear. Thirdly, I have no control over the investigation. Because of the early stage of the work, this paper has rather an explorative character.

The creation of a research center regrouping heterogeneous actors and founded in a context of privatization represents a unique case in the water sector and is the reason why I chose to conduct a single case design. Moreover, this "critical" case will be able to show us how the cooperation between actors of different natures (public-private/industrial-non-profit/French-German) can affect the development of sustainable innovations at local and global levels. This case will be analyzed in the light of the following theoretical framework: first, the sociology of organization should help us to understand the functioning of this Center, the coordination between the actors involved in this cooperation, their vested interests and the resulting potential for creating and exporting new techniques and knowledge in the field of water management. Then, relying on the literature on innovation, I will try to highlight the scope and limits of such cooperation in the innovation process on sustainable water management. I aim to look especially at the tension and/or the complementarity between the various goals of a sustainable innovation (economic, social, ecological) (Tauchmann et al, 2006) and how these diverging goals are managed in a concrete case.

The empirical investigation should be divided into two mains steps. The first, considered as explorative, sets out to trace the constitution of water sectors in France and Germany: This historical description will be a necessary step to better understand the divergent logic and interests of the actors involved in this organization. In a second step, I will describe and analyze the activity of this Center and the context of its foundation. To that aim, I will rely on documents, such as printings from the municipal council, articles in local newspapers, public statements from different associations, articles published in specialized journals, and documents produced by the center. To complete this document analysis, I will conduct semi-structured interviews with various stakeholders, including researchers and engineers involved in this center, former and present members of the senate, managers of the Berlin Water Utility but also with members of private enterprises active in Germany. In addition to the first two methodological tools, I aim to conduct direct observations during public events organized in this centre. This will enable to complement our description of the center's activity and to obtain contacts with potential interview partners.

FIRST EXPLORATION

As the thesis is still at the beginning, the first results are rather assumptions made from a first work of exploration. I will first outline the activity of this Centre and then develop some ideas concerning the scope and limits of such a cooperation.

The activity of the Competence Center

The activity of the centre can be defined along the following lines: Identifying important points of research and development, contributing to the transfer of information and technology in the field of sustainable water management, publicizing results of the research conducted in the centre, offering training, organizing conferences and seminars of experts, creating a network of small and medium-sized businesses through various projects. This Centre aims therefore to become an important intermediary between various actors involved in water management. The various projects carried out by the organization are divided into different areas: information technology in water economics, sustainable management of water resources and technological innovation in the field of waste-water treatment. In total, 30 projects are jointly run by employees of Veolia and the Berlin Water Utility. The purpose of these projects is to contribute to cost reductions and quality improvements, which are two essential issues in water management (KWB, 2001: 4). To that aim, Veolia Water invested 50 million Euros for a period of 10 years following the foundation of the Centre (KWB, 2001: 23).

In the Center, projects are usually conducted by teams made up of engineers and technicians from the Berlin Water Utility and from Veolia. The results of such projects are published by the center, such as through the Schrifftenreihe Kompetenzzentrum Wasserberlin (KWB, 2004: 21;KWB, 2005 : 24). In addition, various symposiums and work groups are organized. For instance, the « Berlin Water Workshop », whose goal is "to promote the information exchange among water experts from the Berlin and Brandenburg region" (KWB, 2004: 26), and has taken place every two months since fall 2004. Finally, the Center offers training possibilities for managers, scholarships for doctoral students, and set up in cooperation with the Technical University the Chair of water management (KWB, 2004: 25-26). Such practices are supposed to contribute to knowledge development through cooperation between various actors and to foster innovations.

The cooperation: scope and limits

At first sight, it could be argued that such a cooperation, between a private multinational, whose main goal is to make a profit, and a municipality, whose purpose should officially be oriented toward the well-being of the population, can impede the activity of such an organization. Moreover, these actors have evolved in distinct national systems making their cooperation potentially more difficult because of the weight of tradition in the development of solutions in the field of water management.

However, the local and international natures of the both shareholders can also complement each other and thus can contribute to the development of innovations. On the one hand, such cooperation enables the implementation of new techniques for water management in Berlin. For instance, the project Mobile IT Service Field Operations introduced a technology enabling employees, during their field work, to have real-time access to all the data they need to solve the problem. This project was developed in France and then exported to Berlin (KWB, 2001: 18). On the other hand, some techniques developed in Berlin could be exported abroad through the commitment of the international player. In such cases, the locality of Berlin can be used as a local laboratory for developing innovations which can be then exported to other cities. For instance, the NASRI project (model for optimizing filtration devices from banks) was developed in Berlin and then implemented in India (KWB, 2005: 13).

This center seems to be strategic for both partners. Veolia could reinforce its legitimacy at the local level by taking part in the economical development of the Region, investing in the research, supporting a network with regional Small and Medium Businesses and making Berlin an international place for water. Therefore, with such an investment, Veolia filled the expectation of local decision-makers and therefore reinforced its contacts at the local level, which in turn could be a strategic advantage in case of a new bid in Berlin. In addition, this Center contributes to the strategy of the Berlin Water Utility to become an international player in the water sector, which seems to be one of the main interests of the local decision makers. Furthermore, such international development does not challenge the respect of public value, as could be the case with investment businesses (see: political debates concerning the activity of Berlin Water International).

CONCLUSION

Despite this first positive description, a deeper field study is needed in order to assess the real effects of this technical-scientific cooperation on both, the various stakeholders involved in the project and on the general development of sustainable innovations in the field of water services. Potential conflicts and vested interests have thus far remained unexplored and still have to be analyzed. In addition to the different logic of these actors, this cooperation can lead to increase the dependence of Berlin toward the private operator. The investment in the cooperation was indeed made for 10 years and it is still unclear how the Center will pursue its activity after the end of the term and how such decisions can influence the cooperation between actors. Moreover, despite the legitimacy that Veolia could gain in investing in this center, a citizen initiative, called "Berlin Roundtable" started in 2006 to struggle against the privatization of the water services in Berlin following a water price rise. How such conflicts can influence the activity of the center remains unclear. Such problems will have to be unraveled during the field research which should take place in 2009.

• Camdessus M., 2003, « Rapport du groupe de travail international sur le financement des infrastructures mondiales d'accès à l'eau », Paris. • Camilleri C. ( 2006), L'eau, un domaine d'excellence pour les entreprises françaises, problèmes économiques, 2.900, p. 29-32. • Drouet D., (1988). L'industrie de l'eau dans le monde, presse de l'E.N.P.C., Paris.

• Gerring, J. 2007. Case study research: principles and practices. Cambridge Univ Pr.

• Goldman, M. ( 2007): How "water for all!" policy became hegemonic: The power of the World Bank and its transnational policy networks, Geoforum, 38 (5), pp. 786-800. • KWB (Kompetenzzentrum Wasser Berlin), (2001). Tätigkeitsbericht des Kompetenzzentums Wasser Berlin. Berlin • KWB, (2002). Tätigkeitsbericht des Kompetenzzentums Wasser Berlin. Berlin.

• KWB, (2004). Tätigkeitsbericht des Kompetenzzentums Wasser Berlin. Berlin.

• KWB, (2005). Tätigkeitsbericht des Kompetenzzentums Wasser Berlin. Berlin.

• Lorrain, D. 1992. "The French model of urban services." West European Politics 15(2): 77-92.

• Tauchmann, H., Hafkesbrink, J., Nisipeanu, P., Thomzik, M., Bäumer, A., Brauer, A., Clausen, H., Drouet, D., Engel, D., Körkemeyer, K., Rothgang, M., Schroll, M., (2006) 11-Blanchet-WWW- YES-2008YES- -edited-2011-06-30.doc -06-30.doc

8th

  World Wide Workshop for Young Environmental Scientists WWW-YES-2008: Urban waters: resource or risks? 13-16 May 2008 4 / 6

  . Innovationen für 8th World Wide Workshop for Young Environmental Scientists WWW-YES-2008: Urban waters: resource or risks? 13-16 May 2008 6 / 6 eine nachhaltige Wasserwirtschaft, Einflussfaktoren und Handlungsbedarf. Physica Verlag, Heidelberg.


L.M.S C M S Silva
email: cmsilva@civil.ist.utl.pt

Castro

L M S S Castro

Continuum damage models with non-conventional finite element formulations

Keywords: hybrid and mixed finite element formulations, non-local damage models, Legendre polynomials

In recent years, some research effort has been devoted to the development of nonconventional finite element models for the analysis of concrete structures. These models use continuum damage mechanics to represent the physically non-linear behavior of this quasi-brittle material. Two alternative approaches proved to be robust and computationally competitive when compared with the classical displacement finite element implementations. The first corresponds to the hybrid-mixed stress model where both the effective stress and the the displacement fields are independently modeled in the domain of each finite element and the displacements are approximated along the static boundary, which is considered to include the interelement edges. The second approach corresponds to a hybrid displacement model. In this case, the displacements in the domain of each element and the tractions along the kinematic boundary are independently approximated. Since it is a displacement model, the inter-element boundaries are now included in the kinematic boundary. In both models, complete sets of orthonormal Legendre polynomials are used to define all approximations required so very effective p-refinement procedures can be implemented. This paper illustrates the numerical performance of these two alternative approaches and compares their efficiency and accuracy with the classical finite element models. For this purpose, a set of numerical tests is presented and discussed.

A c c e p t e d m a n u s c r i p t 1 Introduction

The main objective of the research work reported in this paper is the development of robust and computational competitive numerical techniques for the physically non-linear analysis of concrete structures.

The non-conventional hybrid and mixed finite element formulations used in this work were first developed by Freitas et al [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF] during the nineties. Three main classes are defined, namely the hybrid-mixed, the hybrid and the hybrid-Trefftz formulations. Two models are derived for each formulation, the displacement and the stress models [START_REF] Wang | Elastoplastic structural analysis with hybrid stress elements[END_REF][START_REF] Cismaşiu | The hybrid-Trefftz displacement element for static and dynamic structural analysis problems[END_REF]. All formulations evolve directly from the first principles of mechanics, in particular equilibrium, compatibility and constitutive relations. What distinguishes the three types of formulations is the set of constraints enforced, a priori, on the domain approximations.

In recent years, some of these non-conventional finite element formulations have been extended to non-linear analysis using isotropic damage models [START_REF] Silva | Hybrid-mixed stress model for the nonlinear analysis of concrete structures[END_REF][START_REF] Silva | Hybrid-mixed stress formulation using continuum damage models[END_REF][START_REF] Silva | Hybrid-displacement (trefftz) formulation for softening materials[END_REF][START_REF] Silva | Modelos de Dano em Elementos Finitos Híbridos e Mistos[END_REF][START_REF] Silva | Hybrid and mixed finite element formulations for softening materials[END_REF].

In [START_REF] Silva | Hybrid-mixed stress model for the nonlinear analysis of concrete structures[END_REF][START_REF] Silva | Hybrid-mixed stress model for the nonlinear analysis of concrete structures[END_REF] the hybrid-mixed stress model based on the use of orthonormal Legendre polynomials [START_REF] Pereira | Numerical implementation of a hybridmixed finite element model for reissner-mindlin plates[END_REF] is chosen. The stress and the displacement fields in the domain of each element and the displacements on the static boundary are independently approximated. None of the fundamental relations is enforced a priori and all field equations are enforced in a weighted residual form, ensuring that the discrete numerical model embodies all the relevant properties of the continuum it represents. The Mazars' isotropic model [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non lineaire et à la rupture du béton de structure[END_REF] is adopted and a non-local integral formulation where the damage variable is taken as the nonlocal variable is considered.

In [START_REF] Silva | Hybrid-mixed stress formulation with continuum damage models[END_REF][START_REF] Silva | Hybrid-mixed stress formulation using continuum damage models[END_REF] an improved hybrid-mixed stress model is presented and discussed. The approximation of the stress field in the domain is here replaced by the approximation of the effective stress field. The isotropic damage models presented by Comi and Perego [START_REF] Comi | Nonlocal aspects of nonlocal damage analyses of concrete strucutres[END_REF][START_REF] Comi | A bi-dissapative damage model for concrete with applications to dam engineering[END_REF] are now adopted using a non-local integral model. An alternative technique based on the definition of an explicit enhanced gradient model has also been tested [START_REF] Silva | Modelos de Dano em Elementos Finitos Híbridos e Mistos[END_REF].

The use of hybrid-Trefftz displacement formulations, where the displacements in the domain of each element and the stress field on the kinematic boundary are independently approximated, is reported in [START_REF] Silva | Non-conventional finite element models using continuum damage mechanics[END_REF][START_REF] Silva | Hybrid-displacement (trefftz) formulation for softening materials[END_REF]. The main feature of these models is that the functions used to approximate the displacements are derived from bi-harmonic displacement potentials that solve the Navier equations for a homogeneous elastic material [START_REF] Freitas | Comparative analysis of hybridtrefftz stress and displacements elements[END_REF][START_REF] Cismaşiu | The hybrid-Trefftz displacement element for static and dynamic structural analysis problems[END_REF]. Due to the strategy used in their definition, only while concrete presents a linear elastic behavior the model may be considered as a pure hybrid-Trefftz formulation. When this behavior is no longer valid, it becomes an hybrid-displacement model.

A c c e p t e d m a n u s c r i p t

The use of hybrid-displacement models with Legendre polynomials is reported in [START_REF] Silva | Modelos de Dano em Elementos Finitos Híbridos e Mistos[END_REF][START_REF] Silva | Nonlocal damage theory in hybriddisplacement formulations[END_REF]. An isotropic model with a non-local integral formulation is again adopted.

The main objective of this paper is to compare the numerical efficiency of the two most promising non-conventional finite element formulations implementing damage mechanics models, namely the hybrid-mixed stress model with the independent approximation of the effective stress field and the hybriddisplacement formulation with Legendre polynomials. A set of classical benchmark tests is chosen to illustrate the use of such models and to assess and compare their numerical performance. Comparisons are also made with solutions obtained with the classical finite element formulation, to prove that the proposed models are competitive.

The numerical models are both incremental and iterative and are solved with a modified Newton-Raphson method that uses the secant matrix. The models are implemented in a C program and all tests were run on a Pentium IV Machine, 2996 Mhz, running Linux.

This paper is organised as follows: the formulation of the problem and the adopted damage models are presented in Section 2 and 3. The non-conventional finite element formulations are described in Section 4. The numerical examples are shown in Section 5 and finally, Section 6 summarises the main conclusions and indicates future research work in this field.

Fundamental relations

Consider a domain V limited by the boundary Γ, referred to a cartesian coordinate system. The static boundary Γ σ (or Neumann boundary) and the kinematic boundary Γ u (or Dirichlet boundary) are complementary sub-regions of the boundary Γ, whereon traction-resultants and displacements are respectively prescribed.

The body under analysis is assumed to be homogeneous and isotropic. The model is geometrically linear and only static and monotonic loads are considered. No viscous, thermal or other non-mechanical dissipative effects are taken into account.

The fundamental equilibrium equations may be written in a matrix form as follows:

D σ + b = 0 in V, (1) N σ = t γ on Γ σ ,

A c c e p t e d m a n u s c r i p t

where D is the differential equilibrium operator. The matrix N contains the components of the unit outward normal vector to the static boundary Γ σ . The vector σ lists the independent components of the stress tensor. The vector b represents the body force vector in the domain V and t γ corresponds to the tractions vector on the static boundary Γ σ .

The compatibility equations may be written in the following format:

ε = D * u in V, (2) u = u on Γ u ,
where D * is the differential compatibility operator, adjoint of the differential equilibrium operator D since the model is geometrically linear. The vector ε collects the independent components of the strain tensor and the vector u lists the independent components of the displacement field. The vector u denotes the prescribed displacements on the boundary Γ u .

The constitutive relation depends on the damage model adopted, as detailed in Section 3.

Non-local Damage Model

The mechanical behavior of quasi-brittle materials such as concrete, is characterised by the development of micro-cracks and subsequent evolution to localised macro cracking. The Continuum Damage Mechanics models describe the evolution of the mechanical properties of the continuum as cracking develops. This type of constitutive models are able to describe, with a continuum approach, some of the material properties observed in experiments, such as global softening, stiffness degradation, anisotropy and development of inelastic deformations [START_REF] Simo | Strain-and stress-based continuum damage models -i. formulation[END_REF][START_REF] Mazars | Continuum damage theory -application to concrete[END_REF][START_REF] Laborderie | Phenomenes unilateraux dans un materiau endommageable: modelisation et application a l'analyse de structures en beton[END_REF][START_REF] Frémond | Damage, gradient of damage and principal of virtual power[END_REF][START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF].

In this paper, two different damage models are considered. The first one (Model 1) corresponds to a simple isotropic continuum damage model with only one damage variable used e.g. in the works of [START_REF] Lemaitre | A course on damage mechanics[END_REF][START_REF] Comi | Nonlocal aspects of nonlocal damage analyses of concrete strucutres[END_REF][START_REF] Comi | Symmetric and non-symmetric non-local damage formulations: an assessment of merits[END_REF]. The second model (Model 2) considers two independent damage variables and was introduced by Comi and Perego in [START_REF] Comi | A bi-dissapative damage model for concrete with applications to dam engineering[END_REF][START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF][START_REF] Comi | A non-local model with tension and compression damage mechanics[END_REF]. In both models, plastic strains are neglected. The two adopted models are presented summarily in Tables 1 and2.

The strain softening behavior is well known to produce strain localization with consequent dependence on the data of the finite element model, as for instance the dependence on the mesh and on the degrees of the approximation functions [START_REF] Bažant | Instability, ductility, and size effect in strain-softening concrete[END_REF][START_REF] Jirásek | Modelling of localized damage and fracture in quasibrittle materials[END_REF]. To overcome this problem, several regularisation techniques are proposed in the literature, in particular non-local integral [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF] and gradientenhanced damage formulations [START_REF] Peerlings | Gradient-enhanced damage for quasi-brittle materials[END_REF][START_REF] Comi | Computational modelling of gradient-enhanced damage in quasibrittle materials[END_REF]. Following Pijaudier-Cabot and Bažant

A c c e p t e d m

a n u s c r i p t [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF], the present work assumes a non-local integral analysis for both damage models. When the model with one damage variable is used, the strain energy release rate is adopted as the non-local variable. When using the second model, the strain invariants involved in the definition of the constitutive laws are defined as non-local variables.

As defined in [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF], a generic non-local variable v is computed considering the following weighted average over the whole domain:

v(x) = V W (x, s) v(s) ds ,
where W (x, s) is a weight function taken here as the normalised Gauss function:

W (x, s) = 1 W 0 (x) exp - x -s 2 2 l 2 , W 0 (x) = V exp - x -s 2 2 l 2 ds .
The length l in the previous equation is a geometric length, usually denoted as characteristic length. It works as a localisation limiter and regularises the mathematical problem. According to [START_REF] Bažant | Measurement of characteristic length of nonlocal continuum[END_REF], this length may also be interpreted as a material-dependent parameter related to the width of the fracture process zone. A normalised weight function is chosen because the non-local model should be able to reproduce correctly local uniform fields.

Model 1 has the limitation of considering the same behavior for the material in prevailing tension and compression states, which is not realistic for most of the materials. To overcome this limitation, it is assumed that damage may only appear and develop if the strain tensor trace is positive, tr ε > 0. The constitutive model with the referred assumption is suitable for studying structures subjected mainly to tension stresses and it is competitive due to its simplicity.

Model 2 requires more computational effort than Model 1. This happens because it is not possible to define explicitly the values of d t and d c in terms of the potentials f t and f c . On the other hand, the non-local integral model to be used [START_REF] Comi | A non-local model with tension and compression damage mechanics[END_REF] is more demanding from the computational point of view. The hybrid-mixed stress formulation adopted in this work was for the first time described in [START_REF] Silva | Hybrid-mixed stress formulation using continuum damage models[END_REF]. Compared to the original version of the hybrid-mixed stress formulation [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF] the particularity of the new model is that the approximation of the stress field σ is replaced by the approximation of the effective stress field σ, defined e.g. by Lemaitre in [START_REF] Lemaitre | A course on damage mechanics[END_REF]. The approximations may be expressed as:

σ = S v X in V, (3) u = U v q v in V, u = U γ q γ on Γ σ ,
where the matrices S v , U v and U γ collect the approximation functions and the vectors X, q v and q γ list the associated weights (generalised variables). Since the three fields are approximated independently, it is possible to adopt different degrees of approximation for each one.

Due to the properties presented by Legendre polynomials, it is possible to increase the degree of the approximations without having any problems in terms of numerical stability. This fact enables the implementation of highly efficient p-refinement procedures, as it is possible to define high degree approximations without deteriorating the condition number of the global governing system. In order to ensure the numerical stability of the hybrid-mixed stress model, it is necessary to ensure that the number of stress parameters (number of degrees of freedom associated with the approximation of the static fields) is greater than the number of displacement parameters (number of degrees of freedom associated with the approximation of the kinematic fields). To avoid spurious kinematic modes, it is also important to ensure that the degrees of the Legendre polynomials used in the displacement fields approximation is smaller than the polynomial degrees involved in the effective stress components approximation [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF].

While the concrete is linear elastic, the model proposed coincides with the one described in [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF]. When damage appears, the models are different since the effective stress and the stress field are no longer coincident. In the context of a non-linear analysis with softening, the main advantage of the proposed approach when compared to the one described by Freitas [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF] is that the effective stress field is directly related to the evolution of damage, since it is comparable to the to the strain field, while the stress field is not.

A c c e p t e d m a n u s c r i p t

The generalised strains, e, body forces, Q v , and tractions, Q γ are defined by

e = S t v dV, Q v = U t v b dV, Q γ = U t γ t γ dΓ σ , ( 4
)
in order to ensure the inner product invariance between the pairs of dual discrete variables ( X , e), (q v , Q v ), and (q γ , Q γ ) and the continuum fields they represent.

As demonstrated in [START_REF] Silva | Hybrid-mixed stress formulation with continuum damage models[END_REF][START_REF] Silva | Hybrid-mixed stress formulation using continuum damage models[END_REF][START_REF] Silva | Modelos de Dano em Elementos Finitos Híbridos e Mistos[END_REF], using the definition of the effective stress in the form σ = σ (1d) [START_REF] Lemaitre | A course on damage mechanics[END_REF] and enforcing the fundamental equations (Section 2) on average, in the sense of Galerkin, one obtains the following equilibrium equations for the discrete model:

(A t v -M v ) X = -Q v in V, (5)
(A t γ -M γ ) X = Q γ on Γ σ ,
where the matrices M v , M γ and A v , A γ are defined as follows:

M v = -(D * U v ) t S v d dV + (N * U v ) t S v d dΓ, M γ = U t γ (NS v ) d dΓ σ , (6)
A v = (DS v ) t U v dV, A γ = (NS v ) t U γ dΓ σ .
In a more complex damage model, as the Model 2 described in Section 3, a similar procedure may be followed to define matrices M v and M γ [START_REF] Silva | Modelos de Dano em Elementos Finitos Híbridos e Mistos[END_REF].

The compatibility condition in the discrete model (Equation ( 7)) may be obtained integrating by parts the average enforcement of the compatibility equation in the domain and then replacing in the resulting expression the approximations for the displacements (Equation ( 3)) [START_REF] Silva | Hybrid-mixed stress formulation with continuum damage models[END_REF][START_REF] Silva | Hybrid-mixed stress formulation using continuum damage models[END_REF]:

e = -A v q v + A γ q γ + e, with e = (NS v ) t u dΓ u . ( 7
)
The relation between the independent components of the effective stress tensor and the strain components can be expressed as [START_REF] Silva | Hybrid-mixed stress formulation using continuum damage models[END_REF]:

= F σ , ( 8
)
where F is the symmetric non-singular matrix of elastic constants characterizing a linear reciprocal elastic law.

Introducing the constitutive relation (Section 3) and the generalised strains (Equation ( 4)) in Equation ( 7), we obtain Equation ( 9) that englobes the

A c c e p t e d m a n u s c r i p t

compatibility and the constitutive relations of the discrete model:

X + A v q v -A γ q γ = e, with = S t v F S v dV. ( 9
)
Combining Equation ( 5) and ( 9), one obtains the following solving system for each finite element:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ A v -A γ (A t v -M v ) 0 0 -(A t γ -M γ ) 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ X q v q γ ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ e -Q v -Q γ ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ . ( 10
)
The governing system of the finite element mesh is assembled by direct allocation of the contribution of the elementary systems [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF].

According to Equation ( 10), when a damage mechanism in a finite element is active, the elementary governing system is nonsymmetric, since M v = 0 and M γ = 0. Nevertheless, the governing system remains very sparse and the non-local behavior due to damage appears in the "small dimension" of the system (Figure 1).

Because the non-conventional models adopted in this work use macro-element meshes, it is not possible to control the length of the nonlinear strain localization band through the finite element mesh, as usually happens in a traditional displacement formulation. Consequently, a more refined mesh must be chosen to implement the nonlocal integral model. In this work, the Lobatto points mesh is used for this purpose. Another choice would also be valid. The choice of the Lobatto points mesh is adopted because this mesh is already defined to ensure the accuracy of the numerical integrations. Since the hybrid-mixed stress model requires the knowledge of the damage evolution on the boundary, the Lobatto quadrature rule is used instead of the usual Gauss quadrature rule. In order to capture the strain localization band, it is necessary to ensure that a convenient number of Lobatto control points lie inside the process zone, so the number of control points must be large (at least 10 × 10 per element), as shown in the numerical tests presented in the paper.

The algorithm used in the solution of the non-linear governing system follows a secant Newton-Raphson method. At load step j the iterative algorithm can be described by the following steps:

(1) Initialize the variables by setting 1 and2) in order to define the new values for the damage variable; (d) computation of the secant matrix, A; (e) computation of the residual vector, R; (f) solution of the system A Δsol = -R; (g) update the value for the generalized variables, sol iter = sol iter-1 + Δsol; (h) computation of the new value for the controlling parameter, Error; (4) store the final value for the generalized variables, v (j) = v iter .

v 0 = v (j-1) (2
In the hybrid-mixed stress model, the secant operator A corresponds to the matrix presented in Eq. [START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF]. The solution vector v collects the generalized effective stress parameters, X, and the generalized domain and static boundary displacement variables, q v and q γ . The residual vector R is defined according to Eq.( 10).

Hybrid Displacement Formulation (HD)

Let us express the approximations as:

u = U v q in V, (11)
t = T p on Γ u ,
where the matrices U v and T collect the approximation functions for the displacements in the domain and for the tractions on the boundary Γ u . The vectors q and p list the associated weights. As each field is approximated independently, the approximations may be refined separately.

Enforcing the inner product invariance between discrete and continuum dual quantities, it is possible to obtain the following generalised displacements v and the generalised body forces

Q v : v = T t u dΓ u , Q v = U t v b dV. ( 12
)
In the hybrid displacement formulation, the equilibrium in the domain and the continuity between elements is enforced on average, in a weighted-residual form, using the compatibility equation on Γ u . It can be shown [START_REF] Silva | Hybrid-displacement (trefftz) formulation for softening materials[END_REF] that the first set of equations in the discrete model is:

A c c e p t e d m a n u s c r i p t

Ãq -B p = Q v + Q Γ in V, ( 13
)
where:

à = (D * U v ) t K (D * U v ) dV, (14)
with

K = K (1 -d), B = U t v T dΓ u , Q Γ = U t v t γ dΓ σ ,
In the discrete model, the compatibility equation on the boundary Γ u is given by [START_REF] Silva | Hybrid-displacement (trefftz) formulation for softening materials[END_REF]:

-B t q = -v on Γ u . ( 15
)
The combination of Equation ( 13) and ( 15) leads to the following solving system for each finite element:

⎡ ⎢ ⎣ Ã -B -B t 0 ⎤ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ q p ⎫ ⎪ ⎬ ⎪ ⎭ = ⎧ ⎪ ⎨ ⎪ ⎩ Q v + Q Γ -v ⎫ ⎪ ⎬ ⎪ ⎭ . ( 16
)
Again, the governing system of the finite element mesh is assembled by direct allocation of the contribution of each elementary system [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF].

Because the material is assumed to be non-linear, all the coefficients of the elementary generalised stiffness matrix à are a priori nonzero, see Figure 2(b). Since the dimension of the matrix à is the most important dimension of the governing system, it can become computationally very expensive to store and to solve the non-linear system in the form [START_REF] Freitas | Non-conventional formulations for the finite element method[END_REF]. In order to optimize the numerical performance of the hybrid displacement formulation, two different implementations are considered:

• Implementation 1 (I1): to solve the governing system in the form ( 16);

• Implementation 2 (I2): to update the generalised stiffness matrix K in load step n whenever the number of iterations in load step (n -1) is higher than a certain fixed value. An additional vector on the right hand side of the governing system is introduced, σ 0 :

σ = K ε = K (1 -d 0 ) ε + σ 0 , (17)

A c c e p t e d m a n u s c r i p t

where d 0 corresponds to the damage distribution when the operator K was last updated. It is possible to write:

σ = K ε = K (1 -d) ε = K (1 -d 0 + d 0 -d) ε,
yielding the following definition for vector σ 0

σ 0 = K (d 0 -d) ε,
The equilibrium equation is in this case given by:

à 0 n q -B p = Q v + Q Γ + Q * in V, (18) with
:

à 0 n = (D * U v ) t K (1-d 0 ) (D * U v ) dV ; Q * = (D * U v ) t K ε (d-d 0 ) dV. ( 19
)
Once more, the expressions presented in this paper assume an isotropic model with only one damage variable. If this is not the case, the final expressions are somewhat different and may be found in [START_REF] Silva | Modelos de Dano em Elementos Finitos Híbridos e Mistos[END_REF]. However, all conclusions illustrated are still valid.

For the solution of the non-linear governing system it is adopted an algorithm similar to the one presented for the hybrid-mixed stress model.

Numerical Applications

Three tests are presented: the Hassanzadeh's test, an L-shaped plate and a gravity dam.

The Hassanzadeh's test is used here to compare the performance of the hybridmixed stress model (HMS), the hybrid displacement (HD) model and the classical displacement finite element formulation (FEM).

Initially, the numerical performance of the HMS model is evaluated and compared with the classical FEM, using the numerical results obtained by Claudia Comi, partly published in [START_REF] Comi | A non-local model with tension and compression damage mechanics[END_REF]. In this analysis, the two damage variable model described in Section 3 is chosen. The global load-displacement diagrams and the damage and strain distributions obtained with both models are directly compared. The relative numerical efficiency is qualitatively evaluated using for this purpose the number of iterations required in each step of the loading procedure. It is not useful to compare directly the CPU time since the tests were run using different hardware.

A c c e p t e d m a n u s c r i p t

Then, the numerical performance of the HMS and HD models is assessed assuming a simpler damage model (Model 1 presented in Section 3). The strain field distributions, the CPU time required to complete the analysis and the number of iterations in each step of the loading procedure are compared.

Hassanzadeh's test [19]

The Hassanzadeh's test is illustrated in Figure 3(a). Due to the geometry and to the applied load, only the tension mechanism is activated through the loading history. The structure is analysed as a strain plane problem. In all discretizations a (20 × 20) Lobatto mesh points is used in each element.

HMS versus FEM

The two variable damage model described in Section 3 is used. Following [START_REF] Comi | A non-local model with tension and compression damage mechanics[END_REF] the material parameters are: Table 3 lists the main characteristics of the discretizations associated with each one of the models. The Legendre polynomial degrees adopted to define the approximation for both the stress (S v ) and the displacement (U v ) fields in the domain are given in columns 2 and 3, for the case of Test I and Test II, respectively. The same table lists also the polynomial degrees involved in the definition of the displacement field along the static boundary (U γ ).

a t = 0.
In all discretizations, the approximation is refined near the center of the structure, because the higher stress gradients are located exactly in that region. The FEM meshes use constant strain triangle (CST) elements and correspond to the meshes adopted by Claudia Comi in [START_REF] Comi | A non-local model with tension and compression damage mechanics[END_REF], see Figure 4. For the HMS model two different discretizations are adopted, both using a 7 element mesh (Figure 3(b)), but considering different approximation degrees. Test I refines the approximation mainly at the vertical direction, where Test II uses the same degree of approximation in both directions.

Figure 5 presents the reaction R -prescribed displacement u diagrams, obtained with the numerical simulations and with the experimental test. Both numerical formulations, the HMS and the FEM, correctly model the experimental behavior, except for the experimental bump. On the other hand, the results provided by the different numerical models are almost coincident.

Figure 6 shows that also the damage distribution evolution obtained with both

A c c e p t e d m a n u s c r i p t

numerical models is quite similar.

As shown in Table 3, the use of HMS models leads to larger number of degrees of freedom, N dof . However, the iterative procedure associated with the proposed HMS formulation is less demanding. Figure 7 shows the evolution of the number of iterations required by each incremental load step in terms of the prescribed displacement u. This information is plotted for the HMS Tests I and II and for the FEM hd00d Test. Both formulations use constant local increments and a secant Newton-Raphson technique for the solution of the non-linear governing system. Test hd00d requires 800 load steps until the final prescribed value, u = 0.040 mm, is reached. The test cases run with the HMS model use either 40 or 200 load steps, as indicated in Figure 7. For both formulations it is possible to notice an increase in the number of iterations near the maximum load peak and when the stress dissipation starts to greatly influence the behavior of the structure. In spite of using a larger number of load steps, the curve associated with Test hd00d is always above the curves obtained with the HMS models. When comparing the two tests involving this later model, it is possible to notice a significant decrease in the number of iterations when using 200 load steps instead of 40. This decrease would be even greater if 800 load steps would have been used. This qualitative study indicates that the HMS model is more stable from the numerical point of view than the classical displacement finite element formulations, requiring less number of iterations in each load step to converge1 .

The numerical stability associated with the use of HMS models can be explained by the accuracy of the numerical solution obtained for the strain field, ε, which is used to define and control the material damage evolution and the corresponding mechanical structural behavior. On one hand, in the HMS formulation, the strain field ε is directly computed using the approximation defined for the effective stress fields, σ. When using the classical finite element formulation, the computation of the strain field ε requires the definition of displacement field derivatives. On the other hand, since the damage distribution is continuous, the inter-element continuity of the fields σ/ε is enforced in a weighted residual form in the HMS model through the equilibrium equation on the static boundary. On the contrary, the FEM formulation using the same non-local damage model does not satisfy the inter-element strain continuity. Moreover, the HMS formulation works with macroelements and adopts preferentially a p-refinement instead of an h-refinement, minimizing the number of inter-element boundaries.

A c c e p t e d m a n u s c r i p t

The accuracy of the strain field in each formulation may be confirmed by the analysis of Figures 8 and9, where the evolution of the vertical strain field, ε yy for Tests I and hd00d is shown. This strain component is chosen since it is the one that most influences the non-linear behavior of the structure. The ε yy distribution is quite different in both formulations. The HMS solution shows a smooth variation of this physical quantity in the vertical direction, while the strain field ε yy obtained with the FEM formulation is clearly discontinuous between elements.

HMS versus HD

In this comparative analysis, the model with one damage variable, presented in Section 3, is considered. The corresponding material parameters are given by n = 12, k = 5.8 × 10 -14 MPa, c = 405, l = 1.6 mm, E = 36 GPa and ν = 0.15.

These are the values used by Comi and Perego in [START_REF] Comi | Symmetric and non-symmetric non-local damage formulations: an assessment of merits[END_REF]. All discretizations use the seven element mesh presented in Figure 3(b). The approximation details of each discretization are summarised in Table 4.

According to Figure 10, both models correctly model the global load-displacement diagram. Since the one variable damage model is associated with a larger characteristic length than the model with two damage surfaces [START_REF] Comi | On bifurcation in local and nonlocal materials with tension and compression damage[END_REF][START_REF] Comi | Symmetric and non-symmetric non-local damage formulations: an assessment of merits[END_REF], the global behavior in this case is not so brittle (Figure 7 versus Figure 10).

According to the values presented in Table 4, although the HMS models are associated with a larger number of degrees of freedom when compared to the HD models, the CPU time required by the HMS implementations is always smaller, specially if Implementation I1 (Section 4) is used with the HD formulation. The same happens even if the Implementation I2 (Section 4) is used. This behavior may be justified because the strain filed used in the definition of the damage evolution is obtained more precisely in the HMS model. In this model, the strain field is computed directly from the effective stress field approximation while in the HD formulation the strain field is obtained through the differentiation of the displacement fields modeled in the domain of each element. In terms of vertical strain distribution the final results are quite similar, as illustrated in Figure 11, where the ε yy distribution is presented for all test cases. Only Test C presents a different behavior. However this difference can be overcome by increasing the degree of the approximation in the horizontal direction.

Figure 12 compares the number of iterations required by the convergence of each load step increment, assuming in all test cases a total number of 40 uniform load steps until a maximum prescribed displacement of u = 0.040 mm is reached. When using an implementation I1, the number of iterations per

A c c e p t e d m a n u s c r i p t

load step is similar when using the HMS or the HD models, see Tests A and B versus Test C (I1). As expected, a larger number of iterations is required near the peak of the load-displacement diagram. After that, the number of iterations required by each load increment decreases with the increase of the prescribed displacement. When compared to the case where the bi-dissipative damage model is used, the unloading process is smoother and there is no increase in terms of number of required iterations at the unloading region (Figure 7 versus Figure 12). As CPU times associated with Implementation I1 are unacceptable (see Table 4) the use of the Implementation I2 is highly recommended. In this case, the curve number of iterations versus prescribed displacement presents a very different behavior because the global governing system is updated only when the number of iterations required at the previous load step is greater than a given tolerance.

Since the results of Tests A and B are equivalent in terms of global response, strain field distribution and iterative procedure, one may conclude that the discretization of Test A is sufficiently accurate and that excessive p-refinement associated with a huge number of degrees of freedom (as in Test B) is not needed.

L-shaped plate

Let us consider the L-shaped concrete plate presented in Figure 13. The thickness of the plate is 100mm and an upward vertical displacement, u, at the lowest right corner is prescribed.

The experimental results and the solutions obtained with several numerical simulations are available at the NWD-IALAD website 2 .

The numerical analysis presented here is based on the HMS model. The damage model with one damage variable presented in Section 4 is adopted in this case. The available experimental data are the Young modulus E = 25850 MPa, the Poisson coefficient, ν = 0.18 and the maximum strength in tension, f t = 2.70MP a. The remaining material parameters are defined in order to minimize the differences between the experimental and numerical loadprescribed displacement diagrams. The following values have been assumed:

n = 9

A c c e p t e d m a n u s c r i p t

U ξ,η γ = 7, yielding a total of 1385 degrees of freedom. In each finite element, a (20 × 20) Lobatto mesh points is used to perform all numerical integrations and to compute the non-local variables. Test 2 uses the 10 element mesh shown in Figure 14 and assumes S ξ,η v = 7, U ξ,η v = 6 and U ξ,η γ = 6. This discretization involves a total of 3251 degrees of freedom. In this case a (10 × 10) Lobatto mesh points is used.

The reaction-prescribed displacement u diagrams are presented in Figure 15. It is possible to confirm that both numerical and experimental diagrams are quite similar.

The evolution for the horizontal displacement at the upper left corner of the structure is shown in Figure 16. The numerical response is stiffer than the experimental behavior at the beginning of the loading procedure. This type of behavior is also observed in other numerical simulations and can be justified by the fact that the numerical simulation does not take into account the rotational stiffness of the steel device that embraces the L-shaped concrete structure.

Figure 17 presents the damage evolution and identifies the principal directions of the effective stress tensor. As expected, damage first appears near the reentrant corner with a diagonal orientation and evolves towards an almost horizontal macro crack.

Finally, Figure 18 shows the damage evolution represented in a deformed configuration. The damage appearing in the lower right corner is a local effect induced by the prescribed displacement imposition and does not interfere with the global behavior of the structure.

Koyna dam

The Koyna dam is located in the southern part of India. The geometry of this concrete gravity dam is presented in Figure 19(a). In December 1967, this dam was subjected to an important earthquake. In the last few years, several authors have tried to compute numerically the damage distribution installed due to that earthquake, see for instance [START_REF] Cervera | Seismic evaluation of concrete dams via continuum damage models[END_REF][START_REF] Lee | A plastic-damage concrete model for earthquake analysis of dams[END_REF].

The numerical study presented in this paper aims to identify the highest water lever in the reservoir ensuring dam security conditions considering an initial damage distribution caused directly by the earthquake, the dam self weight and the hydrostatic lateral pressure.

Following the work of [START_REF] Cervera | Seismic evaluation of concrete dams via continuum damage models[END_REF], an initial qualitative damage distribution is taken into account, as illustrated in Figure 20. The initial damage distribution is con- A plane strain behavior is assumed. The finite element mesh presented in Figure 19(b) is used in the analysis. Following [START_REF]ICOLD. V international benchmarkworkshop on numerical analysis of dams, theme a2[END_REF], the rock foundation is included in the finite element analysis. The degrees of approximation are defined by S ξ,η v = 5, U ξ,η v = 4, U ξ,η γ = 3 and a (20 × 20) Lobatto mesh points is defined in each element.

The bi-dissipative damage model presented in Section 4 is used in the analysis. The material parameters for concrete and rock are defined according to [START_REF] Cervera | Seismic evaluation of concrete dams via continuum damage models[END_REF][START_REF] Comi | A bi-dissapative damage model for concrete with applications to dam engineering[END_REF]: The characteristic length is l t = 400 mm. The behavior to compression states is considered to be linear elastic.

•
Figure 21 shows the evolution of the horizontal displacement measured at the top of the dam in terms of parameter q = h water /h dam . When evaluating the load carrying capacity and the safety of a dam against failure, it is a common procedure to consider values greater than unity for parameter q, leading to the so called imminent failure flood analysis.

The damage distribution is presented in Figure 22 for all three cases. It is possible to notice a significant damage evolution near the interface foundationstructure. There is also a slight damage increase at the region where de the geometry of the structure changes.

Interestingly, the comparison of the solutions obtained for the three cases shows that the initial damage value has a very little influence in terms of final collapse load and displacements fields. This conclusion was also obtained by other authors [START_REF] Cervera | Seismic evaluation of concrete dams via continuum damage models[END_REF][START_REF] Comi | A bi-dissapative damage model for concrete with applications to dam engineering[END_REF].

Conclusions

Two promising non-conventional finite element formulations, namely the hybridmixed stress (HMS) and the hybrid displacement (HD) models, are discussed in this paper. The main properties of such models and their numerical per-

A c c e p t e d m a n u s c r i p t

formance are illustrated through the analysis of a set of classical benchmark tests.

The main conclusions can be summarised as follows:

(1) both the HMS and the HD models lead to stable and robust numerical procedures. The HMS models show better convergence rates than the classical displacement FEM implementations. This efficiency can be justified by the accuracy ensured for the effective stress/strain distributions that control the damage evolution behavior; (2) in all performed numerical tests, the quality of the solutions does not depend on the finite element mesh orientation (mesh bias). This behavior results mainly from the use of macroelement meshes associated with the implementation of highly effective p-refinement procedures. As discussed in [START_REF] Jirásek | Modeling of localized inelastic deformation[END_REF], this type of phenomena may influence the quality of the results provided by the classical FEM computations; (3) the HMS model may lead to quasi -equilibrated solutions. Consequently, this model may provide lower bounds for the structure collapse loads.

The physically non-linear analysis of concrete structures is still a challenge and corresponds to an application field where the high performance alternative formulations presented in this paper can be conveniently explored.

The fracture processes strongly influence the structural behavior of concrete structures and may be modeled using continuous or discrete approaches. One strong drawback associated to the application of pure continuous models lies in the fact that fragile rupture is frequently governed by the growth of a dominant crack. Although both approaches continue to deserve the scientific community attention, the current trend is clearly to combine the two methodologies. Therefore, the main expected development corresponds to the implementation of cohesive crack models based on the use of non-conventional finite element formulations. Then, some of the existing alternative schemes used to account for the transition from continuum to discrete approaches [START_REF] Comi | From localized damage to discrete cohesive crack propagation in nonlocal continua[END_REF][START_REF] Simone | Continuous-discontinuous modelling of failure[END_REF][START_REF] Oliver | Continuum approach to material failure in strong discontinuity settings[END_REF] will be studied, implemented, tested and compared.

The other developments expected for the near future have as main goal to make even more competitive the HMS and the HD non-conventional formulations described in this paper. Among several research lines, one can mention the following:

(1) to implement p-refinement procedures based on the use of local functions, such as wavelets or Heaviside functions [START_REF] Silva | Nonlocal damage theory in hybriddisplacement formulations[END_REF]. Using the critical damage concept [START_REF] Comi | Criteria for mesh refinement in nonlocal damage finite element analyses[END_REF], this refinement may be adaptive; (2) the non-conventional finite element models allow the use of local refinement functions. This is one important issue when working with damage models, as the non-linear behavior concentrates usually in a specific zone

A c c e p t e d m a n u s c r i p t

of the structural domain.

(3) to implement an implicit gradient model in order to avoid the definition of an independent spatial averaging [START_REF] Peerlings | Gradient-enhanced damage for quasi-brittle materials[END_REF][START_REF] Peerlings | Some observations on localization in non-local and gradient damage models[END_REF][START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF][START_REF] Engelen | Nonlocal implicit gradientenhanced elastoplasticity for the modelling of softening behaviour[END_REF][START_REF] Simone | Continuous-discontinuous modelling of failure[END_REF]. One of the important advantages of this type of procedure is that the computation of tangent stiffness matrices is simpler than in other types of formulations; (4) to test the numerical performance of the HMS and HD finite element models when using parallel processing techniques; (5) to generalise the existing 2D models for 3D analysis; [START_REF] Comi | Computational modelling of gradient-enhanced damage in quasibrittle materials[END_REF] to generalise the models for cyclic and dynamic analysis.

A c c e p t e d m a n u s c r i p t

Model 1 [START_REF] Lemaitre | A course on damage mechanics[END_REF][START_REF] Comi | Nonlocal aspects of nonlocal damage analyses of concrete strucutres[END_REF][START_REF] Comi | Symmetric and non-symmetric non-local damage formulations: an assessment of merits[END_REF] Damage Variable d (c) for this particular damage model, the internal variable ξ coincides with the damage variable d. Table 3 Hassanzadeh's test [START_REF] Hassanzadeh | Behaviour of fracture process zone in concrete influenced by simultaneous applied normal and shear displacements[END_REF]: comparative study of HMS and FEM models.

Helmholtz free energy density

(a) Ψ = 1 2 (1 -d) ε t E ε + Ψ in (ξ) with Ψ in (ξ) = k (1 -ξ) n i=0 n! i! ln i c 1-ξ State equations (b) σ = ∂Ψ ∂ε = (1 -d) E ε χ = -∂Ψ ∂ξ = -Ψ in (ξ) Y = -∂Ψ ∂d = 1 2 ε t E ε Activation function f (Y -χ) = Y -χ = 1 2 ε t E ε -χ Evolution laws (c) ḋ = ∂f ∂Y γ = γ ξ = -∂f ∂χ γ = γ Kuhn Tucker Conditions f ≤ 0, ḋ ≥ 0, ḋ f = 0 Non-local Variable Y = 1 2 ε t E ε (a)
f t = J 2 -a t I 2 1 + b t r t I 1 -k t r 2 t (1 -α d c ) f c = J 2 + a c I 2 1 + b c r c I 1 -k c r 2 c r i (d i ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1 - 1 -( σe σo ) i d 2 oi (d oi -d i ) 2 se d i < d oi 1 - d i -d oi 1 -d oi

A c c e p t e d m a n u s c r i p t 4

 4 Non-conventional Finite Element Formulations 4.1 Hybrid-mixed Stress Formulation (HMS)

) Error= 10× tol and iter=1 ( 3 )

 3 while Error > tol iter=iter+1 (b) computation of the non-local variable at each Lobatto point; (c) validation of the Kuhn-Tucker conditions (see Tables

  31, b t = 4.4 MPa, k t = 15.1 MPa 2 , c t = 35, d ot = 0.10, (σ e /σ o ) t = 0.80, l t = 1.0 mm, E = 36000 MPa and υ = 0.15.

. 5 ,= 10 ,2

 510 k = 1.1 × 10 -11 MPa, c = 270 and l = 11 mm.A plane stress behavior is considered and the vertical displacement at the lower right corner is prescribed. Two different numerical tests have been performed. Test 1 considers the three element mesh represented in Figure14. The discretization associated with this test uses S ξ,η v Currently, http://nw-ialad.uibk.ac.at/, grupo WP2/TG2, Test n. 1

A c c e p t e d m a n u s c r i p t

  centrated near the structure-foundation interface and near the dam geometric change. Three different damage values are considered in the region with initial diffuse cracking (red region in Figure20): d initial = 0.00 (Test 1), d initial = 0.30 (Test 2) and d initial = 0.60 (Test 3). For the numerical simulations the HMS model is used.

  for concrete: E = 30 GPa, υ = 0.20, ρ = 2630 kg/m 3 , a t = 0.25, b t = 4.4 MPa, k t = 9.14 MPa 2 , c t = 3.3, (σ e /σ o ) t = 0.70 and d ot = 0.0; • for rock: E = 41 GPa, υ = 0.10, ρ = 2700 kg/m 3 , a t = 0.333, b t = 3.3 MPa, k t = 8.6 MPa 2 , c t = 3.3, (σ e /σ o ) t = 0.70 and d ot = 0.0.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Nonzero coefficients in the governing system of one finite element (Equation (10)) (a) without damage and (b) with damage.

Fig. 4 .A c c e p t e d m a n u s c r i p tFig. 5 .A c c e p t e d m a n u s c r i p tFig. 7 .

 457 Fig. 4. Hassanzadeh's test[START_REF] Hassanzadeh | Behaviour of fracture process zone in concrete influenced by simultaneous applied normal and shear displacements[END_REF]: finite element meshes used in the analysis with FEM[START_REF] Comi | A non-local model with tension and compression damage mechanics[END_REF].

Fig. 8 .A c c e p t e d m a n u s c r i p t

 8 Fig. 8. Hassanzadeh's test [19]: ε yy [-] evolution in Test I.

A c c e p t e d m a n u s c r i p tFig. 10 .

 10 Fig. 10. Hassanzadeh's test [19]: reaction R [N]-prescribed displacement u [mm].

Fig. 11 .A c c e p t e d m a n u s c r i p tFig. 12 .

 1112 Fig. 11. Hassanzadeh's test [19]: ε yy strain field distribution[-] when u = 0.040 mm.

A c c e p t e d m a n u s c r i p tFig. 13 .A c c e p t e d m a n u s c r i p tFig. 14 .A c c e p t e d m a n u s c r i p tFig. 15 .A c c e p t e d m a n u s c r i p tFig. 16 .Fig. 20 .A c c e p t e d m a n u s c r i p tFig. 21 .

 131415162021 Fig. 13. L-shaped plate: experimental device (available at the NWD-IALAD website) .

  Dam without initial damage. (c) Initial damage of 0.30. (d) Initial damage of 0.60.

Fig. 22 .

 22 Fig. 22. Koyna dam: damage distribution on the deformed configuration at the final stage of the loading procedure.

  ξ is a scalar internal variable of kinematic nature and variables k, n and c are material parameters. (b) ε, ξ and d are the state variables and σ, χ and Y represent the corresponding associated variables.

A c c e p t e d m a n u s c r i p t Model 2 [ 9 , 2 μ

 292 10,7] Damage Variable d t , d c Helmholtz free energy density (a,b) ψ = 1 2 2μ e : e + k + (tr +ε) 2 + k -(tr -ε) = μ 0 (1d t )(1d c ) k + = k 0 (1d t ) if tr ε ≥ 0 k -= k 0 (1d c ) if tr ε < 0 tr + ε = trε + |trε| 2 , tr -ε = trε -|trε| 2 Constitutive Relation σ = ∂ψ ∂ε = 2μ e + k + (tr + ε) I + k -(tr -ε) IActivation functions(c,d,e)

c i 0 2 s

 02 .75se d i ≥ d oi Kuhn Tucker Conditions f t ≤ 0, ḋt ≥ 0, ḋt f t = 0 f c ≤ 0, ḋc ≥ 0, ḋc f c = 0 Non-local Variable (f ) tr + ε, tr -ε, J ε (a) μ 0 , k 0 representthe initial undamaged shear and bulk moduli, respectively. (b) e denotes the deviatoric part of the small strain tensor. (c) a t , b t , k t , a c , b c , k c and α are non-negative material parameters. (d) I 1 corresponds to the stress tensor first invariant and J 2 = 1 : s. (e) d oi represents the damage value at the peak of the uniaxial diagram.

  As indicated in Figure 3(b). (b) Notation adopted by Claudia Comi.

  As explained in Section 4 and adopting the value 50 as the maximum number of iterations at the previous load step. (b) As defined in Section 4. (c) In a Pentium 4 machine, 2996 Mhz.

Table 1

 1 Continuum damage model with one scalar variable.

Table 2

 2 Continuum damage model with two variables.

Table 4

 4 Hassanzadeh's test[START_REF] Hassanzadeh | Behaviour of fracture process zone in concrete influenced by simultaneous applied normal and shear displacements[END_REF]: comparative study of HMS and HD models.

Please note that the error criteria is probably not the same for both formulations, so the comparison made in this paper is only qualitative. However, it is shown that the three tests lead to identical results in terms of global behavior and damage distribution, so the qualitative comparison is adequate.

L-shaped plate, Test 1: damage variable evolution [-] and effective stress tensor principal directions.

L-shaped plate, Test 1: damage variable evolution [-] on the structure deformed configuration.

Koyna dam: (a) geometry of the structure (dimensions in m) and (b) finite element mesh.

Koyna dam: qualitative initial damage distribution.

Koyna dam: evolution of the horizontal curves q = h water /h dam -displacement measured at the top of the dam [mm].

Koyna dam: damage distribution on the deformed configuration at the final stage of the loading procedure.

Acknowledgements

This work has been partly supported by "Fundação para a Ciência e Tecnologia" through projects "Financiamento Plurianual", PTDC/ECM/71519/2006 and the PhD scholarship SFRH/BD/9050/2002. The guidance of Professors Claudia Comi and Umberto Perego is gratefully acknowledged.

A c c e p t e d m a n u s c r i p t

Nonzero coefficients in the governing system of one finite element (Equation ( 10)) (a) without damage and (b) with damage.

2 Nonzero coefficients in the governing system of one finite element (Equation ( 16)) (a) without damage and (b) with damage.

3

Hassanzadeh's test [START_REF] Hassanzadeh | Behaviour of fracture process zone in concrete influenced by simultaneous applied normal and shear displacements[END_REF]: (a) geometry and (b) adopted finite element mesh.

4

Hassanzadeh's test [START_REF] Hassanzadeh | Behaviour of fracture process zone in concrete influenced by simultaneous applied normal and shear displacements[END_REF]: finite element meshes used in the analysis with FEM [START_REF] Comi | A non-local model with tension and compression damage mechanics[END_REF]. Hassanzadeh's test [START_REF] Hassanzadeh | Behaviour of fracture process zone in concrete influenced by simultaneous applied normal and shear displacements[END_REF]: comparative study of HMS and FEM models.

4

Hassanzadeh's test [START_REF] Hassanzadeh | Behaviour of fracture process zone in concrete influenced by simultaneous applied normal and shear displacements[END_REF]: comparative study of HMS and HD models.