Quang-Huy Nguyen

Not a Member, ASCE Mohammed Hjiaj

ASCE Member

Nonlinear time-dependent behavior of composite steel-concrete beams

Keywords: Composite beam, time effects, nonlinear creep, cracking, mixed formulation

This paper presents a mixed finite element model for the nonlinear time-dependant analysis of composite beams with partial shear connection. The key idea is to consider, as a first approach, a viscoelastic/plastic model for the concrete slab in order to simulate the interaction between the time effects of concrete, such as creep and shrinkage, and the concrete cracking. Creep is taken into account via linear ageing viscoelasticity, while cracking is modeled using an elasto-plastic model with softening. A nonlinear isotropic/kinematic hardening model is adopted for steel behavior and an appropriate nonlinear constitutive relationship is utilized for the shear stud. A consistent time-integration is performed by adopting the Euler backward scheme. Finally, comparisons between the numerical results and experimental data available in the literature are undertaken to validate the accuracy of the model. It is shown that the interaction between cracking and time effects (creep and shrinkage) significantly increases the deflection.

INTRODUCTION

Steel-concrete composite structures are common practice today in bridges and industrial buildings. The advantages of both materials lead to a very economic alternative especially in terms of high bearing capacity. The serviceability of such structures is fundamentally affected by creep and cracking and their interaction, which make the structure behavior becomes prominent nonlinear. Recent publications have highlighted the effect of creep and shrinkage on the overall behavior of composite beams. The authors have adopted in most cases linear ageing viscoelasticity which can be described using either a hereditary integral formulation or a differential formulation with internal variables. The former formulation is not easy to handle and the most straightforward manner to deal with such constitutive law is to approximate the hereditary integral by a discrete sum according to the trapezoidal or the midpoint rule. Therefore, the storage of the entire stress history is required. This method, called the step-by-step or general method, has been proposed by [START_REF] Bažant | Numerical determination of long-range stress history from strain history in concrete[END_REF]. For practical applications, the so-called algebraic methods, which ignore stress history, have been proposed [START_REF] Mcmilan | Method of designing rienforced concrete slabs, discussion by ac janni[END_REF][START_REF] Hansen | Estimating stress relaxation from creep data[END_REF][START_REF] Bažant | Prediction of concrete creep effects using age-adjusted effective modulus mothed[END_REF]. Obviously, the algebraic methods give different levels of accuracy depending on the degree of refinement of the quadrature formulas adopted [START_REF] Dezi | Algebraic methods for creep analysis of continuous composite beams[END_REF]. The problem of the storage of stress history can be elegantly overcome by adopting a differential formulation [START_REF] Bažant | Dirichlet series creep function for aging of concrete[END_REF][START_REF] Bažant | Rate-type law of aging concrete based on maxwell chain[END_REF]. By expanding the kernel of the hereditary integral into a Dirichlet serie, the integral-type stress-strain relation can be transformed into a rate-type stress-strain relation based on Kelvin chain or Maxwell chain spring-dashpot model of ageing viscoelasticity.

In this strategy, the stress or strain history is memorized via a limited number of internal variables. However, the main drawback of differential formulation is the determination of Dirichlet series coefficients which appears to be quite complicated.

It has been shown that linear viscoelasticity with age-dependent constitutive parameter is not able to correctly predict the long term behavior of concrete under high load levels [START_REF] Ngab | Shrinkage and creep of high strength concrete[END_REF]). Indeed, several experimental studies indicate that under high sustained loads, cracks grow and interact with the creep response of the bulk concrete. This interaction may lead to a reduction of the ultimate bearing capacity of structures [START_REF] Bažant | Drying creep of concrete: constitutive model and new experiemnts separating its mechanisms[END_REF]. This is often denoted as nonlinear creep which is a quite complex phenomenon. Several attempts have been made to model the interaction between creep and concrete damage.

2

Nguyen, April 13, 2015

The first models for time-dependent behavior of concrete (creep and relaxation) together with instantaneous non-linear behavior consisted of some kind of superposition in terms of strains between a non-linear explicit stress-strain curve σ = f (ε) (total strain model) and a linear time-dependent model. This approach was adopted by [START_REF] Carol | A model for the non-linear time-dependent behaviour of concrete in compression based on a maxwell chain with exponential algorithm[END_REF] and [START_REF] Fragiacomo | Finite-element model for collapse and long-term analysis of steelconcrete composite beams[END_REF] in which the creep of concrete was taken into account using

Maxwell's generalized rheological model through a step-by-step time increment procedure.

Usually, these types of explicit or total strain approaches are relatively simple and inexpensive compared with other more complex incremental formulations, and allow the engineer to perform the time-dependent static analysis of frames. However, these models seems to be incapable of representing some quite significant experimental features of concrete behavior and more advanced approaches have been devised. One route is to consider the concept of energy activation and rate dependent plasticity (softening) as proposed by [START_REF] Bažant | R-curve modeling of rate and size effects in quasibrittle fracture[END_REF]. In a second approach, a standard rheological model for creep is incorporated into the fictitious crack model in order to simulate time dependency of crack opening [START_REF] Zhou | Time-dependent fracture of concrete: Testing and modelling[END_REF][START_REF] Barpi | Creep and fracture in concrete: A fractional order rate approach[END_REF]. In the last class of models, rheological models are built up by combining either smeared crack representation (softening plasticity) or damage model with ageing linear viscoelasticity. The main advantage of the last class of the model is the ease in the implementation as well as the effectiveness in reproducing key feature of the nonlinear cracking process. [START_REF] Barpi | Creep and fracture in concrete: A fractional order rate approach[END_REF] have enhanced the latter approach using a fractional order rate law. [START_REF] Benboudjema | Interaction between drying, shrinkage, creep and cracking phenomena in concrete[END_REF] proposed a sophisticated hydro-mechanical model accounting for the full coupling between drying, shrinkage, creep and cracking. In their model, a new creep constitutive law, based on microscopic considerations of the role of water, is elaborated and coupled with an improved version of the drying creep model proposed by [START_REF] Bažant | Concrete creep at variable humidity: constitutive law and mechanism[END_REF] in order to account for the interaction between the solid skeleton and the water.

