Le-Hung Tran

Tien Hoang

Gilles Foret

Denis Duhamel

Dinh-Duc Nguyen

Hung Le
email: hungtl@vnu.edu.vn

Tran

Calculation of dynamic responses of railway sleepers on a nonlinear foundation

Keywords: Railway track, railway sleeper, nonlinear foundation, harmonic balance, structural dynamic

The safety of a passing train depends on different factors, of which one of the most important is the behavior of the foundation. Therefore, the effects of the non-linearity of ballast on the dynamic responses of the railway track are a key research interest. In this paper, a new model of railway sleepers posed on a nonlinear foundation has been developed. By coupling the finite element method (FEM) of the sleeper with an analytical model of the periodically supported beam model, the dynamic equation of the sleeper is developed. On the other hand, by considering a periodic series of moving loads, this equation can be transformed to a forced nonlinear oscillation. Iteration procedures have been built to calculate the periodic solution. This method has demonstrated a good convergence of results by comparison with the analytical solution in the linear case. The influence of the nonlinear foundation has been investigated by two examples: cubic-nonlinear and bi-linear foundations. The parametric studies demonstrate that numerical results converge with a small number of iterations.

Introduction

The stability of the railway track is one of the research interests of railway engineers as it is key to the safety of the passing trains. Consequently, many research projects have been undertaken to analyze the dynamic responses of the railway track subjected to moving loads. Firstly, analytical models have been developed for an infinite beam placed on a continuous foundation. Krylov [START_REF] Kriloff | Uber die erzwungenen schwingungen von gleichformigen elastischen staben[END_REF] and Timoshenko [START_REF] Timoshenko | Erzwungene schwingungen prismatischer stäbe[END_REF] are the pioneering researchers in this domain, in which the dynamic stresses are calculated for simply supported beams under moving loads. By using the same method, Fryba [START_REF] Fryba | Vibration of solids and structures under moving loads[END_REF] studied the transverse vibrations of a beam posed on an elastic foundation, in which the beam responses are calculated by the sum of multiple normal modes. The free vibrations of the beam on an elastic foundation has been solved analytically by Timoshenko et al. [START_REF] Timoshenko | Vibration Problems in Engineering[END_REF]. Analytical and numerical methods have been developed for various foundations such as: Winkler, Pasternak, Vlasov or Reissner [START_REF] Wang | Natural frequencies of Timoshenko beams on pasternak foundations[END_REF][START_REF] Adams | Critical speeds and the response of a tensioned beam on an elastic foundation to repetitive moving loads[END_REF][START_REF] Sun | A closed-form solution of a bernoulli-euler beam on a viscoelastic foundation under harmonic line loads[END_REF][START_REF] Chen | Dynamic characteristics of infinite and finite railways to moving loads[END_REF][START_REF] Mallik | Steady-state response of an elastically supported infinite beam to a moving load[END_REF][START_REF] Firat | Dynamic analysis of beams on viscoelastic foundation[END_REF][START_REF] José | Dynamic behaviour of railway tracks on transitions zones[END_REF][START_REF] Ichikawa | Vibration analysis of the continuous beam subjected to a moving mass[END_REF]. In order to take into account the distribution of discrete supports, Mead [START_REF] Mead | Wave propagation in continuous periodic structures: research contributions from southampton[END_REF][START_REF] Mead | Free wave propagation in periodically supported, infinite beams[END_REF] developed the periodically supported beam model under moving loads for the rail. This model type has been also investigated for elastic foundations [START_REF] Metrikine | Vibration of a periodically supported beam on an elastic half-space[END_REF][START_REF] Vostroukhov | Periodically supported beam on a visco-elastic layer as a model for dynamic analysis of a high-speed railway track[END_REF][START_REF] Belotserkovskiy | On the oscillation of infinite periodic beams subjected to a moving concentrated force[END_REF][START_REF] Tran | Calculation of the dynamic responses of a railway track on a non-uniform foundation[END_REF]. A fast analytical method to calculate the dynamic responses of railway sleepers has been presented by Tran et al. [START_REF] Tran | Analytical model of the dynamics of railway sleeper[END_REF][START_REF] Tran | A fast analytic method to calculate the dynamic response of railways sleepers[END_REF]. Most recently, in order to study the influence of non homogeneous foundation on the sleeper responses, Tran et al. [START_REF] Tran | Influence of nonhomogeneous foundations on the dynamic responses of railway sleepers[END_REF] developed a semi-analytical model by coupling a numerical model for the sleeper and a periodically supported beam model for the rail.

The two models which are frequently used to describe the nonlinear behaviors of foundation are cubic-nonlinear and bi-linear foundations. The cubic law describes a foundation reaction which has a cubic dependence on displacement. The bi-linear foundation presents a foundation behavior which has two different linear behaviors: one in compression and one in tension. By taking into account the non-linearity of the foundation on the model of the beam placed on a continuous foundation, the solutions have been calculated by using different techniques. Ding et al. [START_REF] Hu Ding | Convergence of Galerkin truncation for dynamic response of finite beams on nonlinear foundations under a moving load[END_REF][START_REF] Hu Ding | Dynamic response of an infinite Timoshenko beam on a nonlinear viscoelastic foundation to a moving load[END_REF] solved the problem with the help of the Galerkin discretization and a fourth-order Runge-Kutta method. Abdelghany et al. [START_REF] Abdelghany | Dynamic response of non-uniform beam subjected to moving load and resting on non-linear viscoelastic foundation[END_REF], Kargarnovin et al. [START_REF] Kargarnovin | Response of beams on nonlinear viscoelastic foundations to harmonic moving loads[END_REF] and Ansari et al. [START_REF] Ansari | Frequency analysis of finite beams on nonlinear Kelvin-Voight foundation under moving loads[END_REF] used the same method to calculate the responses of the beam under a harmonic load or to calculate the resonance of the beam. Chen et al. [START_REF] Hong Yan Chen | Convergent term of the Galerkin truncation for dynamic response of sandwich beams on nonlinear foundations[END_REF][START_REF] Hong Yan Chen | The Scheme to determine the convergence term of the Galerkin method for dynamic analysis of sandwich plates on nonlinear foundations[END_REF] studied the convergence of Galerkin truncation for the sandwich beam on a nonlinear foundation and a scheme to determine the convergence of this model has been presented. Recently, Ouzizi et al. [START_REF] Ouzizi | Nonlinear dynamics of beams on nonlinear fractional viscoelastic foundation subjected to moving load with variable speed[END_REF] presented a model of the dynamic responses of the beam on a nonlinear frictional viscoelastic foundation with the help of an explicit scheme. The dynamic track responses are computed by using the finite element method [START_REF] Nguyen | Finite element procedures for nonlinear structures in moving coordinates. Part 1: Infinite bar under moving axial loads[END_REF][START_REF] Nguyen | Finite element procedures for nonlinear structures in moving coordinates. part ii: Infinite beam under moving harmonic loads[END_REF][START_REF] Castro Jorge | Finite element dynamics analysis of beams on non-uniform nonlinear viscoelastic foundations under moving loads[END_REF]. A numerical method has been also developed to calculate the solutions of the beam in the case of tensionless foundation [START_REF] Chen | Steady state and stability of a beam on a damped tensionless foundation under a moving load[END_REF][START_REF] Sapountzakis | Nonlinear response of shear deformable beams on tensionless nonlinear viscoelastic foundation under moving loads[END_REF][START_REF] Yang | An explicit periodic nonlinear model for evaluating dynamic response of damaged slab track involving material nonlinearity of damage in high speed railway[END_REF]. The nonlinear responses of beams subjected to moving loads carrying multiple mass-spring-damper attachments has been also studied [START_REF] Samani | Performances of dynamic vibration absorbers for beams subjected to moving loads[END_REF][37][START_REF] Bukhari | Nonlinear vibrations analysis of overhead power lines: a beam with mass-spring-damper-mass systems[END_REF][START_REF] Hoang | Response of a periodically supported beam on a nonlinear foundation subjected to moving loads[END_REF][START_REF] Tran | Calculation of the dynamic responses of rails subjected to moving loads on ballasted railway track[END_REF][START_REF] Tran | Calculation of dynamic responses of a cracked beam on visco-elastic foundation subjected to moving loads, and its application to a railway track model[END_REF]. However, the effect of the non-linearity of the foundation is not easy to take into account in the periodically supported beam model.

In this paper, we present a novel model to calculate the solution of dynamic behavior of a railway sleeper posed on a nonlinear foundation. In Section 2, the dynamic equation of the sleeper is written with the help of a finite element model. When the rails are modeled as a periodically supported beam [START_REF] Hoang | Calculation of force distribution for a periodically supported beam subjected to moving loads[END_REF][START_REF] Hoang | Dynamical response of a timoshenko beams on periodical nonlinear supports subjected to moving forces[END_REF], a relation between the rail displacements and the reaction force in the frequency domain can be obtained and it holds for any foundation behavior. Therefore, the forces applied by the rails on the sleeper can be rewritten in the dynamic equation of the sleepers with the help of the Dirac delta function. By considering that the moving loads are a periodic series, this equation is equivalent to a forced nonlinear oscillator. Then, a numerical method is developed with the help of the harmonic balance method and the iteration procedures for nonlinear oscillators [START_REF] Ali | Introduction to perturbation techniques[END_REF][START_REF] Nayfeh | Nonlinear oscillations[END_REF] in Section 3. The numerical applications are shown in Section 4. In the linear case, the numerical results converge to the analytical solutions [START_REF] Tran | A comparison of beam models for the dynamics of railway sleepers[END_REF]. Thereafter, the dynamic responses of the sleepers are calculated for two types of nonlinear foundation: cubic law and bi-linear foundation. The applications show that the numerical results converge rapidly to the solution of the problem. In addition, parametric studies have been conducted to analyze the influence of nonlinear parameters on the sleeper responses. Finally, conclusions are drawn in Section 5.

Formulations

Let us consider the railway track shown in Fig. 1. In this track, the rails k (with k = [START_REF] Kriloff | Uber die erzwungenen schwingungen von gleichformigen elastischen staben[END_REF][START_REF] Timoshenko | Erzwungene schwingungen prismatischer stäbe[END_REF]) are subjected to moving loads Q (k) j which are characterized by the distance to the first wheel D j (1 ≤ j ≤ K where K is the total number of axles). The mutual distance D j of the moving loads is not restricted to be constant. The rails are periodically supported by the sleepers and each one is separated by a length l. Next, a novel dynamic model of the sleeper is developed by taking into account the non-linearity of the foundation.

Sleeper posed on a foundation

