`Discrete Dynamics in Nature and SocietyVolume 2011, Article ID 415921, 24 pageshttp://dx.doi.org/10.1155/2011/415921`
Research Article

## Stability and Numerical Analysis of the Hébraud-Lequeux Model for Suspensions

1Centro de Investigación Operativa, Universidad Miguel Hernández de Elche, Avenida. Universidad s/n, 03202 Elche (Alicante), Spain
2Departament d'Economia Aplicada, Facultat d'Economia, Universitat de València, Campus dels Tarongers s/n, 46022 València, Spain

Received 27 October 2010; Accepted 16 March 2011

Copyright © 2011 Ángel Giménez et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We study both analytically and numerically the stability of the solutions of the Hébraud-Lequeux equation. This parabolic equation models the evolution for the probability of finding a stress in a mesoscopic block of a concentrated suspension, a non-Newtonian fluid. We prove a new result concerning the stability of the fixed points of the equation, and pose some conjectures about stability, based on numerical evidence.

#### 1. Introduction

Non-Newtonian (or complex) fluids are very common in nature and industry, appearing for instance, in foods, biofluids, personal care products, pharmacology and bioengineering, electronics and optical materials, and energy and plastic production. In fact, one could say that Newtonian (or simple) fluids are rather an exception (if not an idealization), even though they include such a prominent member as water. Attending to their rheologic properties, complex fluids are classified in different categories, including suspensions, colloids, polymer melts, liquid crystals, gels and foams, among others.

In this paper, we will consider numerical approximations for the Hébraud-Lequeux equation for suspensions [1] and, more specifically, for the following model (proposed in [2]): with , and . Furthermore, denotes the characteristic function of the interval and is the Dirac delta function supported at the origin.

We point out the important fact that the function is a probability density function for every fixed , so it must satisfy that and , for any .

Existence of solutions of (1.1) is given in [3]. In [4] a simplified version of this model is studied, in which the term is neglected, proving that in such case the weak solutions of the equation generate a multivalued dynamical system possessing a global attractor with respect to a suitable phase space.

In [5] we considered a lattice dynamical system of the simplified model, obtained by discretizating the derivative and using a Simpson's quadrature for the integral . Also, the Dirac delta function is approximated by a parabolic-type function. Observe that the variable belongs to the whole real line, so the lattice dynamical system consists of an infinite number of ordinary differential equations. We proved that this system defines a continuous semigroup of operators in an appropriate space and, as well, the existence of a global compact connected attractor.

Moreover, by a suitable truncation, finite-dimensional approximations of the lattice system are given, proving that each of the approximating systems possesses a global compact connected attractor and also the upper semicontinuous convergence of these attractors to the one generated by the lattice system.

Finally, an implicit Euler scheme is implemented in the finite-dimensional approximations. Thus, a numerical algorithm approximating the solutions of the original equation is obtained. In [5] the resulting discrete dynamical system is studied, proving the existence of a global compact connected attractor and the upper semicontinuity of the sequence of attractors with respect to the previous finite-dimensional approximations.

In this paper, using this algorithm, we present several numerical approximations of the solutions of the equation.

At the same time, we prove a new result concerning the stability of the fixed points of the equation and give also some conjectures which are supported by the numerical simulations.

More precisely, let us consider the simplified equation (i.e., ). Also, for simplicity we take and . In such a case it is well known [3] that all the probability densities with support in the interval are fixed points. Also, if , then there no more fixed points, and if , then there exists a unique fixed point not contained inside the interval . It is known that this fixed point is asymptotically stable [6].

In this paper, we prove that if , a subset of probability densities with support in the interval is stable.

Furthermore, based on our numerical simulations we make some more conjectures about the stability of fixed points and also about the asymptotic behavior of solutions as time goes to .

First, it is natural to think that when the unique fixed point with support not contained inside the interval if globally asymptotically stable, that is, if the initial data is not a fixed point, then the corresponding solution will converge to this point as time goes to (which would imply that all the other fixed points are unstable). However, it seems that this is not true, as we found numerically a solution converging to a fixed point with support in the interval . Hence, several stable fixed points seem to coexist, although this fact has not been proved yet.