In composite steel-concrete beams with flexible connectors, a slip at the steel-concrete 3 Nguyen, April 13, 2015 interface occurs, which cannot be neglected in the analysis and the design of composite beams. A state of the art on the time-dependent behavior of steel-concrete structures can be founded in [START_REF] Ranzi | State of the art on the time-dependent behaviour of composite steel Űconcrete structures[END_REF]. At the very beginning, time-dependent analysis of composite beams was carried out by assuming full interaction between the concrete slab and the steel beam that is no slip occurs at the steel/concrete interface. This leads to analytical procedures for time-dependent analysis of composite steel-concrete cross-sections in full interaction that have been proposed by [START_REF] Gilbert | Time-dependent analysis of composite steel-concrete sections[END_REF]. Early papers on time-analysis of composite beams taking into account partial shear interaction were published by [START_REF] Bradford | Composite beams with partial interaction under sustained loads[END_REF]; Tarantino andDezi (1992, 1993). [START_REF] Bradford | Composite beams with partial interaction under sustained loads[END_REF] proposed a method using an equivalent transformed section approach for the time-dependent response of simply supported steel-concrete composite beams. [START_REF] Tarantino | Creep effects in composite beams with flexible shear connectors[END_REF] adopted the step-by-step approach to discretize the constitutive concrete model which was combined with the finite-difference method for the space discretization. [START_REF] Virtuoso | Time dependent behaviour of continuous composite beams with flexible connection[END_REF] developed a force-based finite element formulation using the rate-type method to model the time-dependent behavior of the concrete. Recently, closed-form solutions have been proposed in Nguyen et al. (2010a[START_REF] Nguyen | Time-dependent analysis of composite beams with partial interaction based on a time-discrete exact stiffness matrix[END_REF] for the general method, in [START_REF] Faella | Steel and concrete composite beams: "exact" expression of the stiffness matrix and applications[END_REF]; [START_REF] Ranzi | Analytical solutions for the time-dependent behaviour of composite beams with partial interaction[END_REF] beams with partial shear connection taking into account the interaction between the time effects and the nonlinear behavior of concrete (cracking). This model is based on the two fields mixed force-displacement formulation with nonlinear constitutive relationships for each components [START_REF] Zienkiewicz | The finite element method[END_REF][START_REF] Ayoub | Mixed formulation of nonlinear steel-concrete composite beam element[END_REF][START_REF] Nguyen | Analysis of composite beams in the hogging moment regions using a mixed finite element formulation[END_REF].

The mixed formulation is adopted in our work because it combines the advantages of the displacement and force formulations while overcoming most of their limitations, especially for nonlinear analysis of composite beams in partial interaction [START_REF] Spacon | Nonlinear analysis of steel-concrete composite structures: State of the art[END_REF].

Regarding the nonlinear time-dependent behavior of concrete, a viscoelastic/plastic model is proposed in order to take into account the interaction between creep, shrinkage and plasticity (smeared-crack model). The key idea, which is borrowed from [START_REF] Van Zijl | A numerical model for the time-dependent cracking of cimentitious materials[END_REF], is to combine, as a first approach, a viscoelastic creep model with a continuum plasticity model. Creep is formulated via Kelvin chain spring-dashpot model of aging viscoelasticity (differential formulation) proposed by [START_REF] Jirásek | Incremental viscoelastic analysis of composite beams with partial interaction[END_REF]. The concrete cracking is modeled using an 1D elasto-plastic model with softening proposed by the first author [START_REF] Nguyen | Analysis of composite beams in the hogging moment regions using a mixed finite element formulation[END_REF]. A consistent time-integration is performed by adopting the Euler backward scheme. A nonlinear isotropic/kinematic hardening model is adopted for steel behavior [START_REF] Mahnken | Improved implementation of an algorithm for non-linear isotropic/kinematic hardening in elastoplasticity[END_REF]) and an appropriate nonlinear constitutive relationship is utilized for the shear stud [START_REF] Ollgaard | Shear strength of stud connectors in lightweight and normal weight concrete[END_REF][START_REF] Aribert | Modèle calcul élasto-plastique de poutres mixtes à connexion partiell[END_REF]. The proposed FE model is employed to investigate the nonlinear long-term behavior of a two-span continuous composite beam tested by [START_REF] Gilbert | Time-dependent behavior of continuous composite beams[END_REF]. The accuracy of the model is assessed by comparing the predicted deflections against measured deflections. Finally, the contribution of nonlinear creep is evaluated by comparing the deflections predicted by the proposed model against those provided by ageing linear viscoelasticity.

VISCOELASTIC-PLASTIC MODEL FOR CONCRETE

In the proposed viscoelastic-plastic model, it is assumed that total strain rate ε is defined 5 Nguyen, April 13, 2015

as the sum of different contributions

ε = εve + εp + εsh (1)
where εve is viscoelastic strain rate due to concrete creep; εp is strain rate due to plasticity/cracking; and εsh is shrinkage strain rate. The rheological viscoelastic/plastic model is illustrated in Fig. 1. The creep strain rate is evaluated here via an ageing rheological Kelvin chain (with Kelvin elements placed in serial), while the cracking strain is formulated using continuum plasticity theory (smeared cracking model). Further, the shrinkage strain rate can be evaluated by utilizing shrinkage models proposed in design codes (ACI209 1997;[START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures -model b3[END_REF]CEB-FIB 1990).

In this work, concrete creep is modeled using the rate-type formulation with internal variables of the linear viscoelasticity, initially proposed by [START_REF] Bažant | Rate-type law of aging concrete based on maxwell chain[END_REF]. This formulation is based on expanding the creep function J(t, τ ) into a Dirichlet series:

J(t, τ ) ≈ 1 E 0 + m i=1 1 D i (τ ) 1 -exp (- t -τ τ i ) (2)
where τ i=1..m are fixed parameters called the retardation times; E 0 is the asymptotic modulus and D i=1..m are age-dependent moduli which are determined by numerical identification to analytical creep function of CEB-FIB 1990. The details of this identification can be found in [START_REF] Nguyen | Modelling of the nonlinear behaviour of composite beams taking into account time effects[END_REF]. The approximation of the creep function leads to a Kelvin chain of m units with ageing moduli E m (t) and dashpot viscosities η m (t) as shown in Fig. 1. For this chain, the viscoelastic strain can be expressed as

ε ve (t) = σ(t) E 0 + m i=1 ε i (t) (3) 6 Nguyen, April 13, 2015
where ε i (t) is the strain of each Kelvin unit and must satisfy the following rate-type equation:

E i (t) εi (t) + ηi (t) εi (t) + η i (t)ε i (t) = σ(t) (4)
The above differential equation can be solved using the exponential algorithm developed by [START_REF] Bažant | Proccedings of the 1st international conference on structure mechanics in reactor technology, west berlin[END_REF] for ageing concrete creep, which results in the following relationship between the stress increment ∆σ from the time t n to the time t n+1 and the corresponding viscoelastic strain increment ∆ε ve :

∆σ = E ve ∆ε ve + σ n (5)
where E ve is an effective stiffness modulus, which depends only on the current time t n+1 ;

and σ n is a viscous stress which accounts for the whole stress history. Full details of this exponential algorithm can be found in [START_REF] Jirásek | Incremental viscoelastic analysis of composite beams with partial interaction[END_REF].

In the present study, the cracking phenomena is described with a continuum plasticity model in which the plastic strain rate is determined by the plastic flow rule:

εp = λ ∂Φ ∂σ ; λ ≥ 0 and Φ ≤ 0 (6)
where λ is the plastic multiplier and Φ is the yield/craking function. Assuming isotropic hardening/softening, the yield function takes the following form in 1-D:

Φ =      σ -R t (p) in tension -σ -R c (p) in compression (7)
in which R t, c (p) is the hardening/softening functions and p the equivalent plastic strain that describes the hardening/softening behavior. In tension, the following hyperbolic softening 7 Nguyen, April 13, 2015 law is adopted [START_REF] Nguyen | Modelling of the nonlinear behaviour of composite beams taking into account time effects[END_REF]):

R t (p) = f t 1 + f t l eq G t p 2 (8)
where f t is the tension strength; G t the fracture energy; and l eq the equivalent length of the respective FE integration point which is introduced according to the concept of constant fracture energy to avoid the loss of objectivity due to the strain localization [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF][START_REF] Coleman | Localization issues in force-based frame elements[END_REF].

In compression, the expressions of the hardening/softening functions will be determined such as to reproduce the uniaxial stress-strain explicit function given by [START_REF] Krätzig | An elasto-plastic damage model for reinforced concrete with minimum number of material parameters[END_REF] which fits well with experiments. This function is illustrated by three parts in Fig. 2.

Elastic part: below the initial yield/damage stress f cy = f c /3 linear-elastic behavior is assumed. Hardening part: the stress grows until failure strength f c , consequently the tangential stiffness decreases from initial stiffness to zero (horizontal tangent). Softening part: After exceeding the compression strain ε c , localization of damage occurs in this softening region.

The softening branch depends on the fracture energy and on the characteristic length leq to avoid possible ill-posedness problem [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF]. The explicit uniaxial stressstrain expression has the form:

σ (ε) =                            E c ε for ε ≥ -ε cy E ci ε f c + ε ε c 2 1 - E ci ε c f c -2 ε ε c f c for -ε cy > ε ≥ -ε c -1 2 + γ c f c ε c 2f c + γ c ε + γ c 2ε c ε 2 for ε < -ε c (9) 8
Nguyen, April 13, 2015

where

E ci = f 2 c 2E c ε 2 c - f c ε c + 3E c 2 ; ε cy = f c 3E c and γ c = π 2 f c ε c 2 G cl leq -f 2 c 2Ec 2 (10)
with f c being the failure stress and ε c the corresponding strain. E c is the initial modulus of elasticity and G cl the localized crushing energy.

The identification of the hardening and softening functions is performed as follows. According to the Hooke's law, the stress is expressed as

ε = σ E c + ε p (11)
Otherwise, the associated plastic evolution equations reads εp = -ṗ = -λ (12)

Integrating the above relation, one obtains

ε p = -p (13)
From the yield condition, one has

R c = -σ (14)
Inserting Eqs. ( 13) and ( 14) into [START_REF]Time evolution of the most tensile fiber stress at intermediate support[END_REF], one obtains

