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

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.


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 p p t, σ , t ≥ 0, σ ∈ Ê, and D p t ≡ D p t, • α/T 0 |σ|>σ c p t, σ dσ.Furthermore, χ I denotes the characteristic function of the interval I and δ 0 is the Dirac delta function supported at the origin.
We point out the important fact that the function p t, σ is a probability density function for every fixed t, so it must satisfy that p t, σ ≥ 0 and Ê p t, σ dσ 1, for any t ≥ 0. Existence of solutions of 1.1 is given in 3 .In 4 a simplified version of this model is studied, in which the term b t ∂ σ p 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 ∂ 2 σσ p and using a Simpson's quadrature for the integral D p .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., b t 0 .Also, for simplicity we take σ c 1 and T 0 1.In such a case it is well known 3 that all the probability densities with support in the interval −1, 1 are fixed points.Also, if α ≤ 0.5, then there no more fixed points, and if α > 0.5, then there exists a unique fixed point not contained inside the interval −1, 1 .It is known that this fixed point is asymptotically stable 6 .
In this paper, we prove that if α < 0.5, a subset of probability densities with support in the interval −1, 1 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 α > 0.5 the unique fixed point with support not contained inside the interval −1, 1 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 −1, 1 .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 −1, 1 are unstable, since if we take an arbitrary fixed point and a close initial data, usually the numerical solution converges to distant fixed point.

Implicit Finite Difference Scheme
We shall consider the simplified model that results from 1.1 after taking b t ≡ 0. Also, we set σ c 1 for simplicity, although our results can be obtained for an arbitrary σ c > 0. Hence, we consider the equation Here, χ I is the characteristic function on the interval I.
Consider an initial condition p 0 satisfying Then it is proved in 3, Theorem 1.1 that problem 2.1 has a unique solution p t, σ satisfying the following properties for all T > 0: where ν T exists for any T > 0. The last property implies that the average stress τ t Ê σp t, σ dσ belongs to L ∞ 0, T .From δ 0 ∈ H −1 Ê and the above properties, it follows that ∂ t p ∈ L 2 0, T; H −1 Ê .
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 h Δσ.Let −S . . .
σ i i∈ ⊂ Ê be such that σ 0 0 and |σ i 1 − σ i | h.We fix h so as σ −2n 1 −1 and σ 2n 1 1.In our numerical simulations, we will always take p 0 in such a way that its support is contained in −S, S with S hN sufficiently large, and we truncate this unbounded grid for |i| > N. If p 0 has an unbounded support, we can take S large enough in order to achieve S −S p 0 dσ ≈ 1.Similarly, we discretize the time domain with a time step Δt s.The resulting numerical grid is illustrated in Figure 1.
We shall implement an implicit scheme in finite differences by using the following approximations.Let p n h p n i i∈ −N,N denote the vector grid satisfying p n i ≈ p t n , σ i .Since the partial differential equation 2.1 is of evolution type, we approximate ∂ t p at the point t n , σ i by the forward difference In order to approximate ∂ σσ p at the point t n , σ i , we define the operator A N h : Ê 2N 1 → Ê 2N 1 by A N h : 1/ 2h 2 A N , where and we take This choice is necessary for the matrix A N h to satisfy the condition where being β 1, if N is even, and being β 2, if N is odd.This follows from A N h v i 0, for any i.
Observe also that u, v Ê 2N 1 u l 1N S , which together with equality 2.7 is a key fact to prove Theorem 2.1 see below .
We define also the matrix B N h : 1/ 2h B N , where , where the superscript "T" stands for "transpose".
The integral α/T 0 |σ|>1 p t, σ dσ is approximated by which is a slight variation of the Simpson formula.The only difference is that the weight for the terms p 2n 1 and p −2n 1 is equal to zero instead of 1. Observe that we can also set .

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

2.12
This function is continuous, its support is −3/ 2c , 3/ 2c , and it attains the maximum value at σ 0, which is equal to c > 0. Also, the integral Ê δ c σ dσ equals 1 see Figure 2 .Let us take c such that 3/ 2ch n c is even and let n c < 2n 1 so as 3/ 2c < 1, that is, the support of δ c is strictly included in the interval −1, 1 .Then at σ i ih this function looks as

2.14
Denote δ c i δ c σ i .Substituting all the above discretization formulas into the partial differential equation 2.1 , we get the following algebraic system for each integer n ≥ 0:

2.15
Bearing in mind that A N h 1/ 2h 2 A N and denoting by κ s/ 2h 2 s/ 4h 2 , the difference scheme system 2.15 can then be equivalently written in the form which can be expressed in matrix form as where