Fig. 2 presents a finite element model of the sleeper posed on a non linear foundation. In this model, we use the beam element in 2D where each node has 3 degrees of freedom (DOFs) which correspond to the 2 displacements (u x , u z ) and 1 rotation (θ y ). The sleeper is subjected to vertical forces at two nodes: R 1 and R 2 which are the two rail positions. If n n is the number of nodes, the number of DOFs is 3n n . We note that the main dynamic response of the system is in the vertical direction. Each node of the sleeper is related with the foundation by a nonlinear system which describes the behavior of the foundation in two parts: linear and nonlinear. The linear term is modelled by a stiffness k f and damping coefficients ζ f . Meanwhile the nonlinear term of the foundation is written in the function f N L (u, u).

y x z l Q (2) j Q (1) j Fig. 1: Ballasted railway track ζrp Q (2) j krp Q (1) j k f ζ f f N L Middle nodes Right nodes Left nodes Node R 2 Node R 1 x z
Fig. 2: Finite element model of the railway sleeper posed on a nonlinear foundation and subjected to the moving loads

The dynamic equation of the sleeper posed on a foundation is obtained by using a finite element method as follows:

M s ü + K s u + F f (u) = F R (t) (1)
where M s and K s represent respectively the mass and rigidity matrices of the sleeper which are generated with the help of FEM where the notation ˙ denotes the partial derivative with respect to time t. u is a vector of nodal displacements of sleepers and:

u i = u(x i ) = [u x (x i ), u z (x i ), θ z (x i )] T
for the i th node of a sleeper. The vector F R (t) describes a force applied on the sleeper for the two rail positions (see Fig. 2). The vector F f (u) represents the force induced by the foundation on the sleeper and we suppose that this force follows only the vertical direction. For this reason, F f (u) can be calculated as follows:

F f (u) = k f I uz u + ζ f I uz u + I uz f N L (u, u) (2)
The matrix I uz is calculated as follows:

I uz = e uz ⊗ e uz
where the symbol ⊗ denotes the tensor product. The vector e uz which has 3n n components is defined as follows:

e uz = 0 1 0 0 1 0 • • • 0 1 0 T 2.

Rail pad

A rail pad which is located between the rails and the sleeper, can be modeled as a spring-damper system with stiffness k rp and damping coefficient ζ rp . The reaction force of the sleeper to the rail k is expressed as follows:

R k (t) = -ζ rp ẇ(k) r (t) -uR kz (t) -k rp w (k) r (t) -u R kz (t) (3)
where w

(k) r (t) and u R kz (t) are respectively the rail and sleeper vertical displacement at the contact positions between the rail k and sleeper in time domain. So that, we can deduce the expression of the reaction force applied vertically on the beam at two rail positions:

F R (t) = -R 1 (t)e R1 -R 2 (t)e R2 (4)
The vectors e R1 , e R2 are two column vectors which each have 3n n elements. The two vectors are zero everywhere, except at the positions which correspond to the vertical displacement of the two nodes R 1 and R 2 where they have the value 1. Moreover, by using the Fourier transform, Eq. ( 3) can be rewritten in the frequency domain:

Rk (ω) = -k p ŵ(k) r (ω) -ûR kz (ω) (5)
where k p = k rp + iωζ rp is the dynamic stiffness of the rail pad, ω is the angular velocity and i 2 = -1. ŵ(k) r and ûR kz are respectively the rail k and sleeper vertical displacements at the crossing point in the frequency domain. By inserting Eq. ( 23) into Eq. ( 5) (see B), we obtain the following results:

       ŵ(k) r (ω) = k p ûR kz (ω) -Q k k p + K e Rk (ω) = k p k p + K e [K e ûR kz (ω) + Q k ] (6)

Dynamic equation of the railway sleeper

By substituting Eqs. ( 2) and ( 4) into Eq. ( 1), the dynamic equation of the sleeper under a moving load and resting on the nonlinear foundation can be written as follows:

(7) Mü + C u + Ku + f N L (u, u) = k rp w (1) r (t) + ζ rp ẇ(1) r (t) e R1 + k rp w (2) r (t) + ζ rp ẇ(2) r (t) e R2
where w

r (t) denote respectively the displacements of the rails 1 and 2 at the two contact points with the sleeper. M, C, K are the mass, damping and rigidity matrices, which represent the linear part of the model. Beside, the term f N L (u, u) is a vector that describes the non-linear behavior of the model. These matrices are calculated as follows:

     K = K s + k f I uz + k rp I R1 + k rp I R2 C = ζ f I uz + ζ rp I R1 + ζ rp I R2 M = M s
where matrix I R1 is the null matrix with size 3n n ×3n n , except the component which corresponds to the position (R 1 , R 1 ) get the value 1, same for matrix I R2 . In other word, the two matrices are calculated by:

I R1 = e R1 ⊗ e R1 I R2 = e R2 ⊗ e R2
The right side of Eq. ( 7) depends on the rail displacement. By substituting Eq. ( 6) into Eq. ( 7), the dynamic equation of the railway sleeper can be rewritten as a function of the sleeper displacement as follows (see A): ( 8)

Mü + C u + Ku + f N L (u, u) = I R   1 2π +∞ -∞ k 2 p û(ω) k p + K e e iωt dω   -   1 2π +∞ -∞ k p Q 1 (ω) k p + K e e iωt dω   e R1 -   1 2π +∞ -∞ k p Q 2 (ω) k p + K e e iωt dω   e R2
where I R = I R1 + I R2 . This equation describes the dynamic responses of the railway sleeper posed on the nonlinear foundation and it is similar to that of a nonlinear oscillator. The right side of the equation represents the interaction between the rails and the sleeper at the two rail positions. In next step, in order to simplify these terms, we consider only periodic solutions when the moving forces are a periodic series.

Periodic series of moving loads

We consider that the train contains many identical wagons as shown in Fig. 3). The distances D j of each wheel are characterized by: Lw Hw Dw Dw v Fig. 3: Diagram of the periodic series of moving loads

D j =          jH w
for wheel 1 of wagon jH w + D w for wheel 2 of wagon jH w + D w + L w for wheel 3 of wagon jH w + 2D w + L w for wheel 4 of wagon [START_REF] Mallik | Steady-state response of an elastically supported infinite beam to a moving load[END_REF] where D w is the distance of the bogie, H w is the length of the wagon and L w is the distance between wheels 2 and 3 of a wagon (see Fig. 3). This series of moving loads may be used to represent the series of charges for a railway track. In addition, we consider that the loads

Q (k) j
of each wheel j on each rail k are equal:

Q (k) j = Q (10)
By considering infinite periodic series (j ∈ Z), we will use the periodicity of this series to reduce the terms on the right side of Eq. ( 8) and we obtain the following result (see C):

Mü + C u + Ku + f N L (u, u) = I R   +∞ j=-∞ Φ j P j e iωj t   -   +∞ j=-∞ F (1) j e iωj t   e R1 -   +∞ j=-∞ F (2) j e iωj t   e R2 (11)
In the previous equation, we supposed that the solution of the nonlinear problem can admit the same frequencies as the excitation force and the periodical solution of u(t) which can be represented as shown in the Eq. ( 30) (see C). F

(1) j and F

(2) j

represent the two forces applied on the two rails which are normally different in general.

Solution of the problem

In the previous section, the dynamic equation of a railway sleeper has been developed as shown in Eq. [START_REF] José | Dynamic behaviour of railway tracks on transitions zones[END_REF]. We remark that the right side of this equation is the sum of several infinite series. In this section, we present a method to solve this equation with the help of the harmonic balance method and the iterative procedure [START_REF] Mickens | A generalization of the method of harmonic balance[END_REF][START_REF] Mickens | A generalized iteration procedure for calculating approximations to periodic solutions of truly nonlinear oscillators[END_REF][START_REF] Hoang | Frequency dependent iteration method for forced nonlinear oscillators[END_REF] for nonlinear oscillators to develop a numerical method for the dynamic equation of sleepers. Firstly, by using the Fourier development of Eq. ( 11), we obtain the following result:

1 T T /2 -T /2 [Mü + C u + Ku + f N L (u, u)] e -iωj t dt = I R Φ j P j -F (1) j e R1 + F (2) j e R2 (12)
By inserting Eq. ( 30) into the last result, this equation can be reduced as follows:

(13) -D j Φ j + 1 T T /2 -T /2 f N L (u, u)e -iωj t dt = I R Φ j P j -F (1)
j e R1 + F (2)
j e R2
where D j is the dynamic stiffness matrix of the model which is calculated by

D j = ω 2 j M -iω j C -K and:              u(t) = ∞ j=-∞ iω j Φ j e iωj t ü(t) = ∞ j=-∞ -ω 2 j Φ j e iωj t
We remark that Eq. ( 13) is the harmonic balance form of Eq. ( 12). The set of this equation for all j ∈ Z generates a system of equations with regard to Φ j which are to be determined. The railway sleeper is modeled by the finite element method and we note that the index L , M , R , R1 and R2 are respectively the degrees of freedom at the left, middle, right of the sleeper and positions of the forces 1 and 2 (see Fig. 2). In the case of linearity of the foundation, which means that f N L (u, u) = 0, Eq. ( 13) can be explained as follows:

(14)       D j LL D j LR 2 D j LM D j LR 1 D j LR D j R 2 L Dj R 2 R 2 D j R 2 M D j R 2 R 1 D j R 2 R D j M L D j M R 2 D j M M D j M R 1 D j M R D j R 1 L D j R 1 R 2 D j R 1 M Dj R 1 R 1 D j R 1 R D j RL D j RR 2 D j RM D j RR 1 D j RR             Φ L j L Φ L j R 2 Φ L j M Φ L j R 1 Φ L j R       =        0 F (2) j 0 F (1) j 0       
where Φ L j is the solution in the case of a linear foundation and:

Dj R 1 R 1 = D j R 1 R 1 + P j Dj R 2 R 2 = D j R 2 R 2 + P j
The solution of Eq. ( 14) can be calculated easily for all j ∈ Z as follows:

(15)

      Φ L j L Φ L j R 2 Φ L j M Φ L j R 1 Φ L j R       =       D j LL D j LR 2 D j LM D j LR 1 D j LR D j R 2 L Dj R 2 R 2 D j R 2 M D j R 2 R 1 D j R 2 R D j M L D j M R 2 D j M M D j M R 1 D j M R D j R 1 L D j R 1 R 2 D j R 1 M Dj R 1 R 1 D j R 1 R D j RL D j RR 2 D j RM D j RR 1 D j RR       -1        0 F (2) j 0 F (1) j 0       
Here now, we denote the right term of Eq. ( 14) as the vectors F j . For the non-linear foundation, the term f N L (u, u) appears and we will use an iteration procedure to solve this problem. Eq. ( 13) can be rewritten as follows:

D j Φ j = F j + 1 T T /2 -T /2 f N L (u, u)e -iωj t dt (16)
The iteration procedure considers the n first harmonics of the periodic solution, thus, the nodal displacement of sleeper can be determined as:

u nm (t) = n j=-n Φ m j e iωj t (∀m ≥ 1) (17)
By taking Φ 1 j = 0 ∀j, we built a series of vector {Φ m j } such that Φ m j → Φ j when m, n → ∞. The index m is understood as the number of the iteration procedures, which can be defined as follows:

D j Φ m+1 j = F j + Q m j (18)
where Q j represents a vector of non-linear terms which is calculated by:

Q j = 1 T T /2 -T /2 f N L (u, u)e -iωj t dt
The series of vectors {Φ m j } can be determined from Eq. ( 18) as follows:

Φ m+1 j = D -1 j F j + Q m j (19)
The last equation defines recurrent sequences of the vectors {Φ j } with regard to m. In the case of convergence of these sequences for all j ∈ Z when m, n → ∞, Φ m j and Φ m+1 j are replaced by their limit, which we can find once again with Eq. ( 13). It means that these sequences converge to the solution of this equation. Consequently, the approximation of the periodic solution from Eq. ( 17) is found by using the sequence of vectors {Φ j } when m is large enough. In the next section, these sequences will be used to compute the dynamic responses of the railway sleeper.

Numerical applications

A numerical model of the sleeper has been created with the help of ABAQUS software by considering a uniform section beam with the cross-section 0.1927 m× 0.2841 m. The beam length is 2L = 2.41m. The track gauge is 2a = 1.435m. The mesh of the model is generated with 76 nodes and 75 elements. The element type is 2 node-linear beam element (B21) where each node has 3 DOFs (2 displacements and 1 rotation). This element type corresponds to the Timoshenko beam model. We remark that Tran et al. [START_REF] Tran | A comparison of beam models for the dynamics of railway sleepers[END_REF] demonstrated analytically the small difference in sleeper responses (3%) by using the two beam models for the railway sleeper. The material properties of rail pad, sleeper, foundation and parameters of periodic train loads are given in Tab. 1. These parameters will be used for the numerical examples.

Linear foundation

An analytical model of the railway sleeper posed on a Kelvin-Voigt foundation has been developed by Tran et al. [START_REF] Tran | A comparison of beam models for the dynamics of railway sleepers[END_REF]. In this research, the sleeper is modeled with the help of the Timoshenko beam model and its responses are calculated by using the Green's function. When the beam is modeled by FEM, the solutions are given by Eq. [START_REF] Metrikine | Vibration of a periodically supported beam on an elastic half-space[END_REF].

A comparison of sleeper responses of the 2 models is shown in Fig. 4. In this figure, we present a sleeper displacement and sleeper strain at two positions (rail seat and center) of the sleeper in one period of the moving loads, which corresponds to the time where the train moves a distance of a wagon H w . In this example, the displacement at the rail seat is larger than the one at the [START_REF] Tran | A fast analytic method to calculate the dynamic response of railways sleepers[END_REF] middle of the sleeper as demonstrated in Fig. 4a. Besides, Fig. 4b shows that the sleeper is in compression at the rail seat (strain negative), while it is in traction at the center (strain positive). The reference time t = 0 corresponds to the moment when the first wheel passes over the sleeper. In this figure, we can conclude that the two models give the same results, which confirms the validity of the FEM model. Alternatively, by substituting f N L (u, u) = 0 into Eq. ( 18), the non-linear term is avoided: Q m j = 0. The solutions of the problem can be calculated by iteration procedures as follows:

Φ m+1 j = D -1 F j
If the convergence of {Φ j } is satisfied, the iteration procedure converges to the analytical result. Fig. 5 shows the comparison of the sleeper displacement and sleeper strain when the moving load passes the sleeper at the reference time (t = 0). The calculations have been computed with the number of harmonics n = 25. The numerical results agree well with the analytical solutions (blue In this figure, we show the numerical results calculated after 3 iterations: m = 1, 2, 5. We see that, in the case of a linear foundation, the convergence of solutions is obtained after 2 iterations.

Examples of nonlinear foundations

Cubic-nonlinear foundation

In this example, we consider that the nonlinear part of the foundation obeys a cubic law which can be written as follows:

f N Li (u i , ui ) = ε c k c   0 -(u z (x i )) 3 0  
The values ε c = 1 and k c = 440 × 10 7 MNm -3 . Other track parameters are chosen in Tabs. 1 and 2. The numerical results are calculated with a number of harmonics n = 50 and n 0 = 0. As the sleeper responses are the same when the two bogies pass, the numerical results which are shown in this example correspond to a passage of one bogie. Fig. 6 presents the results for different numbers of iterations m. We note that when m ≥ 4, the responses are almost unchanged and the convergence of the iteration procedures is satisfied. Fig. 7 shows the effect of the nonlinear parameter ε c on the sleeper displacement. We have investigated several different values of ε c . The foundation is linear when ε c = 0. The foundation becomes stiffer when this parameter is bigger, thus the amplitude of displacement decreases. This phenomenon is well demonstrated in this figure.

Bi-linear foundation

Several researches demonstrated that the ballast has different linear behaviors in compression and tension. We consider that the constitutive law of the foundation can be described by a stiffness in compression k + f and in tension k - f . In addition, the damping coefficient of the foundation (ζ f ) is considered as unchanged in the two cases. Due to the separation of the linear and nonlinear parts in the model, the nonlinear term can be rewritten as follows:

f N L (u, u) =      k + f -k f u p if u p < 0 k - f -k f u p if u p ≥ 0
where u p presents a component p of the vector displacement u with (p ∈ [1, 3n n ]). In this example, the sleeper displacements are calculated with the two stiffness of foundation: k + f = 440 MNm -1 and k - f = 352 MNm -1 . Fig. 8 shows the numerical results computed by different number of iterations while the number of harmonics is n = 25. The iteration procedures converges when the number is bigger than 8. Fig. 8b shows that the deformed shape of the sleeper at the reference time is almost unchanged when m ≥ 8. Next, the influence of the number of harmonics on the sleeper responses is shown in Fig. 9. When n ≥ 8, we find that the solution has converged. Consequently, this study demonstrates that the lower harmonics are more important than the higher. Put differently, the higher order harmonics can be neglected in this case. Finally, the effect of the nonlinear parameter r = k - f /k + f on the sleeper displacement is shown in Fig. 10. This ratio is calculated by changing the stiffness k - f while maintaining k + f as constant. For each r value, the track responses have been calculated with the same numbers of iterations, but the convergence of solutions is always satisfied. The deformed shape is obtained at the reference time: t = 0, when the first wheel moves over the sleeper. This parameter indicates a foundation with equal stiffness in compression and tension when r = 1, and a foundation that cannot support any tension when r = 0. We see that the sleeper displacements are greater when the foundation has a smaller tension stiffness. Here, we study the influence of the nonlinear parameter on the sleeper strain. In Fig. 11, we plot the sleeper strain at the rail seat in function of the sleeper at and11b). When the nonlinear term of the foundation obeys a cubic law, the sleeper strain at the center is more important than the ones at the rail seat. In the bi-linear case, the relationship between the sleeper strain at the two positions which are shown in Fig. 11 presents clearly two linear behaviors of the foundations.

Conclusions

In this paper, the dynamics response of a railway track posed on a nonlinear foundation has been studied by a semi-analytical model. The iteration method has been developed which permits to obtain fast results. The convergence of the method is also studied with numerical examples (linear and non-linear). Thus, the dynamic responses of the nonlinear track can be calculated quickly.

In addition, this model can be applicable for different nonlinear foundations.

The numerical examples shows that the ratio between the strain at rail seats and center of the sleeper depends on the foundation behavior. This result is significant and can be used to study the foundation in the future work.

The right term of the last equation depends on the rail displacements. It can be rewritten as a function of the sleeper displacement by substituting Eq. ( 6) into Eq. ( 20) as follows:

1 2π +∞ -∞ kp ŵr(ω)e iωt dω = 1 2π +∞ -∞ k 2 p ûR kz (ω) kp + Ke e iωt dω - 1 2π +∞ -∞ kpQe kp + Ke e iωt dω
By combining the previous result and Eq. ( 7), we obtain the following result:

(21)

Mü + C u + Ku + f N L (u, u) =   1 2π +∞ -∞ k 2 p ûR 1z (ω) kp + Ke e iωt dω   e R 1 -   1 2π +∞ -∞ kpQ 1 (ω) kp + Ke e iωt dω   e R 1 +   1 2π +∞ -∞ k 2 p ûR 2z (ω) kp + Ke e iωt dω   e R 2 -   1 2π +∞ -∞ kpQ 2 (ω) kp + Ke e iωt dω   e R 2
where Q 1 (ω) and Q 2 (ω) are the equivalent train loads at the rail 1 and rail 2 respectively (see Appendix B). ûR 1z and ûR 2z are respectively the two displacement of the sleeper at the crossing points with the two rails in the frequency domain. By developing the first and the third terms on the right side of Eq. ( 8), it can be rewritten as follows:

                   1 2π +∞ -∞ k 2 p ûR 1z (ω) kp + Ke e iωt dω   e R 1 = I R 1   1 2π +∞ -∞ k 2 p û(ω) kp + Ke e iωt dω     1 2π +∞ -∞ k 2 p ûR 2z (ω) kp + Ke e iωt dω   e R 2 = I R 2   1 2π +∞ -∞ k 2 p û(ω) kp + Ke e iωt dω  
Finally, by inserting the last result in Eq. ( 21), we have:

Mü + C u + Ku + f N L (u, u) = I R   1 2π +∞ -∞ k 2 p û(ω) kp + Ke e iωt dω   -   1 2π +∞ -∞ kpQ 1 (ω) kp + Ke e iωt dω   e R 1 -   1 2π +∞ -∞ kpQ 2 (ω) kp + Ke e iωt dω   e R 2 (22)
B Periodically supported beam in steady-state

(k) j Q (k) 1 l Rn v D j
-∞ +∞ Fig. 12: Periodically supported beam model

In the frequency domain, Hoang et al. [START_REF] Hoang | Calculation of force distribution for a periodically supported beam subjected to moving loads[END_REF][START_REF] Hoang | Dynamical response of a timoshenko beams on periodical nonlinear supports subjected to moving forces[END_REF] have demonstrated a relation between the reaction force Rk (ω) and the displacement of the rail ŵr(0, ω) in the frequency domain as follows:

Rk (ω) = Ke ŵr(0, ω) + Q k ( 23
)
where: Ke and Q k are the equivalent stiffness and equivalent loads of the system. The two functions are calculated by the parameters of the rail and the train loads as follows:

               Ke(ω) = 4λ 3 r ErIr sin lλr cos lλr -cos ωl v - sinh lλr cosh lλr -cos ωl v -1 Q k (ω) = Ke(ω) vErIr ω v 4 -λ 4 r K j=1 Q (k) j e -iω D j v (24)
where: λr = 4 ρrSrω 2 ErIr . Er, Ir, ρr and Sr are respectively the Young's modulus of rail, the second moment of area of the rail, the density of rail and cross-sectional area of rail. The expressions show that Eq. ( 23) is applicable for any foundation behavior. In this paper, the rail parameters are given in Tab.