ε = - R c E c -p (15) 9
Nguyen, April 13, 2015

By introducing the relations ( 14) and ( 15) into ( 9), one obtains

R c =                            E ci f c R c E c + p - 1 ε 2 c R c E c + p 2 1 + E ci f c - 2 ε c R c E c + p f c for p ≥ p 1 2 + γ c f c ε c 2f c -γ c R c E c + p + γ c 2ε c R c E c + p 2 for p < p ( 16
)
where

p = ε c - f c E c (17)
By solving the equation ( 16) for the hardening/softening function R c , one obtains:

R c =            ξ 1 p 2 + ξ 2 p + ξ 2 3 -ξ 4 p -ξ 3 with p ≥ p 1 3 3 3 i=0 α i p i 2 + 2 i=0 β i p i 3 ± 3 i=0 α i p i + 2 3 E c (ε c -p) with p < p (18)
in which ξ 1..4 , β 0, 1 , α 0..2 and µ 0..3 are the parameters which are determined through the material parameters such as f c , ε c , E c and G cl , as given in Appendix.

Integration of the viscoelastic/plastic model

Numerical implementation of the model requires integrating the rate form of the constitutive relations for a finite time step ∆t = t n+1 -t n . The objective is to compute, for a given set of state variables at the time t n , a shrinkage strain increment ∆ε sh and a total strain increment ∆ε, the respective updated state variables at the end of the time step t n+1 .

An Euler-Backward algorithm is utilized for this purpose since it has the property of being stable [START_REF] Ortiz | Accuracy and stability of integration algorithms for elastoplastic constitutive relations[END_REF][START_REF] Simo | Computational Inelasticity[END_REF]. [START_REF]Time evolution of hogging bending moment at intermediate support[END_REF] Nguyen, April 13, 2015

The numerical integration of the model is carried out by the following procedure:

1. Compute the effective modulus E ve and the viscous stress σ n according to the exponential algorithm [START_REF] Jirásek | Incremental viscoelastic analysis of composite beams with partial interaction[END_REF].

2. Predictor: compute viscoelastic trial stress and test for plastic loading

σ trial n+1 = σ n + E ve ∆ε -∆ε sh (19) Φ trial = Φ σ trial n+1 , p n (20)
If Φ trial ≤ 0 then viscoelastic step: set (•) n+1 = (•) trial n+1 and exit else viscoelastic/plastic step: proceed to step (3). with Newton-Raphson method and update the others variables

p n+1 = p n + ∆λ (22) ε p n+1 = ε p n + sign(σ trial n+1 )∆λ (23)
σ n+1 = σ trial n+1 -E ve sign(σ trial n+1 )∆λ (24)
4. Compute the consistent tangent modulus by

E tg n+1 = E ve 1 - E ve E ve + dR dp (25) 11
Nguyen, April 13, 2015

Qualitative validation of the algorithm

The proposed model is used to simulate a relaxation test. We consider the concrete specimen C30 suggested in the CEB-FIB 1990. The length of the specimen is 200 mm and its cross-section is 100x100 mm 2 . This specimen is subjected to an imposed compressive displacement at its end. Creep and shrinkage functions recommended by the CEB-FIB 1990 are utilized considering that the cement is normal hardening and the humidity (RH) is about 80%.The following material parameters are adopted: As can be seen the total strain is constant over the time because of the imposed displacement (relaxation test). The shrinkage strain is stress-independent and increase in magnitude with time in a monotonic way whereas the creep strain is stress-dependent and varies accordingly. In the present relaxation test, shrinkage will generate tension stresses while creep stresses are compressive. As both the creep and shrinkage strains increase in magnitude but with opposite sign, the yield stress will be reached at a certain time instant. At this stage, plastic deformation will start to develop in order to maintain the stress on the yield surface. However, with the presence of shrinkage strain which is negative, the viscoelastic strain has to vary initially in the contrary way of the shrinkage strain in order to keep the total strain constant. Then the plastic flow is active (at 300 days), the plastic strain starts to increase with time. This leads to a decrease of the viscoelastic strain observed in Fig. 4. The numerical results show that the proposed model seems to be able to represent the cracking phenomena due to shrinkage (Fig. 3).

f c = 38 MPa; f t = 2.

NONLINEAR ISOTROPIC/KINEMATIC HARDENING ELASTO-PLASTIC MODEL

12 Nguyen, April 13, 2015

FOR STEEL

In our work, the von Mises plasticity model with combined isotropic and kinematic hardening is adopted for steel. An extensive description for formulation constitutive rate equations of this model can be found in [START_REF] Doghri | Fully implicit integration and consistent tangent modulus in elastoplasticity[END_REF] and [START_REF] Mahnken | Improved implementation of an algorithm for non-linear isotropic/kinematic hardening in elastoplasticity[END_REF], and we only present here a summary of the resulting set of equations for the 1-D problem:

• Yield surface Φ = |σ -X| -(f y + R) = 0 (26) • Isotropic hardening stress R = R(p) = c 1 -e b p (27)
• Flow rule εp = λ sign(σ -X)

• Equivalent plastic strain evolution ṗ = λ (29)

• Kinematic hardening evolution Ẋ = k 1 εp -k 2 X ṗ ( 30
)
where b, c, k 1 and k 2 are material parameters. The numerical integration algorithm proposed by [START_REF] Mahnken | Improved implementation of an algorithm for non-linear isotropic/kinematic hardening in elastoplasticity[END_REF] is used and implemented in our FE code [START_REF] Nguyen | Modelling of the nonlinear behaviour of composite beams taking into account time effects[END_REF].

CONSTITUTIVE RELATIONSHIP FOR SHEAR CONNECTORS

The constitutive relationship proposed by [START_REF] Ollgaard | Shear strength of stud connectors in lightweight and normal weight concrete[END_REF] is considered for the stud shear connector. The analytical relationship between the shear force Q sc and the slip g of a generic stud is given by

Q sc = Q max 1 -e -c 1 |g| c 2 (31)
where Q max is the ultimate strength of the stud shear connector; c 1 and c 2 are coefficients to be determined from experimental results [START_REF] Aribert | Modèle calcul élasto-plastique de poutres mixtes à connexion partiell[END_REF].

It is worth to mention that the non-hysteretic behavior is adopted for the sake of simplicity. A hysteretic constitutive law of shear connection such as the one in [START_REF] Zona | Finite element response sensitivity analysis of steel-concrete composite beams with deformable shear connection[END_REF] can be easily implemented in the proposed FE model.

13

Nguyen, April 13, 2015

COMPOSITE BEAM ELEMENT FORMULATION

The mixed finite element method is used to formulate the composite beam element [START_REF] Zienkiewicz | The finite element method[END_REF]. The derivation follows the two-field formulation which uses the integral form of equilibrium and compatibility relations to derive the matrix relation between element generalized forces and corresponding displacements. Most of the following formulation is based on the presentation by Ayoub 2001. First of all, we recall in the following the field equations for composite beams with partial shear interaction in a small displacement setting. According to the notations defined in Fig. 5, these equations are summarized as

• Equations of equilibrium ∂D i -∂ sc D i sc -P z = 0 (32)
• Equations of compatibility ∂d ie i = 0 (33)

∂ T sc d i -g i = 0 (34)
• Equations of constitutive law

e i = e i-1 + f i-1 ∆D i (35) D i sc = D i-1 sc + k i-1 sc ∆g i (36)
where

D = [N c N s M c + M s ] T is
∂ =       d dx 0 0 0 d dx 0 0 0 - d 2 dx 2       ; ∂ sc =       1 -1 H d dx       (37)
In the two-field mixed formulation, independent shape functions are used for approximating the displacement and force fields along the element. The two fields are written

d(x) = a(x)q (38) D(x) = b(x)Q + D 0 (x) (39)
where a(x) and b(x) are the displacement and force interpolation matrices, respectively; D 0 (x) is a particular solution accounting for the effects of internal loading on the internal cross-section forces. The mixed beam formulation considered in this study is based on the Hellinger-Reissner variational principle and it is expressed by combining the weak forms of the compatibility equation ( 33) and equilibrium equation ( 32). These weak forms may be expressed as