2.18
Observe that the matrix C depends on the time step s, while M n depends on the steps h and s.It is easy to check that M n is strictly diagonally dominant for all n ∈ AE.In particular, this fact guarantees that M n is nonsingular and the system 2.17 has a unique solution for every n > 0. Recall also the following result.

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 T 0 1.We shall consider two cases depending whether α > 0.5 or α ≤ 0.5 with different initial conditions p 0 , which may or not have support intersecting the interval −1, 1 .We shall not plot cases when the support of p 0 is totally contained in −1, 1 , since in such a case p 0 would be a fixed point of the discrete system 2.15 .Therefore, we will focus on cases such that supp p 0 Case 1 α > 0.5 .In this case all the probability densities whose support is totally contained in −1, 1 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 −1, 1 is given by where D D p f satisfies the equation D √ D α − 1/2 in order for the normalization constraint Ê p f 1 to hold true.Furthermore, this fixed point is stable see 6 .
In the next experiments, we shall take the constants α 0.8 and T 0 1 and we define the spacial and time mesh size Δσ h 0.001 and Δt s 0.1, respectively.Other choices of α > 0.5 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 p 0 with support out of −1, 1 given by

3.2
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 p 0 with support in and out of −1, 1 given by

3.3
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 −1, 1 given by

3.4
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 D p 0 ≈ 0 given by

3.5
The sequence approaches to a fixed point whose support is contained in −1, 1 see Table 4 and Figure 6 .This fact leads to the hypothesis that besides the fixed point 3.1 there are stable fixed points p for 2.1 with support contained in −1, 1 .

3.6
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.
We note that, usually, when we take an arbitrary fixed point with support in the interval −1, 1 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 α < 0.5.
Case 2 α ≤ 0.5 .In this case the probability densities with support totally included in −1, 1 are the unique fixed points.Let us take the constants α 0.3 and T 0 1.We define the spacial  and time mesh size Δσ h 0.001 and Δt s 0.1, 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 p 0 with support out of −1, 1 given by 3.2 .In this case, the sequence converges to a fixed point with support in −1, 1 see Table 6 and Figure 8 .Subcase 2. Consider an initial condition p 0 with support in and out of −1, 1 given by 3.3 .We also obtain a sequence approaching to a fixed point with support in −1, 1 see Table 7 and Figure 9 .Subcase 3. Consider a parabolic initial condition p 0 with support in and out of −1, 1 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 −1, 1 .This fact leads to the hypothesis that if α ≤ 0.5 and is taken an initial condition whose support is not contained in −1, 1 , then the iterates of the discrete system converge to a fixed point with support in −1, 1 .

Stability
Again we take T 0 1.In this section we prove that if α < 0.5, then some fixed points with support totally included in −1, 1 are stable.In fact, due to the numerical simulations, we believe that this result holds true also for α ≥ 0.5.Some examples are also exhibited.
We consider the function  where 0 < a < 1, b > 0, θ σ is smooth on |σ| > a and |θ σ | ≤ 2/ 1 − a 2 for any σ.Sometimes, total or partial derivatives with respect to σ be denoted by subindices like, for example, in θ σσ and p σ below.Lemma 4.1.Denote η 1/ 1 − a 2 and assume that α < 1/2η.Let p • be the unique solution corresponding to the initial data p 0 satisfying 2.2 .Then Proof.Let p • be the unique solution corresponding to the initial data p 0 such that D p 0 > 0. Multiplying 2.1 by θ we have  Hence,

4.4
These calculations are formal, but can be justified by using a cutoff function. Let

4.5
We note that v ∈ H 1 Ê .As D v 0 and Ê f σ dσ 1, we obtain that v is a fixed point.
Denote by then the fixed point v is stable with respect to initial data p 0 satisfying 2.2 in the following sense: for Proof.Let p • be the unique solution corresponding to the initial data p 0 such that D p 0 > 0.
Let us consider now particular examples of stable fixed points.

4.18
where C 0.05, D −0.01, and 4.20 where a 0.001 and c 0.02238, we obtain the numerical results collected in Table 9 see also Figure 11 .
If we take an initial condition closer to v than before, for instance,  where a 0.0001, and c 0.00057414, then we obtain the numerical results collected in Table 10 see also Figure 12 .

4.21
Hence, the solution remains near this fixed point if the initial data is close to it.where the constants c 1 and c 2 are chosen in such a way that p 0 is close enough.For c 1 0.495 and c 2 0.01 we obtain the numerical results given in Table 11 see Figure 13 .For c 1 0.49995 and c 2 0.0001 the numerical results are given in Table 12 see Figure 14 .These numerical results show that the fixed point p f is unstable.

Conclusions
For problem 2.1 with T 0 1, we have proved in Theorem 4.2 that if α < 0.5, then a class of fixed points with support inside the interval −1, 1 is stable.It would be interesting to prove such result for α ≥ 0.5, 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 −1, 1 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 −1, 1 are unstable.
These are interesting open problems to consider for the near future.

Table 10 :
Numerical results: Stable fixed point.