Sofiene Dellagi
email: sofiene.dellagi@univ-lorraine.fr

Anis Chelbi
email: anis.chelbi@planet.tn

Wajdi Trabelsi
email: wajdi.trabelsi@icn-groupe.fr

Joint integrated

Keywords: Integrated production-maintenance strategies, Optimization, Variable production rates, Degradation, and Smoothing production rate

published or not. The documents may come

INTRODUCTION AND LITTERATURE REVIEW

Manufacturing companies must manage several functional areas successfully, such as production, maintenance, quality and marketing. One of the keys to success consists in treating all these services simultaneously. In order to ensure an efficient coordination between them, managers and decision makers have to consider a global systemic approach integrating the interactions between parts or all of those complementary functions. In this perspective, several researchers have investigated ways of integrating maintenance and production after decades during which those two functions had been studied separately. Works on maintenance policies started with Barlow et al. [START_REF] Barlow | Optimum Preventive Maintenance Policies[END_REF] and continued with a huge number of contributions as it can be seen in a survey on maintenance models by Wang [START_REF] Wang | A survey of maintenance policies of deteriorating systems[END_REF]. During the last two decades, several companies have realized that strategies dissociating maintenance and production were ineffective. The need for developing new integrated maintenance-production strategies became evident. Brandolese et al. [START_REF] Brandolese | Production and maintenance integrated planning[END_REF] proposed a strategy for maintaining a multiple machines production system. The planning consists in programming the execution date of each task and the machine that should perform it. They integrated preventive maintenance tasks in the planning placing them as close as possible to the optimal maintenance periods.

Chelbi et al. [START_REF] Chelbi | Analysis of a production/inventory system with randomly failing production unit subjected to a minimum required availability level[END_REF] developed a mathematical model in order to determine both the buffer stock size and the preventive maintenance period for an unreliable production unit which is subject to regular preventive maintenance of random duration. The optimization of the preventive maintenance period and the buffer size was also the problem addressed by Pal et al. [START_REF] Pal | A mathematical model on EPQ for stochastic demand in an imperfect production system[END_REF] in the context of an imperfect production system that may shift to an 'out-of-control' state after a random time. Non-conforming products are reworked incurring an additional cost. They also considered the situation where both the buffer size and the production rate are the decision variables. Rezg et al. [START_REF] Rezg | Optimal strategy of inventory control and preventive maintenance[END_REF] developed an analytical model and a numerical procedure which allow determining a joint optimal inventory control and age-based preventive maintenance policy for a randomly failing production system. Recently, Tambe and Kulkarni [START_REF] Tambe | A superimposition based approach for maintenance and quality planoptimization with production schedule, availability, repair time anddetection time constraints for a single machine[END_REF] developed a selective maintenance and quality control decision optimization framework considering the production schedule of the machine. They derive an optimal maintenance decision, consisting in one of three actions (repair, replace or do-nothing) for the system components along with the optimal sample size, the acceptance number and the time between samples, taking into account the optimal production schedule. The authors used a genetic algorithm in order to solve the problem.

In the same framework, some researchers have taken into account several external constraints. In fact, integrated maintenance-production strategies which take into consideration subcontracting have been studied by Dellagi et al. [START_REF] Dellagi | Preventive maintenance of manufacturing systems under environmental constraints[END_REF]. They developed and optimized a maintenance policy incorporating subcontractor constraints. They demonstrated through a case study, the influence of the subcontractor constraints on the optimal integrated maintenance-production strategy. Dealing also with the same subject, Dahane et al. [START_REF] Dahane | Effects of extension of subcontracting on a production system in a joint maintenance and production context[END_REF] studied analytically the problem of the integration of subcontracting activities and determined the optimal number of subcontracting tasks to be performed during a maintenance cycle. Recently, Nourelfath et al. [START_REF] Nourelfath | An integrated model for production and preventive maintenance planning in multi-state systems[END_REF] dealt also with the same problem of integrating preventive maintenance and production planning, for a production system composed of a set of parallel components. A subsequent work by Nourelfath and Châtelet [START_REF] Nourelfath | Integrating production, inventory and maintenance planning for a parallel system with dependent components[END_REF] assumed the presence of economic dependence and common cause failures in the production system. Mifdal et al. [START_REF] Mifdal | Joint Optimization Approach of Maintenance and Production Planning for a Multiple-Product Manufacturing System[END_REF] considered the same problem in the case of multiple-products manufacturing systems.

Looking at the literature on integrated maintenance-production policies, we noticed that the influence of the production rate on the system degradation in presence of a random demand over a finite planning horizon was rarely addressed in depth. Very few works dealt with this issue. For example, the study of Hajej et al. [START_REF] Hajej | Optimal integrated maintenance/production policy for randomly failing systems with variable failure rate[END_REF] which presents an analytical stochastic optimization model based on the operational age concept, revealed the significant influence of the production rate on the deterioration of the manufacturing system and consequently on the integrated production-maintenance policy. Later, Hajej et al. [START_REF] Hajej | Production/Maintenance Policies Optimization with Operational Age Concept in a Subcontracting Constraint[END_REF] dealt with the same problem by integrating a subcontracting constraint in order to recycle a certain quantity of returned products. In these two studies, the authors started by establishing an optimal production plan which minimizes the total inventory and production cost taking into consideration the subcontractor constraint. Then, using this optimal production plan, they derived an optimal maintenance schedule which minimizes the total maintenance cost. Their approach is sequential. Moreover, they do not consider the influence of having differences between quantities produced in successive periods and the impact of such differences on labour and set-up related costs.

The approach adopted in this study favours the fact of smoothing the production rate between periods in order to avoid the disadvantages of a significant variation (mainly in terms of setup costs). This paper is organized as follows: In the next section we specify the targeted contributions of this work. Then the problem and the working assumptions are described in section 3. Section 4 is dedicated to the working assumptions and the development of the analytical model. A numerical example is presented in section 5 to illustrate the proposed modelling approach and discuss different trade-offs that are considered. Finally, conclusions and potential future work are provided in Section 6.

TARGETED CONTRIBUTIONS

Several industries such as steel production exhibit fluctuation of the production rate between production periods. This yields a great deal of set-up and preparation work like unloading and loading raw materials using special equipment and human resources. Moreover, in the case of production lines, generally buffers with finite capacities are present between machines. The variation of the production rate between periods brings the necessity of changing the buffers' levels involving a certain handling and logistics cost. Hence, it is clear that the production rate fluctuation between periods may have a significant impact on the set-up cost incurred at the beginning of each period. To the best of our knowledge, this issue has not been considered in the literature. The adopted approach in this study favors the fact of smoothing the production rate between periods in order to avoid the disadvantages of a significant variation (mainly in terms of set-up costs). We also propose an integrated (and not sequential) optimization of the maintenance schedule and the production plan in order to derive a near-optimal integrated maintenance-production strategy taking into account the penalty induced by excessive variation of the production rate between periods of the planning horizon. To do so, we express the set-up cost as being proportional to the production rate variation between successive periods of the production plan. Besides, we propose new modelling features in the total integrated cost modelling, particularly when dealing with the shortage and inventory holding costs by considering them separately contrarily to the previous studies mentioned above. The proposed model allows the investigation of different trade-offs. The first one is between the smoothing penalty and the production plan, the inventory, and the preventive maintenance schedule. We also consider the effect of the demand variability on the smoothing penalty, the production and maintenance plans, and on the total expected cost.

PROBLEM STATEMENT

We consider a single machine M subject to degradation and random failures. It produces one type of product whose demand is random and characterized by a Normal probability distribution with known average and standard deviation. The demand is satisfied from a stock (Fig. 1). Production is planned over a finite horizon divided in H production periods of equal length t. It is assumed that the quantity produced at every period is added to the stock at the end of the period and the demand of every period is satisfied at the end of the period. The machine failure rate increases with both time and the production rate. Periodic preventive maintenance actions have to be scheduled to reduce the probability of failure. We are in presence of a production control problem with a state variable, namely the inventory level, together with the following control variables: the production rates for every period and the number of preventive maintenance actions to be performed over the time horizon Ht. The objective is to minimize the total expected cost over the production planning horizon made of H periods. This cost integrates production costs, inventory and shortage costs, preventive and corrective maintenance costs and an additional cost related to a penalty due to the production rate variation between successive periods. As it will be shown in next section, compared to previous works dealing with the same subject, the first contribution of this study consists in the total integrated cost modelling. On one hand, the shortage and inventory holding costs are formulated separately, and on the other hand, the average number of failures is modelled differently considering an increasing failure rate according to both time and production rate.

The second feature consists in the integrated (and not sequential) optimization of the production plan and the maintenance schedule taking into account the influence of the production rate variation on the system degradation. Finally, the third main contribution is the control of the production rate variation between successive periods in order to smooth the production plan by reducing the penalty induced by excessive variation of the production rate over the planning horizon.

ANALYTICAL STUDY

Used notation

The following parameters are used in the mathematical formulation of the model: ∆t : length of a production period H : number of production periods in the planning horizon H.t : length of the finite planning horizon z(p) : produced quantity during period p (p=1,…, H) z max : maximum quantity that can be produced during period p Z={z(1), ….z(p), ..., z(H)}: Production plan over the horizon H.t u(p) : production rate of the system during period p (p=1,…, H); (u(p)=z(p)/ t) u max : maximum production rate of the system (umax=zmax/ t) U={u(1), ….u(p), ..., u(H)}: Set of production rates adopted for each period d(p) : average demand during period p (p=1,…, H) V : variance of demand for all periods p (p=1,…, H) σ : standard deviation of demand for all periods S(p) : inventory level at the end of period p S 0 : the inventory level at the start of period 1 C pr : unit production cost Cs : inventory holding cost of one product unit per time unit. Cl : shortage cost of one product unit. CTS : expected total cost of production, shortage and inventory holding over the time horizon Ht C v : unit set-up cost CM : total average maintenance cost C pm : preventive maintenance action cost C cm : corrective maintenance action cost SC : penalty cost relative to production rate variation mu : monetary units  : probability to satisfy the demand (service level). The cost model to be developed corresponds to a production/maintenance planning problem with constraints related to inventory, production and preventive maintenance.

Woking Assumptions

The following assumptions are considered: -The machine failure rate increases with both time and production rate -Preventive maintenance actions restore the machine to an "as-good-as-new" state.