C Calculation of equivalent load

By substituting Eq. ( 9) into Eq. ( 24) and together with the assumption [START_REF] Firat | Dynamic analysis of beams on viscoelastic foundation[END_REF], a new expression of the equivalent charge Q k (ω) can be expressed as follows:

Q k =
The parameters of Eq. ( 25) are explained in B. Moreover, a property of the Dirac comb [START_REF] Bracewell | The Fourier transform and its applications[END_REF] gives the following result:

∞ j=-∞ e -iω Hw v j = 2π v Hw ∞ j=-∞ δ ω + 2πv Hw j (26)
By introducing this property into Eq. ( 25), the equivalent loads can be expressed as follows: F j e iω j t (28

)
where F j is calculated by: The last equation describes a forced oscillation with the exciting force F j e iω j t with frequency f 0 = v/Hw. We remark that ω j = 2πf j = 2πjf 0 . We suppose that the solution of the nonlinear problem can admit the same frequencies as the excitation force. Therefore, with this assumption, there exists a periodic solution of u(t) which can be represented as follows:

F j = Q ErIrHw   1 + e -iω Dw
u(t) = ∞ j=-∞ Φ j e iω j t (30)
Moreover, by performing the Fourier transform, Eq. ( 30) can be explained in the frequency domain as follows:

û(ω) = 2π ∞ j=-∞ Φ j δ(ω -ω j ) (31)
So that, we can deduce the following result:

1 2π +∞ -∞ k 2 p û kp + Ke e iωt dω = ∞ j=-∞ Φ j P j e iω j t (32)
where P j is calculated by:

P j = k 2 p kp + Ke ω=ω j (33)
When ω j = 0, we have particularly the values: P 0 = krp and F (k) 0 = 2Q (k) Lw Hw . Finally, by combining Eqs. ( 8), ( 28) and [START_REF] Castro Jorge | Finite element dynamics analysis of beams on non-uniform nonlinear viscoelastic foundations under moving loads[END_REF], the dynamic equation of the sleeper can be written by:

Fig. 4 :

 4 Fig. 4: Comparison of the analytical (continuous line) and numerical models (circle)

Fig. 5 :

 5 Fig. 5: Comparison of dynamic responses of sleeper obtained by analytical solution and iteration procedures

Fig. 6 :Fig. 7 :

 67 Fig. 6: Numerical results of cubic-nonlinear foundation computed by different number of iterations

Fig. 8 :

 8 Fig. 8: Numerical results of bi-linear foundation computed by different number of iterations

Fig. 9 :

 9 Fig. 9: Numerical results of bi-linear foundation computed by different number of harmonics

Fig. 10 : 4 . 2 . 3

 10423 Fig. 10: Effect of the nonlinear parameter of the bi-linear foundation on the sleeper displacements

Fig. 11 :

 11 Fig. 11: Influence of the nonlinear parameters on the sleeper strain

Fig. 12

 12 Fig.12presents a periodically supported beam model. In this model, each rail k of the track is modeled by an infinite beam posed on periodic supports. The sleeper spacing is l. The train loads are considered by the concentrated loads Q (k) j . Each load is characterized by its distance D j to the first axle and the train speed v.

Q

QKe vErIr 1 +

 1 e -iω Dw v + e -iω Dw +Lw v

v

  + e -iω Dw +Lw v + e -iω 2Dw +Lw v

( 34 )

 34 Mü + C u + Ku + f N L (u, u) = I R

Table 1 :

 1 Materials properties of the rail pad, sleeper and foundation

	Content	Notation	Value	Unit
	Young's modulus of the sleeper	Es	48	GPa
	Shear modulus of the sleeper	Gs	20	GPa
	Shear coefficient of the sleeper	κs	0.845
	Second moment of the sleeper	Is	1.694×10 -4	m 4
	Sleeper density	ρs	2658	kgm -3
	Stiffness of rail pad	krp	192	MNm -1
	Damping coefficient of rail pad	ζrp	1.96	MNsm -1
	Stiffness of foundation	k f	440	MNm -1
	Damping coefficient of foundation	ζ f	58.8	kNsm -1
	Wagon length	Hw	20	m
	Distance of the wagon bogie wheels Dw	1.8	m
	Distance of wagon inner wheels	Lw	8.5	m

  2.

	Content	Notation	Value	Unit
	Young's modulus of the rail	Er	210	GPa
	Second moment of area of the rail	Ir	3×10 -5	m 4
	Rail density	ρr	7850	kgm -3
	Rail cross-sectional area	Sr	7.69 × 10 -3	m 2
	Sleeper spacing	l	0.6	m
	Moving force rail 1	Q 1	80	kN
	Moving force rail 2	Q 2	80	kN
	Train speed	v	50	ms -1

Table 2 :

 2 Parameters of the periodically supported beam model[START_REF] Hoang | Calculation of force distribution for a periodically supported beam subjected to moving loads[END_REF]

Acknowledgements Tran Le Hung, ID VNU.2021.TTS14, thanks The Development Foundation of Vietnam National University, Hanoi for sponsoring this research.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Data availability statements

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Conflict of interest

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

A Mathematical transformation

By performing the Fourier transform, and then the inverse Fourier transform of the right term of the Eq. ( 7), we can obtain the following result:

kp ŵr(ω)e iωt dω [START_REF] Tran | A fast analytic method to calculate the dynamic response of railways sleepers[END_REF]


P Aimar
email: aimar@chimie.ups-tlse.fr

P Bacchin
email: bacchin@chimie.ups-tlse.fr

Concentrated phases of colloids or nanoparticles: Solid pressure and dynamics of concentration processes

Introduction

The behaviour of colloids concentrated at an interface is strongly affected by surface interactions occurring on the nanoscale between colloids and between the colloids and the surface. For instance, processes such as ultrafiltration, nanofiltration or reverse osmosis which are classically used to purify, eliminate and concentrate colloids or nanoparticles strongly depend on these interfacial phenomena. The level of fouling, its kinetics or even the way colloids build up (porosity, hydraulic resistance or accumulation reversibility) are driven by colloidal properties. It is therefore necessary to establish experimental and theoretical connections between colloidal properties at the local (micro) scale and the efficiency of the concentration process; this knowledge being compulsory for the control of numerous processes dealing with nanoparticles.

In this chapter, our aim is to show how the concept of solid pressure can be a good vehicle to account for particle-particle interactions in the macroscopic description of separation processes like membrane filtration, sedimentation or drying. It will be shown how solid pressure (or to be more precise its variation with the volume fraction) can be related to the way a suspension resists an increase in concentration. We shall also try to illustrate in which way the change of state of the suspension, can be put in relation to the various forms of accumulation on a surface, both from a theoretical and an experimental point of view. The resistance of the dispersion to over-concentration will be defined by osmotic pressure variation and this resistance will be linked to a Péclet number relative to filtration conditions. The importance of a critical volume fraction inducing a phase transition between fluid (dispersed) and solid (condensed) phases will be underlined in order to describe the transition between accumulation phenomena of polarisation concentration and of dense layer formation. The consequence of this transition on the reversibility of accumulated layers will be discussed. The chapter will then investigate how the variation of osmotic pressure with volume fraction can explain the variety of accumulation phenomena occurring at an interface. 1 The concept of solid pressure to depict concentrated colloidal or nanoparticle dispersions

Behaviour of concentrated colloidal dispersion

Dispersions of colloidal or nanoparticles exhibit specific behaviour when concentrated because of surface interactions which become predominant phenomena for submicron particles. These surface interactions are multiple in nature (electrostatic repulsion, Van der Waals attraction or hydrophobicphilic interaction) and have different interaction lengths; thus leading to a complex system when dispersions are concentrated. When investigating properties of concentrated colloids it can be useful to consider the phase diagram [START_REF] Anderson | Insights into phase transition kinetics from colloid science[END_REF] as shown in Fig. 1. Five different phases appear relative to the concentration (x axis) and the degree of destabilisation (y axis):  "gas" phase: diluted dispersion of stable particles having a free and random motion  "liquid" phase: network of stable colloids interacting by repulsion and then moving from and towards equilibrium position (if particles are monodisperse, the network is ordered: colloidal liquid crystal)  aggregate phase: dilute suspension of aggregate  gel phase: network of colloids interacting by attraction characterised by an elastic behaviour  solid phase: solid structure where colloids are in contact (with different compacities according to the fractal dimension of the structure) The degree of destabilisation is dependent on the surface properties (e.g; charge, hydrophobic-philic character) and on the composition of the interstitial liquid (e.g. ionic strength, ionic valency)

Figure 1: Schematised phase diagram of a colloidal dispersion

The boundaries between these areas correspond to phase transitions. The transition from the gas phase to the aggregate phase corresponds to an aggregation that occurs when the dispersion is destabilised (by adding a salt or a coagulant). When concentrating these phases, the creation of the liquid or the gel phase is due to the percolation of surface interactions leading to a network of interacting colloids (by repulsion for the liquid phase and by attraction for the gel). If concentrated again, these phases undergo spinodal decomposition [START_REF] Jönsson | Ultrafiltration of colloidal dispersions -A theoretical model of the concentration polarization phenomena[END_REF][START_REF] Russel | Colloidal Dispersions[END_REF] i.e. an irreversible transition leading to a solid structure where colloids are in contact and interact with strong Van der Waals attractions or chemical bonds (perfect sink of potential energy). In this section, we will consider the concentration of stable colloids or nanoparticles that will lead to the formation of gas, liquid and solid phases. We will focus on the transition between fluid (gas or 3/25 liquid phase corresponding to a dispersed state) and the solid phase, which has important consequences on concentration processes because of the irreversible character of the transition.

Description of concentrated colloidal dispersions

Surface interactions between two particles (in a dilute dispersion) are depicted by DLVO theory. In simple dispersions, this theory allows the relative importance of repulsive and attractive interaction potential energy to be estimated as a function of the separation distance and then allows the stability of the dispersion to be determined. However, the application of DLVO theory to concentrated dispersions (with multi-body particle interaction) and/or for complex particles is tricky and often leads to calculations far from reality.

It is therefore interesting to consider the potential of the concept of solid pressure [START_REF] Bacchin | Modelling of filtration : from the polarised layer to deposit formation and compaction[END_REF]. This notion combines the concept of osmotic pressure for the dispersed state of particles (section 1.2.1) and the concept of pressure yield stress for the solid state of particles (section 1.2.2). It then gives an equation of state [START_REF] Dubois | Equation of state of colloids coated by polyelectrolyte multilayers[END_REF] which is experimentally accessible and which describes the behaviour of the colloidal dispersion over its whole concentration range and its possible phases.

Osmotic pressure in the fluid phase