L δD T ∂d i -e i dx + L δd T ∂D i -∂ sc D i sc -P z dx = 0 ∀ δD , δd (40)
where δD and δd are weighting functions fulfilling the equilibrium and compatibility conditions, respectively. Applying integration by parts for above equation and inserting ( 35), ( 36), ( 38) and (39) into the outcome, one obtains

δQ T G∆q i -F i-1 ∆Q i + q i-1 r + δq T K i-1 sc ∆q i + G T ∆Q i + Q i-1 r -P ext = 0 (41) 15
Nguyen, April 13, 2015

where

G = L b T (∂a) dx (42) F i-1 = L b T f i-1 b dx (43) q i-1 r = Gq i-1 - L b T e i-1 dx (44) K i-1 sc = L (∂ T sc a) T k i-1 sc (∂ T sc a)dx (45) Q i-1 r = G T Q i-1 + L (∂ T sc a) T D i-1 sc dx ( 46
)
P ext = Q ext + L a T P z dx - L (∂a) T D 0 dx (47)
in which Q ext is the vector of nodal applied loads. Since Eq. ( 41) must be hold for any δQ T and δq T , it follows that

   -F i-1 G G T K i-1 sc       ∆Q i ∆q i    =    -q i-1 r P ext -Q i-1 r    (48)
Eq. ( 48) represents in matrix form the two-field mixed formulation for the composite beam with partial interaction. If the first equation in ( 48) is solved for ∆Q i and the result is substituted into the second equation, the following expression results

K i-1 ∆q i = P ext -Q i-1 int ( 49
)
where

K i-1 = G T F i-1 -1 G (50) Q i-1 int = Q i-1 r + G T F i-1 -1 q i-1 r (51)
It should be noted that K i-1 designates the element consistent tangent stiffness at the 16 Nguyen, April 13, 2015

previous Newton-Raphson iteration and Q i-1 int represents the corresponding element resisting force vector. Further, the formulation has been developed by considering the composite beam element depicted in Fig. 6 where the rigid body modes are not included. Since small displacement hypothesis is adopted in this study, rigid body modes can be incorporated with a simple geometric transformation. The element has 7 displacement degrees of freedom and 6 force degrees of freedom as shown in Fig. 6.

State determination algorithm

The implementation of the model in a general purpose finite element analysis program requires the determination of the resisting forces and stiffness matrix that correspond to the displacements of the structural nodes. In a nonlinear structural analysis program each load step corresponds to the application of an external load increment to the structure. The corresponding structural displacement increments are determined and the end displacements are extracted for each element. The process of finding the resisting forces that correspond to the given displacements is known as state determination. The state determination process is made up of two nested stages: (a) the element state determination, in which the element resisting forces for given element displacements are determined; and (b) the structure state determination, in which the element resisting forces are assembled to yield the structure resisting forces. These are then compared with the total applied loads and the difference, if any, yields the unbalanced force vector. In the Newton-Raphson algorithm, the unbalanced forces are then applied to the structure until external loads and internal resisting forces agree within a specified tolerance [START_REF] Spacone | Mixed formulation of nonlinear beam finite element[END_REF]. Note that in the mixed formulation, there are several algorithms depending how the governing equations are linearized. A detailed study of state determination algorithms for nonlinear mixed beam finite elements can be found in [START_REF] Nukala | Variationally consistent state determination algorithms for nonlinear mixed beam finite elements[END_REF]. In this study, a direct element state determination following the proposal of [START_REF] Neuenhofer | Evaluation of nonlinear frame finite element models[END_REF] is used. In this algorithm, the residuals at the section and element levels are pushed to structure level and then they are eliminated during global N-R iterations. One the increment displacement at the element Nguyen, April 13, 2015 level ∆q i are obtained from the i th N-R iteration at the structure level, the element resisting forces are determined by the flowing steps:

• Calculate element force increments and update

∆Q i = F i-1 -1 G∆q i + q i-1 r (52)
Q i = Q i-1 + ∆Q i (53)
• Determine the section force increments (at integration points) and update the section generalized deformations

∆D i = b ∆Q i (54)
e i = e i-1 + f i-1 ∆D i (55)
• Update the slip and determine the interface resisting shear force as well as the consistent tangent shear stiffness through shear connector constitutive law (cf. Section 4)

d i sc = d i-1 dc + (∂ T sc a)∆q i (56)
D i sc = D sc d i sc (57)
k i sc = dD sc dd sc d i sc (58)
• Determine the section stiffness from the section constitutive law through fiber integration and deduce the section flexibility matrix f i (cf. see [START_REF] Nguyen | Analysis of composite beams in the hogging moment regions using a mixed finite element formulation[END_REF]).

• Calculate F i , q i r , K i sc and Q i r by Eqs. ( 43)-( 46) and then the new element stiffness matrix K i and element resisting force vector Q i int are calculated as in ( 50) and ( 51) respectively.

18 Nguyen, April 13, 2015

NUMERICAL EXAMPLE

The purpose of this application is to assess the capability of the proposed model to predict satisfactorily the long-term structural behavior of composite beams with partial interaction. We investigate the long-term deflection of two identical continuous beams (B1 and B2) with two equal spans. The beam B1 was subjected only to its self-weight, i.e. p 0 = 1.92 N/mm, while the beam B2 carried an additional superposed sustained load of 4.75 N/mm, i.e. p 0 = 6.67 N/mm. These continuous beams were tested by [START_REF] Gilbert | Time-dependent behavior of continuous composite beams[END_REF] over a period of 340 days. The dimensions of the tested beams, the loading and the geometric characteristics of the cross-section are shown in Fig. 7. The connector spacing is not given in [START_REF] Gilbert | Time-dependent behavior of continuous composite beams[END_REF] and we assume that the headed studs were uniformly spaced along the beam and the distance between two studs is taken equal to 55 mm. It is worth mentioning that in our model, the connection at the interface slab/joist is modeled by a continuous bond model despite the discrete nature of shear studs. The equivalent continuous connection stiffness is computed by dividing the stiffness of a single shear stud by their spacing along the beam. The concrete compressive strength obtained from the experiments is equal to 27 MPa. Accordingly, we can consider for the computer analysis that the concrete is of grade C30. Shrinkage is assumed to start at 7 days. The well-known B3 model [START_REF] Bažant | Creep and shrinkage prediction model for analysis and design of concrete structures -model b3[END_REF] is selected to predict the creep function and shrinkage strain of concrete.

Table 1 presents the constitutive model parameters which are used for the computer analysis.

Note that, in this table, all symbols are defined as in the corresponding cited reference.

The mid-span deflections obtained by the proposed model are compared against the experimental results and the numerical results obtained by viscoelastic model with cracked analysis [START_REF] Nguyen | Time-dependent analysis of composite beams with partial interaction based on a time-discrete exact stiffness matrix[END_REF] the concrete in hogging zone while in cracked analysis the concrete contribution along 15% of the span length on each side of the internal support is neglected [START_REF] Nguyen | Time-dependent analysis of composite beams with partial interaction based on a time-discrete exact stiffness matrix[END_REF]).

On the other hand, it is well known that the deflection of composite beams is strongly influenced by the shear connection stiffness which is not the same in their two computer analysis.