Secondly, the numerical simulations support the conclusion that every solution of the equation converges to a fixed point as time goes to , in other words, that the omega-limit set of every trajectory is a fixed point. Again, there is still no proof of this fact.

Finally, it is natural to think that most fixed points with support in the interval are unstable, since if we take an arbitrary fixed point and a close initial data, usually the numerical solution converges to distant fixed point.

#### 2. Implicit Finite Difference Scheme

We shall consider the simplified model that results from (1.1) after taking . Also, we set for simplicity, although our results can be obtained for an arbitrary . Hence, we consider the equation with . Here, is the characteristic function on the interval .

Consider an initial condition satisfying Then it is proved in [3, Theorem 1.1] that problem (2.1) has a unique solution satisfying the following properties for all : where exists for any . The last property implies that the average stress belongs to . From and the above properties, it follows that .

Our aim is to approximate the solution to the nonlinear partial differential equation given in (2.1) by solving an implicit discrete problem based on finite differences. We begin by discretizing the unbounded spacial domain via a uniform grid with space step . Let be such that and . We fix so as and . In our numerical simulations, we will always take in such a way that its support is contained in with sufficiently large, and we truncate this unbounded grid for . If has an unbounded support, we can take large enough in order to achieve . Similarly, we discretize the time domain with a time step . The resulting numerical grid is illustrated in Figure 1.

Figure 1: Numerical grid.

Denote by and the usual norm and scalar product in . Let us define the following norms in : for . Obviously, all these norms are equivalent in .

We shall implement an implicit scheme in finite differences by using the following approximations. Let denote the vector grid satisfying . Since the partial differential equation (2.1) is of evolution type, we approximate at the point by the forward difference

In order to approximate at the point , we define the operator by , where and we take . This choice is necessary for the matrix to satisfy the condition where being , if is even, and being , if is odd. This follows from , for any . Observe also that , which together with equality (2.7) is a key fact to prove Theorem 2.1 (see below).

We define also the matrix , where which satisfies , where the superscript “” stands for “transpose”.

The integral is approximated by which is a slight variation of the Simpson formula. The only difference is that the weight for the terms and is equal to zero instead of 1. Observe that we can also set , where

Finally, the -function will be approximated by the following piecewise parabolic function:

This function is continuous, its support is ], and it attains the maximum value at , which is equal to . Also, the integral equals 1 (see Figure 2). Let us take such that is even and let so as , that is, the support of is strictly included in the interval . Then at this function looks as or Denote .

Figure 2: Approximation of Dirac delta function.

Substituting all the above discretization formulas into the partial differential equation (2.1), we get the following algebraic system for each integer :

Bearing in mind that and denoting by , the difference scheme system (2.15) can then be equivalently written in the form which can be expressed in matrix form as where Observe that the matrix depends on the time step , while depends on the steps and . It is easy to check that is strictly diagonally dominant for all . In particular, this fact guarantees that is nonsingular and the system (2.17) has a unique solution for every . Recall also the following result.

Theorem 2.1 (see [5]). Assume that . If and , then and , for all .

#### 3. Numerical Simulation

We present in this section some numerical solutions of the simplified model (2.1) made with the finite difference scheme (2.15). For simplicity, let us consider the case where . We shall consider two cases depending whether or with different initial conditions , which may or not have support intersecting the interval . We shall not plot cases when the support of is totally contained in , since in such a case would be a fixed point of the discrete system (2.15). Therefore, we will focus on cases such that .

Case 1 (). In this case all the probability densities whose support is totally contained in are fixed points of the system. It is known [3] that for the continuous system, the unique fixed point whose support is not totally contained in is given by where satisfies the equation in order for the normalization constraint to hold true. Furthermore, this fixed point is stable (see [6]).
In the next experiments, we shall take the constants and and we define the spacial and time mesh size and , respectively. Other choices of lead to similar results. The numerical simulations will be plotted as a sequence of discrete functions converging to a fixed point. The initial condition and the fixed point will be drawn in black, and between both there are a sequence of plots converging to the fixed point which are drawn in a thinner black. We shall consider different cases depending on the initial conditions, and we shall discuss the results.
Subcase 1. Consider the initial condition with support out of given by The numerical results (see Table 1 and Figure 3) show how the sequence converges to the fixed point (3.1).Subcase 2. Consider the initial condition with support in and out of given by We observe (see Table 2 and Figure 4) the same behavior as in Subcase 1. Subcase 3. Consider a parabolic-type initial condition with support in and out of given by We obtain a similar behavior as in the above cases (see Table 3 and Figure 5).Subcase 4. We start from an initial condition verifying given by The sequence approaches to a fixed point whose support is contained in (see Table 4 and Figure 6).
This fact leads to the hypothesis that besides the fixed point (3.1) there are stable fixed points for (2.1) with support contained in .
Subcase 5. As in Subcase 4, consider an initial condition verifying given by On the contrary, this sequence approaches to the fixed point (3.1) (see Table 5 and Figure 7). In fact, this is the usual behavior of the system in the majority of cases.