The osmotic pressure for colloids has been established via measurements for various dispersions: inorganic nanoparticles [6], phospholipid bilayers [START_REF] Tsao | Osmotic stress measurements of dihexadecyldimethylammonium acetate bilayers as a function of temperature and added salt[END_REF], DNA [START_REF] Podgornik | Molecular Interactions in Lipids, DNA, and DNA-lipid Complexes _Marcel Dekker[END_REF], clay particles such as laponite [START_REF] Michot | diagrams of Wyoming Na-montmorillonite clay. Influence of particle anisotropy[END_REF] or latex particles [START_REF] Bonnet-Gonnet | Osmotic pressure of latex dispersions[END_REF][START_REF] Espinasse | Filtration method characterizing the reversibility of colloidal fouling layers at a membrane surface: Analysis through critical flux and osmotic pressure[END_REF].

A classical theoretical variation in osmotic pressure is presented in Fig. 2 as a function of the volume fraction from computation taking into account the contribution of entropy (Brownian diffusion and excluded volume effects at high volume fraction), repulsive electrostatic repulsion and attractive Van der Waals interactions [START_REF] Jönsson | Ultrafiltration of colloidal dispersions -A theoretical model of the concentration polarization phenomena[END_REF][START_REF] Bacchin | A unifying model for concentration polarization gel-layer formation and particle deposition in cross-flow membrane filtration of colloidal suspensions[END_REF]. The first leg is a linear variation of osmotic pressure with volume fraction (Van't Hoff relationship) and corresponds to the gas phase limited by entropic contribution. When colloids are so concentrated that repulsion between particles occurs, the osmotic pressure increases more severely. This deviation from the linear part is relative to a gas/liquid transition which can, in the case of monodisperse nanoparticles, be linked to an order/disorder transition. When the volume fraction is increased again, the osmotic pressure curve presents an inflexion because of the increasing involvement of Van der Waals attraction which corresponds to a negative contribution in eq 1. A maximum appears when attraction counterbalances repulsion thus leading to instability and to the formation of an aggregated solid phase. This maximum in the osmotic pressure versus volume fraction curve corresponds to spinodal decomposition [START_REF] Russel | Colloidal Dispersions[END_REF] which is linked to an irreversible phase transition between a fluid phase where matter is dispersed to a condensed phase where matter forms a solid phase. Above the maximum, calculation with eq. 1 should give a decrease in osmotic pressure because of Van der Waals interactions. However, this part of the curve (which appears here as a dashed line) can not be experimentally observed : above the spinodal decomposition the particles form a solid state where friction between particles leads to an additional contribution not accounted actually in eq. 1. However, it is no more osmotic pressure (the osmotic pressure is a reversible thermodynamic property whereas friction are dissipative irreversible mechanisms) and then will be introduced in a new section.

Compressive yield pressure of colloids in a solid phase

When spinodal decomposition is reached colloids are in a "solid like" state: particles form a network under attractive interactions. In this state, an important characteristic of the network is the compressive yield stress which represents the resistance of a concentrated suspension to an external compressive field (by analogy with the compressive yield stress which represents the resistance of a concentrated suspension to an external shear field [START_REF] Channell | Shear and compressive rheology of aggregated alumina suspensions[END_REF]. Various experiments have been performed to determine the compressive yield stress [START_REF] Channell | Shear and compressive rheology of aggregated alumina suspensions[END_REF] and this concept has been explicitly incorporated to perform modelling of percolated particle network [START_REF] Buscall | The consolidation of concentrated suspensions. Part 1.-The theory of sedimentation[END_REF][START_REF] Ka Landman | Measurement of the filtration diffusivity D(phi) of a flocculated suspension[END_REF].

For the estimation of compression yield pressure, calculations are classically [START_REF] Landman | Determination of the hindered settling factor for flocculated suspensions[END_REF] made with a power law relationship:

m cp cp crit crit 1 1 1                        (2)
Such an expression describes the solid pressure continuously from the critical osmotic pressure (reached at the critical volume fraction) to an infinite limit for a close packed volume fraction (physically due to Born repulsion). The exponent, m, in this relationship corresponds to the compressibility of the deposit: a small value of m corresponds to a small compressibility: the solid pressure (related to the compression) increasing rapidly with volume fraction. This compressibility can integrate physical properties of particles such as roughness, "deformability", stickiness (degree of VdW attraction, presence of chemical bond).

Solid pressure: definition and example

The solid pressure corresponds to a combination of the osmotic pressure (for the dispersed phase i.e. gas and liquid phase) and the compressive yield stress (for condensed phase i.e. solid phase). It leads to a global equation of state of colloid or nano-particles from the dilute to the concentrated phase which gives a continuous picture of the different phases (gas, liquid and solid) and their associated phase transitions. The solid pressure is dependent on the particle size , the surface physico-chemistry, the composition of the interstitial liquid, the mechanical friction between particles, and so on. It corresponds to the resistance of the dispersion to concentration or compression and is informative on the structuration of particles during a concentration step.

The first and greatest contribution to solid pressure comes from the particle number (Van't Hoff equation) which, for a given volume fraction, is inversely proportional to the particle volume and thus to a 3 . To illustrate the way the solid pressure is correlated to particle size, figure 3 shows two different behaviour patterns for solid pressure as a function of the volume fraction for laponite particles (nanometric clay platelets: diameter of plate 30 nm for 1 nm of thickness at 10-3 M NaCl) [START_REF] Martin | Osmotic Compression and Expansion of Highly Ordered Clay Dispersions[END_REF] and for latex particles (spheres of 120 nm of diameter at 10 -2 M KCl) [START_REF] Espinasse | Filtration method characterizing the reversibility of colloidal fouling layers at a membrane surface: Analysis through critical flux and osmotic pressure[END_REF]. The next most important contribution is from particle-particle interactions and is thus sensitive to parameters that modify the interaction (ionic strength, charge, pH, etc.). A positive contribution is expected for repulsive interactions and a negative one for attractive interactions. Figure 4 illustrates this effect through the measurement of the solid pressure of latex particles at different ionic strengths. An increase in ionic strength leads to a reduction in solid pressure: the decrease in repulsive interaction decreases solid pressure. The solid pressure in a way represents the resistance to the concentration of the dispersion: when repulsive interactions are weak the resistance is low. Solid pressure also depicts phase transition occurring during the concentration of colloidal dispersions. When particles are concentrated, the mean inter-particle distance is reduced. The Van der Waals attractive interaction (with a low distance of interaction) can then become greater than the repulsive interaction (with a large distance of interaction) as depicted through DLVO theory in terms of potential interaction energy to explain dispersion stability. In terms of solid pressure, when attractive interactions become predominant, the osmotic pressure is reduced (the ability of particles to expand in water is reduced). The attractive interactions leads to a change in the curvature of the - function towards low derivative values. Theoretically, an instability is reached when the derivative of the solid pressure with the volume fraction is equal to zero [START_REF] Jönsson | Ultrafiltration of colloidal dispersions -A theoretical model of the concentration polarization phenomena[END_REF] :

ion decomposit spinodal     0  d d (4)
This criterion is related to the definition of the spinodal decomposition [START_REF] Russel | Colloidal Dispersions[END_REF] which, here, relates the irreversible transition of a stable dispersed phase to a condensed phase.

Experimentally, the observation of the transition between stable and unstable dispersion through a maximum in osmotic pressure is difficult to assess. The instability is always preceded by a metastability which leads to a zone of coexistence of fluid and solid phase. However, the sol-gel transition can be assessed when the gel transition occurs over a large range of concentrations (plateau at a volume fraction of around 0.01 for laponite dispersions in fig. 3). Such a plateau is not experimentally accessible for latex dispersions where the transition occurs in a small range of volume fraction. However the transition is distinguishable through experimental analysis of the concentrated phase in terms of the reversibility of packing [START_REF] Espinasse | Filtration method characterizing the reversibility of colloidal fouling layers at a membrane surface: Analysis through critical flux and osmotic pressure[END_REF] or of rheology [START_REF] Bonnet-Gonnet | Osmotic pressure of latex dispersions[END_REF]. Figure 4 shows a fluid-solid transition for latex particles at a volume fraction of around 0.57: the closed symbols correspond to concentrated phase, which cannot be redispersed in solvent after the compression stage (the fluid has undergone an irreversible phase transition). At this point, it is important to note that nanoparticles or colloidal particles present a dual behavior adopting both the behavior of macromolecules (with osmotic pressure) and that of particles (with pressure yield stress). The solid pressure is a key parameter to depict the phase equilibrium and transition occurring in this dispersion in the whole range of concentrations and then to understand the thermodynamics and the interactions in colloidal systems [START_REF] Carrière | Osmotic pressure in colloid science: clay dispersions, catanionics, polyelectrolyte complexes and polyelectrolyte multilayers Colloids and Surfaces A: Physicochem[END_REF]. It represents an equation of state [START_REF] Dubois | Equation of state of colloids coated by polyelectrolyte multilayers[END_REF] containing all the information about colloidal and mechanical interactions.

Figure 5: Schematic representation of solid pressure for different kinds of dispersions and suspensions. Macromolecules exhibit high osmotic pressure. Particles (larger than 10 micrometers) do not exhibit osmotic pressure but resist compression when aggregated (in the solid phase). Nanoparticles and colloids have a dual behaviour and can exhibit a gel behaviour (plateau in osmotic

pressure [START_REF] Mourchid | On viscoelastic, birefringent, and swelling properties of Laponite clay suspensions : revisited phase diagram[END_REF]) at high volume fraction.

These examples illustrate the variety of information that solid pressure can give in a large variety of situations, making the experimental techniques for its determination an indispensable tool for the physical chemist working with colloids or nanoparticles.

Solid pressure and the dynamics of concentration

During a concentration process, the presence of a particles number gradient produces a solid pressure gradient through the system. This solid pressure gradient induces a force on particles which leads to a macroscopic flux of particles from high to low concentration zone. This transport phenomenon can be depicted through a gradient (or collective) diffusion (defined as the collective diffusion of colloids in a concentration gradient in contrast with self-diffusion which acts on a single particle).

When considering a cylinder of section A, a gradient in the solid pressure induces an axial force in the volume, A dx, according to: This force acts on all the particles contained in this volume. The number of molecules is given by:

A dx x x )) ( ) ( (     (5)
p V Adx  ( 6
)
where  is the volume fraction of the dispersion and V p the volume of a particle.

The force acting on one particle is then given by the ratio of relationships ( 5) and ( 6):

dx d V f P     ( 7
)
The dynamics of the concentration process is dependent on solid pressure via this force in a concentration gradient, a velocity or a diffusion coefficient (see box). The dynamics is then consequently a function of: Box 1 : Solid pressure and its consequences in terms of force, velocity and diffusivity This force leads to a particle velocity being:

dx d V m mf u p p      ( 8
)
where m is the mobility of the particle which is dependent on the volume fraction in concentrated dispersions according to:

) ( 0  K m m  (9)
where m 0 is the mobility of particles in dilute dispersion (1/(6a) for a sphere) and K() is the hindered settling coefficient (which corresponds to the ratio of the settling velocity in concentrated to diluted conditions U()/U 0 ). The function K can be usually estimated with the Happel [START_REF] Happel | Viscous flow in multiparticle systems, slow motion of fluids relative to beds of spherical particles[END_REF] or the Batchelor [START_REF] Batchelor | Brownian diffusion of particles with hydrodynamic interactions[END_REF] functions. The velocity of particles induced by the omostic pressure gradient can be considered as diffusive mass transport, N, according to: [START_REF] Bonnet-Gonnet | Osmotic pressure of latex dispersions[END_REF] Comparing this equation to Fick's law,

dx d d d mV u N p p       
    D N
, it is possible to define a diffusion coefficient linked to the variation in osmotic pressure, , with the volume fraction, , by the generalized Stokes Einstein relationship [START_REF] Podgornik | Molecular Interactions in Lipids, DNA, and DNA-lipid Complexes _Marcel Dekker[END_REF]:

       d d K kT V D D p b   (11)
where K() is the hindered settling coefficient, D b is the bulk diffusion coefficient (for infinite dilution) and V p the volume of the colloid. The generalized Stokes Einstein equation accounts for hydrodynamic and thermodynamic interactions through K and  respectively which depend strongly on volume fraction for concentrated dispersions. By using this relationship, the diffusion gradient coefficient corresponding to solid pressure variations can be deduced. The diffusion coefficient is constant for dilute dispersions (gas phase). On increasing the concentration, the diffusion gradient increases rapidly because of the concentration of particles in repulsive interaction (repulsion acts as compressed springs causing extensive collective diffusion of particles). In contrast, for higher volume fractions when attraction becomes important, the diffusion coefficient is reduced. Diffusion is reduced to zero when spinodal decomposition occurs (this nil value is synonymous of irreversibility: the mass cannot return naturally to the bulk by diffusion). In the solid phase, the diffusion coefficient defined by eq. 11 is similar to the concept of filtration diffusivity which plays the role of a diffusion coefficient for determining the time scale of filtration processes in consolidated porous media [START_REF] Ka Landman | Measurement of the filtration diffusivity D(phi) of a flocculated suspension[END_REF]. Such diffusion has been defined from compressional rheology studies according to the definition first given by Landman et al. 1995 [START_REF] Ka Landman | Pressure filtration of flocculated suspensions[END_REF] and developed by Sherwood et al. 1997 [23]:

       d d kr D   1 ( 12
)
where  is the compressive yield stress, k is a convenient pressure scale and r() may be thought of as a hindered settling factor or as a Darcy law constant for flow in porous media (according to the definition given in [START_REF] Ka Landman | Pressure filtration of flocculated suspensions[END_REF]). The generalized Stokes Einstein law (eq. 11) depicts the dynamic transport of colloids from the solid pressure in the solid phase (section 1.2.2.) like the concept of filtration diffusivity (eq. 13).

06/11/2011 10/25

2 Modelling of solid pressure

The osmotic pressure (or solid pressure in the fluid phase) of a dispersion is usually modelled by a power law relationship: [START_REF] Channell | Shear and compressive rheology of aggregated alumina suspensions[END_REF] where the first term proportional to the number of particles (or to the volume fraction) is relative to Van't Hoff relationship and the second term is relative to the increase induced by repulsive interactions between particles (a having positive values for stable repelling particles). For volume fractions over the critical volume fraction, a typical law for compressive yield can be applied according to:

crit if         b p a kT V
crit if                           m cp cp crit crit 1 1 1 (14)
However, relationships 13 and 14 do not depict the fluid-solid transition undergone by the dispersion at the critical volume fraction. This transition is theoretically indicated by a nil derivative of solid pressure with volume fraction (see section 1.2.3.).

To account for the transition, Gaussian fluctuations of the critical solid pressure around the critical point can be considered. The critical solid pressure is distributed around a mean value, crit  and with a standard deviation  with a normal (or Gaussian) distribution characterised by the probability density function -pdf()-, or by the cumulative distribution function -cdf()-. These distribution functions are classically defined by the following equations:

  2 2 2 2 1 ) (    crit crit e pdf crit       (15)         crit crit d pdf cdf ) ( ) ( (16)
The probability density function (pdf) is the density of probability such as the probability that the critical solid pressure will be in the interval [a, b] which is given by the integral of this function between a and b. The pdf function can then give the probability of finding the critical solid pressure between two values. The cumulative distribution function (cdf) is then the probability that the critical solid pressure,  crit , is less than or equal to a given osmotic pressure . The probability of finding the fluid/solid transition occurring is linked to the probability of finding the critical solid pressure between the solid pressure and zero. The probability for this event is:

) 0 ( ) ( ) 0 ( cdf cdf p crit        (17)
This probability is deduced from the solid pressure with a weighable factor leading to a nil derivative of ' at the critical volume fraction (condition for a spinodal decomposition):

) ( ) 0 ( ) ( ' crit pdf cdf cdf        (18)
This negative contribution to the solid pressure is physically linked to the attractive interactions which are the cause for the fluid/solid transition. The parameter  in this model corresponds to the distribution of critical solid pressure around the mean value and can then be related to all the sources of heterogeneity within the colloidal system and to the metastability zone that precedes the spinodal decomposition.

It can be noted, that with such a model, the gradient diffusion coefficient resulting from the derivative of the solid pressure (see box 1) is given by:

06/11/2011 11/25                  ) ( ) ( 1 crit p b pdf pdf d d kT V K D D    (19)
and consequently leads to a nil derivative when = crit .

Such a model for the solid pressure depicts the experimental value of the solid pressure and accounts for the fluid-solid transition. Figure 4 illustrates the model of the experimental solid pressure obtained for latex particles. The model is fitted from a minimisation of the gap between the model and the experimental data by a least squares method 1 using a, b, m and  as variables and with a constraint on the value of the mean critical solid pressure and the critical volume fraction given by experimental determination. The resulting parameters are given in table 1. This model of solid pressure will be used in the next sections to describe the concentration process by filtration (dead end or cross flow), by sedimentation or centrifugation or by drying.

Description of concentration process

It is possible to concentrate a dispersion by: o exerting an external force on particles by sedimentation or centrifugation o extracting the fluid through a filter or a membrane (filtration) or at a gas-liquid interface by evaporation (drying processes) All these processes are controlled by the properties of particles in the dispersion [START_REF] Petsev | Concentrated dispersions of charged colloidal particles: Sedimentation, ultrafiltration and diffusion[END_REF] and in the materials (gel, deposit) created during the concentration process. The solid pressure depicts the behaviour of interacting particles over the whole concentration range (as mentioned in the previous section) and can be used (as presented in next section) to describe the concentration gradient existing during the concentration processes (driven by the filtration rate, gravity or the evaporation rate for filtration, sedimentation and drying respectively).

Filtration

During a filtration, an equilibrium is reached if the permeation drag force, J/m, is balanced by the force for dispersion introduced in eq. 7:

0     m J dx d V P  (20)
This equation can be similarly seen as a flux balance or the application of Darcy law as discussed in box 2. It leads to a relationship where the local differential Péclet number (given for a distance to the membrane) is linked to the solid pressure:

  d K kT V D Jdx p b  ) ( (21) 21
where the first term illustrates the importance of the energy relative to the drag force (induced by the permeate flux J over a distance dx) and the second term concerns the energy relative to the resistance at the over concentration induced by dispersion forces (both terms being normalized by kT). Such an equation can be used to describe the accumulation in the polarised layer near the membrane surface. The relationship to depict the equilibrium during filtration is established with a force balance in eq. 20. However, it has to be noted that the same relationship can be obtained by considering the equilibrium in terms of flux with the classical convective-diffusive flux balance:

0 ) (   dx d D J    ( 22
)
which establishes a balance between a convective flux, J, and a diffusive flux. By writing the diffusion gradient as in box 1 it can easily be demonstrated that eq. 22 can be written as a function of the solid pressure as in eq. 21.

Furthermore, the velocity of the flow through the porous media under a pressure gradient is given by the Darcy law:

dx dp k J    (23)
In a porous medium (here relative to the network of particles accumulated on the membrane), if the flow is sufficiently slow the forces acting on the particles and on the liquid are in equilibrium leading to the sum of the pressure of the liquid and the solid pressure being constant [START_REF] Sherwood | Initial and final stages of compressible filtercake compaction[END_REF]:

        d dp t cons p tan
(24) Such an approach, considering the balance between the permeation drag force and the force due to the pressure loss in the deposit, leads to the following relationship between the deposit permeability k and the hydrodynamic coefficient K () [START_REF] Zick | Stokes flow through periodic arrays of spheres[END_REF]:

      p mV K a k   9 2 2 (25)
According to Eq. 24 and 25, the Darcy law (eq. 23) can then be rewritten as in eq. 20.

At this point, it is interesting to note that very different kinds of balance (drag and dispersion force, convective and diffusive flux, permeation and compression) lead to the same equation. The solid pressure then offers a unique way to depict these equilibria and to account for the properties of the solid (in terms of force of dispersion, diffusivity or compressibility) during concentration processes.

Dead end filtration

During dead end filtration, the accumulated layer is continuously growing. In order to rigorously depict the accumulation, the transient expression for flux balance should be considered:

0       N t  ( 26
)
However, it has been shown [START_REF] Bessiere | Filtration frontale sur membrane: mise en évidence du volume filtre critique pour l'anticipation et le contrôle du colmatage[END_REF] that numerical simulation of the transient equation from eq. 26 leads to the same results in the conventional range of filtration, microfiltration and ultrafiltration as a pseudo-stationary model considering successive equilibrium states over time. In these conditions, the evolution of the accumulated volume of particles (here in m 3 of particles per m 2 of surface area) should be considered with time by assuming that the whole quantity of particles brought by permeation, J, accumulates at the surface: By considering a succession of equilibria with time, eq. 20 can be integrated through this accumulated volume (i.e. from the solid pressure in the bulk,  b , to the solid pressure at the membrane,  m ) according to:

f b t b a V Jdt dx V          0 0 (27
      m b d K kT V D JV Pe p a de ) ( 0  ( 28
)
where the Péclet number here indicates the energy dissipated by friction on all the particles accumulated during the filtration of a volume of liquid equivalent to the particle volume, 6aJV a , over the thermal energy, kT. The numerator of this equation corresponds to the energy relative to the force undergone by the particle near the membrane surface: the particle undergoes the sum of drag forces undergone by all the particles accumulated above the membrane [START_REF] Harmant | Coagulation of colloids in a boundary layer during cross flow filtration[END_REF] considering that the particles are stopped by the membrane and that they transmit their force through repulsive interactions from the top to the bottom layers. The Péclet number increases when filtration conditions are more severe because of an increase in flux, an increase in volume fraction or in filtration time.

Knowing the properties in terms of solid pressure of the colloidal or nanoparticle dispersion, it is then possible to estimate the consequences of the Péclet number (and thus the filtration conditions) on the solid pressure of the materials at the membrane surface. Figure 8 represents the increase in volume fraction at the membrane resulting from an increase of the Péclet number. A more rapid increase in volume fraction at the membrane is observed for the latex dispersion with the higher ionic strength (lower stability). This result is linked to the fact that the increase in ionic strength reduces the repulsive interaction between particles. It leads to a network of particles at the membrane surface presenting less resistance to concentration. Similarly, it should be noted that the Pe number for which a solid phase is formed at the membrane surface is lower for less stable dispersions: these critical Pe number are 32.9 and 13.5 for the latex dispersion at 10 -3 M and 10 -2 M respectively in KCl.

Cross-flow filtration

In cross-flow filtration, when the stationary state is reached, accumulation takes place in a mass boundary layer with a fixed thickness with time at a given position along a membrane. In more complex situations, the thickness can change along the length of the membrane. The accumulation can then be depicted through the equilibrium introduced in eq. 21 integrated through the boundary layer thickness (according to the boundary conditions introduced in figure 7):

      m b d K kT V D J Pe p cf    ) ( 0 ( 29
)
where the first term of this equation is relative to a Peclet number which represents the energy dissipated by friction on a particle during the filtration of a volume of liquid, covering length  over the thermal energy. It has to be noted that the ratio Pe de /Pe cf (given by eq. 28 and 29) is representing the ratio of V a / which is the mean volume fraction in the polarised layer. This Péclet number is linked to the resistance to concentration induced by dispersion forces taking into account the solid pressure (last term of eq. 29). A figure similar to that obtained for Pe de (figure 8) can be plotted for Pe cf . It the same way, the existence can then be noted of a critical Péclet number that will lead to the formation of a solid at the membrane surface. Calculations give Pe cf at 87.8 and 36.1 for the latex dispersion at 10 -3 and 10 -2 M respectively in KCl.

Concentration by sedimentation and centrifugation

As already discussed [START_REF] Davis | An asymptotic description of transient settling and ultrafiltration of colloidal dispersions[END_REF], during sedimentation or centrifugation, the balance between gravity (or centrifugal force) and diffusion lead to an equilibrium which can be described by the following relationship:

0 ) (   dx d D u sed    ( 30
)
where the mass flux induced by the centrifugation velocity, u sed , is counterbalanced by diffusion. As for filtration, this equilibrium can also be written by taking into account the solid pressure through eq. 11 leading to:

0 ) (     dr d V g V p w p p    ( 31
)
where the first two terms represent the net force acting on the particle (the difference between the gravitational force and the buoyancy force:  p and  w being the density of the particle and the water respectively) and the second term is the force resulting from the solid pressure gradient (thermodynamic force). This equilibrium allows the classical Pe number to be introduced for sedimentation according to:

          bot d kT V kT r r g V Pe p bot w p p sed int int ) (    (32)
Likewise, the equilibrium in a centrifugal force field can be depicted:

          bot d kT V kT r r w V Pe p bot w p p centrif int 2 int 2 2 ) (    ( 33
)
where w is the angular rotation speed. r r int r bot

Concentration during the drying process

During a drying process, water is extracted from the dispersion at the liquid/gas interface, selectively concentrating the particles. A gradient of concentration through the thickness of the drying film then appears. The over-concentration at the film-air interface can induce skin formation. Comprehensive modelling [START_REF] Routh | Horizontal drying fronts during solvent evaporation from latex films[END_REF] has shown the importance of the evaporation rate on the processes of film formation.

During film formation, if the particle velocity induced by natural convection is negligible, the particles in the film are only submitted to the force induced by the solid pressure gradient (diffusion transport). The flux of particles in the system can then be written as:

dx d mV N p    (34)
The local volume fraction of solid, , satisfies the mass conservation equation:

0       N t  ( 35
)
By assuming a pseudo-steady state regime, the mass conservation equation in the steady state can be considered to lead to a constant flux of particles:

    dx d mV N p constant ( 36
)
At the free surface, the evaporation of water (where the evaporation rate h is expressed in m 3 /(m 2 .s)) leads to the removal of water while the particles stay in the film. Per volume of liquid extracted, the volume of particles remaining in the film is ) 1 /(    which represents the solid ratio (ratio of solid volume over liquid volume). Assuming a constant evaporation rate, evaporation contributes to an influx of particles at the free surface:



(39) where  0 is the initial volume fraction of particles in the film. Assuming that for the initial stage of drying the volume fraction near the substrate is always close to the initial volume fraction,

b x    0
, integrating eq 38 throughout the volume of accumulated particles in the film leads to:

            i b d K T k V D e e h Pe B p e x e x dry ) ( 1 0 0 0      (40)
where the Péclet number (similar to the one defined in [START_REF] Routh | A process model for latex film formation: limiting regimes individual driving forces[END_REF]) links the drying operating conditions to the dispersion properties. This simplified approach is valid only when the zone where the polarization occurs is thinner than the total film thickness, e, i.e. for the first step of drying (as confirmed by the calculation in section 5.3). When the polarisation layer and the film thickness are of the same magnitude eq. 40 does not apply to the whole thickness: in this particular case the mass flux is zero against the substrate and the volume fraction there increases. ( x=0 >  b ).

Pe number and concentration processes

As demonstrated in the previous section, the concentration processes can be depicted by means of the Péclet number which links the operating conditions of the concentration processes to the properties of the material being concentrated (table 2). By definition, the Péclet number is the ratio between the acting force leading to the concentration processes (J the filtrate flux for filtration, g gravity for sedimentation, w rotation during centrifugation, h the evaporation rate during drying) and the dispersing force (Brownian diffusion). These process operating conditions are linked to the properties of the dispersion which are the solid pressure, , and the sedimentation coefficient, K, integrated from the bulk volume fraction to the interfacial volume fraction (at the membrane surface for filtration, at the bottom of the tube for centrifugation and sedimentation or the gas/liquid interface for drying). It can be noted that for sedimentation and centrifugation equilibrium the sedimentation coefficient, K, does not appear in the Peclet number: during these operations, when equilibrium is reached there is no relative displacement between particles and the liquid phase.

For an ideal solution (i.e. at low concentration when multi-body interactions are negligible), the sedimentation coefficient is around unity (K->1) and the solid pressure, , follows the Van't Hoff law. The integrals can then be simplified to a form involving the constant diffusion coefficient. For instance, in cross-flow filtration this simplification leads to the classical film model where the permeate flux varies with the logarithm of the ratio of the concentration at the membrane to the concentration in the bulk.

Table 2: Péclet numbers describing the mass transfer equilibrium under different concentration processes

Concentration process

Péclet number

Filtration

Dead-end The variation of the Pe number with the solid pressure and the volume fraction is illustrated in figure 11 plotting, in a double graph, the variation of solid pressure with the volume fraction (for the dispersion of latex particles in 10 -3 M KCl in the right plot) and the solid pressure versus Pe number (in the left plot). The left graph can be considered as representing operative functions (which are a general description of mass balance in each processes) and the right graph the equation of state of the dispersed matter (an equilibrium property of the dispersion of colloids or nanoparticles). Such a representation can be used to determine the Péclet number required for compression from an initial concentration to the interfacial concentration. For example, the Pe number needed to increase the concentration from an initial value of 0.001 to an interfacial value of 0.5 can be given by the difference of Pe number for each of these concentrations as represented by the dashed lines in Figure 11. Compression of this type occurs for Péclet numbers of Pe dry = 150, Pe de = 367, Pe cf = 981 and Pe sed or centrif = 19 400. These values represent the energy (in terms of kT) of the acting force required to compress the dispersion to a volume fraction of 0.5 at a given point in the process. It can be seen that for a low volume fraction (ideal dispersion without interaction), the energy to compress dispersion layers becomes greater with sedimentation, centrifugation and cross-flow whereas for a large volume fraction (highly interactive particles) cross flow becomes a more "economic" way to compress layers than dead-end filtration. These values can then help to determine the operation conditions (for example in terms of acting force) relative to the concentration process but also to determine the conditions for which the critical volume fraction is reached locally in the process.

      m b d K kT V D JV Pe p a de ) ( 0  Cross-flow       m b d K kT V D J Pe p cf    ) ( 0 Sedimentation           bot d kT V kT r r g V Pe p bot w p p sed int int ) (    Centrifugation           bot d kT V kT r r w V Pe p bot w p p centrif int 2 int 2 2 ) (    Drying             i b d K T k V D e e h

Critical conditions of concentration processes

The Péclet numbers introduced in the previous section allow us to define critical Péclet numbers that lead to the appearance of a solid at the surface (i.e. when the solid pressure at the membrane reaches the critical solid pressure). These critical Péclet numbers represent the critical friction energy that has to be exerted on particles through the concentration process to form a solid phase (fluid/solid transition in a phase diagram). Before this transition, the energy transferred to the particles by friction is reversibly stored by dispersion (this energy can help the dispersion of accumulated layers if filtration is stopped) whereas, above the transition, the proportion of the energy exceeding the critical Pe number (i.e. Pe-Pe crit ) is dissipated through friction in the deposit. This behavior is the consequence of the irreversibility of the liquid-to-solid transition (spinodal decomposition) which can have various consequences on the way concentration processes work as illustrated in figure 12 and analysed in the next section.

Deposit during filtration

Sediment during centrifugation

Skin formation during drying xg

Critical conditions of filtration

In dead-end filtration, the critical Pe number results from a combination between the permeate flux, J and the accumulated volume of particles, V a (table 2). V a is dependent on the filtered volume, V f (and thus on filtration time), according to equation 27. The critical conditions are linked to the accumulated volume because a particle at the membrane surface supports all the drag that the permeate flux exerts on particles in the polarized layer. Accumulation then leads to the aggregation of the particles near the surface. Figure 13a indicates the values of the pair "permeate flux/filtered volume" for which the solid phase is formed (or not formed) at the membrane surface. No solid phase is formed for conditions of low permeate flux and low filtered volume. These conditions are more restrictive for less stable dispersions (i.e. for 10 -2 M in KCl). In cross-flow filtration, the formation of a solid phase at the membrane surface is linked to the production of the permeate flux and the boundary layer thickness according the critical Pe number, Pe cf (eq. 29). The absence of solid phase is favoured for conditions of low permeate flux and for low boundary layer thickness. Low boundary layer thickness corresponds to conditions of high cross-flow velocity. The classical Leveque relationship has been used to link the boundary layer thickness to the cross flow velocity (by considering a tubular membrane device with a diameter of 6 mm and a length of 1.2 m using the dispersion data given in previous sections). Figure 13b reports the conditions in permeate flux, J and cross-flow velocity for which a solid phase exists on the membrane (conditions of high flux and low cross-flow velocity) or in contrast the conditions for which no solid phase forms.

sub-critical filtration over-critical filtration sub-critical filtration over-critical filtration

dead-end filtration cross-flow filtration These critical conditions allow filtration with no solid phase formation at the membrane (which are classically named sub-critical operating conditions [START_REF] Bacchin | Critical and substainable fluxes : theory, experiments and applications[END_REF]). The absence of solid phase (by definition an irreversible phase) avoids subsequent fouling. In cross-flow filtration, it corresponds to low permeate flux and high cross-flow velocity enabling work in the stationary mode without formation of solid at the membrane. In dead-end filtration, accumulation occurs progressively with filtration time. To avoid solid phase formation at the membrane, the membrane must be periodically rinsed. Studies have shown [START_REF] Bessiere | Low fouling conditions in dead end filtration : evidence for a critical filtered volume and interpretation using critical osmotic pressure[END_REF] that rinsing before reaching the critical filtered volume (i.e. before reaching a critical In actual process applications, however there can be uneven distribution of permeate flux (locally on the membrane) and also, under cross-flow conditions, inhomogeneities in boundary layer thickness because of the progressive build up of material along the membrane surface [START_REF] Bacchin | Distributions of critical fouling conditions: modelling, experiments analysis and consequences for cross flow membrane filtration[END_REF]. Because of these heterogeneities, solid phase can form locally even if the mean conditions are sub-critical.

Critical conditions in sedimentation or centrifugation

For batch sedimentation or centrifugation, eq. 33 has to be solved with global mass balance over the sedimentation or centrifuge cell:

     bot r r bot init r r dr int int   (41)
Eq. 32 and 41 allow the concentration profile to be determined along the cell. Figure 15 represents the evolution of the volume fraction along the centrifuge cell (r int =0.1 m and r bot =0.15 m) for different rotation speeds in rpm. The different rotation speeds: 100, 500 and 1000 rpm correspond to 16 g, 410 g and 1600 g respectively. The volumic mass is taken to be 1300 kg/m 3 for the particles and 1000 kg/m 3 for water. The initial volume fraction of the dispersion is 0.01. It can be seen that the more stable latex particles, dispersed in 10 -3 M KCl, (round symbols) are less compressed than the those dispersed in 10 -2 M in KCl (triangles). The repulsive electrostatic interaction, more pronounced in stable dispersions, leads to greater resistance to compression induced by the external force field. The centrifugal force can overcome the repulsion between particles (which is responsible for the dispersion) when a critical value for the rotation speed is exceeded. This happens when the critical value of the Pe number is reached:

     crit crit d kT V Pe p centrif int  (42)
These critical values are 1.2 10 +5 and 9.6 10 +4 for the latex dispersion at 10 -3 and 10 -2 M respectively in KCl. Figure 16 presents the critical centrifugal rotation speed in rpm relative to the latex particles. The more stable the dispersion, the higher the centrifugal rotation speed required to collapse the particles at the bottom of the cell and to form a sediment. If this speed is not reached, the diffusion induced by repulsive interaction between particles will lead to a redispersion of the particles on withdrawal of the centrifugal force .

Critical conditions in drying

Solving equation 40 gives the vertical concentration gradient in the film during drying. Figures 17 present the evolution of the volume fraction profile with time during drying of a film (initial volume fraction of particles 0.05 and initial thickness 2 mm) of a dispersion of latex in 10 -3 M KCl (Figure 17a) and in 10 -2 M KCl (Figure 17b). These calculations were performed for the first step of the drying process in order to satisfy the assumption made in the simplified model presented in section 3.3 and to avoid taking into account the increase in ionic strength induced by the concentration of the film. The calculations were stopped when the critical volume fraction was reached: the model is based on a constant evaporation rate and does not take account of the reduction in evaporation rate induced by skin formation at the interface. For dispersion of particles exhibiting high repulsive interactions in 10 -3 M KCl, the volume fraction at the interface increases more slowly than for particles with less repulsive interactions (dispersed in 10 -2 M KCl). These results clearly show that the dynamics of concentration by drying is controlled by the solid pressure of the dispersion: stabilized particles resist the concentration process as already mentioned for filtration or centrifugation.

The main consequence on the process is that the formation of a skin at the interface will occur later for stable dispersions (after 1200 s for dispersions in 10 -3 M KCl and after 550 s for dispersions in 10 -2 M KCl). Recent experiments [START_REF] Erkselius | Drying rate variations of latex dispersions due to salt induced skin formation[END_REF] clearly show that operating with destabilized dispersions (under the critical coagulation concentration -ccc -i.e. without aggregation) led to the formation of skin (that dramatically reduced the evaporation rate) more quickly during the drying (salt-induced skin formation). In contrast, if the addition of salt leads to particle aggregation (above the ccc), the consequences on skin formation are the opposite: the aggregation leads to settling of the particles toward the substrate thus limiting skin formation at the gas/liquid interface.

Skin formation occurs when the volume fraction at the interface reaches its critical value: The critical values for the latex dispersions are 14.6 at 10 -3 M KCl and 7.0 at 10 -2 M KCl. These critical conditions result from a critical pair of operating conditions given by the evaporation rate, h, and  0 (e 0 -e) which is the product of the initial concentration with the film thickness reduction (or the volume evaporated per interface area).

Conclusion

The description of the concentration processes in colloidal or nanoparticle dispersions (filtration, centrifugation, drying) has to take into account the presence of surface interactions between entities (electrostatic repulsion, Van der Waals attraction, etc.) as they play a key role among transport phenomena. This source of complexity can be depicted through the concept of solid pressure. Solid pressure can be considered as the equation of state for the solid matter (particles or macromolecules) in the liquid. The solid pressure is considered as a combination of the osmotic pressure (for the dispersed state of matter) and the compressive yield stress (for the condensed state of matter at higher concentrations). This equation of state is at the base of the description of phase transitions in colloidal or nanoparticle dispersions, induced by surface interactions. This chapter shows how the solid pressure can be used to depict the dynamics of the concentration processes. The model demonstrates the existence of critical operating conditions during the concentration process which lead to the formation of a condensed phase on a membrane (fouling during filtration), at the bottom of a centrifugal cell (sediment formation in centrifugation) or at the gas/liquid interface during drying (skin formation). This general approach allowed us to develop a model of the complex interplay between dispersion stability and the dynamic of concentration processes.

Figure 2 :

 2 Figure 2: Classical variation of osmotic pressure with volume fraction for stable colloids

Figure 3 :

 3 Figure 3: Solid pressure versus volume fraction for two stable dispersions with different particle sizes: laponite particles (nanometric clay platelets: diameter of plate 30 nm, thickness 1 nm ) and latex particles (spheres 120 nm in diameter)

Figure 4 :

 4 Figure 4: Solid pressure of latex particles versus volume fraction for different ionic strengths (diamonds and circles at 10 -3 M and 10 -2 M in KCl respectively). Open symbols correspond to fluid (dispersed) phase and closed symbols to solid (condensed) phase. Lines correspond to the fitting of the model described in section 2.

Figure 6 :

 6 Figure 6: Force, f, acting on one particle in a gradient of solid pressure 

  are dependent on: o surface properties (zeta potential), o dispersal properties (ionic strength, ion specificity).

1 aFigure 7 :

 17 Figure 7 : Schematic representation of the filtration under a filtrate flow, J, with the accumulation of matter in a polarized layer of thickness .

Box 2 :

 2 A single law to depict a force balance, a flux balance and the Darcy law

  in KCl (fluid phase) 10-2 M in KCl (fluid phase) 10-3 M in KCl (solid phase) 10-2 M in KCl (solid phase)

Figure 8 :

 8 Figure 8: Volume fraction at the membrane versus the Pe number, Pe de . Calculations were done with equation 28 and for the latex dispersion presented in section 2. The bulk volume fraction is taken at 0.00001.de

Figure 9

 9 Figure 9 : Schematic representation for the sedimentation or centrifugation process.

  During drying, the volume of accumulated particles, V a , (here in m 3 of particles per m 2 of surface area) represents the excess of particles induced by the reduction film thickness from e 0 to e:

  Figure 10 : Schematic representation of the drying processus.

Figure 11 :

 11 Figure 11: Evolutions (in the right sector) of the solid pressure, , with the volume fraction, , for the dispersion of latex particle in KCl 10 -3 M and (in the left sector) of the solid pressure versus the Pe number for different concentration processes ( b for this calculation is taken at 10 -4 )

Figure 12 :

 12 Figure 12: Critical conditions of concentration a) during filtration with formation of a deposit b) during centrifugation with formation of sediment c) during drying with formation of a skin

Figures 13 a

 13 Figures 13 a and b: Representation of operation conditions leading to the presence or to the absence of solid phase at the membrane for dead-end a) and cross-flow b) filtration. Conditions are dependent on the permeate flux, J and the filtered volume V f in dead-end filtration and on the permeate flux and the cross-flow velocity in cross-flow filtration.