For example, the shear connection secant stiffness given by the exponential equation ( 31) is 1317 MPa for 1 mm of slip. This value is indeed much larger than the one taken in cracked analysis, ie 400 MPa (see [START_REF] Nguyen | Time-dependent analysis of composite beams with partial interaction based on a time-discrete exact stiffness matrix[END_REF].

In order to investigate the time effect of concrete in the plastic range, we consider that the beam is instantly subjected at 7 days to a distributed load p 0 . The loading is supposed to be constant during an analysis of 70 years. Figs. 9a and9b show the evolutions of the mid-span deflection in function of loading for t = 7 days and in function of time for p 0 = 20, 30 and 40 kN/m, respectively. Again, we observe that the deflection increases more or less significant over time. Further, the results indicated that the time effect is more important when the beam is in the nonlinear range (plasticized and cracked). For instant, with p 0 = 20 kN/m the beam is still in linear range (see Fig. 9a) and the creep and shrinkage make increase the deflection about 8 mm while when the beam is plasticized and cracked (with p 0 = 40 kN/m) the time deflection is about 26 mm. Fig. 10 presents the time evolution of the bending moment at intermediate support. As can be seen, the bending moment distribution is less important when the beam is more loaded. It can be explained by the fact that under p 0 = 20 kN/m the beam is not almost plasticized as under p 0 = 40 kN/m. Therefore the shrinkage makes it totally plasticized/cracked. As a result, there is an important bending moment redistribution compared to the case of 40 kN/m as we can see in Fig. 10.

The time evolution of the most tensile fiber stress is shown in Fig. 11. As can be seen, the instantaneous stress under service load p 0 = 6.67 kN/m (real load level of the test [START_REF] Gilbert | Time-dependent behavior of continuous composite beams[END_REF]) is remain smaller than the tensile stress resistance f t . The time effect makes 20 Nguyen, April 13, 2015

the tensile stress go up to f t at t = 10 days. This fiber is then "plasticized/cracked" therefore the softening curve of stress over time is obtained (see black line in Fig. 11). However, for the three other load levels one can observe that the tensile stress decreases over time.

That means that the slab cross-section at intermediate support is instantly cracked. The most compressed fiber stress at mid-span versus time curves are plotted in Fig. 12. Because of statically indeterminate beam, one can observe the stress relaxations for all load levels.

Further, it can be noted that the slab cross-section at mid-span partially compressed at t = 7 days can become totally tensed under time effect. It can be seen from Fig. 12 for the case p 0 = 6.67 kN/m. Fairly good agreement between the predicted deflections and the measured deflections has been found. The predicted deflections were also compared against the ones obtained by the cracked analysis in [START_REF] Nguyen | Time-dependent analysis of composite beams with partial interaction based on a time-discrete exact stiffness matrix[END_REF]). Good agreement has been observed. The numerical results have been indicated that the deflection can be significantly larger when the nonlinear behavior of concrete is combined with creep and shrinkage. Further, the time effect can produce an important bending moment redistribution in the statically indeterminate beam under serviceability loads. This redistribution becomes less significant if the beam is plasticized/cracked before the acting of time effect. However, further experimental research needs to be conducted to confirm this.

CONCLUSIONS

21

Nguyen, April 13, 2015

ξ 1 = E 2 c E 2 ci ε 2 c -4E 2 c E ci ε c f c + 4E 2 c f 2 c 4 (E ci ε c -2f c + f 2 c /(E c ε c )) 2 (59) ξ 2 = E c E ci ε c f c (E c ε c + E ci ε c -2f c ) -2E c f 2 c (E c ε c + f c ) 2 (E ci ε c -2f c + f 2 c /(E c ε c )) 2 (60) ξ 3 = E c ε 2 c f c (E c -E ci ) 2 (E c E ci ε 2 c -2E c ε c f c + f 2 c ) (61) ξ 4 = E c (E c E ci ε 2 c -2E c ε c f c + 2f 2 c ) 2 (E c E ci ε 2 c -2E c ε c f c + f 2 c ) (62) α 0 = 27E 2 c ε c γ c - 18E 3 c ε 2 c f c γ c -E 3 c ε 3 c (63) α 1 = 18E 3 c ε c f c γ c + 3E 3 c ε 2 c (64) α 2 = -3E 3 c ε c (65)
α 3 = E 3 c ( 66
)
β 0 = 6E 2 c ε c f c γ c -E 2 c ε 2 c ( 67
)
β 1 = 2E 2 c ε c (68) β 2 = -E 2 c (69) 22
Nguyen, April 13, 2015

f (MPa) c E (MPa) ε c (%) c G (kN/m) , eq c l (mm) t f (MPa) t G (N/m)

  3. Corrector: Determine ∆λ by solving Φ σ trialn+1 -E ve sign(σ trial n+1 )∆λ, p n + ∆λ = 0 (21)

  9 MPa; E c = 33500 MPa; ε c = 0.22%; G c = 30 kN/m; G t = 0.06 kN/m; and l eq = 21 mm. The comparison of the time evolution of the stress predicted by the proposed model against simulation with the linear viscoelastic model are presented in Fig. 3. It can be observed that, under creep and shrinkage effects, the stress relaxes from a negative value (compression) until it reaches the tensile strength at 300 days. At this stage, the plastic flow (cracking) begin which results in a decrease of the stress. This drop is necessary as the stress need to remain inside the yield surface. Fig. 4 depicts the time evolution of the strains.

  section generalized force vector; D sc the bond force per unit length; P z = [0 0 p z ] T with p z the external uniformly distributed load; d = [u c u s v] T the displacement vector; e = [ε c ε s κ] T section generalized deformation vector; g the relative slip between the concrete slab and the steel beam; superscript i indicates the current 14 Nguyen, April 13, 2015 iteration of the Newton-Raphson (N-R) iteration loop; f the flexibility matrix of the crosssection derived from uniaxial constitutive models for both steel and concrete. ∆ denotes the increments of the corresponding quantities; the operator ∂ and ∂ sc are given by

  in Fig.8for both beams. It can be seen that the numerical results are in reasonable agreement with the experimental data. Further, the deflection predicted by the proposed model is smaller than the one of cracked analysis. It can be explained on one hand by the fact that the proposed model takes into account the stiffness of 19Nguyen, April 13, 2015

  In this article, a FE model for the nonlinear time-dependant behavior of composite beams with partial shear connection has been proposed. Two fields mixed force-displacement formulation has been used. A viscoelastic/plastic has been developed for the time-dependent behavior of concrete. This model takes into account the interaction between creep, shrinkage and cracking of concrete. Creep and shrinkage are formulated via Kelvin chain springdashpot model of aging viscoelasticity. The concrete cracking is modeled using an 1D elasto-plastic model with softening proposed by the authors. The proposed FE model has been applied to two-span continuous composite beam to investigate the long term behavior.

  FIG. 2. Stress-strain diagram for uniaxial behavior

FIG. 5 .FIG. 9 .FIG. 10 .FIG. 11 .

 591011 FIG. 4. Time evolution of strains

TABLE 1 .

 1 Input values of the constitutive models

	c f ′	h	w	Creep and shrinkage B3 model c / a c m n	/ v s	E	(28)	Ciment Curing
	(MPa)	(%)	(kg/m 3 )	(kg/m 3 )				(mm)	(MPa)	type
	38	50	192	500	2.96	0.5	0.1	70	29182	2	moist
				Viscoelastic-plastic concrete model
	c


Haythem Gammoudi
email: haythemgammoudi@hotmail.fr

Fatma Belkhiria

Saloua Helali

Ben Assaker

Ibtissem Gammoudi

Fabien Morote

Ahmed Souissi

Mosbah Mokhtar Karyaoui

Mosbah Amlouk

Touria Cohen-Bouhacina

Radhouane Chtourou

Chemically grafted of single-walled carbon nanotubes onto a functionalized silicon surface

Keywords:

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

Introduction

Durant the last decade, the carbon nanotubes (CNTs) have emerged as new materials with a variety of potential applications [START_REF] Walters | Elastic strain of freely suspended single-wall carbon nanotube ropes[END_REF] [START_REF] Kobayashi | Extremely intense Raman signals from single-walled carbon nanotubes suspended between Si nanopillars[END_REF], and [START_REF] Tae | Fabrication of single walled carbon nanotube three-dimensional networks inside the pores of a porous silicon structure[END_REF]. This is due to their remarkable mechanical, thermal, optical and electrical properties.

However, many of these interesting and unique properties can only be realized once the CNTs are integrated into more complex assemblies: CNTs/Nano composites applications [START_REF] Valentinia | Modification of fluorinated single-walled carbon nanotubes with aminosilane molecules[END_REF] [START_REF] Sander | Individual single-wall carbon nanotubes as quantum wires[END_REF], and [START_REF] Collins | Extreme oxygen sensitivity of electronic properties of carbon nanotubes[END_REF], particularly in the field emission, solar cells, electronics domains and enhanced functional materials [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF] and [8]. The use of SWCNTs in the majority of these fields requires control over their placement and coverage on a surface, because such surfaces on CNTs play an important role to enhance their performance [START_REF] Yu | Application of polymer interlayers in silicon carbon nanotube hetero junction solar cells[END_REF]. Also, it is reported that pSi seems an ideal biomaterial due to its biodegradability, large internal surface area and biocompatibility, which leads to many applications such as: in biosensor design, tissue engineering and drug delivery [START_REF] Sweetman | Dual silane surface functionalization for the selective attachment of human neuronal cells to porous silicon[END_REF].

