Raouf Boucekkine
email: raouf.boucekkine@univ-amu.fr

Fouad El Ouardighi
email: elouardighi@essec.fr

OPTIMAL GROWTH WITH POLLUTING WASTE AND RECYCLING

Keywords: Capital accumulation, sustainability, waste, recycling. JEL Classification. Q57, C61

We study an optimal AK-like model of capital accumulation and growth in the presence of a negative environmental externality in the tradition of Stokey (1998). Both production and consumption activities generate polluting waste. The economy exerts a recycling effort to reduce the stock of waste. Recycling also generates income, which is fully devoted to capital accumulation. The whole problem amounts to choosing the optimal control paths for consumption and recycling to maximize a social welfare function that notably includes the waste stock and disutility from the recycling effort. We provide a mathematical analysis of both the asymptotic behavior of the optimal trajectories and the shape of transition dynamics. Numerical exercises are performed to illustrate the analysis and to highlight some of the economic implications of the model.

Introduction

Since the Meadows report (1972), economic research on the limits of the contemporaneous growth regime and the design of optimal sustainable policies has been central to many agendas. While the very first attempts along these investigation lines have been directed at the issue of non-renewable resources (see for example, [START_REF] Stiglitz | Growth with exhaustible natural resources: Efficient and optimal growth paths[END_REF][START_REF] Hartwick | Intergenerational equity and investing rents from exhaustible resources[END_REF], substantial efforts have been devoted to studying the impact of pollution on growth and social welfare, especially since the '90s. A fundamental theoretical contribution to this topic is that of [START_REF] Stokey | Are there limits to growth[END_REF]. Within various optimal growth settings, Stokey studied the implications of pollution externalities for optimal capital accumulation (and therefore for optimal growth). A remarkable outcome of this study is the analysis of the AK economy with and without pollution externalities. In the latter case, the economy optimally starts on exponentially growing balanced growth paths with strictly positive growth rates. When pollution produces negative welfare losses, the economy no longer follows these virtuous paths; rather, it converges to optimal steady states (therefore, with zero growth). In other words, pollution drastically limits (optimal) growth. Several authors have extended Stokey's framework to account for more technological or ecological ingredients; [START_REF] Boucekkine | Technological vs. ecological switch and the environmental Kuznets curve[END_REF] is one of the most recent extensions. When pollution is irreversible (that is, when the environmental absorption capacity irreversibly declines above a certain pollution stock level), these authors show that the optimal relationship between income and pollution can take a much richer set of forms, although optimal exponential growth invariably vanishes under pollution externalities. This paper incorporates polluting waste and recycling into this type of optimal relationship analysis. More precisely, we are interested in polluting waste for which natural absorption takes extremely long. A characteristic example of such pollution is that of plastic waste, that is, bags, bottles, etc., which can be absorbed by the environment only after four to six centuries. Because of this extremely long time, no natural amenities can be reasonably assumed, and recycling efforts are required to avoid massive accumulation with harmful consequences in the long run. In the case of plastic waste, it seems that insufficient efforts have been deployed to prevent such a scenario, as evidenced by the emergence of several plastic vortexes in the oceans [START_REF] Kaiser | The dirt on ocean garbage patches[END_REF]. [START_REF] Jambeck | Plastic waste inputs from land into the ocean[END_REF] calculated that 275 million metric tons of plastic was generated in 192 coastal countries in 2010, with 4.8 to 12.7 million metric tons entering the ocean. In spite of the significant development of recycling and energy recovery activities, post-consumer plastic waste predominantly goes to landfill (PlasticsEurope, 2015). [START_REF] Jambeck | Plastic waste inputs from land into the ocean[END_REF] consider that without waste management infrastructure improvements, the cumulative quantity of plastic waste available to enter the ocean from land is predicted to reach 80 million metric tons by 2025.

Several authors have already modeled waste generation and recycling activities within economic frameworks. Of course, these models are largely found in the industrial organization literature (see for example, [START_REF] Martin | Monopoly power and the recycling of raw materials[END_REF][START_REF] Grant | Recycling and market power: A more general model and re-evaluation of the evidence[END_REF]. Nonetheless, waste and recycling are increasingly being examined from a more macroeconomic perspective. The recent macroeconomic literature contains studies of the impact of recycling on aggregate fluctuations (see De Beir et al., 2010, for example). There have also been several attempts to incorporate waste and recycling in a sustainability analysis, as we are doing in our paper, most often in decentralized equilibrium frameworks with technological progress (see for example, the recent paper by [START_REF] Fagnart | Quantitative versus qualitative growth with recyclable resource[END_REF].

Below we provide a first-best central planner analysis in line with Stokey's seminal paper incorporating waste (both from consumption and production) and recycling. Perhaps the work that is most closely related to ours is that of [START_REF] Lusky | A model of recycling and pollution control[END_REF]. Lusky also solved a central planner problem regarding waste and recycling. However, our problem differs from his along three essential dimensions. First, Lusky has strictly concave production functions with labor as a unique input, whereas capital accumulation is an essential feature of our model and we use an AK production technology as in Stokey. Second, recycling is modeled differently in [START_REF] Lusky | A model of recycling and pollution control[END_REF]: recycling produces a consumable good (so it increases consumption possibilities), while in our model recycling output goes to capital accumulation given that our focus is sustainable growth via capital accumulation, to which recycling contributes. Third, the social welfare functions are not the same. In particular, while in both papers waste produces negative externalities, recycling increases instantaneous utility in [START_REF] Lusky | A model of recycling and pollution control[END_REF] via the recycled consumption good whereas the recycling effort supposes a strictly concave welfare loss in our set-up.

The problem we consider is an infinite time horizon problem with two control variables (consumption and recycling effort) and two states (capital and stock of waste). The state equations are linear mainly due to the linearity of the production function and of waste generation processes. Using a version of the maximum principle, we can extract the (necessary and sufficient) optimality conditions, and study the asymptotic properties of the optimal paths. The economic implications of this analysis are then evaluated in light of the sustainability literature à la Stokey.

The paper is organized as follows. Section 2 gives the specifications of our central planner problem with waste and recycling. Section 3 analyzes the mathematical properties of the model. In Section 4, we perform some quantitative exercises and we bring out the main lessons we can draw from the point of view of sustainability. Section 5 concludes the paper.

Model

Assume that the social planner owns a stock of productive capital that provides a continuous flow of revenue, () aK t , where ( ) 0

Kt  denotes the capital stock at time t , and 0 a  , the constant marginal unit of revenue generated by the productive capital stock. The production function is linear in the stock of capital, as in [START_REF] Stokey | Are there limits to growth[END_REF] and [START_REF] Boucekkine | Technological vs. ecological switch and the environmental Kuznets curve[END_REF]  being the marginal wasting impact of the revenue generation process and the current consumption, respectively.

Though most of the polluting waste is generally related to productive processes [START_REF] Klassen | Plant-level environmental management orientation: The influence of management views and plant characteristics[END_REF], it is not clear whether the consumption process generates more waste than the revenue generation process. In the case of the plastics industry, polluting waste is predominant in either process depending on the market segment and the polymer type (PlasticsEurope, 2015).

In this regard, we assume that    . By construction (no capital accumulation), waste is generated only by consumption in [START_REF] Lusky | A model of recycling and pollution control[END_REF]. The same is assumed in the macroeconomic literature in line with De [START_REF] De Beir | Life cycle of products and cycles[END_REF], which relies heavily on the related industrial organization (e.g., [START_REF] Martin | Monopoly power and the recycling of raw materials[END_REF]: waste produced from consumption is used one period ahead as an input in the recycling sector. In this paper, waste can come from both consumption and production (or capital utilization), and some of the essential properties of optimal paths may depend on whether    .

To reduce the stock of polluting waste, denoted by ( ) 0 Wt  , the social planner may invest in recycling efforts ( ) 0 vt  over time. We assume that the waste generating processes and recycling operations are mutually independent so that the recycling efforts are nonproportional to the waste emissions. This assumption allows for unbounded recycling efforts, i.e., ( ) ( )

v t w t  
, to account for the possibility of reduction of past waste emissions. The environmental absorption capacity of polluting waste is approximated by zero. This approximation is consistent with the extremely long time needed for natural absorption of plastic waste. In economic terms, it implies that the social planner cannot benefit from any natural abatement of pollution waste1 .

Finally, a fixed proportion of recycled waste is supposed to generate additional revenues, () vt  , and therefore to positively influence the capital accumulation process, 10   being the marginal proportion of the recycled waste that adds to capital accumulation. An illustration of this assumption is related to the plastics industry, where 60 million tons of plastics diverted from landfills are equivalent to over 60 billion euros (PlasticsEurope, 2015).

In more elaborated industrial organization models of recycling, the recycling sector may produce profits (as in [START_REF] Martin | Monopoly power and the recycling of raw materials[END_REF], which are later redistributed to the owners. In our central planner setting, the idea is pretty much the same: recycling not only decreases the level of polluting waste, but it also generates an income, which contributes to capital accumulation.

Both functions of recycling help alleviate the sustainability problem faced by the economy.

Only abatement plays this role in [START_REF] Stokey | Are there limits to growth[END_REF] for example. If 0   , the income generation channel of recycling is shut down, and we are closer to Stokey's framework regarding the impact of pollution control instruments.

Based on these assumptions, the endogenous capital accumulation process is described as follows:

( ) ( ) ( ) K aK t c t v t     , 0 (0) 0 KK 
where a positive difference between the total revenues from capital and recycled waste, and current consumption results in investment in productive capital, while a negative difference leads to disinvestment. The initial endowment in productive capital is given by 0 0 K  .

The dynamics of polluting waste are given by:

( ) ( ) ( ) W aK t c t v t     , 0 (0) 0 PP 
where the recycling efforts are such that 0 W   , 0 t  , for a given initial stock of waste, 0 0 W  . Regarding the objective function of the social planner, we make the following assumptions.

At each period, the instantaneous social utility is given as the difference between the utility drawn from current consumption and the costs incurred from the stock of waste and the recycling efforts, respectively. The instantaneous utility from current consumption is a concave function, that is, ln ( ) ct . In addition, the stock of waste entails negative externalities such as environmental pollution and destruction of the biomass (e.g., [START_REF] Barnes | Invasions by marine life on plastic debris[END_REF]

( ) ( ) ( ) K aK t c t v t     , 0 (0) 0 KK  (2) ( ) ( ) ( ) W aK t c t v t     , 0 (0) 0 WW  (3)
As mentioned in the introduction, our social welfare function differs from Lusky's in that it includes a welfare loss due to the recycling effort and has no additional (recycled) consumption term. Note that we consider a pollution (via aggregate waste) negative externality while the macroeconomic literature of recycling and fluctuations does not (see for example De [START_REF] De Beir | Life cycle of products and cycles[END_REF]. The same comparison holds with the ecological sustainability literature (see [START_REF] Fagnart | Quantitative versus qualitative growth with recyclable resource[END_REF].

We now come to the mathematical resolution of the optimal control problem considered. We limit our presentation to the case of an interior solution.

Analysis

Skipping the time index for convenience, the current-value Hamiltonian is:

    22 ln 22 eW v H = c aK c v aK c v              (4) where () t   and () t   are costate variables, 1, 2 j 
, that evolve according to:

  r a a       (5) r eW   (6)
Necessary conditions for optimality are:

11 0 c H = c c           (7) 0 v H = v v            (8)
Because the stock of polluting waste has a negative marginal influence on the social planner's objective function, its implicit price should be non-positive, i.e., 0

  . Along with 0   ,
this should result in strictly positive consumption and recycling effort respectively in ( 7) and ( 8).

The Legendre-Clebsch condition of concavity of the Hamiltonian with respect to the control variables is satisfied, as the Hessian: Plugging the respective expressions of c and v from ( 7) and ( 8) in (4) results in the maximized Hamiltonian:

    2 2 0 1 ln 1 2 2 eW H = aK                
from which the Hessian matrix: is negative semi-definite. This ensures that the necessary conditions are also sufficient for optimality. 

Plugging the value of c and v from ( 7) and ( 8) in ( 2) and (3), respectively, the equations:

  1 K aK           (9) W aK            (10)
along with ( 5) and ( 6), form the canonical system in the state-costate space.

We now prove that as in Stokey (1998), the dynamic system will not converge to a balanced growth path despite the assumed impacts of recycling on the stock of pollution and on income generation. We first show the existence of steady states and then assess stability.

Proposition 1. In the case of a farsighted social planner, i.e., ra  , and 1

  , the steady state is unique and given by:

  T S S S S K W c v                           1 1 1 1 1 1 1 1 1 a r a r a r r a r a e a r                                                            (11)
where

  ar         
, and the superscript ' S ' stands for steady state. Otherwise, the steady state is zero.

Proof. Equating the RHS of ( 9)-( 10)-( 5)-( 6) to 0 and solving by identification and substitution, we get:

    11 S a ar                     11 S ar ar            
and S K and S W as given in (11). Note that ra  implies that

    1 r a a         ,
which allows for a feasible steady state stock of waste. Conversely, for any ra  , the steady state stock of waste is not feasible. Plugging the above expressions in ( 7) and ( 8), respectively, and simplifying, yields S c and S v in (11). From (11), it can be shown that the limiting transversality conditions are satisfied for the saddle-paths because: This ensures the uniqueness of the globally optimal solution.  From the canonical system ( 5)-( 6)-( 9)-( 10), the isoclines of 0 K  and 0 W  are given by:

    1 lim ( ) ( ) lim 0 1 rt rt tt e e t K t                         2 lim ( ) (
      0 1 S S S K a r eW r a r K ra a r a eW           (18)       0 1 1 S S S W a r eW r a r K a r a r eW                (19)
Using ( 18)-( 19) and ( 11), Figure 1 describes the sensitivity of the state variables at the steady state to the parameters. From Figure 1, the marginal wasting impact of both the revenue generation process ( ) and the consumption process (  ) on the capital stock is negative. In contrast, the influence of these processes on the stock of waste is different because an increase in the marginal wasting impact of the revenue generation (consumption) process results in less (more) waste stock. In contrast, an increase in the marginal impact of the recycling effort on the capital accumulation process ( ) reduces both the capital stock and the waste stock. Note that the influence of the discounting rate and the marginal revenue coefficient on the stock of capital and of waste is similar to that of the marginal polluting impact of revenue generation ( ) and the marginal impact of the recycling effort on the capital accumulation process ( ). The marginal revenue coefficient ( a ) has a similar impact on the stock of capital and the stock of waste as the marginal wasting impact of the consumption process (  ). Finally, a greater cost coefficient of the waste stock ( e ) lowers the waste stock and does not affect the capital stock, and vice- versa. We now move to the stability analysis and investigate the structure of the associated stable manifolds.

Proposition 2. The steady state exhibits a (local) two-dimensional stable manifold if the social planner is moderately farsighted (i.e.,

 

1 ra  

), and a one-dimensional stable manifold otherwise.

Proof.

To analyze the stability of the steady state, we compute the Jacobian matrix of the canonical system ( 9)-( 10)-( 5)-( 6), that is:

        2 2 2 2 2 2 1 0 0 1 00 0 0 a J a r a a e r                                                 
Given that  and  are evaluated at their steady state value, we compute the determinant:

            2 1 1 2 1 1 J ae a r ae a r                               
which has a positive value for a moderately farsighted social planner (

 

1 ra   ), and negative otherwise. As shown in [START_REF] Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF], a negative determinant is a necessary and sufficient condition for the Jacobian matrix to have one negative eigenvalue and either three positive eigenvalues or one positive eigenvalue and two with positive real parts. In terms of dynamic behavior, this corresponds to the case of a one-dimensional stable manifold.

In the case of a moderately farsighted social planner (

 

1 ra   ), we use Dockner's formula [START_REF] Dockner | Local stability in optimal control problems with two state variables[END_REF] to determine the sum of the principal minors of J of order 2 minus the squared discounting rate, that is:

            2 2 2 11 11 ar a a r e a a r e                                 
The necessary and sufficient conditions that ensure that two eigenvalues have negative real parts and two have positive real parts, which corresponds to the case of a two-dimensional stable manifold, are 0 J  and 0  . The sign of  is negative, which implies that a two- dimensional stable manifold (saddle-point) exists in the case of a moderately farsighted social planner.  According to Proposition 2, if the social planner is not moderately farsighted, the zero steady state cannot be reached from some or all initial states. Conversely, if the social planner is moderately farsighted, a saddle-point, either zero or positive, exists that can be reached. We may dig deeper in the analysis and draw further properties of the optimal paths, in particular about the existence of oscillatory transitions to the steady states. This gives a rough idea of the ability of the model to generate non-monotonic optimal trajectories and fluctuations, a central issue in the recent recycling-related macroeconomic literature [START_REF] De Beir | Life cycle of products and cycles[END_REF].

Proposition 3. Assuming a farsighted social planner (i.e., ra  ), for any given a , e , and 1   , there exists a threshold 0   such that for any   , the convergence to the steady state is oscillatory.

Proof. To determine whether the optimal path is monotonic or follows cyclical motions, we compute the expression [START_REF] Dockner | Local stability in optimal control problems with two state variables[END_REF]:

            2 2 2 11 4 1 8 1 1 ar J a a r e ae a r                                               
A positive (negative) sign of  indicates that convergence to the saddle-point is monotonic (spiraling) near the steady state. Because the sign of  is ambiguous, a limit value analysis highlights the role played by a , e ,  ,  and  in the sign of  for a given ra  (Tab. 1).

The results suggest that  is generally positive, which implies that convergence to the saddle-point is monotonic near the steady state in general. However, given that 0 lim 0

  
and lim

     , it can be shown that 0     .     2 2 2 0 12 lim a e             lim a       2 2 0 lim e a a r     lim e       2 2 0 lim a a r        1 lim        2 0 lim 2 a a r e           lim              2 2 0 1 lim lim 8 6 2 a a r ae e r r a e a a                         
Tab. 1. Limit value analysis Therefore, we conclude that there exists a threshold 0   such that for any   , we have 0  .  According to Proposition 3, if the social planner is farsighted, convergence to the locally stable steady state is either monotonic or oscillatory, depending on the magnitude of the wasting impact of the revenue generation process.

Proposition 4. The saddle paths of the control and state variables are:

3 1 1 3 2 4 1 () t t SS ct B B e B B e
(20)

3 1 1 3 2 4 () t t SS v t B B e B B e (21) 3 1 56 () t t S K t K B e B e (22) 3 1 78 () t t S W t W B e B e ( 23
)
where 18 ,..., BB are constants of integration and 13 ,0 . Proof. The linear approximation of the canonical system ( 5)-( 6)-( 9)-( 10) around the steady state is: Dockner's formula (1985), the four eigenvalues associated with the Jacobian matrix of the canonical system are:

  r a a       r eW         2 2 SS SS K aK                        2 2 SS SS W aK                      Using
2 2 2 3 4 2 12 1 1 1 1 4 1 2 4 2 2 2 4 2 2 ar r r r r J a a r e
and S  and S  are given in the proof of Proposition 1. As expected, two eigenvalues, 1 and 3 , have a negative sign and two have a positive sign, 2 and 4 . Choosing the negative roots for convergence, the time paths of the costate and state variables are written as: These equations involve 10 unknowns (i.e., 18 ,..., BB , (0), ( 0) ) that can be solved with the 10 following equations, which are drawn from the above expressions and the linearized versions of ( 5)-( 6)-( 9)-(10):

12 (0) S BB 34 (0) S BB 0 5 6 S K K B B 0 7 8 S W W B B         1 1 3 2 3 4 0 SS r a r a B r a B a B B                       1 3 3 4 7 8 0 S S r r B r B e W B B                      1 2 3 4 1 5 3 6 1 2 3 4 2 0 S S S S S S S B B B B aK a B a B B B B B                                            1 2 3 4 5 6 1 2 3 4 1 7 3 8 2 0 S S S S S S S B B B B a K B B B B B B B B                                               2 2 2 2 2 1 1 3 2 3 4 7 8 2 0 S S S r a r a B r a B a a r B B e W B B                                               2 2 2 2 2 1 3 3 4 1 7 3 8 0 SS r r B r B e rW r B r B                 
This system of equations can be solved numerically. Finally, using ( 7)-( 8) yields ( 20)-( 23).

Numerical example

We now give a numerical example to suggest the economic insight that could be gained from our model. To this end, we use the following parameter values (Table 2). (Table 3).

The main relationship to be investigated is that between the optimal stocks of capital and pollution to uncover a possible environmental Kuznets curve as in [START_REF] Stokey | Are there limits to growth[END_REF] or [START_REF] Boucekkine | Technological vs. ecological switch and the environmental Kuznets curve[END_REF]. Regarding this relationship, the results of our experiments are reported in Figure 2 below. Several remarkable features can be deduced from these figures. The first result is that the relationship between capital (or income) and pollution is non-monotonic, if recycling generates additional income (that is 0   ). This is true when production is more wasteful than consumption (Figure 2a) and in the opposite case (Figure 2).

In both cases, we observe that the stock of pollution decreases while capital rises initially, but in the last stage of convergence to the steady state both stocks go down. In contrast, also in both cases, the relationship is permanently monotonic if recycling does not generate income (that is, 0   ): the pollution stock decreases to its steady state value whereas capital increases to its corresponding stationary value. No turning point is observed in such cases.

2.a. Production more wasteful than consumption 2.b. Consumption more wasteful than production Fig. 2. Phase-planes of stock of capital and stock of waste It is worth pointing out two interesting connected outcomes at this stage. First, the kind of non-monotonicity we get is not of the environmental Kuznets curve kind, which is therefore not the rule in our optimal recycling model. Second, this non-monotonicity only arises when recycling generates additional income. To understand these remarkable properties, it is useful to have a look at the optimal control time paths. They are reported in Figure 3.

In all the parametric configurations, optimal consumption starts at a relatively low level, in contrast to recycling. We are therefore in a situation where both controls act as substitutes. In our calibrated economy, priority is given to pollution control in the short run to decrease the stock of polluting waste as quickly as possible. In the transition to the steady states, consumption increases steadily to its corresponding stationary value when recycling has no additional income ( 0.1   ) or rises then decreases to the steady state value when recycling does generate income ( 0   ). Whatever the parameterization, recycling effort always starts at a high value and follows an unambiguous monotonically decreasing path to its steady state value. With these elements in hand, one can rationalize the optimal relationship between income and pollution we have uncovered. When recycling generates income, capital increases markedly in the initial stage of transitional dynamics because recycling is highest at this stage, and income from recycling goes entirely to capital accumulation. This significant increment in the stock of capital ends up boosting production and therefore consumption (notably with respect to the case where recycling does not generate income). In the medium term, consumption is quite high while recycling drops sharply. In the last stage of the transitional dynamics, the consumption level is still markedly higher than in the case 0   , while income from recycling is much lower than in the initial stage of the transition dynamics (because the amount recycled in also much lower). The conjunction of the two latter forces pushes capital down in the ultimate adjustment stage given the law of motion of capital (2), generating the turning point observed and described above.

Therefore, the non-monotonicity observed is mainly due to the timing of optimal recycling, which massively takes place in the initial periods. The welfare implications of such timing can be seen in Figure 4, where the social welfare function is computed for different time horizons.

The initial sharp drop in social welfare is due to the initial intense recycling period.

Conclusion

In this paper, we have studied the sustainability of a Stokey-inspired AK model, one that considers a negative environmental externality that arises both from production and consumption. Instead of the typical abatement technologies, we present novel recycling modelling where recycling also generates income that is fully devoted to capital accumulation.

We have studied the qualitative properties of the resulting optimal control problem, notably in terms of optimal asymptotic states, stability and transition. We have also worked out a numerical example and got some highly interesting economic results in comparison with the seminal framework of Stokey, both in terms of the optimal pace of recycling and the relationship between income and pollution.

In particular, the role played by recycling as an income generator is crucial in the sense that it gives rise to a contraction of both the consumption and capital stock in the long run after an expansion phase. Whether polluting waste is predominantly due to production or consumption, when recycling generates additional income, greater consumption and lower capital stock are obtained in the long run compared with the situation when recycling does not create additional income. In parallel, when recycling generates additional income, greater recycling effort and lower stock of waste are resulted in the long run than when recycling has no additional income.

Of course, this is just a preliminary investigation; further analyses involving control-state constraints and alternative specifications within the same class of models along with alternative calibrations are needed to corroborate and complement this study.

  negative definite. This guarantees a maximum of the Hamiltonian.Lemma 1. The necessary conditions are sufficient for optimality.

Fig. 1 .

 1 Fig. 1. Sensitivity of the steady state to the parameter values

   ,  and  , we show that various structurally different solutions are possible. The induced steady state values are reported in the following table

  effort Fig. 3. Time paths of control variables (logarithmic time scale)

Fig. 4 .

 4 Fig. 4. Overall social utility

  . .

	The flow of revenue allows for a certain consumption level, ( ) 0 ct  . The revenue generation
	process and the consumption decisions are both assumed to generate polluting waste,
	( ) 0 wt  , where ( ) w t	( ) aK t  	( ) c t	, ,0

  . These negative externalities are valued as an increasing convex function of the stock of waste, that

	is,	2 ( ) 2 eW t	,	0 e  . Lastly, the recycling effort generates an increasing quadratic cost,
	denoted by	2 ( ) 2 fv t	,	0 f  . Without loss of generality, we set	1 f  .
	Denoting the discounting rate by	0 r  , and assuming an infinite planning horizon, the social
	planner's optimal control problem is:
		U		22 ( ) 22 ( ) eW t v t    0 ln ( ) d rt e c t t      	(1)
	subject to:

This assumption is optimistic because accumulation of polluting waste can lead to negative environmental absorption capacity that might create additional negative externalities (see ElOuardighi et al.,

2014)


Laurent Chupin
email: laurent.chupin@math.univ-bpclermont.fr

Sébastien Martin
email: sebastien.martin@parisdescartes.fr

Stationary Oldroyd model with diffusive stress: mathematical analysis of the model and vanishing diffusion process

Keywords: viscoelastic fluid, Oldroyd model, diffusive stress MSC. 76D03, 76A10, 35Q35

In this paper we address the stationary Oldroyd model with diffusive stress in two ways: first we present the mathematical analysis of the model. Second, as the link between the diffusive model and the standard one is questionable from the mathematical point of view, we discuss, by means of numerical simulations, the behaviour of the model with respect to vanishing diffusion. In particular, numerical results in 2D or 3D suggest that the solution of the diffusive model converges to the solution of the non-diffusive model at order 1 in the L 2 -norm and H 1 -norm.

Introduction 2 Introduction

The Oldroyd model describes the behaviour of a viscoelastic fluid. Its principle is built upon a description of the shear stress that interpolates between a purely viscous contribution and a purely elastic contribution. Let Ω be a bounded open set in R 3 and f be a given vector function on Ω. We are looking for a vector function u : Ω → R 3 , a scalar function p : Ω → R and a symmetric tensor function σ : Ω → R 3×3 , representing the velocity, the pressure and the elastic extra-stress of the fluid satisfying the following set of equations:

           Re (u • ∇u) -(1 -r)∆u + ∇p = div σ + f, in Ω, div u = 0, in Ω, We (u • ∇σ + g a (∇u, σ)) + σ -D ∆σ = 2rD(u),
in Ω, u = 0, on ∂Ω, D ∂ n σ = 0, on ∂Ω.

(1)

Here, Re denotes the Reynolds number that quantifies the inertial effects in the fluid flow; We is the Weissenberg number related to a relaxation time that characterizes the elasticity of the fluid; D is a diffusive parameter of the elastic stress; r ∈ [0, 1] is an interpolation parameter : cases r ∈ (0, 1) are often referred as Jeffreys models whereas case r = 1 is referred as the Maxwell model. Besides, D(u) (resp. W(u)) denotes the symmetric (resp. skew-symmetric) part of the velocity gradient. The function g a , -1 ≤ a ≤ 1, is a bilinear mapping related to the total derivative, in which the parameter a interpolates between the so-called upper-convected model (a = 1) and lowerconvected model (a = -1). Note that the case a = 0 is known as the so-called corotational model. The function g a is defined as

g a (∇u, σ) = W(u) • σ -σ • W(u) + a(D(u) • σ + σ • D(u)).
Remark 1 In standard derivations of Oldroyd model from kinetic models for dilute polymers, the diffusive term D∆σ is routinely omitted, on the grounds that it is several orders of magnitude smaller than the other terms in the equation. It physically corresponds to a centre-of-mass diffusion term in the dumbell models and it is in the range of about 10 -9 to 10 -7 when the macroscopic length-scale of the domain is of order 1, see [START_REF] Barrett | Existence of global weak solutions to some regularized kinetic models for dilute polymers[END_REF].

The Oldroyd model also has a transient version which, notably, has been extensively analysed in many ways: the transient version of the standard Oldroyd model (without diffusion: D = 0) has been the subject of intensive studies: in particular P.-L. Lions & N. Masmoudi [START_REF] Lions | Global solutions for some Oldroyd models of non-Newtonian flows[END_REF] proved a global existence result of weak solutions for any data in the corotational case only. In the general case, C. Guillopé & J.-C. Saut [START_REF] Guillopé | Existence results for the flow of viscoelastic fluids with a differential constitutive law[END_REF] proved the existence and uniqueness of local strong solutions; besides, if the fluid is not too elastic and if the data are sufficiently small, then solutions are global. Then L. Molinet & R. Talhouk [START_REF] Molinet | On the global and periodic regular flows of viscoelastic fluids with a differential constitutive law[END_REF] proved that the smallness assumption on the elasticity of the fluid could be relaxed. The transient version of the diffusive Oldroyd model (with diffusion: D > 0) has been studied by a few authors: recently P. Constantin & M. Kliegl [START_REF] Constantin | Note on global regularity for twodimensional Oldroyd-B fluids with diffusive stress[END_REF] proved the existence of global strong solutions in 2D for the Cauchy problem of the diffusive Oldroyd-B model (i.e. a = +1) and uniqueness of the solution among a class of strong solutions. Notice also that other regularizations of the standard Oldroyd model have been studied, see in particular [START_REF] Barrett | Existence of global weak solutions to some regularized kinetic models for dilute polymers[END_REF].

Let us now discuss the stationary case. To our knowledge, the only available result is due to M. Renardy [START_REF] Renardy | Existence of slow steady flows of viscoelastic fluids with differential constitutive equations[END_REF] and focuses on the standard model (without diffusion) only: existence and uniqueness of strong solutions is proved under the assumption of small regular data. The method used by M. Renardy is based on a reformulation of the Oldroyd model as a "Newtonian generalized" fluid: the contribution of the stress div σ is expressed as an implicit function of the velocity field u and, then, an iterative scheme is built upon this fully nonlinear system. When dealing with stress diffusion, the mathematical analysis of the model can be approached with a completely different framework, as the diffusive contribution in the stress equation drastically changes the mathematical properties of the system. So far, the method that we present in the present paper is based on a classical weak formulation and then on energy estimates. Let us underline the main differences with the standard Oldroyd model studied by M. Renardy: first, the diffusive model makes it possible to handle with irregular data and boundary conditions ; second, in the corotational case, the smallness of the data is not needed anymore.

Let us give some motivation for the study of the vanishing diffusion process which consists in investigating the asymptotic regime D → 0 in the diffusive Oldroyd model. As a matter of fact, the vanishing viscosity process describes the transition between a (second order) diffusive model towards a (first order) transport model. Thus, boundary conditions have to be discussed. For this, let us draw some parallel with the behaviour of scalar transport models, as the Oldroyd equation is a convection equation. In the framework of scalar conservation laws, a classical approach to deal with boundary conditions is to introduce the so-called vanishing viscosity method. But the analysis highly depends on the type of boundary conditions that are considered. In the case of Dirichlet boundary conditions, introducing artificial diffusion leads, when passing to the limit, to a problem in which boundary conditions cannot be assumed to be satisfied pointwise because of the degeneracy of the operator. Nevertheless, Bardos, LeRoux & Nedelec [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF] exhibited a condition that expresses the way the boundary conditions are satisfied at the limit. In particular, it includes the fact that when the flux is monotone, the boundary condition is active when information enters the domain whereas it is not active when it goes out. In the end, Bardos, LeRoux & Nedelec provide a so-called weak entropy formulation, still based of the vanishing viscosity method, which leads to a well-posed problem: at the limit (i.e. without diffusion but with Dirichlet boundary conditions that have to be understood in a weak sense), the problem admits a unique solution. When dealing with Neumann boundary conditions, well-posedness in a similar set-ting is questionable: does the vanishing viscosity method may lead to the selection of a unique solution? To our knowledge, this is an open question. In the context of the diffusive Oldroyd model, as the order of mag nitude of D is about 10 -7 , the vanishing diffusion process has to be addressed. Besides, at least for the Oldroyd model, classical solutions do exist in the non-diffusive model [START_REF] Renardy | Existence of slow steady flows of viscoelastic fluids with differential constitutive equations[END_REF] with suitable boundary conditions. Dealing with a system of conservation law (such as the Oldroyd model) is much more difficult than dealing with a scalar conservation law but numerical simulations may play a role in the analysis of the asymptotic regime D → 0 by providing ideas on the behaviour of the diffusive model with respect to the standard one.

The present paper is composed of three sections. In Section 2 we give the weak formulation of the problem, the mathematical framework and the main results: existence, uniqueness, regularity. Section 3 is devoted to the numerical simulation of the models: an algorithm is proposed and benchmark tests are considered in order to investigate the vanishing viscosity process.

Weak formulation of the stationary problem and main theoretical results

For the sake of physical relevance, the mathematical analysis presented in this section is presented in dimension 3. However, results are valid in dimension 2.

The variational formulation of Problem (1) is written

                     Find (u, σ) ∈ V × W such that, for all (v, τ ) ∈ V × W , Re Ω (u • ∇u) • v + (1 -r) Ω ∇u : ∇v + Ω σ : D(v) = f, v , We Ω (u • ∇σ + g a (∇u, σ)) : τ + Ω σ : τ + D Ω ∇σ : ∇τ = 2r Ω D(u) : τ , (2)
where the functional spaces V and W are defined by

V = {u ∈ H 1 0 (Ω) 3 ; div u = 0}, W = {σ ∈ H 1 (Ω) 3×3 ; σ = T σ},
and where •, • denotes the duality bracket between H -1 and H 1 0 . We will use the following norms:

u 2 V = Ω ∇u 2 , σ 2 W = Ω σ 2 + ∇σ 2 .
The main theorem which is proved in this article concerns an existence result for Problem [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF]. In its general form, it requires some assumptions on the data. For this, we introduce the following constants:

C (I) := 8|a|C 2 Ω We f H -1 min(1 -r, D) 2 , C (II) := √ 2r min(1 -r, D) 4|a|C 2 Ω We 1 -1 -C (I) ,
where C Ω is a constant which only depends on the domain Ω.

Remark 2

1. Note that C (II) is defined provided C (I) ≤ 1.
2. For a = 0, we have by continuity C (I) = 0 and

C (II) = √ 2r f H -1 min(1 -r, D) . 3. Constant C Ω is related to the Sobolev injection W ⊂ L 4 (Ω) 3×3 : ∀σ ∈ W, σ L 4 ≤ C Ω σ W .
4. The only physical parameter that is not involved in the definition of constants C (I) and C (II) is the Reynolds number Re. In other words, the existence result that is further described does not depend on the value of the Reynolds number.

Theorem 1 Let Ω be a Lipschitz bounded open set in R 3 and f ∈ H -1 (Ω) 3 . Let Re ≥ 0, We ≥ 0, 0 < r < 1, -1 ≤ a ≤ 1 and D > 0.

• Existence. If C (I) ≤ 1 then, there exists a solution (u, σ) of problem (2) which satisfies

2r u 2 V + σ 2 W ≤ C 2 (II) . (3)
Moreover there exists p ∈ L 2 (Ω) such that (u, p) satisfies (1) 1,2,3 in the sense of distributions.

• Uniqueness. Problem (2) admits at most one solution if one of the following conditions is satisfied:

a) f H -1 is small enough; b)
Re and We are small enough.

• Regularity. If Ω is of class C ∞ and if each component of f belongs to C ∞ (Ω) then each component of any solution of (2) belongs to C ∞ (Ω)
and the considered solution satisfies (1) in a classical sense.

Let us mention two corollaries.

Corollary 2 Problem (2) admits a unique solution if one of the following conditions is satisfied:

a) f H -1 is small enough; b)
Re and We are small enough.

Corollary 3 If a = 0, there exists a solution for all data.

The proof is decomposed into four parts. In section 2.1 we show the existence for the weak formulation (2) using a Galerkin approximations and compactness results to perform the limit. The existence of a pressure is obtained by De Rham theory. In section 2.2, we prove the uniqueness of the solution and, in section 2.3, we investigate the regularity of the weak solutions and prove that, if the data are regular, so is the solution.

Proof of existence

As a preliminary, V × W is endowed with the scalar product (•,

•) V ×W de- fined by ((u, σ), (v, τ )) V ×W = 2r(u, v) V + (σ, τ ) W .
As V and W are separable Hilbert spaces, we consider a countable orthonormal basis (w k ) k∈N in the space V , and a countable orthonormal basis (s k ) k∈N in the space W . We use the notation V k := span(w 1 , ..., w k ) and W k := span(s 1 , ..., s k ). For each fixed integer k ∈ N, we would like to define an approximate solution (u k , σ k ) of ( 2) by

u k = k i=0 α i,k w i , σ k = k i=0 β i,k s i , satisfying the variational problem Re Ω (u k • ∇u k ) • v k + (1 -r) Ω ∇u k : ∇v k + Ω σ k : D(v k ) = f, v k , We Ω (u k • ∇σ k + g a (∇u k , σ k )) : τ k + Ω σ k : τ k + D Ω ∇σ k : ∇τ k = 2r Ω τ k : D(u k ), (4)
for all (v k , τ k ) ∈ V k ×W k . Equations (4) form a system of nonlinear equations for α 1,k , ..., α k,k , β 1,k , ..., β k,k , and the existence of a solution of this system is not obvious. We use the following lemma:

Lemma 1 Let X be a finite dimensional Hilbert space with scalar product (•, •) and norm • , and let P be a continuous mapping from X into itself such that

∃R > 0; ∀ξ ∈ X ( ξ = R =⇒ (P (ξ), ξ) ≥ 0). (5)
Then there exists ξ ∈ X, ξ ≤ R, such that P (ξ) = 0.

The related proof, based on the Brouwer fixed point theorem, can be found in [17, p.166]. We only note that the result which is proved in [START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF] corresponds to the case where the inegality (P (ξ), ξ) ≥ 0 in the assertion ( 5) is a strict inequality. The case of a large inequality holds too, with the same proof.

We apply this lemma to prove the existence of (u k , σ k ) as follows: let X be the space defined as X := V k × W k , endowed with the scalar product inherited from V × W . Let P k the mapping from X into itself defined by, for all ((u, σ),

(v, τ )) ∈ X 2 , (P k (u, σ), (v, τ )) X = 2r Re Ω (u • ∇u) • v + (1 -r) Ω ∇u : ∇v + Ω σ : D(v) -f, v + Ω σ : τ + D Ω ∇σ : ∇τ + We Ω (u • ∇σ + g a (∇u, σ)) : τ -2r Ω D(u) : τ .
The continuity of P k is obvious. Let us show that condition (5) holds. Denoting ξ = (u, σ) ∈ X we have

(P k (ξ), ξ) X = 2r(1 -r) Ω ∇u 2 + Ω σ 2 + D Ω ∇σ 2 + We Ω g a (∇u, σ) : σ -2r f, u .
From the Hölder inequality and the Sobolev injection W ⊂ L 4 (Ω) 9 , there exists a positive constant C Ω which only depends on the domain, such that

σ L 4 ≤ C Ω σ W . Thus, the contribution g a satisfied Ω g a (∇u, σ) : σ ≤ 2|a| Ω ∇u σ 2 ≤ 2|a| ∇u L 2 σ 2 L 4 , which yields Ω g a (∇u, σ) : σ ≤ 2|a|C 2 Ω u V σ 2 W . (6)
We deduce that

(P k (ξ), ξ) X ≥ 2r(1 -r) u 2 V + min(1, D) σ 2 W -2|a|C 2 Ω We u V σ 2 W -2r f H -1 u V .
Using the definition of the norm ξ 2 X = 2r u 2 V + σ 2 W , we obtain

(P k (ξ), ξ) X ≥ min(1 -r, D) ξ 2 X - √ 2|a|C 2 Ω We √ r ξ 3 X - √ 2r f H -1 ξ X .
We write this inequality as

(P k (ξ), ξ) X ≥ ξ X (-α ξ 2 X + β ξ X -γ).
• If a = 0 (i.e. α > 0), we deduce that (P k (ξ), ξ) X may be positive for some value of ξ if the discriminant β 2 -4αγ is nonnegative; in this case we have

(P k (ξ), ξ) X ≥ 0 ⇐⇒ β -β 2 -4αγ 2α ≤ ξ X ≤ β + β 2 -4αγ 2α .
In particular, by denoting

C (II) = β -β 2 -4αγ 2α ,
we obtain

ξ X = C (II) =⇒ (P k (ξ), ξ) X ≥ 0.
We also note that C (II) > 0 (except for the case f = 0 (i.e. γ = 0) where we can use instead of C (II) any constant C such that 0 < C < β α ).

• If a = 0 (i.e. α = 0), we easily have

(P k (ξ), ξ) X ≥ 0 ⇐⇒ ξ X ≥ γ β .
We then use

C (II) = γ β .
As for the case a = 0 we have C (II) > 0 (except for f = 0 (i.e. γ = 0) where we can use instead of C (II) any positive constant).

Finally, assumptions given by Lemma 1 are satisfied taking R = C (II) as soon as β 2 -4αγ ≥ 0, which is equivalent to the condition C (I) ≤ 1 introduced in Theorem 1. By Lemma 1, we deduce that, for any k ∈ N, Equations ( 4) admit a solution (u k , σ k ) which satisfies

2r u k 2 V + σ k 2 W ≤ C 2 (II) . (7)
This estimate [START_REF] Constantin | Note on global regularity for twodimensional Oldroyd-B fluids with diffusive stress[END_REF] implies that the sequence (u k , σ k ) remains bounded in V × W . Thus there exists some (u, σ) ∈ V × W and a subsequence (still denoted by k) such that (u k , σ k ) ⇀ (u, σ) for the weak topology of V × W , as k → +∞. As H 1 ⊂ L 2 with compact injection, (u k , σ k ) → (u, σ) for the strong topology of L 2 (Ω) 3 × L 2 (Ω) 9 . The consequences are twofold:

• We can pass to the limit in all the terms of Equations ( 4), and deduce that (u, σ) is a solution of (2) ;

• Due to the usual property of the weak limit, we have

2r u 2 V + σ 2 W ≤ lim inf k→+∞ 2r u k 2 V + σ k 2 W ≤ C 2 (II) .
This concludes the proof of existence of a solution (u, σ) of problem (2) satisfying Eq. (3). Now let us prove the existence of a pressure field associated to the incompressibility condition. For a solution (u, σ) of problem (2), we have, for all v ∈ D(Ω) such that div v = 0,

Re(u • ∇u) -(1 -r)∆u -div σ -f, v = 0.
By De Rham theorem, there exists a pressure p ∈ D ′ (Ω) such that Eq. ( 1) holds in D ′ (Ω). The regularity of ∆u, u • ∇u, div σ and f implies that the pressure p is more regular. For instance, by Sobolev embeddings, we have

u ∈ V ⊂ L 6 (Ω) 3
and since ∇u ∈ L 2 (Ω) 3×3 , then the convective term u•∇u belongs to L 3/2 (Ω). Thus, we have

∇p = ∆u + div σ + f -Re(u • ∇u) ∈ H -1 (Ω).
This implies that p ∈ L 2 (Ω), see [17, p.14].

Proof of uniqueness

Let (u 1 , σ 1 ) and (u 2 , σ 2 ) be two solutions of (2) and introduce the difference (u, σ) = (u 2 , σ 2 )-(u 1 , σ 1 ). By subtraction we obtain, for all (v, τ

) ∈ V ×W , Re Ω (u • ∇u 2 + u 1 • ∇u) • v + (1 -r) Ω ∇u : ∇v + Ω σ : D(v) = 0, We Ω (u • ∇σ 2 + u 1 • ∇σ + g a (∇u, σ 2 ) + g a (∇u 1 , σ)) : τ + Ω σ : τ + D Ω ∇σ : ∇τ = 2r Ω τ : D(u).
Taking (v, τ ) = (2r u, σ) and adding the two equations, we obtain

2r(1 -r) u 2 V + min(1, D) σ 2 W ≤ -2rRe Ω (u • ∇u 2 ) • u -We Ω (u • ∇σ 2 + g a (∇u, σ 2 ) + g a (∇u 1 , σ)) : σ.
By definition of the bilinear function g a , we have

Ω g a (∇u, σ 2 ) : σ = (a + 1) Ω (∇u • σ 2 ) : σ + (a -1) Ω (σ 2 • ∇u) : σ.
Using the Hölder inequality and the Sobolev injection W ⊂ L 4 (Ω) 9 we deduce that

Ω g a (∇u, σ 2 ) : σ ≤ 2C 2 Ω u V σ 2 W σ W .
Estimating the other terms as in Equation ( 6), we obtain

2r(1 -r) u 2 V + min(1, D) σ 2 W ≤ 2rC 2 Ω Re u 2 V u 2 V + 3C 2 Ω We σ 2 W u V σ W + 2|a|C 2 Ω We u 1 V σ 2 W .
As the solutions (u 1 , σ 1 ) and (u 2 , σ 2 ) satisfy Equation (3), we deduce that

2r(1 -r) u 2 V + min(1, D) σ 2 W ≤ C 2 Ω C (II) √ 2r 2rRe u 2 V + 3We √ 2r u V σ W + 2|a|We σ 2 W . (8)
From Young inequality we deduce

3We √ 2r u V σ W ≤ 2rRe u 2 V + 9We 2 4Re σ 2 W .
Consequently, Estimate (8) reads 2rA u 2 V + B σ 2 W ≤ 0 where A and B are given by

A = (1 -r) - 4ReC 2 Ω C (II) √ 2r , B = min(1, D) - C 2 Ω C (II) √ 2r 9We 2 4Re + 2|a|We .
Using the fact that C (II) tends to 0 as f H -1 tends to 0, we deduce that for f H -1 small enough, the coefficients A and B are positive. In the same way, if Re is small enough then A > 0 and, if We is small, then B > 0. These inequalities imply u V = σ W = 0, which means (u 1 , σ 1 ) = (u 2 , σ 2 ).

Strong solutions

As previously, we notice that u ∈ L 6 (Ω) 3 , so that u • ∇u ∈ L 3/2 (Ω) 3 . We also have div σ ∈ L 2 (Ω) 3 and if f is regular then the regularity of the Stokes problem, see [START_REF] Cattabriga | Su un problema al contorno relativo al sistema di equazioni di Stokes[END_REF], implies that

u ∈ W 2,3/2 (Ω) 3 , p ∈ W 1,3/2 (Ω).
In the same way, u • ∇σ and g a (∇u, σ) belong to L 3/2 (Ω) 3×3 . The regularity of the Laplace problem implies that

σ ∈ W 2,3/2 (Ω) 3×3 .
Repeating such a process, we find that u ∈ W s,3/2 (Ω) 3 , p ∈ W s-1,3/2 (Ω) and σ ∈ W s,3/2 (Ω) 3×3 for any s ∈ N. The proof is concluded.

Numerical investigation of the vanishing diffusion process

Numerical simulations have been performed in order to illustrate the influence of the diffusive parameter on the stationary solution. After giving some comprehensive details on the numerical method, we will present two simulations:

• The flow past a cylinder in a plane channel is a benchmark case from the literature. The vanishing diffusion process is faced with its formal limit. Let us recall that the non-diffusive model investigated by M. Renardy [START_REF] Renardy | Existence of slow steady flows of viscoelastic fluids with differential constitutive equations[END_REF] deals with Dirichlet boundary conditions for the velocity field with an additional assumption: the normal component of the velocity is zero. Thus, no fluid enters or exits the domain, avoiding possible compatibility problems with boundary conditions on the stress tensor, which do not need to be imposed in this particular situation. In the diffusive model, there is no such limitation on the Dirichlet boundary condition. The benchmark case allows us to investigate the convergence of the solution with diffusion to the solution without diffusion, relaxing the assumption on the boundary condition.

• The driven cavity is considered. In this problem, the normal component of the velocity is zero but the discontinuity of the prescribed velocity field at the corners of the domain lead to a loss of regularity of the solution. Therefore, strong solutions do not exist and the framework given by M. Renardy is not valid anymore. Nevertheless, the finite element solution is well-defined (for both the diffusive case and the non-diffusive case) and, thus, we investigate the convergence of the solution with diffusion to the solution without diffusion, despite the lack of regularity of the solutions.

The numerical simulations described in Section 3 have been led in the corotationnal case, i.e. a = 0. The Johnson-Segalman issue is addressed in the concluding remarks in Section 4.

Algorithm

The numerical method is based on the transient version of the Oldroyd model: stationary states are attained as the long-time solution of the evolutive model. If the solution converges in time, then a stationary state is obtained through this process. Note also that for a = 0 existence of a stationary solution is ensured and, if the data are sufficiently small, the stationary solution is unique. As a consequence, in a non-uniqueness regime, if the algorithm converges, then it is likely to select a stationary solution among all the possible ones.

In order to describe the application of the characteristics method for the particular convective terms in the viscous and elastic equations, we consider the artificial time dependent velocity field u(t, x) = u and an artificial time dependence in all involved functions, by introducing the notation φ(t, x) = φ(x). By convenience, we immediatly drop he overscripts •. Then, we can introduce the following material (or total derivative) operator associated to the velocity field u:

D Dt = ∂ ∂t + u • ∇ . (9)
In this way, the initial equations take the following form.

• the momentum equation can be written as

Re Du Dt -(1 -r)∆u + ∇p = div σ + f ;
• the incompressibility equation is unchanged : div u = 0 ;

• the constitutive equation can be written as We Dσ Dt + g a (∇u, σ) + σ -D ∆σ = 2rD(u).

Next, by approximating the total derivative by the characteristics method, the previous equations can be approximated by

• momentum equation: Re u n+1 -u n • κ n ∆t -(1 -r)∆u n+1 + ∇p n+1 = div σ n+1 + f n+1 ;
• incompressibility equation: div u n+1 = 0 ;

• constitutive equation:

We

σ n+1 -σ n • κ n ∆t + g a (∇u n+1 , σ n+1 ) +σ n+1 -D∆σ n+1 = 2rD(u n+1 ),
where the index n + 1 denotes the approximation at the artificial time t n+1 of the introduced time dependent functions and ∆t denotes an artificial time step. Moreover, the function κ n is defined at each spatial point of the domain by κ n (x) = κ(t n+1 , x; t n ), that denotes the position at time t n of the point placed in x at time t n+1 and moving along the integral path (characteristic curve) defined by the velocity field u, so that κ n (x) can be obtained from the solution of the final value ODE problem:

dκ dτ (t n+1 , x; τ ) = u τ, κ(t n+1 , x; τ ) , κ(t n+1 , x; t n+1 ) = x.
In order to pose an approximated problem (which is actually dependent on the parameter ∆t), The semi-discretized problem is completed with the boundary conditions now written in terms of u n+1 and σ n+1 , instead of u and σ and, also, some initial condition for both the velocity field and the elastic constraint. In this way, a stationary solution of the Oldroyd model with diffusive stress is obtained by letting t tend to ∞.

Before describing in details the effective problems to solve, let us highlight the fact that the discretized problem is strongly coupled. As we are only interested in the long time asymptotic solution, a very simple trick for easier computations consists in splitting the coupling between the velocity field and the elastic stress. In other words, we consider an explicit term in the 3.2 Flow past a cylinder in a plane channel in 2D

Numerical profiles

As a benchmark case, we consider a viscoelastic flow past a cylinder in a plane channel. The cylinder ensures the fact that classical Poiseuille profiles are not a solution of the system in the viscoelastic case. This case was studied experimentally in [START_REF] Verhelst | Visco-elastic flow past circular cylinders mounted in a channel: experimental measurements of velocity and drag[END_REF], and computationally in [START_REF] Alves | The flow of viscoelastic fluids past a cylinder: finite-volume high-resolution methods[END_REF][START_REF] Fan | Galerkin/least-square finiteelement methods for steady viscoelastic flows[END_REF] and recently in [START_REF] Goyal | Direct simulations of spherical particles sedimenting in viscoelastic fluids[END_REF].

The channel-with-cylinder geometry is built upon a channel with length L = 40R, height H = 4R. Introducing Cartesian ccordinates, the channel is thus described by x ∈ (0, L), z ∈ (-2R, 2R). The cylinder has a radius R and its center of mass is located at x = 10R and z = 0, preserving geometric symmetry with respect to the line {z = 0}. No-slip boundary conditions are imposed at {z = ±2R} and a parabolic (Poiseuille) velocity profile is imposed at the inlet {x = 0} and at the outlet {x = L}: it is characterized either by the maximal velocity value U max or the input/output flux Q. In all the computations, R has been set to 0.25 and the flux has been set to Q = 1.0 which corresponds to a maximal velocity U max = 2Q/(3R) = 8/3 at the inlet/outlet. Considering the benchmark case studied in [START_REF] Goyal | Direct simulations of spherical particles sedimenting in viscoelastic fluids[END_REF], numerical simulations have been performed with Re = 0.67, We = 1.0, r = 0.6, a = 0 and different values of D. Note that the focus on the corotational model (i.e. a = 0) guarantees the existence of a stationary solution for both the diffusive and standard models. The study of the Johnson-Segalman (i.e. a = 0) model will be discussed in concluding remarks.

Figure 1 represents the solution of the stationary Oldroyd model without diffusion, near the cylinder, computed with an unstructured mesh made of 136879 triangles and 69342 vertices. This solution is the reference solution when investigating the behaviour of the model with respect to the diffusive parameter.

Error analysis

Figure 2 illustrates the influence of the diffusion parameter D on the velocity field and the elastic stress tensor, as it tends to 0. The benchmark case D 1). Rates of convergence are numerically derived (see Figure 2). Linear regressions have been addressed for small values of D (see Figure 2). It shows that the error in the L 2 -norm and the H 1 -norm behaves at order 1 for both the velocity and the elastic constraint.

u (D) -u (0) L 2 u (D) -u (0) H 1 σ (D) -σ (0) L 2 σ (D) -σ (0) H 1 1e-04 3.2e-3 1.
We point out the fact that the previous observations related to the order of approximation are still valid in the non-inertial regime (e.g. Re = 0.0, r = 0.5, We = 1.0) or inertial regime (e.g. Re = 10.0, r = 0.5, We = 1.0).

Let us highlight the fact that one could expect that difficulties arise in the vinicity of boundaries identified as inlets, i.e. where the fluid enters the domain. Indeed Dirichlet boundary conditions on the velocity field, at the inlet, enforce the information on the stress tensor to propagate inside the fluid domain because of the transport process. As the stress tensor is transported from the inlet into the fluid domain, Neumann boundary conditions on the elastic stress tensor do not seem to be the most relevant choice in terms of numerical stability. However, numerical simulations do not reveal any instability or boundary layer phenomenon. Let us notice that other boundaries (outlets corresponding to velocity fields expressing the fact that the fluid goes outside the domain, and lateral boundaries corresponding to zero velocity field) are not identified as a source of problems, which is also confirmed by numerical simulations.

Driven cavity in 2D

Numerical profiles

The driven cavity problem is considered. The fluid domain is a unit square. Boundary conditions on the velocity field are prescribed: u = 0 is imposed on lateral and bottom boundaries whereas u = (-1, 0) is imposed on the top boundary. The interest of this problem relies on the loss of regularity of the velocity field, even in the case of the linear Stokes problem: the solution does not belong to H 1 (Ω) anymore, because of the discontinuity of the Dirichlet boundary conditions located at the corners of the geometry [START_REF] Cai | An error estimate for two-dimensional Stokes driven cavity flow[END_REF]. Nevertheless, the finite element solution is well defined. Therefore, we observe the behaviour of the elastic tensor with respect to the loss of regularity.

Numerical simulations have been performed with Re = 10.0, We = 1.0, r = 0.5, a = 0 and different values of D. Simulations evidence several facts: Figure 3 reveals a large vortex whose typical size corresponds to the size of the domain: this is due to the inertial effects. The loss of regularity is observed on the pressure profiles: singularities are revealed at the top corners of the domain, at which extremal values are reached. Besides, the influence of the diffusion on the tensor profiles is important, playing the role of a regularizing effect, as it is expected, see Figures 4 and5. Furthermore, velocity profiles are little affected by the diffusion parameter D, as long as it is small enough (less than 10 -3 ).

In the end, the driven cavity test suggests that the loss of regularity does not affect the convergence of the solution towards the solution of the Oldroyd problem without diffusion.

Error analysis

The numerical tests have been led for various values of the Reynolds number Re, Weissenberg number We and retardation parameter r, see Figures 678910111213. In all cases, the following estimates hold:

σ (D) -σ (0) L 2 = O(D), σ (D) -σ (0) H 1 = O(D), u (D) -u (0) L 2 = O(D), u (D) -u (0) H 1 = O(D)
. Note that even if the velocity fields do not belong to H 1 , the lack of regularity is localized at the corner of the domain. In particular, the loss of regularity is the same for both u (D) and u (0) so that the difference between the two velocity fields is likely to belong to H 1 . This suggests that, despite the lack of regularity of the solution that prevents us from falling into the scope of the work of Renardy [START_REF] Renardy | Existence of slow steady flows of viscoelastic fluids with differential constitutive equations[END_REF], the Oldroyd model is a rigorous limit of the diffusive Oldroyd model as the vanishing diffusive solution converges to the non-diffusive solution in L 2 and even in H 1 .

Driven cavity in 3D

Numerical profiles

The driven cavity problem is considered. The fluid domain is a unit cube. Boundary conditions on the velocity field are prescribed: u = 0 is imposed on lateral and bottom boundaries whereas

u(x, y, 1) = 4 x (1 -x) y (1 -y) e x
is imposed on the top boundary. Numerical simulations have been performed with a structured mesh, size 30 × 30 × 30, in the inertial regime Re = 4.0, elastic regime We = 1.0 and r = 0.5 and different values of D.

Streamlines associated to the velocity field are shown on Figures 14-16.

Error analysis

As in the 2D case, error analysis in 3D by means of numerical simulations shows that the solution of the diffusive problem is close to the solution of the non-diffusive problem at order 1 in the L 2 -norm or in the H 1 -norm. More precisely, the behaviour of the error with respect to D can be described as follows, see Figures 17 and18:

• Figure 17 or 18, zone I. For values of D less than 10 -9 , the computed error is persistant as the machine precision, when computing the error, has been attained.

• Figure 17 or 18, zone II. The convergence result associated to the asymptotic regime D → 0 is ensured by the behaviour of the error corresponding to small values of D, namely D ∼ 10 -9 -10 -3 . In this range, the magnitude of the order of convergence for both the velocity field and the elastic tensor is identified as being equal to 1.

• Figure 17 or 18, zone III. For values of D greater that 10 -3 , which lack physical relevance anyway, the corresponding solution cannot be assumed to be close to the solution of the non-diffusive problem as the diffusive effect on the elastic tensor is predominant over the memory effect.

Hence, the computations suggest that the convergence order with respect to D does not depend on the dimension of the problem.

Concluding remarks

As stated in Theorem 1, existence of a stationary solution of the diffusive model is guaranteed if

C (I) := 8|a|C 2 Ω We f H -1 min(1 -r, D) 2 ≤ 1. (10)
Uniqueness of the solution is guaranteed if the source term is small enough or if the Weissenberg and Reynolds numbers are small enough. It is also noted that examples of nonuniqueness exist. They are even the source of wellknown phenomenon like shear banding, see [START_REF] Olmsted | Johnson-Segalman model with a diffusion term in cylindrical Couette flow[END_REF]. For example in [START_REF] Guillopé | Global existence and one-dimensional nonlinear stability of shearing motions of viscoelastic fluids of Oldroyd type[END_REF] or in [START_REF] Boyer | Numerical study of viscoelastic mixtures through a Cahn-Hilliard flow model[END_REF] it is shown that for r > 8/9 (and D = 0), there may be three solutions for the stationary Oldroyd problem, forming in the case of a Couette flow a number of shear bands. Besides, as presented in Section 3, convergence in strong topologies in 2D and 3D has been illustrated on benchmark tests, suggesting that the formal asymptotics related to the vanishing diffusion process may be valid at least in some realistic physical situations. However, the numerical experiments presented in Section 3 have been led with the corotational case only, i. e. a = 0, and for moderate Weissenberg numbers. Therefore, it is natural to question the validity of this conjecture for different values of a and W e, which is discussed in the two next subsections.

On the Johnson-Segalman model

Before addressing the vanishing viscosity process for the Johnson-Segalman model, we can question whether condition expressed by Eq. ( 10) is technical or necessary when a = 0, since the asymptotics only makes sense if the stationary solution is defined for any value of D. The Johnson-Segalman

On the high Weissenberg regime

The high Weissenberg regime does not fall into the range of validity of the study, in particular when considering condition expressed by Eq. ( 10) that ensures existence of a solution. Therefore, at this stage, well-posedness of the problem with large values of the Weissenberg number is an open problem. However, from the numerical point of view, simulations with the algorithm proposed in this paper failed. We argue two possible directions that may justify this observation: first, ill-posedness of the problem could explain the numerical instabilities that are observed in the simulations ; second, the splitting procedure of the algorithm may be the cause of the instabilities as it tends to segregate the two steps of the computations when the Weissenberg number is large, leading to unphysical oscillations of the solution of the transient model and, in the end, numerical instabilities which prevent us from defining a stationary state.

⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ +0.8 ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ +0.5 ⊠ ⊠ ⊠ ⊠ +0.1 0.0 -0.1 -0.5 ⊠ ⊠ ⊠ ⊠ -0.8 ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ -1.0 ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠

Figure 1 :

 1 Figure 1: Solution of the stationary Oldroyd model near the cylinder, with Re = 0.67, r = 0.6, a = 0, We = 1.0, D = 0.0. (A): norm of the velocity field, (B): pressure p, (C): elongational component of the elastic tensor σ 11 -σ 22 , (D): shear component of the elastic tensor σ 12 .

Figure 2 :

 2 Figure 2: Flow past a cylinder in a plane channel in 2D: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 3 : 5 - 5 σ 11 -σ 22 5 - 5 σ 11 -σ 22 5 - 5 σ

 35511551155 Figure 3: Solution of the stationary Oldroyd model for the driven cavity problem with Re = 10.0, r = 0.5, a = 0, We = 1.0, D = 0.0: (t.) Norm of the velocity field, (b.) pressure field.

Figure 4 :Figure 5 :

 45 Figure 4: Solution of the stationary Oldroyd model for the driven cavity problem with Re = 10.0, We = 1.0, r = 0.5, a = 0. Elongational component of the elastic tensor σ 11 -σ 22 . (t.): D = 0.000, (m.): D = 0.001, (b.): D = 0.100. 29

Figure 6 :

 6 Figure6: Lid-driven cavity in 2D with Re = 0.0, We = 0.5, r = 0.2, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 7 :

 7 Figure7: Lid-driven cavity in 2D with Re = 0.0, We = 0.5, r = 0.5, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 8 :

 8 Figure8: Lid-driven cavity in 2D with Re = 0.0, We = 1.0, r = 0.2, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 9 :

 9 Figure9: Lid-driven cavity in 2D with Re = 0.0, We = 1.0, r = 0.5, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 10 :

 10 Figure10: Lid-driven cavity in 2D with Re = 5.0, We = 0.5, r = 0.2, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 11 :

 11 Figure11: Lid-driven cavity in 2D with Re = 5.0, We = 0.5, r = 0.5, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 12 :

 12 Figure12: Lid-driven cavity in 2D with Re = 5.0, We = 1.0, r = 0.2, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

Figure 13 :Figure 14 :

 1314 Figure13: Lid-driven cavity in 2D with Re = 5.0, We = 1.0, r = 0.5, a = 0: influence of the diffusive term on the error on the velocity (l.) and the elastic constraint (r.).

2 Figure 15 : 12 -Figure 16 :

 2151216 Figure 15: Solution of the stationary Oldroyd model for the 3D driven cavity problem with Re = 4.0, We = 1.0, r = 0.5, a = 0, D = 0.0: cutting planes sampling the elongational components of the elastic tensor (t.) σ 11 -σ 22 , (m.) σ 11 -σ 33 and (b.) σ 22 -σ 33 . 36

Table 1 :

 1 Flow past a cylinder in a plane channel in 2D: computation of the error with respect to the diffusive parameter for the velocity and the elastic constraint, in the L 2 -norm and H 1 -norm.

	6e-1	2.2e-1	6.6e+1

evidences the convergence of the solution corresponding to the diffusive model to the solution of the non-diffusive model, in classical L 2 -and H 1 -norms (see Table

Table 2 :

 2 Simulations of the Johnson-Segalman model with diffusive constraint for the 2D benchmarks: (P) stands for the benchmark related to the flow past a cylinder in a plane channel and (C) stands for the driven cavity problem. Symbol means that a stationary state has been captured whereas symbol ⊠ means that the algorithm has failed in capturing a stationary state.

source term of the viscous equation instead of an explicit one: div σ n+1 is replaced by div σ n . Besides, in order to use classical approximation spaces, the incompressibility condition is treated by duality. Thus, the effective problems to solve, at each time step t n+1 , now reads

with the following operators

and the source terms defined as

Now we have to deal with the spatial discretization of the weak formulation. For this, we consider a triangular mesh of the domain Ω. The first two equations may be analyzed as a Stokes system which is, therefore, treated with P + 1 -P 1 or P 2 -P 1 finite elements (for the velocity and pressure fields), in order to ensure the inf-sup condition at the discrete level. When dealing with a diffusive term D > 0, the Oldroyd equation is treated with P 1 finite elements which naturally approximates H 1 (Ω). For comparison needs, when computing the solution without diffusion D = 0, the Oldroyd equation is treated with P 0 finite elements which naturally approximates L 2 (Ω) corresponding to the regularity of the stress tensor. model with diffusive constraint was particularly studied in [START_REF] Olmsted | Johnson-Segalman model with a diffusion term in cylindrical Couette flow[END_REF] in a cylindrical Couette flow. But in the most general case, i.e. without assumption on the geometry and boundary condition, does existence of a solution depend on a choice of parameters that should belong to a limited range of admissible values in order to satisfy Eq. ( 10) (or a sharper condition)? Or can the condition expressed in Eq. ( 10) be weakened or even purely forgotten in order to ensure existence of a sationary solution? Supplementary computations have been led and proven to be helpful in the understanding of the diffusive Johnson-Segalman model:

• Table 2 tends to show that Eq. ( 10) is certainly necessary: too large values of |a| need to be counterbalanced by large values of D for the algorithm to succeed in defining a stationary state.

• In both benchmark problems, computations with a = ±0.5 and a = ±1.0 (with large values of D) reveal velocity profiles which are quite different from those observed with smaller values of a.

• The case a = ±0.1 is somehow intriguing as the solution with diffusive constraint is defined for every value of D. It suggests that, for small values of |a| Eq. ( 10) may not be necessary. Alternatively the discretization process (in time and space) of the Oldroyd system may induce some relaxation of the condition on the resulting finite-dimensional Oldroyd system.

We would like to underline that these numerical observations are mostly compatible with theoretical results related to the standard Oldroyd models [START_REF] Guillopé | Existence results for the flow of viscoelastic fluids with a differential constitutive law[END_REF][START_REF] Lions | Global solutions for some Oldroyd models of non-Newtonian flows[END_REF][START_REF] Molinet | On the global and periodic regular flows of viscoelastic fluids with a differential constitutive law[END_REF]]: global-in-time solutions exist in the case a = 0 without any assumption and in the case a = 0 under smallness assumption on the data ; relaxing this smallness assumption in the Johnson-Segalman model drastically weakens the existence result, as the existence of local-in-time solutions only can be guaranteed. These results, related to the transient version of the standard model, shed a light on both theoretical and numerical results: in our algorithm, stationary states are obtained by letting time tend to +∞, which is possible only if a global-in-time solution of the corresponding transient problem exists. In the vanishing diffusion limit, condition expressed by Eq. ( 10) may express the degeneracy of the behaviour of solutions with respect to time in the Johnson-Segalman models.