Table 1: Numerical results: Subcase 1
Table 2: Numerical results: Subcase 2.
Table 3: Numerical results: Subcase 3.
Table 4: Numerical results: Subcase 4.
Table 5: Numerical results: Subcase 5.
Figure 3: Graphics sequence: Subcase 1.
Figure 4: Graphics sequence: Subcase 2.
Figure 5: Graphics sequence: Subcase 3.
Figure 6: Graphics sequence: Subcase 4.
Figure 7: Graphics sequence: Subcase 5.

We note that, usually, when we take an arbitrary fixed point with support in the interval and an initial data close to it, then the numerical solution converges as time goes to a fixed point situated far away. Hence, we think that most of the fixed points of such type are unstable. However, we prove in Section 2 that some of them are stable if .

Case 2 (). In this case the probability densities with support totally included in are the unique fixed points. Let us take the constants and . We define the spacial and time mesh size and , respectively. As in Case 1, the numerical simulations will be plotted by a sequence of plots converging to a fixed point. The initial condition and the fix point are both drawn in a thicker black. Let us take the same first three initial conditions as in Case 1.Subcase 1. Consider an initial condition with support out of given by (3.2). In this case, the sequence converges to a fixed point with support in (see Table 6 and Figure 8).Subcase 2. Consider an initial condition with support in and out of given by (3.3). We also obtain a sequence approaching to a fixed point with support in (see Table 7 and Figure 9).Subcase 3. Consider a parabolic initial condition with support in and out of given by (3.4). A similar behavior is obtained (see Table 8 and Figure 10).
Observe that in the three cases we obtain a sequence converging to different fixed points, although all of them have support . This fact leads to the hypothesis that if and is taken an initial condition whose support is not contained in , then the iterates of the discrete system converge to a fixed point with support in .

Table 6: Numerical results: Subcase 1.
Table 7: Numerical results: Subcase 2.
Table 8: Numerical results: Subcase 3.
Figure 8: Graphics sequence: Subcase 1.
Figure 9: Graphics sequence: Subcase 2.
Figure 10: Graphics sequence: Subcase 3.

#### 4. Stability

Again we take . In this section we prove that if , then some fixed points with support totally included in are stable. In fact, due to the numerical simulations, we believe that this result holds true also for . Some examples are also exhibited.

We consider the function where , is smooth on and for any . Sometimes, total or partial derivatives with respect to be denoted by subindices like, for example, in and below.

Lemma 4.1. Denote and assume that . Let be the unique solution corresponding to the initial data satisfying (2.2). Then

Proof. Let be the unique solution corresponding to the initial data such that . Multiplying (2.1) by we have Hence, These calculations are formal, but can be justified by using a cutoff function.

Let be such that , , , , and . We define by We note that . As and , we obtain that is a fixed point.

Denote by the norm in .

Theorem 4.2. Let . If then the fixed point is stable with respect to initial data satisfying (2.2) in the following sense: for any there exists such that if and , then

Proof. Let be the unique solution corresponding to the initial data such that . Denote . Multiplying (2.1) by z, we have
We use that Denote , . Thus
Now so
Hence Denote . Using (4.6), (4.13), and the inequality , for any (an interval of ) and (see [7, page 129]), we have It is known [4] that for any there exists such that if . Hence by Lemma 4.1, we have
We note that . Hence, for any there exists such that implies From here the stability follows.

Let us consider now particular examples of stable fixed points.