One particular approach for implementing these properties is the formation of two scale structured surfaces, micron scale structured surfaces that also possesses nanometer scale structure and the attachment of SWCNTs in this specific surface.

The first tries to attach SWCNTs in specific surface was in 2000 with the dicyclohexylcarbodiimide (DCC), mediated assembly of CNTs onto a surface when Liu et al. [START_REF] Liu | Organizing single-walled carbon nanotubes on gold using a wet chemical self-assembling technique[END_REF] attached shortened SWCNTs to gold. In such work, the amine group forms a cysteamine mole cule which is used to create an amide linkage to the carboxylic acid on SWCNTs. The remaining thiol group on cysteamine was then used to attach to the gold substrate. The same authors used the AFM to visualize the SWCNTs Au surfaces where they observed needle like protrusions which were described as SWCNTs perpen dicularly to the substrate [START_REF] Liu | Organizing single-walled carbon nanotubes on gold using a wet chemical self-assembling technique[END_REF]. This initial report was followed up in 2001 by another study in which they described the alignment of carboxylic acid functionalized SWCNTs on silver. Also, they re ported that the negatively charged carboxylic acid groups were electrostatically attracted to the positively charged silver surface [START_REF] Wu | Chemical alignment of oxidatively shortened single-walled carbon nanotubes on silver surface[END_REF]. Since that time, many attempts have been published, with CNTs assembled on gold using various amine terminated thiol monolayers with the predominant application of electrochemical sensing [START_REF] Diao | Vertically aligned single-walled carbon nanotubes by chemical assembly methodology properties and applications[END_REF], [START_REF] Gooding | Protein electrochemistry using aligned carbon nanotube arrays[END_REF].

Nevertheless, carbon nanotubes were not chemically attached to silicon until 2006 when Yu et al. [START_REF] Yu | Preparation and characterization of an aligned carbon nanotube array on the silicon (100) surface[END_REF] attached shortened SWCNTs to silicon by first hydrogenating a silicon (100) surface and then reacting with ethyl undecylenate to reach an ethyl undeca noate self assembled monolayer. The ester terminal of the mono layer was reduced to an alcohol whereupon the shortened carboxylated SWCNTs were covalently attached using DCC coupling. It was found from AFM observations that vertically aligned SWCNTs were attached to the surface. Cyclic voltammetry of the surfaces revealed that the SWCNTs to silicon electronic connection was not strong, indicating that the self assembled Monolayer (SAM) hindered electron transduction.

In 2007, Yu et al. [START_REF] Yu | Direct attachment of well-aligned single-walled carbon nanotube architectures to silicon (100) surfaces: a simple approach for device assembly[END_REF] refined chiefly the approach by hydrox ylating silicon (100) surfaces and then chemically attaching short ened SWCNTs using DCC coupling. The SWCNTs arrays formed were shown to be oriented vertically and cyclic voltammetry of the surfaces showed that the surface was suitable for electronic ap plications. In the same period, Flavel et al. [START_REF] Flavel | Patterned attachment of carbon nanotubes to silane modified silicon[END_REF] attached SWCNTs to an amine terminated silane layer (APTES) on silicon. This method produced well aligned SWCNT arrays using appropriate reaction between an amine and carboxylic acid to produce a stable amide linkage. Vertically aligned CNTs arrays produced by wet chemical assembly have been tested as chemical sensors [START_REF] Yu | Direct attachment of well-aligned single-walled carbon nanotube architectures to silicon (100) surfaces: a simple approach for device assembly[END_REF], biological sensors [START_REF] Gooding | Protein electrochemistry using aligned carbon nanotube arrays[END_REF], and photovoltaic devices [START_REF] Bissett | Photocurrent response from vertically aligned single-walled carbon nanotube arrays[END_REF] and [START_REF] Yu | Ruthenium porphyrin functionalized single-walled carbon nanotube arrays-a step toward light harvesting antenna and multibit information storage[END_REF] where they have shown promise as a simple method of producing functional surface bound CNTs arrays.

In 2010, Cameron et al. [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF] reported both the chemical attach ment of functionalized SWCNTs to amino silane modified pSi and the formation of silane patterns on pSi. In this context, their results show that pSi pore sizes and porosity increase across the substrate, SWCNTs coverage decreases concurrently, and their observations confirm that with decreasing APTES density on the surface, the coverage of the attached SWCNTs also decreases. In 2012, the same authors used pSi functionalized with APTES to immobilized carboxylated SWCNTs [8]. They found that the CNTs surface coverage noticeably altered the field emission properties of the surface. Also, their work covered the fact that there was an optimal coverage at which point there is enough CNTs to produce a large amount of current proving that CNTs surface coverage has also been found to have a profound effect on solar cell performance, cell adhesion and thin film transistors. By 2013, the same researchers found that both the chemical attachment and patterning of SWCNTs on pSi exhibit the ability of SWCNTs to capture human neuroblastoma cells [START_REF] Shearer | Nanostructured biointerfaces created from carbon nanotube patterned porous silicon films[END_REF].

The novelty in this work is indeed to modify the silicon topog raphy by a PPS structure [START_REF] Yooa | Black surface structures for crystalline silicon solar cells[END_REF] and [START_REF] Shia | Refined Nano-textured surface coupled with SiNx layer on the improved photovoltaic properties of multicrystalline silicon solar cells[END_REF] to increase the surface contacted with the carbon nanotube; and then attached chemically SWCNTs, on PPS surface and into the nano pores of PPS structure.

Indeed, the problem that persists is the implantation of carbon nanotubes on PPS surface and into the nano pores of PPS structure, which imposes to investigate the pore and pyramid size adapted to the size of SWCNTs.

Micron scale structured surfaces can being done by using anisotropic etching generally undertaken by using sodium hy droxide (NaOH), a pyramidal form of silicon was obtained during this chemical attack. Furthermore, a chemical etching time opti mization is necessary to obtain a good distribution of all pyramids along the surface of Si [START_REF] Seidel | Anisotropic etching of crystalline silicon in alkaline solutions[END_REF] [START_REF] Juvonen | High Efficiency single crystalline silicon solar cells[END_REF], and [START_REF] Fan | Differences in etching characteristics of TMAH and KOH on preparing inverted[END_REF].

Nanometer scale structure surfaces can being done by using electrochemical etching in HF solutions diluted with ethanol [START_REF] Xiao | Reflectivity of porous-pyramids structured silicon surface[END_REF] [27], and [START_REF] Marrero | Effect of porous silicon stain etched on large area alkaline textured crystalline silicon solar cells[END_REF] to forming the nano porous coating on wafers of Si pyramids pre textured.

In the present work, an efficient method is provided to control the pore size using a circular electrode that covers the entire sample surface to ensure the homogeneity of these pores and their sizes.

For the implantation of carbon nanotubes on the PPS, we chose the silanization of SWCNTs functionalized by carboxylic acid eCOOH with a coupling agent APTES to enhance the interfacial adhesion between nanotubes and the matrix of the PPS surface and to help to the dispersion and homogeneous implanting of the SWCNTs on the same surface [START_REF] Bethune | Cobalt catalysed growth of carbon nanotubes with single-atomic-layer walls[END_REF] and [START_REF] Popov | Carbon Nanotubes properties and application[END_REF]. This silanization method of SWCNTs allows us to control both the time and the suitable quantities used as well as concentration [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF].

Experimental details

Materials

Single wall carbon nanotubes (SWCNTs) (Fig. 1), containing carboxylic acid functionality were purchased from Sigma Aldrich. The diameter range of the SWCNTs was 4e5 nm, the length was 0.5e1.5 mm and the purity is above 90%. APTES (Fig. 1): with purity of 99% (Aldrich) was used for attached functionalized SWCNTs with a carboxylic acid on a PPS structure (Fig. 1). The following reagents were used without further purification.

Substrates preparation

P type (1 0 0) oriented silicon wafers were employed for the experiment. The resistivity of the samples was between 0.01 and 0.02 U cm. The size of samples was 20 Â 20 ±1 mm. First, the wafers were firstly washed by acetone, ethanol, deionized water and iso propanol. After that, these wafers were immersed in NaOH (1 M) solution at 85 C during 6 min to remove Si damage [START_REF] Huang | Key technique for texturing a uniform pyramid structure with a layer of silicon nitride on monocrystalline silicon wafer[END_REF] and [START_REF] Xi | Investigation of texturization for monocrystalline silicon solar cells with different kinds of alkaline[END_REF].