-The planning time horizon starts with a new or as good as new machine. However, it doesn't end necessarily with a perfect preventive maintenance action. The equipment could be replaced by a different one or maintained differently after the end of the planning period. -The quantity z(p) produced in period p is added to the inventory at the end of the period and then the demand d(p) is satisfied. -Preventive maintenance and repair actions durations are negligible.

-The units short are lost (no backorders) and a shortage cost is incurred.

Formulation of the integrated production-maintenance total expected cost

We develop in what follows the expressions of the three major components of the total expected cost incurred over the planning horizon. These costs are: the production and inventory cost, the maintenance cost, and the set-up cost (smoothing penalty). Recall that the decision variables are the production rates for every period and the number of preventive maintenance actions to be performed over the planning horizon.

Production, shortage and inventory costs formulation

The constraints related to inventory balance, service level satisfaction, and maximum production quantity per period, are expressed as follows:

Inventory balance: 1 0,.... ) 1 ( ) 1 ( ) ( ) 1 (    H , = p + p d + p z + p S = + p S (1) With: ) 1 ( ) 1 ( ) 1 ( 0 d z + S = S 
Service level satisfaction:

  1 ... 0, 0 ) 1 ( Prob      H , = p p S  (2)
Production capacity:

H = p z p z 1,..., ) ( 0 max    (3)
Accordingly, the inventory evolves as shown in the example below in Fig. 2:

Figure 2: Inventory evolution over the planning horizon Ht Thus, the expression of the total cost of production, shortage and inventory over the finite horizon Ht is given by:

              0 0 1 1 1 1 1 H H H pr S p S p p p p CTS Cs t S p Cl S p C z p                                           (4)
With: 1 Our approach (Eq. 4) makes this important distinction between these unitary costs that are not the same in practice.

The constraint relative to the service level is stochastic. It is possible to transform this constraint to a deterministic one (Hajej et al. [START_REF] Hajej | Optimal integrated maintenance/production policy for randomly failing systems with variable failure rate[END_REF]). In fact, it provides the minimal quantities to produce at each period p in order to satisfy the service level.

      1 ,..., 0 ) 1 ( 0 ) 1 ( Prob ), (          H p z p z p S p S   (5) Where    ), ( p S z
represents a minimum cumulative production quantity expressed as follows:

      1 0,... ) ( ) 1 ( ) ( ) ( 1 1 1 ), (       H , p p S + p d + V = z + p + p p S    (6)
With: V p+1 : variance of demand at period p+1;  p+1 : cumulative Gaussian distribution function with mean d(p+1) and finite variance V p+1 .   -1 p+1 : inverse Gaussian distribution function

The maintenance cost formulation

We consider a periodic preventive maintenance (PM) policy with minimal repair at failure. The interval [0, Ht] is partitioned into N equal parts each of length kt plus a possible remaining duration shorter than kt. Perfect preventive maintenance actions are performed periodically at times i× kt, i=1,…,N. For a given production plan Z (Z={z(1), ..., z(p), ..., z(H)}), and a given PM period kt, the total maintenance cost is given by:

      k U, C + k N C = k U, CM cm pm    (7)
N(k) represents the integer part of (H/k), which corresponds to the number of PM actions. k is the decision variable through which the maintenance policy can be defined. (U,k) stands for the average number of failures (minimal repairs) over the horizon Ht.

The lemma below provides the expression of the average number of failures:

Lemma   ) 1 ( ) ( 1)) ( ( ) ( ) ( ) 1 ( ) ( ) ( )) 1 ( ( ) ( ) ( ) ( 1 ) ( ) ( 2 max ) ( ) 1 ( 0 max ) ( 0 max ) ( 1 1 ) ) ( ( 0 max ) ( 1 1 ) ( 0 2 max ) 1 ( 0 max 0 max 1 > k k N H k k N H = i n k N i Δt n k iN Δt n k N = k k N H Δt n k N k N = j k = i n j i Δt n ij Δt n j Δt Δt λ u u i k + dt t λ u u + dt t λ u u + dt t λ u u + Δt Δt λ u u i k + dt t λ u u + dt t λ u u = k U,                                                                                   (8)
u ij represents the production rate during period i after the j th PM action.

Proof

Recall that PM actions are performed every k production periods, which means every kt restoring the system to as good as new state. The planning horizon Ht is divided into (N(k) +1) periods, each one is noted j, j{0,1,2…. ,(N(k))}. Hence, j=0 corresponds to the first PM interval [0, kt], the second PM period [kt, 2kt] corresponds to period j=1, etc. Note that the last PM action is performed at instant (N(k).k.t), and the interval (N(k)-1).k.t), (N(k).k.t) corresponds to j = (N(k)-1). To illustrate this decomposition, Fig. 3 Average number of failures during a PM period j Recall that the system failure rate evolves according to time and to the production rate which varies from one period to another. In order to take into account the influence of the production rate variation on the system degradation, we use Cox model [START_REF] Cox | Regression models and life-tables[END_REF]. Formally, the failure rate is expressed according to the production rate during each period and the nominal failure rate.

The following expression is proposed (Hajej et al. [START_REF] Hajej | Optimal integrated maintenance/production policy for randomly failing systems with variable failure rate[END_REF]):

          , p n t u p t g u p    
Where λ p (t,u(p)) represents the instantaneous failure rate function at period p as a function of the production rate u(p). λ n (t) is the failure rate for nominal conditions which corresponds to the failure rate when the production system is operating with its maximum production rate, and

      max up g u p u  .
Although this linear relationship between the production rate and the failure rate does not necessarily hold for any type of systems, it is adopted here like in (Hajej et al. [START_REF] Hajej | Optimal integrated maintenance/production policy for randomly failing systems with variable failure rate[END_REF], [START_REF] Hajej | Production/Maintenance Policies Optimization with Operational Age Concept in a Subcontracting Constraint[END_REF], [START_REF] Hajej | Joint optimization of maintenance and production policies with subcontracting and product returns[END_REF]) as an approximation.

Let u ij represent the production rate during period i after the j th PM action.

      ) ( ,..., 0 , ,..., 1 ) , ( k N k j i with u ij 
Hence, the failure rate during period i after the j th PM action is given by:

      ) ( ,..., 0 , ,..., 1 ) , ( ) ( ) , ( max k N k j i u u t u t ij n ij ij     
The evolution of the production system's failure rate within any preventive maintenance cycle period j = h is presented in Fig. 4 In the same way, S 2j , corresponding to the average number of failures during the second production period after the j th PM action, is given by:

    2 1 2 max max 0 t j j j n n u u S t dt t t u u                 
Generalizing this to all production periods after the j th PM actions, we obtain:

            1 max max 0 1 for 2,3... t i j ij ij n n u u S t dt k i t t i k u u                     
Hence, the average number of failures during any period j is expressed as follows:

  ) ( ,..., 0 ) ( )) 1 ( ( ) ( ) ( 2 max ) 1 ( 0 max 0 max 1 k N j t t u u i k dt t u u dt t u u k i n j i t n ij t n j                                   
Average number of failures over the period 0, (N(k).k.t  Considering all periods j (j=0 to N(k)-1), the total average number of failures over the interval 0, (N(k).k.t  is expressed as follows:

                                         1 ) ( 0 2 max ) 1 ( 0 max 0 max 1 ) ( )) 1 ( ( ) ( ) ( k N j k i n j i t n ij t n j t t u u i k dt t u u dt t u u   

Average number of failures over the interval (N(k).k.t, HΔt

Looking at the remaining period N(k).k.t, Ht after the last PM action is performed during the planning horizon, one can distinguish three possible cases:

 Case 1: this period doesn't exist. This means that the last preventive maintenance action coincides with H. t, the end of the planning horizon.

Formally: H-(N(k).k) = 0  Case 2: we have only one production period of length t, in the interval N(k).k.t, Ht.

Formally: H-(N(k).k) = 1  Case 3: we have more than one production period in the interval N(k).k.t, Ht

Formally:

H-(N(k)*k) > 1
Taking into account these three possible situations, the average number of failures in the interval N(k).k.t, Ht is given by:

  ) 1 ( ) ( )) 1 ( ( ) ( ) ( ) 1 ( ) ( 1 ) ( ) ( 2 max ) ( ) 1 ( 0 max ) ( 0 max ) ( 1 1 ) ) ( ( 0 max ) ( 1                                                                k k N H k k N H i n k N i t n k iN t n k N k k N H t n k N t t u u i k dt t u u dt t u u dt t u u    
The first term of this sum corresponds to case 2, and the second and third terms correspond to case 3. Obviously, in the first case (case 1), the number of failures is equal to zero.

Average number of failures over the planning horizon H.t

Summing up the obtained average number of failures for all periods, we obtain: max

) ( ) 1 ( 0 max ) ( 0 max ) ( 1 1 ) ) ( ( 0 max ) ( 1 1 ) ( 0 2 max ) 1 ( 0 max 0 max 1                                                                                                               k k N H k k N H i n k N i t n k iN t n k N k k N H t n k N k N j k i n j i t n ij t n j t t u u i k dt t u u dt t u u dt t u u t t u u i k dt t u u dt t u u k         U

End of the proof

Penalty cost relative to the production rate variation: variable set-up cost

As mentioned earlier, one of the objectives of this study is to smooth the production rate between periods. To do so, we consider a variable set-up cost, which is proportional to the difference of the quantity produced from one period to the next one. Indeed, reducing or increasing the quantity to produce may induce in many situations labour related costs (overtime, or hiring/firing costs) and more or less complexity of the set-up. We define C v as unit set-up cost. It corresponds to the direct cost of resources needed to adapt the production unit to the production quantity change between successive periods. Hence, the penalty cost relative to production rate variation is given by:

    1 1 1 H v p SC C z p z p        (9)

The total integrated maintenance/ production/inventory cost taking into account the set up cost

Using the developed mathematical expressions above in this section, the total cost optimization problem can be formulated as follows:

) ( ) , ( SC CM CTS Min k   Z (10) Subject to: 1 0,..., ) 1 ( 1) ( ) ( ) 1 (    H = p + p d + p z + p S = + p S   1 ... 0, 0 ) 1 ( Prob      H , = p p S  H = p z p z 1,..., ) ( 0 max   

Solution method

Formally, solving the constrained minimization problem stated above (Eq. 10), yields the determination of optimal production quantity for every period {z(1)*, z(2)*, ..., z(p)*, …, z(H)*} and the optimal period (k*t) of PM actions to be performed over the finite horizon Ht. Due to the obvious complexity of the model, we developed an iterative search procedure to find a near-optimal solution. It is presented below in four steps.

 Step 1: Use Eqs. ( 5) and ( 6) to determine the minimum cumulative production quantity to produce in each period in order to satisfy the service level constraint given by Eq. ( 2).

 Step 2: Perform a simple enumeration by varying the production quantity for each period p (using a relatively small step) between the minimum value obtained in step 1 and the maximum one z max given in constraint (3) in order to consider all possible production plans.

 Step 3: For every production plan established in step 2, vary k, number of production periods between preventive maintenance actions, between 1 and H, and calculate each time the corresponding total cost (Eq. 10).

 MATHEMATICA ® software was used to perform the calculations and obtain a near-optimal solution for all instances of the problem considered in this paper.

NUMERICAL EXAMPLE AND TRADE-OFFS ANALYSIS

Reference example

A numerical example is presented in this section in order to illustrate the use of the analytical model developed in the previous sections. We assume that the finite planning horizon is made of 36 periods of one month each (t = 1 month). The other input data are summarized in Table 1 below ('mu' stands for 'monetary units):

Table 1: Numerical data

As for the demand, the variance is supposed to be equal for all periods: V = 4.52. The average demand per period is presented in Solving the problem using the procedure described in section 3, we obtain the following nearoptimal maintenance/production plan presented in Table 3 N*: Optimal number of preventive maintenance actions over Ht k*: Optimal number of production periods between preventive maintenance actions Table 3: Near-optimal maintenance/production plan The total cost related to this plan is equal to : 5.24429×10 7 mu.

Sensitivity analyses

 Effect of C v on the production plan, the inventory, and the preventive maintenance schedule

In the perspective of smoothing the production plan through the production rate control between periods of the planning horizon, we studied the impact of the variation of the unit set-up cost between successive periods, Cv, on the obtained solution.

Considering the same input data, we obtained the following results presented in Table 4: Effect of the variation of the unit set-up cost Fig. 5 below shows the trade-off between the unit-setup cost and the average production quantity variation between successive periods, the optimal number of preventive maintenance actions, and the average end of month inventory. First, one can easily notice that, as expected, when the unit set-up cost between successive periods increases, the average variation of the production quantity, between successive periods, decreases. This clearly illustrates the effect of production smoothing.

Furthermore, it is interesting to notice that the decrease of the average variation of the production quantity between successive periods yields an increase of the quantity to be produced in some periods, and consequently a higher average end of month inventory. Indeed, we know that the quantity to be produced in each period should be between a certain minimum given by the service level constraint (Eqs. 5 and 6) and a given maximum corresponding to the capacity of the machine (Eq. 3). These limits are considered in the first two steps of the solving procedure. Hence, as the smoothing penalty increases, reducing the average variation of the production quantity between successive periods can be achieved either by increasing or by decreasing successive production quantities between the two limits. This is achieved using the model considering simultaneously the trade-off between inventory and shortage costs. The near-optimal solutions, for this specific example, display an increase of the quantities to be produced and therefore an increase of the average end of month inventory. In the same time, the increase of the production rate through the production periods affects the system degradation and induces the need for a more frequent preventive maintenance as it can also be seen in Fig. 5.

The above observations illustrate the fact that the developed model allows the simultaneous consideration of the different trade-offs between the smoothing penalty, the production plan, the inventory, and the preventive maintenance schedule.

We performed several other numerical experiments with different input parameters settings to check the smoothing effect. We present below a sample of the obtained results corresponding to a situation involving a higher shortage cost. While keeping the original combination of input parameters, we increased the shortage cost Cl from 1 mu to 5 mu. A priori, for such a situation, the quantities to be produced should be higher throughout the production periods since shortages are more costly while production and inventory holding costs remain the same.

The following near-optimal maintenance/production plan has been obtained ( The total cost related to this near-optimal plan is equal to : 7.38×10 8 mu.

Compared to the situation with a lower shortage cost (Cl = 1 mu), one can easily observe comparing Tables 3 and5 that, in total, more products are produced in the second situation (a total of 15782, compared to a total of 14186 items). Moreover, this increase of total production over the same horizon yielded an increase of the number of preventive maintenance actions that have to be performed (from 6 to 7 actions). This is due to the fact that the machine failure rate increases with the production rate.

Again, in this situation, we studied the impact of the smoothing penalty on the optimal strategy and on the inventory.

The following results were obtained ( One can easily notice the same trade-offs as in the case of (Cl=1 mu) discussed above. The increase of Cv induces a decrease of the average production quantity variation between successive periods, the need for more preventive maintenance (higher N*) and the increase of inventory.

 Effect of demand variability on the production and maintenance plans, and on the total expected cost While keeping the original set of input parameters (Table 1) as well as the monthly average demand (Table 2), we ran experiments varying the demand variance. The following nearoptimal maintenance/production plans have been obtained (Tables 7, 8 and9 One can clearly see that a demand with a higher variability causes the system to produce more and generate higher levels of inventory. Hence, because of higher production rates yielding a faster degradation of the machine, more frequent preventive maintenance is needed (N* increases). We also note that higher variability of demand yields an increase of the total expected cost.

In fact, as the demand variance increases, obviously, the minimum quantity determined in

Step 1 increases in order to satisfy the required service level. The system is then forced to produce more per period without exceeding the machine capacity. Here comes the trade-off between the smoothing penalty, the inventory, and the maintenance policy. Going further with a relatively high variance (V= 20), the obtained near-optimal integrated policy consists in producing at the maximum production rate for all periods. In this situation, smoothing is perfectly achieved (set-up cost, SC, is null).

CONCLUSION

In this paper we considered a manufacturing system made of a single machine that is subject to random failures. In order to reduce the failure frequency, a periodic preventive maintenance policy with minimal repair at failure is adopted. The production system must satisfy a random product demand over a finite planning horizon with a given required service level. The finite horizon is divided into equal periods. A mathematical model has been developed in order to derive a joint near-optimal production plan and preventive maintenance schedule minimizing the total cost of production, maintenance and inventory, taking into account constraints related to the production system capacity, the service level to be satisfied and the inventory balance.

Compared to previous works, a different cost modelling approach has been proposed. Furthermore, an integrated (not sequential) optimization of the production plan and the maintenance schedule has been proposed. It takes into account the influence of the production rate variation on the production system failure rate. Finally, expressing the set-up cost as being proportional to the production rate variation between successive periods, allowed reducing the penalty induced by excessive variation of the production rate over the planning horizon.

The proposed modelling approach has been illustrated through a numerical example. A sensitivity study with respect to the set-up cost has been performed in order to show the impact of the production plan smoothing on the obtained solution.

The trade-off between the smoothing penalty and the production plan, the inventory, and the preventive maintenance schedule has also been discussed. We finally studied the effect of the demand variability on the smoothing penalty, the production and maintenance plans, and on the total expected cost.

Extensions to this work are under consideration. One of them consists in the investigation of more realistic situations in which preventive maintenance actions are imperfect and have nonnegligible durations. Another extension of this work would consist in going beyond the limitation imposed by the linear relationship between the production rate and the failure rate of the system. More general expressions should be investigated empirically and analytically.

Figure 1 :

 1 Figure 1: Problem statement

  f(t) : probability density function associated with the time to failure of M F(t) : probability distribution function associated with the time to failure of M R(t) : reliability function: 1-F(t)  n (t) : nominal failure rate corresponding to the maximum production rate (depends only on time) λ(t,u(p)) : machine failure rate function during period p (p= 1,…, H) (depends on time and on the production rate during period p) tu : time units

Figure 3 :

 3 Figure 3: Example of PM periods over a planning horizon Ht =13 time units

Figure 4 :

 4 Figure 4: Evolution of the production system failure rate within a preventive maintenance cycle (period j = h)

Step 4 :

 4 Determine the values of the decision variables ((z(1)*, z(2)*, ... z(p)*, …, z(H)*), k*) which yielded the minimal total cost.

F

  (t): Weibull distribution with : Scale parameter : β = 16 & Shape parameter : = 3 Corrective maintenance cost C cm = 3000 mu Preventive maintenance cost C pm = 500 mu Inventory data Shortage cost per unit Cl = 1 mu/unit short Inventory holding cost per unit per time unit Cs = 0.65 mu/unit/tu C pr = 7 mu/unit Maximum production rate Maximum production quantity per period u max = 500 units per month z max = 500 units Nominal set-up cost C v = 5 mu

Figure 5 :

 5 Figure 5: Evolution of the average production quantity variation between successive periods, and the average end of month inventory according to the unit set-up cost C v

  below:

	λ(t,u)
	PM		PM
		(u ih /u max )λ n (t)
			S ih
	S 1h	S 2h	t	t
	hkt	J=h	(h+1)kt

Table 2

 2

	below:

Table 2 :

 2 Monthly mean demand

Table 4 :

 4

	C v : Unit set-up cost between successive periods	Average production quantity variation between successive periods	N*	k*	Average end of month inventory
	2	143.71	4	9	156.94
	5	126.17	6	6	464.56
	7	75.31	7	5	955.58
	8	74.57	7	5	1200.53
	10	58.34	9	4	1221.47

Table 5

 5

	):

Table 5 :

 5 Near-optimal maintenance/production plan with Cl = 5 mu and Cv = 5 mu

Table 6

 6

	):

Table 6 :

 6 Effect of the variation of the unit-setup cost with Cl = 5 mu

Table 7 :

 7 Near-optimal maintenance/production plan with demand variance V=8

	):

Table 8 :

 8 Near-optimal maintenance/production plan with demand variance V=10

	Near-optimal production plan over the finite horizon HΔt = 36 months	N*	k*
	Jan Feb Mar Ap May Jun Jul Aug Sep Oct Nov Dec
	500 500 500 500 500 500 500 500 500 500 500 500	12	3
	500 500 500 500 500 500 500 500 500 500 500 500
	500 500 500 500 500 500 500 500 500 500 500 500

Table 9 :

 9 Near-optimal maintenance/production plan with demand variance V =20We summarize the obtained results V in table 10 below:

	V: Demand Variance	N*	k*	Average production quantity variation between successive periods	Average end of month inventory	Total Cost (mu)
	4.52	6	6	126.17	464.56	5.25×10 7
	8	7	5	90.8	987.33	5.77×10 7
	10	7	5	27.25	1741.83	6.07×10 7
	20	12	3	0	1980.25	6.57×10 7

Table 10 :

 10 Effect of the variation of the demand


J Mathiaud
email: julien.mathiaud@cea.fr

L Mieussens
email: luc.mieussens@math.u-bordeaux.fr

A Fokker-Planck model of the Boltzmann equation with correct Prandtl number for polyatomic gases

Keywords: Fokker-Planck model, Ellipsoidal-Statistical model, polyatomic gas, H-theorem, Prandtl number

We propose an extension of the Fokker-Planck model of the Boltzmann equation to get a correct Prandtl number in the Compressible Navier-Stokes asymptotics for polyatomic gases. This is obtained by replacing the diffusion coefficient (which is the equilibrium temperature) by a non diagonal temperature tensor, like the Ellipsoidal-Statistical model is obtained from the Bathnagar-Gross-Krook model of the Boltzmann equation, and by adding a diffusion term for the internal energy. Our model is proved to satisfy the properties of conservation and a H-theorem. A Chapman-Enskog analysis shows how to compute the transport coefficients of our model. Some numerical tests are performed to illustrate that a correct Prandtl number can be obtained.

Introduction

Numerical simulations of rarefied gas flows are fundamental tools to study the behavior of a gas in a system which length is of same order of magnitude as the mean free path of the gas molecules. For instance, these simulations are used in aerodynamics to estimate the heat flux at the wall of a re-entry space vehicle at high altitudes; numerical simulations are also used to estimate the attenuation of a micro-accelerometer by the surrounding gas in a micro-electro-mechanical system; a last example is when one wants to estimate the pumping speed or compression rate of a turbomolecular pump.

Depending on the problem, the flow can be in non-collisional, rarefied, or transitional regime, and several simulation methods exist. The most common method is the Direct Simulation Monte Carlo method (DSMC) proposed by Bird [3]: it is a stochastic method that simulate the flow with macro particles that mimic the real particles with transport and collisions. It works well for rarefied flows, but can be less efficient when the flow is close to continuous regime (even if there are other version of DSMC that are designed to work well in this regime: see for instance [START_REF] Radtke | Low-noise monte carlo simulation of the variable hard-sphere gas[END_REF][START_REF] Dimarco | Numerical methods for kinetic equations[END_REF]). For non-collisional regimes, the Test Particle Monte Carlo (TPMC) method is very efficient and is a common tool for vacuum pumps [START_REF] Davis | Monte carlo calculation of molecular flow rates through a cylindrical elbow and pipes of other shapes[END_REF]. For transition regimes, there exist mature deterministic solvers 1 based on a direct approximation of the Boltzmann equation that can be used [START_REF] Titarev | Efficient deterministic modelling of three-dimensional rarefied gas flows[END_REF][START_REF] Baranger | Locally refined discrete velocity grids for stationary rarefied flow simulations[END_REF][START_REF] Chigullapalli | Unsteady 3d rarefied flow solver based on boltzmann-esbgk model kinetic equations[END_REF][START_REF] Chen | A unified gas kinetic scheme with moving mesh and velocity space adaptation[END_REF][START_REF] Kolobov | Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement[END_REF][START_REF] Wu | Oscillatory rarefied gas flow inside rectangular cavities[END_REF]. We refer to [START_REF] Mieussens | A survey of deterministic solvers for rarefied flows[END_REF] or further references on these solvers. Among these solvers, some of them are based on model kinetic equations (BGK, ES-BGK, Shakhov): while these models keep the accuracy of the Boltzmann equation in the transition regime, they are much simpler to be discretized and give very efficient solvers [START_REF] Li | Gas-kinetic numerical studies of three-dimensional complex flows on spacecraft re-entry[END_REF][START_REF] Chigullapalli | Unsteady 3d rarefied flow solver based on boltzmann-esbgk model kinetic equations[END_REF][START_REF] Titarev | Efficient deterministic modelling of three-dimensional rarefied gas flows[END_REF][START_REF] Chen | A unified gas kinetic scheme with moving mesh and velocity space adaptation[END_REF][START_REF] Baranger | Locally refined discrete velocity grids for stationary rarefied flow simulations[END_REF]. These models simply replace the Boltzmann collision operator by a relaxation operator toward some equilibrium distribution, but still retain the elementary properties of the Boltzmann operator (collisional invariants, Maxwellian equilibrium). Two of them (BGK and ES-BGK) also satisfy the important entropy dissipation property (the so-called "H-theorem") [START_REF] Andriès | The Gaussian-BGK model of Boltzmann equation with small Prandtl number[END_REF].

Very recently, in a series of paper [START_REF] Jenny | A solution algorithm for the fluid dynamic equations based on a stochastic model for molecular motion[END_REF][START_REF] Gorji | Fokker-Planck model for computational studies of monatomic rarefied gas flows[END_REF][START_REF] Gorji | A Kinetic Model for Gas Mixtures Based on a Fokker-Planck Equation[END_REF][START_REF] Gorji | A Fokker-Planck based kinetic model for diatomic rarefied gas flows[END_REF], Jenny et al. proposed a very different and innovative approach: they proposed to use a different model equation known as the Fokker-Planck model to design a rarefied flow solver. In this model, the collisions are taken into account by a diffusion process in the velocity space. Like the model equations mentioned above, this equation also satisfy the main properties of the Boltzmann equation, even the H-theorem. Instead of using a direct discretization of this equation, the authors used the equivalent stochastic interpretation of this equation (the Langevin equations for the position and velocity of particles) that are discretized by a standard stochastic ordinary differential equation numerical scheme. This approach turned out to be very efficient, in particular in the transition regime, since it is shown to be insensitive to the number of simulated particles, as opposed to the standard DSMC. However, this Fokker-Planck model is parametrized by a single parameter, like the BGK model, and hence cannot give the correct transport coefficients in near equilibrium regimes: this is often stated by showing that the Prandtl number-which is the ratio between the viscosity to the heat transfer coefficient-has an incorrect value. The authors have proposed a modified model that allow to fit the correct value of the Prandtl number [START_REF] Gorji | Fokker-Planck model for computational studies of monatomic rarefied gas flows[END_REF], and then have extended it to more complex flows (multi-species [START_REF] Gorji | A Kinetic Model for Gas Mixtures Based on a Fokker-Planck Equation[END_REF] and diatomic [START_REF] Gorji | A Fokker-Planck based kinetic model for diatomic rarefied gas flows[END_REF]). However, even if the results obtained with this approach seem very accurate, it is not clear that these models still satisfy the H-theorem.

In a previous paper [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF], we have proposed an alternative modification of the Fokker-Planck equation for which we were able to prove the H-theorem and to obtain the correct Prandtl number. Our approach follows more or less the ideas used to transform the BGK model into the ES-BGK model, and our new model is called the ES-Fokker-Planck model (ES-FP), for monoatomic gases. In the present paper, we continue our work and extend our ES-FP model to polyatomic gases. To do so, we follow the approach used in [START_REF] Andriès | The Gaussian-BGK model of Boltzmann equation with small Prandtl number[END_REF] to derive a ES-BGK model for polyatomic gases. Again, we are able to prove that this model satisfies the H-theorem and that it gives the correct transport coefficients in its hydrodynamic limit. We also show numerical experiments that confirm our analysis for some space homogenous problems.

The outline of this paper is the following. In section 2, our model and its properties are proved, including the H-theorem. A Chapman-Enskog analysis is made in section 3 to compute the transport coefficients obtained with this model. In section 4, we show that our model can include a variable parameter that allow to fit the correct Prandtl number and the viscosity parameter for polyatomic gases. Our numerical tests are shown in section 5, and conclusions and perspectives are presented in section 6.

2 The Fokker-Planck model for polyatomic gases

The Boltzmann equation for polyatomic gases

We consider a gas described by the mass density of particles f (t, x, v, I) that at time t have the position x, the velocity v and an internal energy parameter I (note that both position x and velocity v are vectors and that I is a scalar). The internal energy of a particle is equal to ε(I) = I 2 δ , δ being linked to the number of degrees of freedom of the polyatomic gas (δ = 2 for diatomic gases). The corresponding macroscopic quantities are (ρ, ρu, E) = (1, v, 1 2 |v| 2 ) + ε(I) f , where ρ, ρu, and E are the mass, momentum, and energy densities, and φ = R 3 φ(v, I) dvdI for any velocity dependent function. The temperature T of the gas is defined by relation E = 1 2 ρ|u| 2 + 3+δ 2 ρRT , where R is the gas constant, and the pressure is p = ρRT . We can also define a translational energy E tr , a translationnal temperature T tr , an internal energy E int and an internal energy temperature T int through:

E tr = 3 2 ρRT tr = 1 2 |v -u| 2 f (1)
E int = δ 2 ρRT int = ε(I)f , (2)
so that:

T = 3 3 + δ T tr + δ 3 + δ T int (3)
The evolution of the gas can be governed by the following Boltzmann equation (if we do not take into account exchanges in internal energy)

∂ t f + v • ∇ x f = Q(f, f ), (4)
where Q(f, f ) is for instance the Borgnakke-Larsen kernel of collision described in ( [START_REF] Desvillettes | Sur un modèle de type borgnakke-larsen conduisant à des lois d'énergie non linéaires en température pour les gaz parfaits polyatomiques[END_REF][START_REF] Desvillettes | A kinetic model allowing to obtain the energy law of polytropic gases in the presence of chemical reactions[END_REF]) which allows reversible exchanges of energy between internal energy and translationnal energy of the particles.

It is well known that this operator conserves the mass, momentum, and energy, and that the local entropy H(f ) = f log f is locally non-increasing. This means that the effect of this operator is to make the distribution f relax towards its own local Maxwellian distribution, which is defined by

M p (f )(v, I) = ρΛ δ (2π) 3/2 (RT ) (3+δ)/2 exp - |v -u| 2 2RT - I 2/δ RT , (5)
where

Λ δ = +∞ 0 exp(-I 2/δ )dI -1
. For the sequel, it is useful to define another macroscopic quantity, which is not conserved: the "temperature" tensor, defined by

Θ := 1 ρ (v -u) ⊗ (v -u)f . (6)
In an equilibrium state (that is to say when f = M (f )), Θ reduces to the isotropic tensor RT I.

2.2

The ES-Fokker-Planck model for polyatomic gases

Some facts about the ES-Fokker Planck model for monoatomic gases

In [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF] we extended the standard Fokker-Planck model for the Boltzmann equation by modifying its collision kernel D(f ) through:

∂ t f + v • ∇ x f = D(f ) = 1 τ ∇ v • (v -u)f + Π∇ v f , (7)
Π = (1 -ν)RT Id + νΘ.
with τ a characteristic relaxation time and ν a parameter allowing to recover the correct Prandtl number for monoatomic gases when ν = -1 2 . Now we remind the results proved in [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF].

Proposition 2.1. The tensor Π is symmetric positive definite for every tensor Θ if, and only if,

- RT λ max -RT < ν < RT RT -λ min , (8)
where λ max and λ min are the (positive) maximum and minimum eigenvalues of Θ. Moreover Π is positive definite independently of Θ as long as: -

1 2 < ν < 1.
Proposition 2.2. The operator D satisfies:

D(f ) = 1 τ ∇ v • ΠG(f )∇ v f G(f ) = 1 τ ∇ v • Πf ∇ v log f G(f ) , (9)
where G(f ) is the anisotropic Gaussian defined by

G(f ) = ρ det(2πΠ) exp - (v -u)Π -1 (v -u) 2 , ( 10
)
which has the same 5 first moments as f and Π for temperature tensor:

(1, v, 1 2 |v| 2 )G(f ) = (ρ, ρu, E), (v -u) ⊗ (v -u)G(f ) = ρΠ.
If ν satisfies (8) and ν < 1, it conserves mass, momentum, energy and entropy is dissipated:

(1, v, 1 2 |v| 2 )D(f ) = 0, D(f ) log f ≤ 0, ( 11
)
Moreover as expected the equilibrium property holds:

D(f ) = 0 ⇔ f = G(f ) ⇔ f = M (f ).

The ES-Fokker-Planck model for polyatomic gases

Our extension of the previous model to polyatomic gases is based on two ideas. First, the variation of the internal energy of the particles due to collisions are taken into account by an additional diffusion term (with respect to the internal energy variable) in the collision operator. Then, the transfers between translational and internal energies are taken into account by an additional parameter θ, like in the construction of the ESBGK model for polyatomic gases proposed in [START_REF] Andriès | The Gaussian-BGK model of Boltzmann equation with small Prandtl number[END_REF]. This leads to our ES-Fokker-Planck model for polyatomic gases:

∂ t f + v • ∇ x f = D(f ), (12)
where the collision operator is defined by

D(f ) = 1 τ ∇ v • (v -u)f + Π∇ v f + ∂ I (δf I + δ 2 2 RT rel I 2-2 δ ∂ I f ) , ( 13
)
where τ is the relaxation time, and Π is the following combination between the temperature tensor Θ, the translational temperature tensor RT tr Id and its equilibrium value RT Id:

Π = (1 -θ) ((1 -ν)RT tr Id + νΘ) + θRT Id, (14)
where the coefficients ν and θ are some free parameters, and the temperature T rel is defined by

T rel = (1 -θ)T int + θT. ( 15
)
This model is the natural extension of the monoatomic ESFP model to get the polyatomic anisotropic Gaussian as equilibrium for polyatomic gases, as it was done for the ESBGK model for polyatyomic gases (see [START_REF] Andriès | The Gaussian-BGK model of Boltzmann equation with small Prandtl number[END_REF]). Now, we give a few comments about parameters θ and δ. First, it can be seen that θ governs the reversible exchanges between internal and translational energies. Indeed, for a space homogeneous problem (we neglect the space dependence of the problem), it can easily be seen that the rate at which the energy is transferred from internal energy to translational energy is 3+δ

2δθ τ E int . The corresponding relaxation time then is 3+δ 2δθ τ , and a relaxation collision number on the collision rate time-scale τ (see [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]) can be defined by Z = 3+δ 2δθ which is clearly parametrized by θ. Consequently, it is clear that when θ = 0 we no longer have exchanges between internal and translational energies. Note that this model is compatible with the monoatomic model: indeed, when δ = 0 we immediatly recover the monoatomic model since in this case T tr = T .

For the following, it is usefull to rewrite model ( 12) a bit differently. We define the generalized tensor Ω by

Ω = Π 0 0 δ 2 2 RT rel I 2-2 δ , (16)
so that the collision operator of (12) reads

D(f ) = 1 τ ∇ v,I • v -u δI f + Ω∇ v,I f ,
where ∇ v,I denotes the derivative operator with respect to both v and I variables. We can now define the generalized Gaussian G p (f ) for polyatomic gases as:

G p (f ) = ρΛ δ det(2πΠ)(RT rel ) δ/2 exp - 1 2 v -u δI T Ω -1 v -u δI . (17)
One can note that contrary to the case of monoatomic gases the tensor Ω depends on the variables of the phase space through the internal parameter. The generalized Maxwellian M p (f ) of equilibrium for polyatomic gases (5) can also be defined as:

M p (f ) = ρΛ δ (2π) 3/2 (RT ) 3+δ 2 exp   - v-u √ 2RT I 1/δ √ RT T v-u √ 2RT I 1/δ √ RT   (18)
We now look for the properties of the collision kernel D(f ) as it was done in the monoatomic case.

Proposition 2.3. The operator D has two other equivalent formulations:

D(f ) = 1 τ ∇ v,I • ΩG p (f )∇ v,I f G p (f ) , (19)
and

D(f ) = 1 τ ∇ v,I • Ωf ∇ v,I log f G p (f ) , (20)
Moreover G p has the same 5 first moments as f

(1, v, 1 2 |v| 2 + I 2/δ )G p (f ) = (ρ, ρu, ρE),
and has its translationnal temperature tensor equal to

(v -u) ⊗ (v -u)G p (f ) /ρ = Π.
Proof. Relations ( 19) and ( 20) are simple consequences of the relations

∇ v G p (f ) = -Π -1 (v - u)G p (f ) and ∂ I G p (f ) = -2I 2/δ-1 δRT rel G p (f ) .
The moments of G p (f ) are given by standard Gaussian integrals.

Proposition 2.4 (Condition of definite positiveness of Π).

Assume 0 < θ < 1, the tensor Π is symmetric positive definite if, and only if,

- RT tr + θ 1-θ RT λ max -RT tr < ν < RT tr + θ 1-θ RT RT tr -λ min , (21)
for every tensor Θ ( λ max and λ min being its (positive) maximum and minimum eigenvalues). Moreover Π is positive definite independently of the eigenvalues of Θ as long as:

- 1 2 < ν < 1, (22)
Proof. Since Θ is symmetric and positive definite, it is diagonalizable and its eigenvalues are positive. Moreover, because of the definition of Π, both Θ and Π are diagonalizable in the same orthogonal basis. Then there exists an invertible orthogonal matrix P such that

P ΘP T =   λ 1 0 0 0 λ 2 0 0 0 λ 3   ,
and

P ΠP T =   γ 1 0 0 0 γ 2 0 0 0 γ 3   .
with

γ i = (1 -θ) ((1 -ν)RT tr + νλ i ) + θRT .
Clearly Π is positive definite if, and only if, γ i > 0 for every i, which leads to:

- RT tr + θ 1-θ RT λ max -RT tr < ν < RT tr + θ 1-θ RT RT tr -λ min ,
which is [START_REF] Wu | Oscillatory rarefied gas flow inside rectangular cavities[END_REF]. Note that RT tr = trace(Θ)/3 is the mean value of the eigenvalues λ i of Θ. Consequently, when λ max gets close to the RT tr , then λ min gets close to RT tr too-which happends when f is close to a Maxwellian equilibrium-and hence our analysis shows that all the γ i are positive, indendently of ν. This is consistent with ( 21) wich shows that ν can take any value.

The sufficient condition( 22) is obtained as follows. At worse λ min is equal to zero and θ is equal to 0 so that ν < RT tr RT tr -0 = 1.

At worse λ max is equal to trace(Θ) = 3RT tr and θ is equal to 0 so that

ν > - RT tr 3RT tr -RT tr = - 1 2 .
Now, we state that D has the same conservation and entropy properties as the Boltzmann collision operator Q. Proposition 2.5. We assume ν satisfies [START_REF] Wu | Oscillatory rarefied gas flow inside rectangular cavities[END_REF] and that ν < 1. The operator D conserves the mass, momentum, and energy:

(1, v, 1 2 |v| 2 + I 2/δ )D(f ) = 0, (23)
it satisfies the dissipation of the entropy:

D(f ) log f ≤ 0,
and we have the equilibrium property:

D(f ) = 0 ⇔ f = G p (f ) ⇔ f = M p (f ).
Proof of the conservation property. We take any function φ(v) and by integration by parts, we get

D(f )φ = - 1 τ ((v -u)f + Π∇ v f ) • ∇ v φ -(δf I + δ 2 2 RT rel I 2-2 δ ∂ I f )∂ I φ . With φ(v) = 1, we directly get D(f ) = 0. With φ(v) = v i , we get D(f )v i = - 1 τ (v i -u i )f + Π i,j ∂ v j f = 0, since v i f = ρu i = ρ v i f and ∂ v j f = 0. For φ(v, I) = 1 2 |v| 2 + I 2/δ , we get D(f ) 1 2 |v| 2 + I 2/δ = - 1 τ ((v i -u i )f + Π i,j ∂ v j f )v i -(δf I + δ 2 2 RT rel I 2-2 δ ∂ I f ) 2 δ I 2 δ -1 = - 1 τ (3ρRT tr -ρtrace(Π)) - 1 τ (ρδRT int -ρδRT rel ) .
We note that ( 14) implies trace(Π) = (1 -θ) (3RT tr ) + θ(3RT ). Then we use the definitions ( 3) and ( 15) of T and T rel to obtain that the last integral is zero.

Proof of the entropy property. This proof, which rather long, is based on two main ideas. First, as noticed in the proof of proposition 2.4, Θ and Π are symmetric and can be diagonalized in the same orthogonal basis. Like it is done in [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF], we will now suppose that we work in this basis so that:

Θ =   λ 1 0 0 0 λ 2 0 0 0 λ 3   , and Π =   γ 1 0 0 0 γ 2 0 0 0 γ 3   .
with

γ i = (1 -θ) ((1 -ν)RT tr + νλ i ) + θRT .
Moreover, we will make use of the well known Jensen inequality: for every convex function x → C(x) and every mean µ, we have C(µ(a 1 , a 2 , . . . , a n )) ≤ µ(C(a 1 ), C(a 2 ), . . . , C(a n )), for every real a 1 , a 2 , . . . , a n , with equality if, and only if, all the a i are equal or C is linear. This inequality will be used twice. First, since 3 i=1 λ i = trace(Θ) = 3RT tr , we have RT tr = 1 3 3 i=1 λ i = µ(λ 1 , λ 2 , λ 3 ), where µ is the arithmetic mean of the λ i . Consequently, since the function C(x) = 1/x is convex, the Jensen inequality induces C(T T r ) ≤ µ(C(λ 1 ), C(λ 2 ), C(λ 3 )), which reads

1 RT tr ≤ 1 3 3 i=1 1 λ i . (24)
Second, the definition (3) of the equilibrium temperature reads T = µ(T tr , T int ), where µ(a, b) = (3a + δb)/(3 + δ) is a weighted average. Consequently, the Jensen inequality implies

C(T ) ≤ µ(C(T tr ), C(T int )), which reads 1 T ≤ 1 3 + δ ( 3 T tr + δ T int ). ( 25
)
Now we compute D(f ) log f with an integration by parts to get

D(f ) log(f ) = - 1 τ 3 i=1 ((v i -u i ) f + γ i ∂ v i f ) ∂ v i f f -δf I + δ 2 2 RT rel I 2-2 δ ∂ I f ∂ I f f = - 1 τ 3 i=1 γ i λ i ((v i -u i ) f + λ i ∂ v i f ) ∂ v i f f + 3 i=1 1 - γ i λ i ∂ v i f (v i -u i ) - 1 τ δf I + δ 2 2 RT rel I 2-2 δ ∂ I f ∂ I f f = - 1 τ 3 i=1 γ i λ i ((v i -u i ) f + λ i ∂ v i f ) ∂ v i f f + (v i -u i ) λ i - (v i -u i ) λ i + 1 τ 3 i=1 1 - γ i λ i f - 1 τ δf I + δ 2 2 RT rel I 2-2 δ ∂ I f ∂ I f f = - 1 τ 3 i=1 γ i λ i 1 f λ i ((v i -u i )f + λ i ∂ v i f ) 2 + 1 τ 3 i=1 γ i λ i ((v i -u i ) f + λ i ∂ v i f ) (v i -u i ) λ i + 1 τ 3 i=1 1 - γ i λ i f - 1 τ δf I + δ 2 2 RT rel I 2-2 δ ∂ I f ∂ I f f = (A) + (B) + (C) + (D)
The first term (A) is clearly non positive. By using another integration by parts and the definition of the λ i , it is easy to see that (B) = 0. For (C), we use the definition of γ i to find

(C) = 1 τ 3 i=1 1 - (1 -θ) ((1 -ν)RT tr + νλ i ) + θRT λ i f .
Then we expand the right-hand side

(C) = 3 τ 1 -(1 -θ)ν -((1 -θ)(1 -ν)RT tr + θRT ) 1 3 i 1 λ i , (26)
and we use inequality (24) to find

(C) ≤ 3θ τ 1 - T T tr f . (27)
For (D), we try again to make appear a non positive term. We find

(D) = - 1 τ δf I + δ 2 2 RT rel I 2-2 δ ∂ I f ∂ I f f = - 1 τ T rel T int δf I + δ 2 2 RT int I 2-2 δ ∂ I f ∂ I f f + 2 δ I 2 δ -1 RT int - 2 δ I 2 δ -1 RT int - 1 τ 1 - T rel T int δ I∂ I f = - 1 τ T rel T int 1 f δ 2 2 RT int I 2-2 δ δf I + δ 2 2 RT int I 2-2 δ ∂ I f 2 + 1 τ T rel T int δf I + δ 2 2 RT int I 2-2 δ ∂ I f 2 δ I 2 δ -1 RT int + 1 τ 1 - T rel T int δ f = (D1) + (D2) + (D3) (28)
Then (D1) is clearly non positive. For (D2), an integration by parts and the definition (2) of T int implies (D2) = 0. This shows that (D) ≤ (D3) = 1 τ 1 -T rel T int δ f . Finally, we collect all the inequalities obtained for (A), (B), (C), and (D) to find

D(f ) log(f ) ≤ 3θ τ 1 - T T tr f + θ τ δ 1 - T T int f = (3 + δ)θ τ f 1 -T 1 3 + δ ( 3 T tr + δ T int ) ≤ 0, (29)
from the Jensen inequality (25).

Proof of the equilibrium property. First, we assume that f = M p (f ), as defined by [START_REF] Chigullapalli | Unsteady 3d rarefied flow solver based on boltzmann-esbgk model kinetic equations[END_REF]. Then it can easily be seen that ( 1) and (2) imply T tr = T int = T , while (6) implies Θ = RT Id. Using ( 14) and ( 15) and the previous relations yields T rel = T and Π = RT Id. This shows that G p (f ), as defined by [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF], satisfies G p (f ) = M p (f ), and hence f = G p (f ). Finally, the formulation (19) of the collision operator implies D(f ) = 0. Now, we assume that D(f ) = 0, and we want to prove that f = M p (f ). This is where the entropy inequality is used: indeed, D(f ) = 0 implies D(f ) log f = 0, which is used as follows. First, this relation with inequality (29) imply 1 3+δ ( 3 Ttr + δ T int ) ≤ 1 T . This means that the Jensen inequality (25) is in fact an equality, and hence T tr = T int = T . Now, we remind that D(f ) log f = (A) + ((C) + (D3)) + (D1), where (A), ((C) + (D3)), and (D1) are non positive (see the proof of the entropy inequality). Then D(f ) log f = 0 implies in particular (A) = 0, (C) + (D3) = 0 and (D1) = 0. Note that (A) = 0 and (D1) = 0 imply

∂ v i f f + (v i -u i ) RT = 0 and ∂ I f f + 2 δ I 2 δ -1 RT = 0.
This shows that f = α exp(i

(v i -u i ) 2 2λ i -I 2 δ
RT ). Moreover, (C) + (D3) = 0 simplifies to (C) = 0, since the equality of relative and internal temperatures implies (D3) = 0. This also implies (C) =

3 τ (1 -(1 -θ)ν)(1 -RT 1 3 i 1 λ i ). Consequently, (C) = 0 implies 1 3 i 1 λ i = 1
RT : this means that Jensen inequality (24) is in fact an equality, and hence

λ 1 = λ 2 = λ 3 = RT . Finally, we have obtained f = α exp(-|v-u| 2 2RT -I 2 δ
RT ), and hence f = M p (f ), which concludes the proof.

Equation on internal energy

While the formulation of our model with the internal energy parameter I is useful to prove the H-theorem, using the internal energy ε(I) = I 2/δ turns out to be more convenient for numerical simulations (see section 5.1). If we denote by g(t, x, v, ε) the corresponding distribution function (that is to say, such that g dε = f dI), we get the following equation

∂ t g + v • ∇ x g = D(g), (30)
where the collision operator is defined by

D(g) = 1 τ ∇ v • (v -u)g + Π∇ v g + ∇ ε 2(ε - δ 2 RT rel )g + ∇ ε 2(RT rel εg) . (31)
The equilibrium and its Gaussian extension are now

M(g) = ρΛ g (2π) 3/2 ε δ-2 2 (RT ) δ 2 exp - 1 2 (v -u) 2 + ε RT and G(g) = ρΛ g det(2πΠ) ε δ-2 2 (RT rel ) δ 2 exp - 1 2 (v -u) t Π -1 (v -u) - ε RT rel with Λ g = R ε δ-2 2 exp(-ε)dε -1 .
When one deals with diatomic gases (δ = 2) one recovers the same equation as before since ε(I) is equal to I. The main difference with the ε-formulation of the problem is that the entropy inequality is transformed into a relative entropy inequality, that is

D g (g) log g ε δ-2 2 ≤ 0,
due to the change of variables.

Chapman-Enskog analysis

From ( 12) and ( 23), it can easily be shown that the moments of f satisfy the following conservation laws

∂ t ρ + ∇ x • ρu = 0, ∂ t ρu + ∇ x • (ρu ⊗ u) + ∇ x • Σ(f ) = 0, ∂ t E + ∇ x • (Eu + Σ(f )u + q(f )) = 0, (32)
where Σ(f ) and q(f ) denote the stress tensor and the heat flux, defined by

Σ(f ) = (v -u) ⊗ (v -u)f q(f ) = ( 1 2 |v -u| 2 + I 2/δ )(v -u)f . (33)
The Chapman-Enskog procedure consists in looking for an approximation of Σ(f ) and q(f ) up to second order with respect to the Knudsen number ε, which is defined below, as functions of ρ, u, T and their gradients.

To do so, we now write our model in a non-dimensional form. Assume we have some reference values of length x * , pressure p * , and temperature T * . With these reference values, we can derive reference values for all the other quantities: mass density

ρ * = p * /RT * , velocity v * = √ RT * , time t * = x * /v * , distribution function f * = ρ * /(RT * ) 3/2 .
We also assume we have a reference value for the relaxation time τ * . By using the non-dimensional variables w = w/w * (where w stands for any variables of the problem), our model can be written

∂ t f + v • ∇ x f = 1 ε D(f ), (34)
where ε = v * τ * x * is the Knudsen number. Note that since we always work with the non-dimensional variables from now on, these variables are not written with the ' in (34).

Note that an important consequence of the use of these non-dimensional variables is that RT has to be replaced by T in every expressions given before. Namely, now Π is defined by

Π = (1 -θ)((1 -ν)T tr Id + νΘ) + θT Id, (35)
the temperatures are now defined by

E tr = 3 2 ρT tr , E int = δ 2 ρT int , and E = 1 2 ρ|u| 2 + 3+δ 2 ρT, (36)
and the Maxwellian of f now is

M p (f ) = ρΛ δ (2π) 3/2 (T ) (3+δ)/2 exp - |v -u| 2 2T - I 2/δ T .
Now, it is standard to look for the deviation of f from its own local equilibrium, that is to say to set f = M p (f )(1 + εg). However, this requires the linearization of the collision operator D around M p (f ), which is not very easy. At the contrary, it will be shown that it is much simpler to look for the deviation of f from the anisotropic Gaussian G p (f ) defined in [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF]. Since it can easily be seen that M p (f ) and G p (f ) are close up to O(ε) terms, this expansion is sufficient to get the Navier-Stokes equations. We shall then prove the following proposition: Proposition 3.1. The moments of the solution of the kinetic model (34) satisfy, up to O(ε 2 ), the Navier-Stokes equations

∂ t ρ + ∇ • ρu = 0, ∂ t ρu + ∇ • (ρu ⊗ u) + ∇p = -∇ • σ, ∂ t E + ∇ • (E + p)u = -∇ • q -∇ • (σu), ( 37
)
where the shear stress tensor and the heat flux are given by

σ = -µ ∇u + (∇u) T -α∇ • u , and q = -κ∇ • T, ( 38
)
with the following values of the viscosity and heat transfer coefficients (in dimensional variables)

µ = τ p 2(1 -(1 -θ)ν) , α = (γ -1) - (1 -ν)(1 -θ) θ 5 3 -γ and κ = 5 + δ 6 τ pR, (39)
and γ = δ+5 δ+3 . Moreover, the corresponding Prandtl number is

Pr = 3 2(1 -(1 -θ)ν) . ( 40
)
In the first steps of the proof, we state some approximation properties.

Proposition 3.2. We write f as f = G p (f )(1 + εg). Then we have:

(1, v, 1 2 |v| 2 + I 2 δ )G p (f )g = 0. ( 41
)
Moreover, if we assume that the deviation g is an O(1) with respect to ε, then we have

Π = T Id + O(ε), T tr = T + O(ε), T int = T + O(ε), T rel = T + O(ε). ( 42
)
and

G p (f ) = M p (f ) + O(ε). (43)
Proof. We have already noticed that (1, v,

1 2 |v| 2 )G p (f ) = (1, v, 1 2 |v| 2 + I 2 δ )f (see proposi- tion 2.
3), which gives (41). Relations (42) and ( 43) are obtained by simple Taylor expansions of ( 35) and ( 17) with respect to ε. Now, we show that Σ(f ) and q(f ) can be approximated up to O(ε 2 ) by using the previous decomposition.

Proposition 3.3. We have

Σ(f ) = pId + ε 1 -(1 -θ)ν Σ(M p (f )g) + ρ (1 -ν)(1 -θ)(T tr -T ) 1 -(1 -θ)ν Id + O(ε 2 ), ( 44
)
θρ(T tr -T ) = 1 3 ε trace(Σ(M p (f )g)) + O(ε 2 ) (45) q(f ) = εq(M p (f )g) + O(ε 2 ). ( 46
)
Proof. By using f = G p (f )(1 + εg), we have

Σ(f ) = Σ(G p (f )) + εΣ(G p (f )g) = ρΠ + εΣ(G p (f )g). (47)
But since by definition Σ(f ) = ρΘ(f ), (35) also implies:

ρΠ = (1 -θ)((1 -ν)ρT tr Id + νΣ(f )) + θρT Id.
Using this relation in (47 Consequently, the Navier-Stokes equations can be obtained provided that Σ(M p (f )g) and q(M p (f )g) can be approximated, up to O(ε). This can be done by looking for an approximation of the deviation g itself: it is obtained by using the decomposition f = G p (f )(1 + εg) into (34), which gives

) gives Σ(f ) = (1-θ)((1-ν)ρT tr Id+νΣ(f ))+θρT Id+εΣ(G p (f )g), which, with (43)
∂ t G p (f ) + v • ∇ x G p (f ) + O(ε) = 1 ε D(G p (f )(1 + εg)). First, we state what the expansion of D(G p (f )(1 + εg)) is.
Proposition 3.4. We have

D(G p (f )(1 + εg)) = ε 1 τ M p (f )Lg + O(ε 2 ),
where L is the linear operator

Lg = 1 M p (f ) ∇ v • (T M p (f )∇ v g) + ∂ I δ 2 2 T I 2-2 δ M p (f )∂ I g . (48)
Proof.

Formulation (19) of D(f ) gives D(G p (f )(1 + εg)) = 1 τ ∇ v,I • ΩG p (f )∇ v,I G p (f )(1 + εg) G p (f ) . ( 49
) A direct computation gives ∇ v,I Gp(f )(1+εg) Gp(f )
= ε∇ v,I g. Using this relation and the definition ( 16) of Ω gives:

D(G p (f )(1 + εg)) = ε 1 τ ∇ v • (ΠG p (f )∇ v g) + ∂ I δ 2 2 T rel I 2-2 δ G p (f )∂ I g .
Then the final result is obtained by using relations (42) and (43) to approximate Π, G p , and T rel up to O(ε) by T Id, M p (f ) and T , respectively.

This proposition shows us that g satisfies

M p (f )Lg = τ (∂ t G p (f ) + v • G p (f )) + O(ε)
, and then using (43) gives

Lg = τ M p (f ) (∂ t M p (f ) + v • ∇ x M p (f )) + O(ε). (50)
We now have to solve this equation approximately to obtain an approximation of g. First, the right-hand side of (50) is expanded, and the time derivatives of ρ, u, and T are approximated up to O(ε) by their space gradients by using (32), (44), and (46). Then we get

1 M p (f ) (∂ t M p (f ) + v • ∇ x M p (f )) = A(V, J) • ∇T √ T + B(V, J) : ∇u + O(ε), where V = (v -u)/ √ T , J = I T δ/2 and A(V, J) = |V | 2 2 + J 2/δ - 5 + δ 2 V and B(V, J) = V ⊗ V - |V | 2 + 2J 2/δ 3 + δ Id.
Then (50) reads

Lg = τ A(V, J) • ∇T √ T + B(V, J) : ∇u + O(ε). (51)
This equation can be solved by showing that A and B are eigenvectors of L. Indeed, we have the following proposition.

Proposition 3.5. The components of A and B satisfy LA i = -3A i and LB i,j = -2B i,j .

Proof. By using the change of variables V = (v -u)/ √ T and J = I T δ/2 , the definition (48) reduces to

Lg = 1 M 0 ∇ V • (M 0 (V, J)∇ V g) + 1 M 0 ∂ J M 0 (V, J) δ 2 2 J 2-2 δ ∂ J g + O(ε), (52)
where

M 0 (V, J) = Λ δ (2π) 3 2 exp - |V | 2 2 -J 2/δ
. The direct computation of LA i and LB i,j is easily obtained from (52) and is left to the reader.

This property shows that we can look for an approximate solution of (51) as a linear combination of the components of A and B. We find

g = - τ 3 ∇T √ T • A(V, J) - τ 2 ∇u : B(V, J) + O(ε).
Then we just have to insert this expression into Σ(M p (f )g) and q(M p (f )g) to get, after calculation of Gaussian integrals,

Σ(M p (f )g) = - τ p 2 (∇u + (∇u) T -(γ -1)(∇ • u)Id) + O(ε), q(M p (f )g) = - τ p(5 + δ) 6 ∇T + O(ε).
Consequently, we take ν in the larger interval [-

RTtr+ θ 1-θ RT λmax-RTtr , 1]
. This means that ν is a non constant value, since T , T tr , and λ max depend on time and space. This allows us to hope for a Prandtl number smaller than 1, while the proposition 2.5 ensures that the model is well defined. However, it is not obvious that this interval is large enough to adjust any small value of the Prandtl number. Our idea is that the Prandtl number makes sense when the flow is close to the equilibrium, that is to say when f is close to its own local Maxwellian M p (f ). In such case, Θ = RT I + O(ε), and all the eigenvalues of Θ are close to RT up to O(ε). This implies λ max = RT (1 + Cε) and T tr = T (1 + Dε), and since by definition λ max > RT tr , then C > D. Consequently, the smallest possible value of ν behaves like -1 (1-θ)(C-D)ε which can be as small as possible for small ε. This shows that ν can take any arbitrary value between -∞ and 1 when ε is small enough, and hence the Prandtl number can be made arbitrarily small.

Case of a diatomic gas: air

For air, the Prandtl number is estimated to Pr ≈ 0.71 . Since θ = 0.2 is generally used in literature [START_REF] Andriès | The Gaussian-BGK model of Boltzmann equation with small Prandtl number[END_REF] for air, one gets that we need to use ν air ≈ -1.3908.

In other words, by defining ν as a non constant value that satisfies [START_REF] Wu | Oscillatory rarefied gas flow inside rectangular cavities[END_REF] and is lower than 1, we can adjust the correct Prandtl number provided that ν can be set to ν air near the equilibrium regime.

This analysis suggests a very simple definition of ν : we propose to use the smallest negative ν such that:

• Π remains strictly definite positive,

• ν ≥ ν air .
This leads to the following definition:

ν = max ν air , - RT tr + θ 1-θ RT λ max -RT tr . ( 54
)
Note that in most realistic cases ν is equal to ν air . Indeed, ν = ν air implies λ max > 1.82RT tr for θ = 0.2 because T = 3 5 T tr + 2 5 T int ≥ 3 5 T tr , which can only happen in case of a highly non equilibrium flow with strong directional non isotropy: such a case is very specific and is usually not observed in aerodynamical flows, for instance.

Numerical tests

Due to the diffusion operator used in our model, we find it attractive to use a stochastic simulation method (like in [START_REF] Jenny | A solution algorithm for the fluid dynamic equations based on a stochastic model for molecular motion[END_REF], see also [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF]). In this paper, our study is restricted to space homogeneous problems, since this is sufficient to illustrate the main property of our model. Further studies will be done in inhomogeneous cases and comparisons with full Boltzmann equation and/or ESBGK model will be provided in a forthcoming paper. In this section, we use our model written in energy variable (see section 2.3):

∂ t g = 1 τ ∇ v • (v -u)g + Π∇ v g + ∂ ε 2(ε - δ 2 RT rel )g + ∂ ε 2(RT rel εg) . (55)
In the test cases presented here, all numerical weights ω i are equal, and we just have to define the dynamics of the numerical particles. To solve diffusive problems, it is well-known that using Brownian motion is a good way to proceed (see [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF]): the corresponding stochastic ordinary differential equation is called the Ornstein-Uhlenbeck process that reads

dV i (t) = - dt τ (V i (t) -u) + A v dB v (t) (62) dε i (t) = - 2dt τ ε i (t) - δ 2 RT rel + 2 RT rel ε i (t) dB ε (t) (63)
for each 1 ≤ i ≤ N . The quantities dB v (t) and dB ε (t) respectively are three dimensional and one dimensional Brownian processes. The matrix A v has to satisfy A v A v T = Π. Several choices are available for A v : we find it simpler to use the Choleski decomposition. For the time discretization we use a backward Euler method for velocities. The complete algorithm for the velocity part of the algorithm is explained and validated in [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF].

We now focus on the internal energy equation. Since the diffusion coefficient depends on the square root of the internal energy, this latter variable must remain positive, and hence a simple Euler algorithm cannot work in this case. We will adapt the Milshtein scheme (see [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF]) that reads on our problem:

ε n+1 i =   ε n i + dt τ RT rel B ε 2 + (δ -1)RT rel dt τ   / 1 + 2dt τ , (64)
with B ε a random number chosen from a standard normal law. The global scheme on velocities and internal energies preserves the mass but does not preserve the momentum and the energy, like most Monte-Carlo methods. However, these quantities are preserved in a statistical way. Moreover at the end of each time step, the distribution is renormalized to keep the mean velocity constant and to decrease the error in Θ. Now we focus on the validation of our numerical method. We do not re-validate the velocity scheme alone since its was done in [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF]. We are going to validate the energy part which is new and also the global numerical scheme to see if we can recover the correct Prandtl number.

Validation of the stochastic simulation of the internal energy equation

First, we remind that in an homogeneous case, the temperature T does not change and that the internal temperature varies according to the following law:

d dt T int = - 2 τ (T int -T rel ) = 2 τ θ(T -T int ) (65)
so that the internal energy density is E int (t) = δ 2 RT int (0) exp(-2θ t τ ). We use 1 million particles for validation and the following parameters:

• the initial internal energies are defined by 1000 2 (2s + 1) where s follows a uniform law in [0, 1];

• θ is set to 1 so that T rel = T , and δ 2 RT is equal to 7.500.000;

The internal energies are defined by 40.000s + 40.000 where s follows a uniform law between [0, 1]. The timestep is set to 0.02. The final simulation time is 10τ , where the relaxation time is set to τ = 1 to catpure the relaxation towards equilibrium for velocities and internal energies (whose characteristic time is 5τ ).

On figure 7 we can observe after some time the convergence of the directional tranlational temperatures T 11 , T 22 , T 33 (diagonal components of Θ). towards the same translational temperature T tr , and later the convergence of all these temperatures with the internal energy temperature T int to the equilibrium temperature T .

At convergence we recover the Maxwellian distribution on velocities (figure 8) and the exponential one on internal energies (figure 9).

On figure 10 we observe the convergence of the heat flux q and the quantity T 11 -T tr towards zero. Due to the numerical noise induced by the random method, one has to consider only the beginning of the curves for a relevant comparison with the theory. Equation (57) shows that the profile of log (|q/q(0)|) should have a slope -3 since τ = 1. Numerically we obtain -3.0029. For the temperature, equation (60) shows that the slope of log |T 11 (t) -T tr (t)|/|T 11 (0) -T tr (0)| should be 2(1 -ν(1 -θ)) = 4.2253: numerically, we obtain 4.2091. Therefore the results are in very good agreement with the theory.

Second test case

In this second test case we want to test the correction (54) of ν. To do so we only change the law of the first component of the velocity to 5000s 4 -1000 where s follows a uniform law in [0, 1], which induces now a very nonisotropic tensor Θ. We keep the same numerical parameters except for final time which is set to 1, since after that ν is not corrected and we recover the previous test case. We still have a convergence of the temperatures as before (figure [START_REF] Gorji | A Fokker-Planck based kinetic model for diatomic rarefied gas flows[END_REF]): since we stop the computations at t = 1 relaxation between internal energy and translationnal energy is still in progress. Nonetheless since our initial date is anisotropic in velocity, the correction on ν is activated so that the Prandtl is not correct (we are far from equilibrium...) as it can be seen in figure [START_REF] Gorji | Fokker-Planck model for computational studies of monatomic rarefied gas flows[END_REF]. As a consequence initially the log of the difference (T 11 -T tr ) is not a linear function of time. But as soon as ν is equal to ν air we go back to the first test case and it works (figure 13).

Conclusion

In this paper, we have proposed an extension of the ES-Fokker-Planck model of the Boltzmann equation to polyatomic gases. This model satisfies the properties of conservation of the Boltzmann equation, and it has been defined so as to allow for correct transport coefficients. This has been illustrated by numerical tests for homogeneous problems. Moreover, it has been proved that this new model satisfies the H-theorem. An application of our model to non-homogeneous flows will be proposed in a future work.

  , yields (44). Taking the trace of the previous relation gives (45). The approximation (46) is deduced from the decomposition f = G p (f )(1 + εg), the even parity of G p (f ) with respect to v, and relation (43).

Figure 1 :

 1 Figure 1: Time evolution of the mean internal energy (δ = 6).

Figure 2 :

 2 Figure 2: Time evolution of the relative error on the mean internal energy (δ = 6).

Figure 3 :

 3 Figure 3: Theoretical internal energy distribution at equilibrium vs numerical distribution at final time (δ = 6): the equilibrium state is reached.

Figure 4 :Figure 5 :

 45 Figure 4: Time evolution of the mean internal energy (δ = 2).

Figure 6 :

 6 Figure 6: Theoretical equilibrium vs numerical equilibrium (δ = 2).

Figure 7 :

 7 Figure 7: Convergence of the directional translational temperatures and the internal temperature to their equilibrium value.

Figure 8 :Figure 9 :

 89 Figure 8: Maxwellian equilibrium for the distribution of the first component of the velocities, at time t = 10s: number of numerical particles on the y axis as a function of the V x component

Figure 10 :

 10 Figure 10: Convergence to the equilibrium of the heat flux and a directional temperature.

Figure 11 :

 11 Figure 11: Convergence of the directional translational temperatures and the internal temperature to their equilibrium value.

Figure 12 :

 12 Figure 12: Time evolution of the parameter ν and of the numerical Prandtl number (estimated by (61)).

Acknowledgments. This study has been carried out in the frame of "the Investments for the future" Programme IdEx Bordeaux -CPU (ANR-10-IDEX-03-02).

Now we use proposition 3.3 and we come back to the dimensional variables to get Σ(f ) = pI -µ ∇u + (∇u) T -α∇ • u , and q = -κ∇ • T, up to O(ε 2 ), with the following values of the viscosity and heat transfer coefficients

-γ , and κ = 5 + δ 6 τ pR.

Using these relations into the conservation laws (32) proves that ρ, ρu, and E solve the Navier-Stokes equations (37) up to O(ε 2 ) with transport coefficients given by (39) and that the Prandtl number is

4

The ES-Fokker Planck model with a non constant ν

This short section shows how we adjust the Prandtl number in all cases, following a procedure proposed for monoatomic gases [START_REF] Mathiaud | A Fokker-Planck model of the Boltzmann equation with correct Prandtl number[END_REF].

Some limits of the model with a constant ν

In the previous section, we have found that the Prandtl number obtained for the polyatomic ES-Fokker Planck model is

and can be adjusted to various values by choosing a corresponding value of the parameter ν. Moreover, we have seen that the model is well defined (that is to say that the tensor Π is positive definite for every f ) if, and only if, -1 2 < ν < 1. This last condition leads to the following limitations for the Prandtl number:

1 < P r < +∞.

However, the correct Prandtl number for polyatomic gases, which depends on δ, is smaller than 1, and hence cannot be obtained with our model. In the next section, we show that this analysis, which is based on the inequality ( 22) is too restrictive, and that there is a simple way to adjust the correct Prandtl number.

Recovering the correct Prandtl number

General case

The previous analysis relies on inequality (22) of proposition 2.4 that does not take into account the distribution f itself: this inequality ensures the positive definiteness of Π independently of f . However, proposition 2.4 also indicates that Π is positive definite for a given f if ν satisfies [START_REF] Wu | Oscillatory rarefied gas flow inside rectangular cavities[END_REF]. This inequality is less restrictive, since it depends on f via the temperature T and the extremal eigenvalues of Θ.

Note that here, we use dimensional variables only. Consequently, the parameter ε is here the molecular internal energy, and should not be confused with the Knudsen number used in section 3.

Our goal is to show numerically that we are able to capture the correct Prandtl number for polyatomic gases. This numerical illustration is based on the following remarks. First, the density, velocity, and temperature of f are constant in time: this is due to the conservation properties of the collision operator. Second, the heat flux q defined for g as q := ( 1 2 |v -

so that

Finally, the tensor Θ satisfies the relaxation equation

and its trace satisfies

which shows that Θ tends to RT I for large times, and hence ν tends to the constant value ν air (see ( 54)). Then for large times t > s, ν can be assumed to be constant, and (58) can be solved to give

Using equations ( 60) and (57), one gets

for i = 1, 2, 3, so that we are able to recover the Prandtl number by looking at the long time values of Θ(t) and q(t).

A stochastic numerical method

The distribution function is approximated with N numerical particles so that

where ω i is the numerical weight of numerical particle i and δ V i (t) is the Dirac function at the particle velocity V i (t) and particle internal energy ε i (t). Moreover α is defined through the constant density

• the initial value of δ 2 RT int is set to 2.500.000;

• τ is chosen to be equal to 1;

• the final time is equal to 5 which is sufficient to reach the equilibrium;

• the time step is equal to 0.005 so that we do 1000 iterations.

We will test two values of δ to validate the model.

Polyatomic case: δ = 6

Using equation ( 65) we can compare the behavior of the first moment with the theoretical result.

There is a very good agreement between our results and the theory: we recover the exponential convergence in time of the mean internal energy in figure 1, and the relative error is of order 10 -3 when we compare the analytical solution with the numerical solution of equation 65 in figure 2. Moreover the comparison between the theoretical equilibrium which is Λ g ε δ-2 2 (RT ) -δ 2 exp(-ε/RT ) is very good (figure 3).

Diatomic case: δ = 2

The conclusions are the same for δ = 2:

• good agreement for exponential convergence on internal energy (figure 4),

• relative error is very small (figure 5),

• the correct equilibrium is obtained in figure 6.

Numerical results for the polyatomic model

We now present two test cases to validate the ES Fokker-Planck for polyatomic gases. In the first one the flow is close to isotropic so that the ν correction (54) is not active, while in the second one it will be activated in during an initial transient phase.

First test case: a diatomic gas

In this test case, we use air whose Prandtl number is approximately 0.71. We take the standard value θ = 0.2, and the Prandtl number is adjusted with ν air = -1.3908. The initial data is not too anisotropic so that the correction on ν is not active.

We use 2 million particles. We choose three independent laws for the three components of velocity of the numerical particles:

• the first component is equal to 1000s 4 -200 where s follows a uniform law between [0, 1],

• the second component of velocity is equal to 1000s -500 where s follows a uniform law between [0, 1],

• the second component of velocity is equal to 1000s -500 where s follows a uniform law between [0, 1].