Example 4.3. Take and where , and , . Here and .
This function satisfies , , and so that (4.6) holds. Hence, by Theorem 4.2, is stable in the above sense.
Taking an initial condition close to , for instance, where , we obtain the numerical results collected in Table 9 (see also Figure 11).
If we take an initial condition closer to than before, for instance, where , then we obtain the numerical results collected in Table 10 (see also Figure 12).

Table 9: Numerical results: Stable fixes point.
Table 10: Numerical results: Stable fixed point.
Figure 11: Graphics sequence: Example 4.3.
Figure 12: Graphics sequence: Example 4.3.

Hence, the solution remains near this fixed point if the initial data is close to it.

Example 4.4. Take and where and , . Here and .
This function satisfies , and so that (4.6) holds.
Hence, by Theorem 4.2, is stable.

Example 4.5. Take and where , and , . Here , and .
This function satisfies , and so that (4.6) holds. Hence, by Theorem 4.2   is stable.
We can see that for a given more than one stable point exists. Finally, we shall give an example of an unstable fixed point for . Let us consider the fixed point given by and the initial condition where the constants and are chosen in such a way that is close enough. For and we obtain the numerical results given in Table 11 (see Figure 13). For and the numerical results are given in Table 12 (see Figure 14). These numerical results show that the fixed point is unstable.

Table 11: Numerical results: Unstable fixed point.
Table 12: Numerical results: Unstable fixed point.
Figure 13: Graphics sequence: Unstable fixed point.
Figure 14: Graphics sequence: Unstable fixed point.

#### 5. Conclusions

For problem (2.1) with , we have proved in Theorem 4.2 that if , then a class of fixed points with support inside the interval is stable. It would be interesting to prove such result for , as the numerical simulations suggest that this is true. It follows from such hypothesis that the unique fixed point with support not contained in the interval is not globally asymptotically stable, that is, not all solutions which are not fixed points converge to it as time goes to , although it seems that most of solutions behave in this way.

On the other hand, we have obtained in this paper numerical simulations of the solutions leading to the following conjectures.(1)Every solution of the equation converges to a fixed point as time goes to , that is, the omega-limit set of every trajectory is a fixed point.(2)Most fixed points with support in the interval are unstable.

These are interesting open problems to consider for the near future.

#### Acknowledgments

This work was partially supported by Spanish Ministerio de Ciencia e Innovación, Projects MTM2008-00088 and MTM2009-11820, the Consejería de Innovación, Ciencia y Empresa (Junta de Andalucía), Proyecto de Excelencia P07-FQM-02468, and the Consejería de Cultura y Educación (Comunidad Autónoma de Murcia), Grant 08667/PI/08.

#### References

1. R. G. Larson, The Structure and Rheology of Complex Fluids, Oxford University Press, New York, NY, USA, 1999.
2. P. Hébraud and F. Lequeux, “Mode-coupling theory for the pasty rheology of soft glassy materials,” Physical Review Letters, vol. 81, no. 14, pp. 2934–2937, 1998.
3. E. Cancès, I. Catto, and Y. Gati, “Mathematical analysis of a nonlinear parabolic equation arising in the modelling of non-Newtonian flows,” SIAM Journal on Mathematical Analysis, vol. 37, no. 1, pp. 60–82, 2005.
4. J. M. Amigó, I. Catto, Á. Giménez, and J. Valero, “Attractors for a non-linear parabolic equation modelling suspension flows,” Discrete and Continuous Dynamical Systems. Series B. A Journal Bridging Mathematics and Sciences, vol. 11, no. 2, pp. 205–231, 2009.
5. J. M. Amigó, Á. Giménez, F. Morillas, and J. Valero, “Attractors for a lattice dynamical system generated by non-Newtonian fluids modeling suspensions,” International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, vol. 20, no. 9, pp. 2681–2700, 2010.
6. E. Cancès and C. Le Bris, “Convergence to equilibrium of a multiscale model for suspensions,” Discrete and Continuous Dynamical Systems. Series B. A Journal Bridging Mathematics and Sciences, vol. 6, no. 3, pp. 449–470, 2006.
7. H. Brezis, Análisis Funcional, Alianza Editorial, Madrid, Spain, 1984.