Pyramidal silicon structures were formed after chemical etching; the silicon wafers were rinsed with plentiful deionized water and dried by blowing Nitrogen gas. Second, porous pyramids structured silicon layers were prepared by electrochemical anod ization of heavily electrolyte solution of HF (25%) and EtOH, respectively. The anodization cell is a Teflon beaker where the c Si wafer acts as the anode and the cathode is formed by a platinum wire [3] [26], and [START_REF] Searson | The fabrication of porous silicon structures[END_REF]. Finally, a computer controlled current source was used to fabricate porous pyramids structured silicon from selected current density and time data, as shown in Table 1.

After the electrochemical attacks, all the samples were rinsed with ethanol and dried in the ambient air for the formation of an oxide layer SiO x [START_REF] Aggarwal | Porous silicon surface stability: a comparative study of thermal oxidation techniques[END_REF] [35] [START_REF] Peden | Thermally grown Si 3 N 4 thin films on Si (100) Surface and interfacial composition[END_REF], and [START_REF] Thogersen | An experimental study of charge distribution in crystalline and amorphous Si nanoclusters in thin silica films[END_REF].

Dispersion and implantation steps of SWCNTs were done on three steps. First step consists on the dispersion of SWCNTs in a 25 ml volumetric flask, carbon nanotubes (1 mg) were mixed with Dim ethylformamide (DMF, 10 ml). It was then sealed with a Teflon cap and dispersed via ultra sonication for 30 min until the mixture was visually homogeneous (Fig. 2).

Second step concerns the use of the coupling agent APTES which was prepared with different concentrations of 1% and 5% in the mixtures with ethanol [START_REF] Terracciano | Optical characterization of aminosilane-modified silicon dioxide surface for biosensing[END_REF]; this solution was made using a 100 ml Teflon beaker. Each sample of PPS structure emerged from 8 h in the various concentrations of APTES as shown in Table 1, then rinsing with ethanol to remove the APTES that has not been attached to the surface of PPS structure. This treatment leads to the formation of APTES network, which progressively cover the inside of the pores and the entire porous surface.

In third step, the solution of dispersed SWCNTs was deposited by spray method on different PPS structure substrates with different concentrations of APTES. The obtained samples were heated for 1 h at 120 C for drying and tailings disposal. The fabrication process of all the samples step is schematically shown in Graphical abstract (Fig. 3).

Characterization techniques

Fourier Transform infrared spectra were recorded in trans mittance mode at room temperature from 350 to 4000 cm À1 with a Fourier Transform infrared spectrometer (Perkin Elmer Spectrum 1000). The morphology of the samples was characterized with HITACHI S 4800 scanning electron microscopy (SEM) and Bioscope II, Bruker, Atomic Force Microscope (AFM). AFM is a high resolution scanning probe surface analysis technique, based on the measure ment of interaction forces between a fine tip and the sample sur face, these forces depending on the tip sample distance [START_REF] Nony | Dissipation induced by attractive interaction in dynamic force microscopy: contribution of adsorbed water layers[END_REF] and [START_REF] Simon | Study of two grafting methods for obtaining a 3-Aminopropyltriethoxysilane monolayer on silica surface[END_REF]. AFM images were recorded with a BioscopeII AFM setup (Veeco Brucker, Santa CA) equipped with a G scanner (maximum XYZ scan range of 150 mm  150 mm  12 mm). Samples were scanned in tapping mode using PPP NCL silicon probes (NANO SENSORS™) with a spring constant of about 35 N m À1 and a cor responding measured resonance frequency of about 165 kHz. All scans were done with scan rates between 0.3 Hz and 1 Hz (ac cording to the scan size and the scanning mode). AFM data were processed using the Nanoscope (version 7.30, Veeco).

Results and discussion

SEM morphological study

Fig. 3 shows the SEM results obtained on realized structure. Due to anisotropic etching of mono crystalline silicon with NaOH so lution, the whole wafer surface was homogeneously covered (85% coverage) by pyramidal structure of silicon (PS) (Fig. 3a).

Furthermore, the sizes of the pyramids are not uniform (2 mme7 mm), the silicon surface with an average pyramid size of 4.5 mm and the facets of the pyramids are all smooth.

This microstructuring of silicon provides an increase in the surface area with keeping the same size of samples, which make it an ideal biomaterial in terms of acceleration of their biodegrad ability [START_REF] Shearer | Nanostructured biointerfaces created from carbon nanotube patterned porous silicon films[END_REF].

Then the silicon wafers where electrochemically etched to form multilayered PPS on the pyramids. It is generally assumed that pore initiation occurs at surface active sites defects or irregularities. The properties of PPS layer (porosity, thickness, pore diameter and microstructure) mainly depend on the operating parameters of layer formation, including the HF concentration ratio, the presence of surface active compound additives, the duration and the tem perature as well as the wafer type and its resistivity [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF]. The in vestigations have been focused on visualization techniques, because the optical, electrical and physical properties of PPS layer depend strongly on its microstructure.

Fig. 3b, d displays the optical micro and nano imagery of PPS samples with low and high magnification. It can be seen that the surface of Si wafer with PPS shows black color (Fig. 3b) comparing to Fig. 3a, which indicates much lower reflectance of the silicon surface. When the magnifications are low, only pyramids could be observed (Fig. 3b). As the magnification increases, one can see that the porosification was already fabricated on every facet of pyra mids, which are distributed uniformly around the whole substrate, (Fig. 3d). The diameters of the PPS pores are identical and their sizes are around 20 mm as shown in the zoom of Fig. 3d (for the sake of clarity). So, it can be noted that the PPS multilayer is nano structured, and the pyramids remain unchanged after the electro chemical etching. Thus, they are clearly observed in the plan view (Fig. 3b). The obtained porosity was about 80% which is in good agreement with the result found elsewhere [START_REF] Shearer | Nanostructured biointerfaces created from carbon nanotube patterned porous silicon films[END_REF].

AFM morphological study

We used AFM to confirm and also to supplement the informa tion obtained by various other techniques such those of electron microscopy.

As shown in Fig. 4, the physico chemical treatment of silicon with APTES has been effective. Indeed, the silanized sample has a relatively rougher surface than that of the untreated substrate [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF] and [START_REF] Daniel | Aligned carbon nanotube thin films from liquid crystal polyelectrolyte inks; chemically grafted carbon nanotube surface coverage gradients[END_REF], the grain size is slightly larger and also note the pres ence of many aggregates. The average roughness is about 14 nm for the PPS structure (Fig. 4b) and 3 nm for the silanized surface (Fig. 4d).

Similarly, in the case of the study of the Chemical fixation of SWCNTs to an amine terminated silane (APTES), AFM results confirm those obtained by SEM. So, AFM imaging was completed for three samples at different concentrations of APTES, 0%, 1% and 5% to determine their influence on the bonding strength and the surface coverage. As shown in Fig. 5, the deposition of nanotubes for untreated surface is not very dense (Fig. 5 a,b,c). The contrast of the corresponding height image (Fig. 5 a,c) is not well resolved because the deposit nanotubes are not very dense [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF] and [START_REF] Daniel | Aligned carbon nanotube thin films from liquid crystal polyelectrolyte inks; chemically grafted carbon nanotube surface coverage gradients[END_REF], and more, the grain size of the silicon surface are of the same order of magnitude as those of the nanotubes. The "derived signal image" (Fig. 5b) highlights the sub structures such as nanotubes or grains of the reference substrate and gives also a "3D view" of the topography.

As shown in Fig. 5 e, f, g one can make the same observation for the surface silanized with a low APTES concentration (1%). In this case, the nanotubes are well deposited but with no preferred orientation.

However, for intermediate concentration of APTES (between 1 and 5%), the surface can be completely covered can be well covered but the nanotubes are not well attached. But, once we begin to scan the surface, the tip gets polluted, possibly by clinging to the tip SWCNTs which is causing duplication of nanotubes observed in the corresponding AFM images and that complicates the character ization (results not shown here). On the other hand, the contrast is striking for the concentration of APTES of 5%, the SWCNTs density is high and the surface is completely covered by well deposited entangled nanotubes as observed in Fig. 5 g,h, i So, it is obvious that treatment with silane optimizes the recovery by the nanotubes and ensures their good adhesion. This finding regarding the implanta tion of CNTs is consistent with the one reported by Cameron et al. [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF].