Figure 14 :

 14 Figure 14 : Application of critical filtration conditions to the determination of sub-critical filtration in dead-end and cross-flow filtration modes avoiding the formation of solid phase at the membrane.

  filtration time) allows successive filtration and rinse steps without solid phase formation at the membrane.

Figure 14 :

 14 Figure 14: Particle concentration profiles during centrifugation of latex particles.

Figure 16 :

 16 Figure 16: Critical value for the rotation speed during centrifugation leading to the formation of a sediment.

Figures 17 a

 17 Figures 17 a and b: Particle volume fraction versus the film thickness (symbols) as a function of time fordispersions of latex particles . We also report the increase of the volume fraction (line) at the gas/liquid interface as the film thickness decreases during drying.

Table 1 :

 1 Parameters of the model (eq. in section 2) for solid pressure plotted in figure4

		latex 10 -3	M	latex 10 -2	M
	R	4.00E-08	4.00E-08
	Vp	2.6808E-22	2.6808E-22
	a	8.04E+05	1.39E+07
	b	6.67		13.1
	m	0.2738564	0.25317859
		5000		5000

24/25

7 References

Skin formation

No skin formation

Figure 18: Representation of the critical pair of drying operating condition leading to skin formation. Skin formation is favoured for high evaporation rates, h, and high values of the term  0 (e 0 -e) (i.e. for example for high initial concentrations in the film). Destabilisation of the dispersion favours the operating conditions where skin formation occurs. Figure 18 presents the range of operating conditions where skin formation will occur at the interface. It can be seen that the conditions where skin formation will not occur is favoured:  when drying more stable dispersions (particles in 10 -3 M KCl);  when drying a dilute dispersion (low initial volume fraction,  0 );  when operating at low evaporation rates;  when the volume of water evaporated per interface area (e 0 -e) does not exceed a critical value. The approach proposed here allows these drying regimes to be quantified from the properties of the dispersion in terms of solid pressure.