We conclude that the increase in surface coverage is tributary to an increase in reactive surface sites and the formation of the more favorable amide bond, while the entangled is most likely due to the less hydrophilic APTES functionalized surface [START_REF] Shearer | Chemically grafted carbon nanotube surface coverage gradients[END_REF] [10], [START_REF] Shearer | Nanostructured biointerfaces created from carbon nanotube patterned porous silicon films[END_REF] [41] [START_REF] Shearer | Preparation and characterisation of vertically aligned single-walled carbon nanotube arrays on porous silicon[END_REF], and [START_REF] Constantopoulos | Carbon nanotubes anchored to silicon for device fabrication[END_REF].

FTIR structural study

Fig. 7 shows FTIR spectra of the realized structures. To increase the coverage of SWCNTs on the silicon surface, a PPS surface is synthesized in which an oxide layer SiO 2 was formed during the drying time in ambient air (Fig. 6a), because the instability of SieH bonds, a layer of APTES was deposited to produce an amine terminated monolayer (Figs. 3, 4b and6) that reacts with the car boxylic acid attached at SWCNTs (Figs. 3 and6) that forms (or leading to) an amide layer (Figs. 3, 4c and6). In the case of Fig. 7a, the PPS surface displays characteristic SieH vibrational modes at [START_REF] Mawhinney | FTIR study of the oxidation of porous silicon[END_REF], and [START_REF] Sweetman | Dual silane surface functionalization for the selective attachment of human neuronal cells to porous silicon[END_REF], A weak oxidation peak (silanol) assigned to 974 cm À1 is attributed to SieOH groups which facili tates the fixation of APTEs on the PPS surface [START_REF] Torabi | Study of physical, chemical and morphological alterations of smectite clay upon activation and functionalization via the acid treatment world[END_REF] and a peak at 874 cm À1 is attributed to O3SieH wagging mode.

Once the APTES surface has been filed (Fig. 7b), we observe the degradation of SieH x bonds due to their instability and a slight increase in the oxidation of SieOeSi during the time. Moreover, a broad peak appearing at 900 cm À1 is attributed to NeH 2 scissoring instead of the low oxide peak SiOeH at 974 cm À1 , and the one located at 1550 cm À1 corresponds to NeH 2 bending [START_REF] Stuart | Infrared spectroscopy: fundamentals and applications[END_REF]. In the same line, three small primary amine peaks assigned at 790, 1250 and 1650 cm À1 are related to the urethane bond within the APTES [START_REF] Michael | Exciton-mediated hydrosilylation on photoluminescent nanocrystalline silicon[END_REF] [50], and [START_REF] Rosu | Cascaval, IR-change and yellowing of polyurethane as a result of UV irradiation[END_REF] On the other hand, there is also a broad vibra tional peak around 3010 cm À1 attributed to the terminal CH 3 group of the APTES. Also, it is found that the silicon oxide (SieOeSi) stretching (bending) vibrational peak at 1150 cm À1 still exists but with low intensity.

We deduce then, that after the filing of APTES there is the appearance of the amine groups in the PPS surface that promotes binding of SWCNTs. The final FTIR spectrum displayed in Fig. 7c corresponds to SWCNTs deposited in PPS with APTES surface. Upon SWCNTs attachment (Fig. 7c), the amine peaks at around 1250, 1550 and1650 cm À1 disappear, CeH stretching peaks remain and two new peaks appear at 1200 cm À1 and 2945 cm À1 which are attributed to the carboxylic acid groups (O]OH) attached on the surface of SWCNTs (Fig. 1). The large bands at 790 and 850 cm À1 with low intensity is also indicates that the chemical functional amine (NH 2 ) groups associated with each APTES are reserved after deposition of SWCNTs which returns that the amount of amine does not react with all the acid carboxylic. Two new bands appeared at 1200 cm À1 and 3400 cm À1 corresponding to the amide bond thus confirms the reaction of amine with carboxylic acid group.

The FTIR spectra support the formation of an amine terminated monolayer and the creation of an amide linkage between the APTES and functionalized SWCNTs (Fig. 4), that confirms the successful implantation and dispersion of each SWCNTs into the PPS APTES surface. The higher coverage of SWCNTs observed, (Fig. 3e), is due to the favorable and stable amide linkage as well as the fact that the APTES layer does not include any germinal groups further increasing the stability of the attachment, and the APTES mediated attachment was considered the best method to produce SWCNTs PPS substrate.

Conclusions

In this work, we successfully constructed a composite structure combining pyramids, nano porosification, APTES and carbon nanotubes has been manufactured by both chemical and electro chemical etchings. Nano pores are formed by changing the time and the current of the electrochemical etching of a Si substrate structured pyramids. The optimal condition for the manufacture of this structure is obtained with a current of 60 mA and 35 s etching time. Single walled carbon nanotubes functionalized by carboxylic acid was silanized and implanted in two stages: First step consists of the surface preparation by direct immersion of PPS structure in APTES solution at room temperature and ambient pressure. Second step leads to the implantation of SWCNTs by spray method on PPS structure. Validation of PPS structure created and the chemical attachment of SWCNTs into this structure by APTES molecules were confirmed by AFM, SEM and FTIR.

Further investigation on i (v) characteristics as well as the thermal properties of such structures is in progress. Finally, this new combined structure may be of interest for various sensitivity applications such as gas sensors, waste water, biosensors and so on also, this study paves the way for possible use of this structure as chemical sensor functionalized by biological entity.

Fig. 1 .

 1 Fig. 1. Schematic procedure: (step1) Formation of pyramids on silicon surface; (step2) Porousification of pyramids on silicon surface; (step3) Deposition of APTES in the surface of silicon porous structured pyramidal; (step4) Deposition of carbon nanotubes on silicon porous structured pyramidal surface with APTES. Substrates preparation.

Fig. 3 .

 3 Fig. 3. SEM Images of: (a) Pyramidal silicon structure (PS), (b) Porous pyramidal silicon structure (PPS), (c) Porous pyramidal silicon structure with APTES (APPS), (d) Porous pyramidal silicon structure with a low-scale view and a large-scale view (PPS), (e) Porous pyramidal silicon structure with APTES and SWCNTs with a low-scale view and a largescale view (CNAPPS).

Fig. 4 .

 4 Fig. 4. AFM Height image of: (a,b) Porous silicon structure (PPS); (c,d) Porous silicon structure with APTES (APPS).

Fig. 5 .

 5 Fig. 5. AFM image of SWCNTs deposited on: (a,b,c) Structured porous silicon non functionalized with APTES; (d, e, f) Structured porous silicon functionalized with APTES -Low concentration (1%); (g, h, i) Structured porous silicon functionalized with APTES -High concentration (5%); (a, c, d, f, g, i) Height images; (b, e, h) derived signal.

Fig. 6 .

 6 Fig. 6. Schematic representation of (a) oxidized Si (b) oxidized Si after amine immobilization and (c) oxidized Si after amine immobilization and carboxylic acid attachment.

Fig. 7 .

 7 Fig. 7. FTIR spectra of: (a) Pores silicon structured pyramidal (PPS); (b) Pores silicon structured pyramidal with APTES (APPS), (c) Porous pyramidal silicon with APTES and SWCNTs structure (CNAPPS).

Table 1

 1 Preparation parameters of the porous pyramidal silicon structures and concentration of APTES and SWCNTs.

	Samples	Chemical etching (min)	Electrochemical etching		Concentration of APTES (%)	Concentration of SWCNTs (M)
			Current (mA)	Time (S)
	PS	6
	PPS	6	60	35
	APPS1	6	60	35	1
	APPS5	6	60	35	5
	CNAPPS1	6	60	35	1	1
	CNAPPS5	6	60	35	5	1

Fig. 2. Dispersion of SWCNTs in a volumetric flask.

Acknowledgments

The authors gratefully acknowledge Dr. Guy Garry for giving us access to the Transmission Electron Microscope at Thales Research and Technology France, Palaiseau, France.