Dimitris K Kardaras
email: kardaras@aueb.gr

Stavros Kaperonis

Stavroula Barbounaki
email: sbarbounaki@yahoo.gr

Ilias Petrounias
email: ilias.petrounias@manchester.ac.uk

Kostas Bithas
email: kbithas@eesd.gr

An Approach to Modelling User Interests using TF-IDF and Fuzzy Sets Qualitative Comparative Analysis

Keywords: User interests, Fuzzy Sets Qualitative Comparative Analysis, TF, IDF

Modelling and understanding user interests are particularly important tasks for designing services and building systems for customized solutions in web personalization and recommender systems. User generated content (UGC) constitutes a significant source of information for capturing user interests. This paper, suggests an approach to user profiling that analyses the Term Frequency (TF) and the Inverse Document Frequency (IDF) of selected tourism services by utilising the Fuzzy set Qualitative Comparative Analysis (FsQCA). It analyses a sample of customer reviews that are collected from tourism web sites. This paper considers the amount of money that customers spent during their hotel stay, as the outcome set in the FsQCA analysis. The results produce causal combinations of services that are necessary and sufficient for building customer interests models that best lead to the outcome and argue for the applicability of the FsQCA in modelling user interests.

Introduction

Recommender systems RC utilise techniques spreading from statistics, to AI and machine learning in order to capture user interests, build user and products/services profiles and suggest the most appropriate products or services to them. RC draw on several methods for developing user references models, with user-generated-content (UGC) to represent a source with rich customer information [START_REF] Martínez-Garcia | Profile of business and leisure travelers on low cost carriers in Europe[END_REF][START_REF] Baka | The becoming of user-generated reviews: Looking at the past to understand the future of managing reputation in the travel sector[END_REF]. Since social media platforms allow users to exchange experience, feedbacks, opinions, complaints, etc., they provide significant information for capturing and understanding user interests [START_REF] Hunecke | Attitude-based target groups to reduce the ecological impact of daily mobility behavior[END_REF]. Web personalisation is another area where user profiling is necessary for developing customised web interfaces, supporting personalised search [START_REF] Zhang | A hybrid fuzzy-based personalized recommender system for telecom products/services[END_REF] that allow users to retrieve search results according to their personal needs.

User Profiling in Tourism

Building user interests models has also been the focus of e-tourism research studies. Drawing on behavioural, socio-economic and demographic data analysis several researchers shed light into understanding people's travel behaviour [START_REF] Hunecke | Attitude-based target groups to reduce the ecological impact of daily mobility behavior[END_REF]. Indeed, surveys on travellers' preferences have shown that the travel selection process is complex depending among others, on personality and mood related factors, service quality issues, the Word-Of-Mouth (WOM) and the eWOM. Customers often express their experience by publishing their reviews. Sentiment analysis of user reviews provides the means for capturing and modelling users' preferences, emotions and attitudes, thus refining market segregation by grouping customers with similar needs and incentives and predicting customers' travel behaviour more precisely [START_REF] Wedel | Market segmentation: Conceptual and methodological foundations[END_REF].

Collantes and Mokhtarian [START_REF] Collantes | Subjective assessments of personal mobility: What makes the difference between a little and a lot?[END_REF] claim that a variety of personality factors such as: personality traits, travel-related behaviours, lifestyle characteristics, and travel trends, determine the subjective assessment of travelling and tourism services. Other researchers have noticed that travel behaviour is influenced by travel experiences and feelings [START_REF] Handy | Driving by choice or necessity?[END_REF][START_REF] Sheller | The new mobilities paradigm[END_REF]. It is also argued that it is important to analyse human behaviour characteristics in order to understand how customers react to alternative transport policies [START_REF] Schade | Acceptability of urban transport pricing strategies[END_REF]. Other travel research studies have analysed environmental factors that influence travel and tourism. Stradling and Anable [START_REF] Stradling | Individual transport patterns[END_REF], argue that environmental characteristics, such as workplace, shops and site topography affect travel choices. Several approaches have been proposed for building user interests models. Kim and Chan [START_REF] Kim | Learning implicit user interest hierarchy for context in personalization[END_REF], have proposed a hierarchical model for representing user interests. The user profile is constructing by analysing documents that users have visited on the web. The documents' analysis yields a list of user interests, which subsequently are grouped upon their similarity on the hierarchical interests' model. It is argued that there exist four classes of information contexts that need to be specified when attempting to understand user interests [START_REF] Joung | A user model for personalization services[END_REF]. The general information class that refers to personal characteristics such as name, contact details, demographics of the user.

The event class represents user's activities. The preference class refers to user's interests. The social network class explains user's connections and interactions with other users. The preference class is usually discovered by analysing various sources such as relevant documents that the user has published [START_REF] Joung | A user model for personalization services[END_REF][START_REF] Bakalov | A Hybrid Approach to Identifying User Interests in Web Portals[END_REF]. Several representational approaches have been proposed for representing user interests. Most frequently though there are three different formats namely: keywords, semantic networks and concept-based representations [START_REF] Gauch | User profiles for personalized information access[END_REF][START_REF] Michlmayr | Learning user profiles from tagging data and leveraging them for personal (ized) information access[END_REF]. Keywords representing domains of interests are associated with weights indicating the strength of user interests for a particular topic. Polysemy and Synonymy are problems associated with keywords. Semantic networks, address these problems, by representing keywords with nodes that are connected with each other, including co-occurrences. Concept-based representations resemble semantic networks in structure but they differ in having nodes to represent abstract topics rather than keywords [START_REF] Gauch | User profiles for personalized information access[END_REF][START_REF] Michlmayr | Learning user profiles from tagging data and leveraging them for personal (ized) information access[END_REF]. User profiles can be used in various ways such as: during personalised information retrieval, that is when a system detects relevant documents and information according to users' interests, during re-evaluating the relevance of documents taking into consideration what documents a user has retrieved and during query processing, when a user query can be modified based on user interests [START_REF] Saleheen | UIWGViz: An architecture of user interest-based web graph vizualization[END_REF].

It is argued that filtering and clustering techniques are very useful in reducing the number of concepts that are found on the web in order to be used in formulating user profiles. However, [START_REF] Saleheen | UIWGViz: An architecture of user interest-based web graph vizualization[END_REF], argues that these techniques lack effectiveness for they produce the same structure of interests for users with different needs. Research show that while many systems produce and use user profiles, e.g. in web personalisation, recommender systems there exists no definite procedure for deriving user interests [START_REF] Saleheen | UIWGViz: An architecture of user interest-based web graph vizualization[END_REF][START_REF] Magnini | Improving user modelling with content-based techniques[END_REF][START_REF] Lehmann | Interactive visualization for opportunistic exploration of large document collections[END_REF][START_REF] Bastian | Gephi: an open source software for exploring and manipulating networks[END_REF]. This paper addresses the need for investigating alternative ways of developing user interests' models and suggests the analysis of the TF-IDF with the FsQCA.

Methodology

The aim of the paper is to identify the causal combinations that are necessary and sufficient to represent customer interests. This paper utilises the FsQCA in order to analyse the TF and IDF of UGC and produce causal combinations that best lead to an outcome. The FsQCA is particularly important for investigating intertwined relationships between multiple factors that affect a dependent variable or contribute to the realisation of certain outcome [START_REF] Chari | Alternative pathways to utilizing customer knowledge: A fuzzy-set qualitative comparative analysis[END_REF]. The FsQCA analyses the sets of relationships among causes. In FsQCA variables are modelled as sets. The FsQCA models allow a detailed analysis of how alternative conditions of causes combine and contribute to high membership scores of the outcome [START_REF] Rihoux | Configurational comparative methods: Qualitative comparative analysis (QCA) and related techniques[END_REF]. FsQCA may detect multiple paths, i.e. alternative causal combinations that can lead to high levels of the same outcome [START_REF] Chari | Alternative pathways to utilizing customer knowledge: A fuzzy-set qualitative comparative analysis[END_REF][START_REF] Skarmeas | Examining the role of CSR skepticism using fuzzy-set qualitative comparative analysis[END_REF]. Data in this paper is collected from customer reviews published on hotel web sites. Causal combinations may be represented by tourism services terms such as room, view, cleanliness, etc., in the set of selected documents. The outcome set in this paper, is represented by the large amount of money spent by the customer. Other outcome sets can also be considered. Thus, this paper aims to identify the combinations of customer hotel services interests that best reflect customer's spending. A sample of the data collected is analysed in this paper. The steps of the methodology are shown below:

1. Select documents published by user ) ( i u .

2. Identify the terms that will constitute the causal combinations and specify the term that will represent the outcome set. 3. Calculate the (TF) and the (IDF) for each identified term. The fuzzy union, is defined as

) , max( ) ( B A B A      , (2)
The fuzzy intersection is defined as

) , min( ) ( B A B A      (3)
and the fuzzy complement is calculated as

A A     1 (4) 6.
Calculate the consistency and the coverage of the solutions using formulas [START_REF] Baka | The becoming of user-generated reviews: Looking at the past to understand the future of managing reputation in the travel sector[END_REF] and ( 3) respectively.



  X Y X Y X y Consistenc ) , min( ) (  [24] (5)
   Y Y X Coverage ) , min( [24] (6)
where ) (X is the membership degree of each causal combination and ) (Y is the membership degree of the outcome set. 7. Identify best combinations, by selecting the combinations that exhibit a consistently rate above a threshold (in this paper is set at 0.8) and the highest possible coverage. Simplify solutions into the final set of causal combinations.

The final causal combinations indicate the hotel services that customers who spend large amount of money consider as the most important.

Data Analysis: Illustrative Example

This paper analyses reviews collected from five (5) hotel customers. Then, for simplicity reasons, five (5) terms representing hotel services are selected from the total set of terms identified in the reviews. The outcome set large amount of money spent (LMSp) by each user during his/her hotel stay is represented as triangular fuzzy numbers (TFN). The membership function

) (x f A of TFN ) , , ( ~b m a A
can be calculated according to the following equation [START_REF] Lin | A fuzzy-based decision-making procedure for data warehouse system selection[END_REF]:

                      otherwise b m b x m b m b x a m m x a a m a x x f A , 0 , , , , (7)
where a, m, b are real numbers. The linguistic scales which are used and their corresponding TFNs adopted in this study are shown in table 1 ) ( i u . Then, the weights for each term result from using formula [START_REF] Martínez-Garcia | Profile of business and leisure travelers on low cost carriers in Europe[END_REF]. The results are shown in Table 2.

1 0 0 0 0 0 2 0 0 0 0 1 3 0 0 0 1 0 4 0 0 0 1 1 5 0 0 1 0 0 6 0 0 1 0 1 7 0 0 1 1 0 8 0 0 1 1 1 9 0 1 0 0 0 10 0 1 0 0 1 11 0 1 0 1 0 12 0 1 0 1 1 13 0 1 1 0 0 14 0 1 1 0 1 15 0 1 1 1 0 16 0 1 1 1 1 17 1 0 0 0 0
The cells in the truth table take the value (1) or (0) representing true or false. Thus, permutation number 3 is read (Quietness=false, Sea View=false, Staff Friendliness=false, Cultural Activities=true, Restaurant=false). Next the membership degrees for all combination for each user are calculated drawing on the fuzzy sets operations theory. Table 4 shows the membership degrees for the first 17 combinations.

 =  (Quietness=false  Sea View=false  Staff Friendliness=false  Cultural Activities=true  Restaurant=false) =  (not (Quietness), not (Sea View), not (Staff Friendliness), Cultural Activities, not (Restaurant)).
The  (Quietness=false) =  ((1- (Quietness)) = (1-0.3)=0.7. Similar calculations are performed for all terms thus, 3 C  =min(0.7; 0.5; 0.6; 0.3)=0.3. After all membership degrees are calculated the consistency and coverage degrees are determined. Table 5 shows the results for the first 17 combinations.



 X (0.3+0.1+0.5+0.3+0.3)=1.5.

Therefore the consistency for combination number 3=0.733.

Regarding the coverage, by applying formula (6),   ) , min( Y X

1.5 and   Y 2.9 thus coverage=0.37.

According to FsQCA the best causal combinations should exhibit as high as possible consistency and coverage. However, the higher the consistency is the lower the coverage. Assuming a threshold value of 0.8 for the consistency firstly and then the higher possible coverage, the analysis results into two causal combinations; the combinations number 12 and 16 extracted from Table 3, are shown in Table 6. A closer look at the combinations reveals that "quietness" is not within the customers interests at all. It is not a necessary service. Thus, restructuring the causal combination the analysis results that customers who spend a large amount of money, show interest in ➢ (Sea View) AND (Staff friendliness) AND (Cultural activities) AND (Restaurant) OR ➢ (Sea View) AND (Cultural activities) AND (Restaurant).

In order to simplify the causal combinations, the "staff friendliness" could be omitted for it does not appear on both combinations.

Conclusions-Future Research

This study suggests that the FsQCA can be used for modelling users' interests. Data selected from customer reviews is analysed by utilising the TF and the IDF. The application of the FsQCA results into useful insights that can be used to understand customer priorities and build customer profiles. Future research can focus on examining the applicability of the FsQCA to handle multiple outcome sets and to specify terms' priorities. When applying the FsQCA method in large data sets with a long list of factors, the truth table and the set of possible causal combinations can become cumbersome to analyse. Thus, future research can focus on combining the FsQCA analysis with other techniques that will be used in pruning the size of the truth table and reduce the causal combinations to manageable size.

  Apply the FsQCA and produce User Interests causal combinations. a. Produce the truth table of all possible permutations of the terms considered.Each permutation is a possible causal combination. b. Calculate membership degrees for each combination. Its calculation is performed drawing on the fuzzy sets operations theory. Assume two fuzzy

	tk W 	tk TF	*	log(	tk i N d	)	[23]	(1)
	4. Calculate the weight of each term where, tk W , represents the weight of term ) ) ( k t using the following formula: ( k t , tk TF , is the term frequency for term ) ( k t , i N , is the total number of documents published by user ) ( i u and tk d , represents the number of documents that contain term ) ( k t . 5.

sets

A ~and B ~then:

Table 1 .

 1 . Linguistic scales and corresponding TFNs for Large Amount of Money-Spent fuzzy sets

	Linguistic scale	Triangular fuzzy scale	Mean of fuzzy numbers
	Very High	(0.75, 1.00,	1.00)	1.00
	High	(0.50, 0.75,	1.00)	0.75
	Medium	(0.25, 0.50,	0.75)	0.50
	Low	(0.00, 0.25,	0.50)	0.25
	Very Low	(0.00, 0.00,	0.25)	0.00
	The linguistic scales represent indicate to what extent a customer is included to the set
	of those who spend large amount of money during their hotel stay. The TF and IDF
	scores (step 3) are calculated by using the KNIME tool, for all documents published
	by each user

Table 2 .

 2 The term weights and the membership degree for money spent for each Customer  .Table3shows part of the truth

	Large Amount Spent membership		Terms' Weights Customer	tk W based on TF-IDF for each
	degree Outcome Set (Y)	Customer	Quietness	Sea View	Staff Friendli ness	Cultural Activities	Resta urant
	0.50	1	0.30	0.50	0.40	0.70	0.70
	0.70	2	0.30	0.70	0.60	0.70	0.90

table . Table 3 .

 .3 The truth table (part of) show all possible permutations of the terms

	Causal Combination	Quietness	Sea View	Staff Friendliness	Cultural Activities	Restaurant

Table 4 .

 4 Membership degrees for combinations for each customer

	Causal	Customer	Customer	Customer	Customer	Customer
	Combination	1	2	3		4	5
	1	0.3	0.1	0.4		0.3	0.3
	2	0.3	0.3	0.4		0.3	0.3
	3	0.3	0.1	0.5		0.3	0.3
	4	0.5	0.3	0.5		0.3	0.3
	5	0.3	0.1	0.2		0.3	0.3
	6	0.3	0.3	0.2		0.3	0.3
	7	0.3	0.1	0.2		0.3	0.3
	8	0.4	0.3	0.2		0.3	0.3
	9	0.3	0.1	0.3		0.3	0.3
	10	0.3	0.3	0.3		0.5	0.3
	11	0.3	0.1	0.3		0.3	0.3
	12	0.5	0.4	0.3		0.5	0.4
	13	0.3	0.1	0.2		0.3	0.3
	14	0.3	0.3	0.2		0.4	0.3
	15	0.3	0.1	0.2		0.3	0.3
	16	0.4	0.6	0.2		0.4	0.6
	17	0.3	0.1	0.1		0.3	0.3
	The membership degree of combination number 3	3  , for customer-1, see framed C
	cell in table 4, is calculated as follows by using formulas (3) and (4):
	Consider combination number 3 membership degree	C
						3

Table 5 .

 5 Causal combinations' Consistency and Coverage

	Causal Combination	Consistency	Coverage
	1				0.785714286	0.379310345
	2				0.8125	0.448275862
	3				0.733333333	0.379310345
	4				0.789473684	0.517241379
	5				0.916666667	0.379310345
	6				0.928571429	0.448275862
	7				0.916666667	0.379310345
	8				0.933333333	0.482758621
	9				0.846153846	0.379310345
	10				0.882352941	0.517241379
	11				0.846153846	0.379310345
	12				0.904761905	0.655172414
	13				0.916666667	0.379310345
	14				0.933333333	0.482758621
	15				0.916666667	0.379310345
	16				0.954545455	0.724137931
	17				1	0.379310345
	The consistency for combination number 3 is calculated, by applying formula (5) as
	follows: Consider the outcome column (Y) shown in Table 2 and the membership
	degrees (X) of combination number 3, for all users as shown in Table 4. Then,
		, min( Y X	)	
	min{min(0.3;0.5)+min(0.1;0.7)+min(0.5;0.1)+min(0.3;0.7)+min(0.3;0.9)=
	min(0.3+0.1+0.1+0.3+0.3)=1.1.

Table 6 .

 6 The two necessary and sufficient causal combinations

	Causal	Quietness Sea	Staff	Cultural	Restaurant
	Combination		View	Friendliness	Activities
	16	0	1	1	1	1
	12	0	1	0	1	1


Di Zhou

Yuxin Zhao
email: zhaoyuxin@hrbeu.edu.cn

Chang Liu
email: liuchang407@hrbeu.edu.cn

Yanlong Liu
email: yanlong_liu@hrbeu.edu.cn

Speedup of Network Training Process by Eliminating the Overshoots of Outputs

Keywords: Output Overshoots, Local Linearization, Weighted Mean, Artificial Neural Network

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Artificial neural networks (ANNs) have developed rapidly and have been used wildly in many fields, such as speech recognition, computer vision, games and classification of skin cancer [START_REF] Hinton | Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups[END_REF][START_REF] Gould | Darwin: A framework for machine learning and computer vision research and development[END_REF][START_REF] Silver | Mastering the game of go with deep neural networks and tree search[END_REF][START_REF] Esteva | Dermatologist-level classification of skin cancer with deep neural networks[END_REF]. So far, the studies of ANNs can be divided into three aspects. The first aspect is the topology of network. The numbers of nodes and layers play a significant role in applications of networks. ANNs have developed multiple branch types, such as convolution neural network, recurrent neural network, extreme learning machine [START_REF] Mahendran | Visualizing deep convolutional neural networks using natural pre-images[END_REF][START_REF] Weninger | Introducing currennt: The munich opensource cuda recurrent neural network toolkit[END_REF][START_REF] Huang | Extreme learning machine: theory and applications[END_REF]. However, there is still no substantial instruction on the topology of network. The second aspect is to combine the neural network with other optimization algorithms. Since lots of optimization algorithms exist, this aspect is mainly used for special applications. In [START_REF] Lin | Identification and prediction of dynamic systems using an interactively recurrent self-evolving fuzzy neural network[END_REF] and [START_REF] Yu | A short-term load forecasting model of natural gas based on optimized genetic algorithm and improved bp neural network[END_REF], ANNs are combined with fuzzy and genetic theory to identify nonlinear system and a short-term load forecasting of natural gas. Network training algorithms are the third aspect and most popular aspect.

Training algorithms play a decisive role in the accuracy of ANNs. In general, a large number of historical data have to be used while training network. Therefore, how to handle big data efficiently and extract features accurately are the key to training network. In [START_REF] Rumelhart | Learning representations by backpropagating errors[END_REF], Rumelhart proposed learning representations by back-propagating errors. BP algorithm has good convergence property, however, it may lead to the problem of gradient vanishing. For the disadvantages of BP algorithm, ReLU is proposed to replace the sigmoid activation function [START_REF] Dahl | Improving deep neural networks for lvcsr using rectified linear units and dropout[END_REF]. In [START_REF] Srivastava | Dropout: a simple way to prevent neural networks from overfitting[END_REF], dropout of nodes is proposed to prevent ANNs from overfitting. In addition, the parallel and online sequential learning algorithms are introduced to quicken the training process [START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF][START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF].

Almost all existing algorithms focus on gradient while few studies take the overshoots of outputs into consideration. An eliminating overshoot method (EOM) is proposed to improve BP algorithm. An overshoot definition is given on the basis of traditional BP (TBP) algorithm. Meanwhile, local linearization and weighted mean algorithms are adopted to train ANNs. This paper is organized as follows. In Sect.2, we briefly describe the overshoot definition and its mathematical method. Sect.3 introduces the numerical simulation. Finally, the conclusion and the future work are summarized in Sect.4.

The Overshoot Definition and Mathematical Method

In this study, the OD combines the outputs with gradient decent. The OD is given as all outputs must be on the same side of the excepted output with the output at the first iteration while training network. The purpose of training network is to adjust the parameters on the basis of a given set of samples to minimize the errors between the actual and expected outputs. Not only gradients but also outputs have to be considered while training networks. A cost function is defined to characterize the errors. At present, there are many kinds of cost functions. The Euclidean metric is used as the cost function in this study, which is defined as follows:

1 ( ) ( ) 2 T E O D O D      (1)
O is the actual output vector, D is the expected response output. T denotes transpose operation. E is the total instantaneous error energy.

The weight and threshold matrixes of the th i layer are defined as follows:

1 1 , 1,1 ,1,2 ,1, ,2,1 ,2,2 ,2, ,1,1 ,1,1 ,1,1 i
i i i i N i i i N i i i i w w w w w w W w w w                  (2) 1 2 i i i i i N b b B b         (3)
,, i j k w denotes the weight value between the th j node of the th i layer and the th k node of the ( 1) th i  layer.

j i b denotes the threshold value of the th j node in the th i layer.

i N denotes the number of the nodes in the th i layer.

The output of the th i layer is defined as follows:

1 2 i i i i i N o o O o         (4)
According to the equation ( 2)-( 4), the output of the ( 1) th i  layer can be expressed as follows: 1 ( , , )

i i i i O F W B O   (5)
F is the action function.

The correction i W  and i B  are defined as follows:

1 1 1 ,1,1 ,1,2 ,1, ,2,1 ,2,2 ,2, , ,1 , ,2 , , i i
i i i i i i i N i i i N i i i N i N i N N E E E w w w E E E dE w w w W dW E E E w w w                                                     (6) 1 2 i i i i i i N E b E dE b B dB E b                                      (7)
In equation [START_REF] Huang | Extreme learning machine: theory and applications[END_REF],

i dE dB
is the jacobian matrix of E with respect to the threshold vector i B denoted as B i J . The jacobian matrix can be calculated layer by layer.

1 1 L i i B i i L L i i dE dE dO dO dO J dB dO dO dO dB      , (8)
Where 1, 2, , iL 

.

In equation( 8), L dE dO denotes the jacobian matrix of E with respect to the actual output L O . L denotes the number of all layers in network.

1 1 2 2 12 L L L L L L L L L NN L L L L N dE E E E O d O d O d dO O O O                        . ( 9
)
L i d denotes the expected response.

i O . 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2 2 1 1 2 1 1 1 1 2 i i i i i i i i i i i i L i i i i i i i L i i i i L L L i i i L o o o o o o o o o dO o o o dO o o o o o o                                                         ( 10
) i i dO dB denotes the jacobian matrix of i O with respect to i B . 1 1 1 1 2 2 2 2 1 2 1 2 i i i i i i i i i i i i N i i i i i i i N i i i i N N N i i i N O O O b b b O O O dO b b b dB O O O b b b                                            (11)
According to the equations ( 7)- [START_REF] Dahl | Improving deep neural networks for lvcsr using rectified linear units and dropout[END_REF], the initial correction i B  can be obtained.

The expression of the initial correction i W  can be derived with the same method above. In this study, the sigmoid function is chosen as the action function. The output locates in the interval (0,1). It can be obtained from the equation ( 8) that the correction is determined by the errors of outputs directly.

Assuming that the overshoot arises at the th n iteration, it indicates the corrections exceed the range. Accordingly, the OD should be adopted to eliminate the overshoots for achieving the optimal solution as soon as possible. Since the traditional line search algorithm is time consuming, local linearization and weighted mean methods are used to adjust the parameters to eliminate overshoots in this paper. The local linearization method is based on the relationship between cost function and network parameters. When overshoots occur, the actual outputs of the latest two iterations locate on both sides of the excepted output, which are close to the expected output. In this situation, the error changes little, which is in line with the local linearization condition. The network parameters are regarded as the independent variables, while the cost function is the corresponding variable. Linear equations are established by the latest two iterations, which are defined as follows:

1

i i i n n n n i i i n n n n E E B B B E E B B B                   (12) i n B  is the gradient of E with respect to i B at the th n iteration. i n B  1 1 () ()
is the updated network parameters.

n E is the total error at the th n iteration.

The updated network parameters can be obtained by solving the equations when the gradients are nonzero. According to equation ( 12), the updated parameters can be divided into two cases. For the first case, if the cost function is a convex function ,))

i i i i n n n n B B B B  .
The weighted mean method is adopted to deal with this situation. The specific process is defined as follows:

1 1 11 | | | | | | | | | | | | ii ii nn nn i i i i n n n n BB BB B B B B               (13)
Similarly to the first case, the updated parameters are accepted if the updated parameters satisfy the overshoot definition, or this step will continue to circulate until overshoots are eliminated.

Numerical Simulation and Result Analysis

In this section, a numerical simulation is conducted to verify the effectiveness and practicability of the EOM algorithm. A three layers feed-forward ANN is designed.

The sigmoid function is chosen as the action function. The numbers of neurons in input, hidden and output layer are 22, 30 and 58. The input and output samples are selected from the solutions of the downwelling atmosphere monochromatic radiative transfer model equation (MonoRTM). Before training the network, all samples are normalized.

In order to avoid the denominator is 0, the network parameter is accepted when the absolute values of its gradients in the two iterations satisfy the overshoot definition and are both below 0.000001. This operation will also speed up the training process.

Similarly to the TBP methods, the OD is also introduced to train the network layer by layer. However, contrary to the direction the error propagates, the OD propagates forwardly. The parameters of the next layer are adjusted on the basis of no overshoot on the previous layers. The TBP algorithm is simulated to compare with the proposed method. The structure of network and initial parameters of the two methods are set to the same. The maximal number of iterations and the acceptable error are set to 50 and 0.000001.

The total instant error of all samples at each iteration is recorded and shown in Fig. 1. The blue and red curves represent the total error of the TBP and the EOM method respectively. As can be seen from the figure, the errors of the two methods almost overlap in the initial stage. This is mainly because no overshoot occurs at this stage. The network parameters meet the OD at each iteration which implies that they are not corrected again.

However, the error of the EOM method is reduced at a slower speed than the TBP method in the middle stage. The experimental data shows that errors of the most sample are still large, but overshoots appear at some sample points. In order to satisfy the OD, the parameters leading to overshoots are corrected again, resulting in that the error cannot reduce in the direction of the steepest descent. Therefore, the blue curve is lower than the red curve in the figure. However, the gradient of the blue curve slows down gradually while the red curve almost shows a linear decline trend.

In the final stage, the errors of the samples are small. Overshoots in many sample points lead to a slow speed of training network with traditional methods. However, the EOM method takes great advantage of the good suppression of overshoots and minimizes the error quickly. As shown in Fig. 1, it takes 14 iterations to achieve the optimal solution for the TBP method, however, only 11 for the EOM method. It can be manifested from Fig. 2, no overshoot appears in the initial stage. It indicates the parameters satisfy the OD. The blue and red curves overlap each other which is consistent with the result in Fig. 1. The correctness of the results in Fig. 1 can be confirmed in this stage.

The blue curve implies the overshoot appears when training network with traditional BP method in the middle stage. In order to meet the OD, the parameters have to be corrected again by the local linearization and weighted mean methods. The red curve depicts the elimination of overshoots of outputs while training neural network. In addition, the EOM method takes 8 iterations to approximate the excepted output which is half of that of the TBP method. Furthermore, the output curve of the EOM method is smoother than that of TBP method.

In order to illustrate the effectiveness of the proposed algorithm to eliminate overshoot intuitively, the training trajectories of 58 output samples of the two methods are shown in Fig. 3 and Fig. 4. The blue points in Fig. 3 and green points in Fig. 4 represent the trajectories of the TBP and EOM method respectively. The red points in both figures are the excepted outputs. Overshoots appear at 90 percent of samples in Fig. 3. And overshoots are obvious especially when the initial error is around 0.15. However, no overshoot appears in Fig. 4 which is consistent with the results in Fig. 2. No overshoot for all samples indicates the proposed method is universal. And the OD is universal for other neural network model.

Conclusion

In this paper, an improved training method and the OD are proposed for eliminating the overshoots of outputs. The outputs and gradients of the network are taken into consider. The training process of neural network is derived in matrix form and the specific expressions are given. Local linearization and weighted mean methods are introduced to optimize training algorithm. In the final stage of training network, the proposed method eliminates all overshoots perfectly. The advantage of eliminating overshoot makes it possible to reach the optimal earlier than the traditional method.

The simulation demonstrates effectiveness of the EOM method and its universality for all samples. Therefore, the EOM method is superior to the traditional method. In addition, a drawback of the EOM method is that it trains network at a slow speed in the middle stage. It is also the direction of our future work.



  denotes the jacobian matrix of the output of the ( 1) th i  layer 1 i O  with respect to the output of the th i layer

.

  Then the overshoot of the updated parameters is judged according to the OD. If overshoots are eliminated, the updated parameters are accepted and continue the next iteration. Otherwise, the above step has to be circulated until the overshoots are eliminated. For the second case, the cost function is a non-

Fig. 1 .

 1 Fig. 1. The total errors of outputs of the EOM and BP methods. The red and blue lines represent the outputs of EOM and BP methods, respectively.

Fig. 2 .

 2 Fig. 2. The outputs of single sample of the EOM and BP methods. The red and blue lines represent the outputs of EOM and BP methods, respectively. The green line is the excepted output

Fig. 3 .

 3 Fig. 3. The output trajectories of the TBP method.

Fig. 4 .

 4 Fig. 4. The output trajectories of the EOM method.

Acknowledgements

This paper is funded by the International Exchange Program of Harbin Engineering University for Innovation-oriented Talents Cultivation, National Natural Science Foundation of China (Grant Nos. 41676088) and the National Key Research and Development Project of China (2017YFC1404100).


Benoît Trouette
email: benoit.trouette@u-pem.fr

Lattice Boltzmann simulations of a time-dependent natural convection problem

Keywords: Thermal lattice Boltzmann method, double distribution function, time-dependent natural convection, Richardson extrapolation, convergence study

A two-dimensional double Multiple-Relaxation-Time thermal lattice Boltzmann method is used to simulate natural convection flows in differentially heated cavities. The buoyancy effects are considered under the Boussinesq assumption. Flow and temperature fields are respectively solved with nine and five discrete velocities models. Boundary conditions are implemented with the classical bounce-back or a "on-node" approach. The latter uses popular Zou and He and Counter-Slip formulations. This paper evaluates the differences between the two implementations for steady and time-dependent flows as well as the space and time convergence orders.

Introduction

The Lattice Boltzmann Method (LBM), derived from the lattice gas automata [START_REF] Mcnamara | Use of the Boltzmann equation to simulate lattice-gas automata[END_REF], has been developed as an alternative numerical scheme for solving the incompressible Navier-Stokes equations. It has demonstrated its ability to simulate several physical systems [START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF]. Its straightforward implementation, natural parallelism and easy boundary condition treatment [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] make it very efficient and accurate for hydrodynamical flows. However, for thermal flows, its efficiency is a pending issue. Indeed, thermal lattice Boltzmann methods can be split into three main classes. The first one relies on increasing the number of discrete velocities in order to make the original athermal LBM able to solve correctly the temperature field [START_REF] Alexander | Lattice Boltzmann thermohydrodynamics[END_REF]. The main disadvantages of this approach are the loss of the cellular automata transport scheme and a fixed Prandtl number. The second one is the hybrid approach. Mass and momentum are solved with a lattice Boltzmann model while the convection-diffusion equation for temperature is solved with a classic macroscopic solver like Finite-Volume or Finite-Element methods [START_REF] Lallemand | Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions[END_REF][START_REF] Mezrhab | Hybrid lattice-Boltzmann finite-difference simulation of convective flows[END_REF][START_REF] Jami | Lattice Boltzmann method applied to the laminar natural convection in an enclosure with a heat-generating cylinder conducting body[END_REF]. Counterparts are an harder parallelization and the resolution of a linear system. Finally, the third one, the double distribution function formulation uses two evolution equations: one for the mass and momentum conservation and one for the temperature [START_REF] Guo | A coupled lattice bgk model for the boussinesq equations[END_REF][START_REF] Orazio | Computational Science -ICC 2003, chapter Boundary Conditions for Thermal Lattice Boltzmann Simulations[END_REF][START_REF] Mezrhab | Double MRT thermal lattice Boltzmann method for simulating convective flows[END_REF]. This approach preserves the natural advantages of lattice methods: the local formulation and the explicit time evolution scheme. Besides, supplementary distribution functions can be added to take into account additional effects like magnetic field [START_REF] Ashorynejad | Magnetic field effects on natural convection flow of a nanofluid in a horizontal cylindrical annulus using lattice Boltzmann method[END_REF]. Mass and momentum equations can be solved with the generalized, or Multiple Relaxation Times (MRT), formulation [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability[END_REF][START_REF] Humières | Multiple-Relaxation-Time Lattice Boltzmann Models in Three Dimensions[END_REF]. Two approaches stand out for the thermal part: the internal energy formulation [START_REF] He | A novel thermal model for the lattice Boltzmann method in incompressible limit[END_REF][START_REF] Peng | Simplified thermal lattice Boltzmann model for incompressible thermal flows[END_REF] and lattice Boltzmann schemes designed for convection-diffusion equations [START_REF] Suga | Numerical schemes obtained from lattice Boltzmann equations for advection diffusion equations[END_REF]. The former might be suitable for non straight boundary condition [START_REF] Yan | Numerical simulation of heat transfer and fluid flow past a rotating isothermal cylinder -A LBM approach[END_REF] and fluid-solid conjugate heat transfer [START_REF] Wang | A lattice Boltzmann algorithm for fluid-solid conjugate heat transfer[END_REF][START_REF] Meng | Lattice Botlzmann simulations of conjugate hear transfer in high-frequency oscillating flows[END_REF]. The latter offers more flexibility in the discrete velocities set [START_REF] Huang | Numerical study of lattice Boltzmann methods for a convection-diffusion equation coupled with Navier-Stokes equations[END_REF] and recent developments show that evolution equation can be solved with a MRT collision operator [START_REF] Mezrhab | Double MRT thermal lattice Boltzmann method for simulating convective flows[END_REF][START_REF] Moussaoui | A computation of flow and heat transfer past three heated cylinders in a vee shape by a double distribution MRT thermal lattice Boltzmann model[END_REF][START_REF] Mohamad | Lattice Boltzmann Methods, Fundamentals, Applications, with Computer Codes[END_REF][START_REF] Wang | Lattice Boltzmann simulation of thermal convective flows in two dimensions[END_REF]. In many cases, lattice Boltzmann method is used for stationary flows and do not enjoy the benefits of the inherent unsteady formulation. Indeed, the physical transient flow has to be simulated and calculations of accurate solutions would be very expensive compared to stationary solvers. Considering this background, this work aims to analyze results and convergence orders of the double MRT method for a timedependent natural convection flow configuration. Reference values are extracted from benchmark solutions on differentially heated cavities [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high rayleigh number[END_REF][START_REF] Xin | An extended Chebychev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF]. Proper comparisons with theoretical second-order accuracy both in space and time [START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF] are also made thanks to Richardson Extrapolation. Because of numerical instability, this analysis can hardly be handled using the classical LBM method with the Bhatnagar-Gross-Krook (BGK) collision term. This paper is organized as follows. Section 2 presents the general configuration of the differentially heated cavity under Boussinesq approximation. In section 3, the numerical method is described: a double MRT lattice Boltzmann method using respectively nine and five velocities models for the fluid flow and the energy equation is detailed. The obtained results for two configurations are presented in Section 4. Comparisons with reference and convergence orders are discussed. Finally, the last section is dedicated to the concluding remarks.

Problem description

The configuration studied is the natural convection in a two-dimensional cavity heated differentially on vertical side walls. The configuration is illustrated in figure 1 where W is the width and H the height of the cavity. The cavity aspect ratio is A = H/W . The gravity vector is directed in the negative y-coordinate direction. The Boussinesq hypothesis is used and only small temperature variation from the mean temperature are admitted.

Hot

             ∇ • u = 0 ∂u ∂t + u • ∇u = -∇p + P r Ra ∇ 2 u + θe y ∂θ ∂t + u • ∇θ = 1 RaP r ∇ 2 θ (1)
where u = (u, v), p and θ are respectively the velocity, pressure and temperature fields and e y is the unit vector in the ydirection. This dimensionless system was obtained using the characteristic length W , buoyancy velocity scale U = √ gβW ∆T , time scale W/U and pressure scale ρU 2 . Here, ρ is the mass density, g the gravitational acceleration and β the coefficient of thermal expansion. The dimensionless temperature is defined as follows:

θ = T -T ref T h -T c with T ref = T h + T c 2 (2)
and T h and T c are respectively the prescribed temperatures of hot and cold walls. The Rayleigh and Prandtl numbers are control parameters of the problem and are written:

Ra = gβ∆T W 3 να and P r = ν α ( 3
)
where α is the thermal diffusivity, ν the kinetic viscosity, and ∆T = T h -T c the temperature difference between the hot and the cold walls. For all calculations, the Prandtl number is P r = 0.71 (air at STP) and the aspect ratio A and the Rayleigh number Ra are free parameters. The cavity boundary conditions for the velocity field are a no-slip conditions, i.e. u = v = 0 on all walls. On the left and right walls, the prescribed temperature boundary conditions is written:

θ(x = 0, y) = + 1 2 and θ(x = 1, y) = - 1 2 (4)
Along the horizontal bottom and top walls, the zero-flux condition reads ∂θ ∂y y=0 = 0 and ∂θ ∂y y=A = 0

For all calculations, solution is initialized with an isothermal fluid at rest, .i.e.

u(u, v) = 0 and θ = 0 (6) Two configurations will be studied. The first one, the classical square cavity (A = 1) with Ra = 10 6 , is used as a validation. In the second case, an aspect ratio A = 8 and a Rayleigh number Ra = 3.4 × 10 5 , lead to a time-dependent periodic flow, also referenced as a benchmark solution. In all cases, the Prandtl number is set at 0.71.

Numerical methods

We consider the double Multiple-Relaxation-Times thermal lattice Boltzmann (MRT-TLB) model introduced in [START_REF] Mezrhab | Double MRT thermal lattice Boltzmann method for simulating convective flows[END_REF][START_REF] Moussaoui | A computation of flow and heat transfer past three heated cylinders in a vee shape by a double distribution MRT thermal lattice Boltzmann model[END_REF][START_REF] Wang | Lattice Boltzmann simulation of thermal convective flows in two dimensions[END_REF]. Two sets of distribution functions are used. The evolution equations for both mass and momentum and thermal distributions will be described as well as the coupling. Boundary conditions treatment will also be presented.

The general principle of lattice Boltzmann methods is the following: the phase space is discretized into a regular lattice, or mesh, and into a finite set of N symmetrical discrete velocities, e i=1,...,N . Lattice nodes are x j=1,...,Nx×Ny where N x and N y stand respectively for the number of nodes in x-and y-directions. A set of distribution functions, associated with discrete velocities, is defined on each node. Their evolutions are governed by a free transport, associated with the discrete time step δ t , along the discrete velocities. Then, a collision process makes the distribution functions relax through an equilibrium.

During collisions, macroscopic quantities such as mass and momentum are conserved.

Multiple relaxation time method for mass and momentum

The time evolution equation for the mass and momentum conservations can be written as follows:

f (x j + e δ t , t + δ t ) -f (x j , t) = -Q(f (x j , t) -f eq (x j , t)) + F (x j , t) (7)
where Q is the collision matrix and F (x j , t) an external force. The collision matrix Q can be written in the general form:

Q = M -1 SM (8)
and the evolution equation (7) becomes f (x j + e δ t , t + δ t ) -f (x j , t) = -M -1 S(m(x j , t) -m eq (x j , t))

+ F (x j , t). (9)
The quantities, f (x j , t), f eq (x j , t), m(x j , t), m eq (x j , t) and F (x j , t) are N -tuple vectors and the superscript eq stands for the equilibrium values. The general notation for a N -tuple vector φ expressed at time t is:

φ(x j , t) = (φ 1 (x j , t), . . . , φ N (x j , t)) T φ(x j + e δ t , t) = (φ 1 (x j + e δ t , t), . . . , φ N (x j + e δ t , t)) T [START_REF] Orazio | Computational Science -ICC 2003, chapter Boundary Conditions for Thermal Lattice Boltzmann Simulations[END_REF] where T denotes the transpose matrix. The velocity distribution function f (x j , t) is expressed in the velocity space V = R N while its corresponding moment m(x j , t) is in the moment space M = R N . With the above expressions, the mapping between discrete velocity space V and moment space M is achieved by the transformation matrix M which maps the vector f (x j , t) to the vector m(x j , t):

m(x j , t) = M f (x j , t) and f (x j , t) = M -1 m(x j , t) (11)
The operator M is a square matrix and its rank is equal to the number of discrete velocities N . With the N = 9 velocities model (d2q9), discrete velocities are:

e i =    (0, 0) i = 1 (1,
0)e, (0, 1)e , (-1, 0)e , (0, -1)e i = 2 -5 (1, 1)e, (-1, 1)e, (-1, -1)e, (1, -1)e i = 6 -9 [START_REF] Ginzburg | Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium[END_REF] where e = δ x /δ t , and δ x and δ t are the lattice spacing and discrete time steps. On the lattice, CFL consideration generally leads to δ x = δ t = 1. To determine the transformation matrix M , the ordering of the moments must be prescribed first. Here, the choice is: m(x j , t) ≡ m = (m 1 , m 2 , m 3 , m 4 , m 5 , m 6 , m 7 , m 8 , m 9 ) T = (ρ, j x , j y , e, p xx , p xy , q x , q y , ) T

where ρ is the mass density, j = (j x , j y ) = ρ(u, v) = ρu is the impulsion and u is the flow velocity, e, p xx and p xy are secondorder moments corresponding to energy, and two off-diagonal components of the stress tensor, respectively. q x and q y are the third-order moments corresponding to x-and y-components of the energy flux and is the fourth-order moment of energy square. With this ordering of the moments, the transformation matrix can be constructed thanks to the Gram-Schmidt orthogonalization [START_REF] Lallemand | Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions[END_REF]:

M =               1 1 1 1 1 1 1 1 1 -4 -1 -1 -1 -1 2 2 2 2 4 -2 -2 -2 -2 1 1 1 1 0 1 0 -1 0 1 -1 -1 1 0 -2 0 2 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 0 0 -2 0 2 1 1 -1 -1 0 1 -1 1 -1 0 0 0 0 0 0 0 0 0 1 -1 1 -1               (14)
M being an orthogonal matrix, its inverse can be computed according to the formula:

M -1 = M T (M M T ) -1 .
In the case of an isothermal fluid, the conserved variables are only the density ρ and the momentum j. The other moments are nonconserved moments and as suggested in Ref. [START_REF] Lallemand | Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions[END_REF], they relax linearly towards their equilibrium values. For the d2q9 model, the non-conserved moments (i = 4 -9) are:

m eq 4 = e eq = -2ρ + 3ρ 0 u • u m eq 5 = p eq xx = ρ 0 (u 2 -v 2 ) m eq 6 = p eq xy = ρ 0 uv m eq 7 = q eq x = -ρ 0 u m eq 8 = q eq y = -ρ 0 v m eq 9 = eq = ρ -3ρ 0 u • u (15)
where the constant ρ 0 is the mean density in the system. With the above equilibrium moments, the sound speed of the lattice is c s = 1/ √ 3. As mentioned, the collision process makes the moments relax towards their equilibrium values as follows:

m c = m -S(m -m eq ) ( 16
)
where superscript c denotes the post-collision state and S is the diagonal relaxation matrix. The post-collision vector f c is then carried out as

f c = f -Q(f -f eq ). ( 17
)
The diagonal matrix S is given by S = diag(1, 1, 1, s e , s ν , s ν , s q , s q , s )

where s i ∈ (0, 2) for the non-conserved moments. The choice of the relaxation parameters s i can be determined by a linear stability analysis [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability[END_REF]. In this work, as suggested by [START_REF] Luo | Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations[END_REF] we will use

s ν = s e = s = 2 6ν + 1
and

s q = 8 2 -s ν 8 -s ν (19)
That leads to the Two-Relaxation Time (TRT) model [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions[END_REF]. Let's note that, based on the product Λ = 1 sν -1 2 1 sq -1 2 , the current choice gives Λ = 3/16 and leads to correct positions of walls when bounce-back boundary conditions are used. Other choices are possible. For example, Λ = 1/4 leads to a better stability of the LBM and provides steady-state solutions dependant only on the equilibrium function [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF][START_REF] Ginzburg | Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium[END_REF]. Values Λ = 1/6 and 1/12 respectively remove high-order diffusion and advection errors. The choice s i = 1/τ BGK recovers the single relaxation time or Bhatnagar-Gross-Krook (BGK) model [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]. Thanks to this MRT formulation, the different models for the collision process can be tested only by changing the values of the relaxation coefficients. The macroscopic fluid variables density ρ and velocity u are obtained from the moments of the distribution functions as follows:

ρ = 9 i=1 f i and ρ 0 u = 9 i=1 e i f i (20)
The corresponding form of the equilibrium distribution is given by the formula [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability[END_REF]:

f eq i = w i ρ 1 + 3(e i • u) + 9 2 (e i • u) 2 - 3 2 u 2 (21)
where ω 1 = 0, ω 2-4 = 1/9 and ω 5-9 = 1/36 are lattice constants. This particular form of the equilibrium distribution function is related to the incompressible formulation of the Navier-Stokes equations [START_REF] He | Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation[END_REF]. This expression will be useful when formulating the boundary conditions.

Multiple relaxation time method for thermal problem

The same way, the energy conservation is modelled by an evolution equation. The distribution functions, noted g, obey to the following equation:

g(x j + e δ t , t + δ t ) -g(x j , t) = -Q(g(x j , t) -g eq (x j , t)) (22)
The corresponding lattice has now five discrete velocities, and reads:

e i = (0, 0) i = 1 (1, 0)e, (0, 1)e, (-1, 0)e, (0, -1)e i = 2 -5 (23)
Like in the previous section, Q represents the collision operator Q = M -1 SM , and the transformation matrix is given by

M =       1 1 1 1 1 0 1 0 -1 0 0 0 1 0 -1 -4 1 1 1 1 0 1 -1 1 -1       . ( 24
)
Like the dynamic part, this matrix is invertible and orthogonal. The temperature θ is the only conserved quantity and can be computed by:

θ = N i=1 g i (25)
The equilibrium moments, m eq corresponding to the distribution functions g, can be written as:

m eq 1 = θ, m eq 2 = uθ, m eq 3 = vθ, m eq 4 = aθ, m eq 5 = 0 ( 26
)
where u = (u, v) is the macroscopic computed velocity and a is a constant. The diagonal relaxation matrix is given by:

S = diag(1, s α , s α , s e , s ν ) (27)
The choice of s i is discussed in details in [START_REF] Wang | Lattice Boltzmann simulation of thermal convective flows in two dimensions[END_REF]. We take here

1 s e - 1 2 = 1 s ν - 1 2 = 1 6 and 1 s α - 1 2 = √ 3 6 (28)
Theses parameters lead to the thermal diffusivity

α = √ 3(4 + a) 60 ( 29
)
where the constant a must be maintained a < 1 in order to avoid numerical instability of the d2q5 model. The value of a will be determined with the physical problem parameters and will be discussed in section 3.5.

The corresponding form of the equilibrium for the distribution functions g is given by the formula

g eq i = ω i θ (1 + 3e i • u) (30)
where ω 1 = 0 and ω 2-5 = 1/4 are d2q5 lattice constants.

External forcing

In a natural convection problem under Boussinesq hypothesis, the external forcing represents the buoyancy effect. In the coordinate system, this force is given by

F (x j , t) = F = f e y with f = -ρ 0 βgθ∆T (31)
The projection of F on the velocity space related to mass and momentum conservation is given by:

F i = -3ω i f e i • e y c 2 (32)
A brief review of force term implementations and their consequences in lattice Boltzmann model can be found in [START_REF] Mohamad | A critical evaluation of force term in lattice Boltzmann method, natural convection problem[END_REF].

Boundary conditions

Two popular kinds of boundary conditions (BCs) are used to impose the macroscopic conditions on velocity and thermal fields described in Sec. 2. The first implementation uses the bounce-back, or anti bounce-back, conditions for both velocity and thermal fields (BB-BCs). With this method, the position of physical walls is half a lattice spacing beyond the last fluid node. The second formulation uses Zou and He [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF] implementation for the velocity field and a Counter-Slip approach for thermal field [START_REF] Inamuro | A non-slip boundary condition for lattice Boltzmann simulations[END_REF][START_REF] D'orazio | Application to natural convection enclosed flows of a lattice Boltzmann BGK model coupled with a general purpose thermal boundary condition[END_REF]. This particular choice is noted ZHCS-BCs in the following. A review of common implementation of boundary conditions for advection-diffusion models can be found in [START_REF] Huang | Numerical study of lattice Boltzmann methods for a convection-diffusion equation coupled with Navier-Stokes equations[END_REF]. Both methods are then briefly described. Two remarks are in order here: 1) With BB-BCs, physical walls are not on nodes but halfway beyond the last fluid node. Then ,the use of a multiple relaxation time model (MRT or TRT) is strongly recommended [START_REF] Wang | Lattice Boltzmann simulation of thermal convective flows in two dimensions[END_REF]. Indeed, with the BGK model, the boundary location depends on the single relaxation time τ BGK . 2) As a consequence, with N fluid nodes, the domain length in lattice units is N when BB-BCs are used and N -1 with ZHCS-BCs.

Bounce-back approach

As mentioned, when the bounce-back scheme is applied, the effective boundaries in the x-direction are in x = 1/2 and x = N x +1/2. No-slip boundary conditions are then realized by the bounce-back scheme. The incoming unknown distribution function f j (x f , t + δ t ) is equal to the outgoing post-collision distribution function f c i (x f , t):

f i (x f , t + δ t ) = f c i (x f , t) (33)
where x f is a fluid node adjacent to a boundary and the incoming distribution function f i corresponds to e i = -e i .

For the temperature field, a bounce-back like scheme can also be applied [START_REF] Wang | Lattice Boltzmann simulation of thermal convective flows in two dimensions[END_REF]. For a wall at a fixed temperature θ w , the following boundary condition is used:

g i (x f , t + δ t ) = -g c i (x f , t) + 2 √ 3 α θ w (34)
An adiabatic wall can be realized with the anti bounce-back condition:

g i (x f , t + δ t ) = g c i (x f , t) (35)

Zou-He and Counter-Slip boundary conditions

For the dynamic part, the no-slip condition for the velocity field is imposed by Zou and He formulation [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF]. This approach relies on the bounce-back rule for non-equilibrium part of the distribution functions. For thermal counterpart, fixed temperature or zero-flux are set with a Counter-Slip approach [START_REF] Orazio | Computational Science -ICC 2003, chapter Boundary Conditions for Thermal Lattice Boltzmann Simulations[END_REF][START_REF] D'orazio | Application to natural convection enclosed flows of a lattice Boltzmann BGK model coupled with a general purpose thermal boundary condition[END_REF]. The incoming unknown distributions are assumed to be at the equilibrium. Both approaches are consistent with the second-order accuracy [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF]. Derivation of boundary conditions is not presented here. The implementation and the particular treatment for corners can be found in the following references [START_REF] Inamuro | A non-slip boundary condition for lattice Boltzmann simulations[END_REF][START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF][START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF][START_REF] Silva | Corner boundary conditions in lattice Boltzmann method[END_REF]. Let us underline that this "on-node" approach allows to impose more sophisticated boundary conditions like pressure boundary conditions for the dynamic part and non-zero fluxes for the thermal part. Furthermore, this formulation is more portable and can also be used with the internal energy formulation method that uses d2q9 model to simulate the thermal effects [START_REF] Peng | Simplified thermal lattice Boltzmann model for incompressible thermal flows[END_REF][START_REF] Dixit | Simulation of high Rayleigh number natural convection in a square cavity using the lattice Boltzmann method[END_REF].

It is important to remind that the choice of boundary conditions defines the size of the computational domain and, as a consequence, the physical space step.

Problem set-up

Setting up a problem in lattice Boltzmann method is not always obvious and some cautions must be taken before each simulation.

The results are generally presented in dimensionless units. Then, we need to convert values obtained on the lattice (in lattice units). The procedure presented here resumes the general procedure given in Ref. [START_REF] Latt | Choice of units in lattice Boltzmann simulations[END_REF]. In dimensionless formulation, conservation of mass, momentum and energy is given by system (1). Using a Chapmann-Enskog procedure [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF], evolution equations ( 7) and ( 22) lead to the following approached system:

           ∇ L • u L = 0 ∂u L ∂t L + u L • ∇ L u L = -∇ L p L + ν L ∇ 2 L u L + f L e y ∂θ L ∂t L + u L • ∇ L θ L = α L ∇ 2 L θ L (36)
In this section, values with subscript L denotes lattice values while classic notation stands for dimensionless ones. The following scalings are used to express system (1) on the lattice

x = δ x x L , t = δ t t L , u = δ x δ t u L , p = ρ 0 δ 2 x δ 2 t p L (37)
The lattice parameters can then be identified and related to dimensionless parameters Ra, P r, δ x and δ t :

ν L = P r Ra δ t δ 2 x , α L = 1 RaP r δ t δ 2 x and f L = δ 2 t δ x f (38)
However, in lattice Boltzmann simulations, it is rather convenient to work with the Mach number Ma defined by Ma = U/c s where U is a reference velocity and c s is the sound speed. Using √ gβ∆T L as a reference velocity, we obtain on the lattice the following relation:

ν L = P r Ra Ma √ 3 1 δ x , α L = ν L P r and δ t = Ma √ 3 δ x (39)
The parameters are now Ra, P r, δ x = 1/N x and the Mach number Ma. The time step varies linearly with Ma as shown in the previous relation. All parameters are then known and it is convenient to chose Ma < 0.3 to make sure the flow stays in the incompressible regime. In practical, with the MRT formulation, the coefficient a is first set through the relation a = 20Ma/δ x √ RaP r -4. In order to avoid numerical instability of the d2q5 model, the value of a must be kept such as a < 1 [START_REF] Ginzburg | Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium[END_REF][START_REF] Wang | Lattice Boltzmann simulation of thermal convective flows in two dimensions[END_REF].

Results and discussion

System (1) is solved with the double MRT thermal lattice Boltzmann formulation for two configurations. In a first part we present a classical validation test case: the square cavity with Ra = 10 6 . In a second part, we study a periodic flow obtained for the 8:1 aspect ratio cavity and Ra = 3.4 × 10 5 . Convergence orders and extrapolated values are then computed by the Richardson extrapolation. In both cases the Prandtl number is P r = 0.71.

Richardson Extrapolation

In order to check the numerical convergence and the order of the numerical scheme, the Richardson Extrapolation (RE) is used. Detailed principles and assumptions can be found in [START_REF] Nicolas | Benchmark Solution for a Three-Dimensional Mixed-Convection Flow, Part 2: Analysis of Richardson Extrapolation in the Presence of a Singularity[END_REF]. RE is constructed with Taylor expansion and then supposes sufficiently smooth exact solution f exact . Small space (or time) step must be small enough so that the discrete solution f hi is in the convergence region. The general principle of RE consists of computing numerical solutions of the discretized problem on N different uniform grids of size h i , with h 1 the coarsest grid and h N the finest one. The solution for the mesh h i is written in the form

f hi = f exact + C α h α i + O(h α+1 i ) ( 40
)
where C α is a coefficient. Thus, by using three grids (N = 3) such as h 1 /h 2 = h 2 /h 3 , the approximations α, Cα and f extra are given by

α = ln[(f h1 -f h2 )/(f h2 -f h3 )] ln(h 1 /h 2 ) Cα = f h2 -f h3 h α 2 -h α 3 f extra = f h3 -Cα h α 3 (41)
with Cα = C α + O(h N -1 ) and f extra = f exact + Cα h α+1 N . As a consequence, the approximation f extra will be closer to the asymptotic solution f exact with decreasing h N and increasing α. In practical, the use of RE under good conditions increases of one order the numerical solution. In this paper, RE is used with N = 3 and N = 4 such as h 1 /h 2 = h 2 /h 3 and h 1 /h 2 = h 3 /h 4 respectively (see Ref. [START_REF] Nicolas | Benchmark Solution for a Three-Dimensional Mixed-Convection Flow, Part 2: Analysis of Richardson Extrapolation in the Presence of a Singularity[END_REF] for the N = 4 case). Since RE is based on Taylor expansion, the same principle can be applied to temporal convergence.

Steady flow

In this section, the aspect ratio is fixed to the unity. The Rayleigh number is Ra = 10 6 . This configuration has been extensively studied with various numerical methods [START_REF] De Vahl | Natural convection in air square cavity: a benchmark numerical solution[END_REF][START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high rayleigh number[END_REF][START_REF] Guo | A coupled lattice bgk model for the boussinesq equations[END_REF]. Here we will use the spectral results of Le Quéré [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high rayleigh number[END_REF] as a comparison for the quantities of interest. To our knowledge, Le Quéré provides the most accurate values in the literature. For this benchmark, we need to compute the following quantities: the maximum horizontal velocity at mid-height u max = max(u(x, y = A/2)), the maximum vertical velocity at mid-width v max = max(v(x = A/2, y)) and their respective positions y max and x max . These values are then scaled with the diffusive velocity α/W . The Nusselt number Nu and the averaged Nusselt number Nu are calculated at the hot wall:

Nu(y) = - 1 θ h -θ c ∂θ ∂x wall (42)
and

Nu = 1 H H 0 Nu(y) dy (43)
Integral calculations are performed with a second-order accuracy Newton-Cotes quadrature (trapezoid or mid-point rule).

Positions of maximum and minimum Nusselt values at the hot wall are also estimated. The steady flow is assumed converged when the following criteria is reached:

9 i=1 Nx×Ny j=1 ||f i (x j , t + 100 δ t ) -f i (x j , t)|| 1 9 i=1 Nx×Ny j=1 ||f i (x j , t + 100 δ t )|| 1 ≤ 10 -9 (44)
where || • || 1 denotes the L 1 norm. The Mach number is set to Ma = 5 × 10 -2 . It has been shown [START_REF] Wang | Lattice Boltzmann simulation of thermal convective flows in two dimensions[END_REF] that there is no influence of the Mach number in this stationary problem. Tables 1 and2 shows the numerical results for BB-BCs and ZHCS-BCs respectively. At first, using finer grids makes the values converge through the expected ones. On the finer grid, N x × N y = 301 2 , the maximum difference is less than 2% for the BB-BCs and less than 1% for ZHCS-BCs. However, this difference is not significant. We can notice that some quantities do not converge through Le Quéré ones. For example, the extrapolated value x of the maximum vertical velocity is slightly greater than the reference. Thanks to the RE, extrapolated values and convergence order are computed. The convergence orders are quite acceptable for both boundary conditions sets. The second-order convergence (here with space step) is obtained for most of the quantities of interest. The use of three and four points shows that the convergence region seems to be larger for the ZHCS-BCs formulation.

N 2: Convergence of values with P r = 0.71, Ra = 10 6 and Ma = 5 × 10 -2 for ZHCS-BCs. The reference values of Le Quéré [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high rayleigh number[END_REF] are also presented as the results of RE on three (N = 3) and four (N = 4) meshes. Values in parenthesis are relative differences with reference or convergence order.

We would like to underline that the use of a single time relaxation model for the temperature (BGK) does not change results in a significant way.

To sum up, these observations are in good agreement with references. Both boundary conditions approaches lead to the expected results. The relative difference between the two implementations is less than 0.5%. Furthermore, the correct convergence orders validate our implementation and allows us to simulate a time-dependent natural convection case.

Periodic flow

In this section, the 8:1 aspect ratio differentially heated cavity is simulated. As presented in [START_REF] Xin | An extended Chebychev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF][START_REF] Christon | Computational predictability of time-dependent natural convection flows in enclosures (including a benchmark solution)[END_REF], with the Rayleigh number Ra = 3.4 × 10 5 , the flow admits only one time-dependent skew-symmetric solution. This configuration has been studied by multiple authors and stands for a benchmark for time-dependent natural convection. A synthesis can be found in [START_REF] Christon | Computational predictability of time-dependent natural convection flows in enclosures (including a benchmark solution)[END_REF]. Two probes are set in the cavity. Their respective coordinates are (x 1 , y 1 ) = (0.1810, 7.3700) and (x 2 , y 2 ) = (0.8190, 7.3700).

The quantities under study are mean values of horizontal velocity u and temperature θ at probe position 1, mean value of pressure difference between points 1 and 2 ∆p 12 , mean value of the Nusselt number Nu at hot wall and their respective perturbations u , θ , ∆p 12 , Nu . The pressure is obtained through the relation p = ρc s . For a quantity φ, the mean value φ and the oscillatory component φ are computed as:

φ = 1 τ t+τ t φ(x, t) dt (45)
and

φ = max(φ(x, t)) -min(φ(x, t)) (46)
where τ represents the period (based on temperature evolution at probe 1) for which the average is computed.

The dimensionless temperature evolution is given in figure 2 while iso-lines of the temperature field over one period are presented in figure 3. The temperature field is stratified in the center while oscillations are visible at the top and the bottom of the cavity. At first sight, a single frequency time-dependent flow is observed. The periodic state is reached for dimensionless times greater than 600. Quantities of interest are extracted during a period τ for times greater than 1000. Tables 3 and4 present the evolution quantities of interest for different meshes with Ma = 0.1. The reference values of Xin and Le Quéré [START_REF] Xin | An extended Chebychev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF] for this problem are also given as the results of the RE on three and four meshes. First, the mean values approach the reference for both boundary conditions sets except the pressure difference ∆p 12 in the case of BB-BCs (Tab. 3). Unlike most quantities, the relative difference on ∆p 12 obtained on coarse grids is important (78%). Even for the finest grids, it remains greater than 20%. For this particular variable, the RE produces nearly a first accuracy space convergence order ( α = 0.74 and 0.82 respectively with three and four meshes). With ZHCS-BBs (Tab. 4), Finally, the convergence orders are discussed (Tab. 5). With space step, RE shows that theoretical predicted orders are better recovered with ZHCS-BCs approach. Like in the pressure field case, space convergence orders for the mean velocity u 1 and temperature θ 1 are close to the unity with the BB-BCs approach. With ZHCS-BBs, the second-order accuracy is obtained for most quantities, mean value and oscillatory component. An explanation could be that when bounce-back conditions are used, the grid is not totally regular: the wall is half a space-step beyond the last fluid node. Hence, the condition of regular grid for the RE is not fulfilled when BB-BCs are applied. With time step, analysis on mean values provides a second-order accuracy except for the Nusselt number. We suppose this quantity outside the convergence region. The theoretical order is also recovered for oscillatory data. It is important to notice that the choice of the boundary condition formulation does not change the temporal convergence order. We notice that simulations were also performed with the simplified d2q9 internal energy thermal model for temperature introduced in [START_REF] Peng | Simplified thermal lattice Boltzmann model for incompressible thermal flows[END_REF]. The results on finer meshes are quite close to those obtained with the double MRT approach. But, as it is based on a single relaxation time, simulations are unstable for small relaxation rates, typically less than 0.55, and a rigorous convergence study is not possible as it requires too important meshes. We should emphasize that this problem is very expensive to simulate and a parallel implementation is highly recommended. As mentioned in Sec. 3.5, the time step varies linearly with both space step and Mach number. For guidance, time calculations are given in table 5 as a function of the Mach number. Calculation times could probably be reduced using a mesh refinement algorithm developed for non uniform grids [START_REF] He | Some Progress in Lattice Boltzmann Method. Part I. Nonuniform Mesh Grids[END_REF][START_REF] Shu | Simulation of natural convection in square cavity by Taylor series expansion-and least squares-based lattice Boltzmann method[END_REF]. In this paper, regular meshes are uses as required for the Richardson extrapolation. [START_REF] Xin | An extended Chebychev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF] are also presented as the results of RE on three (N = 3) meshes. Values in parenthesis are relative differences with reference or convergence order. For guidance, the total cpu time in hours on AMD 6276 computers is given.

N x × N y u

Ma u

To summarize this section, predicted results and convergence orders were obtained for a periodic thermal flow. As expected, the Mach number is a crucial parameter in time-dependent problem. When the pressure field is a relevant quantity, ZHCS-BCs approach gives better results and should be used.

Conclusion

In this paper, the double multiple-relaxation-time collision model for the lattice Boltzmann equation has been implemented to simulate two natural convection problems under the Boussinesq assumption. Two kinds of popular boundary conditions have been tested. The bounce-back approach where walls are halfway beyond the last fluid node, and a "on-node" approach constructed with Zou and He and Counter-Slip formulations. The Richardson extrapolation has systematically been used to compute an extrapolated solution and convergence orders from the numerical results. The results are in good agreement with reference benchmark studies. The double MRT method is numerically more stable than double BGK or MRT-BGK formulations and thus allows rigorous convergence studies. With the finest space and time steps, less than 1% relative difference error is observed for most of the quantities of interest. The non-negligible error made on the pressure field and the small convergence order obtained with bounce-back conditions are corrected with the "on-node" formulation. The latter provides a better framework to study numerical properties. It should be used when the pressure field is a relevant quantity. The time-dependent flow study highlights the significant effect of the Mach number, or time step, on the oscillatory data. The expected second-order time accuracy has been found for both mean and oscillatory values and makes the multiple-relaxation-time thermal lattice Boltzmann model suitable for time-dependent thermal flows.

Figure 1 :

 1 Figure 1: Differentially heated cavity with insulated horizontal walls and constant temperature vertical walls. The aspect ratio is A = H/W . Points 1 and 2, respectively in (x 1 , y 1 ) and (x 2 , y 2 ), are time history points. The dimensionless formulation of incompressible Navier-Stokes and energy equation coupled with Boussinesq hypothesis for

Figure 2 :Figure 3 :

 23 Figure 2: Dimensionless temperature θ as a function of time at probe (x 1 , y 1 ) = (0.1810, 7.3700)

Table 1 :

 1 x × N y u max Convergence of values with P r = 0.71, Ra = 10 6 and Ma = 5 × 10 -2 for BB-BCs. The reference values of Le Quéré[START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high rayleigh number[END_REF] are also presented as the results of RE on three (N = 3) and four (N = 4) meshes. Values in parenthesis are relative differences with reference or convergence order.

		y	v max	x
	75 2 64.7722(0.10%) 0.8484(0.18%)	220.907(0.16%) 0.03791(0.24%)
	100 2 64.7999(0.05%) 0.8491(0.10%)	220.772(0.10%) 0.03786(0.37%)
	150 2 64.8216(0.02%) 0.8496(0.04%)	220.658(0.04%) 0.03780(0.53%)
	200 2 64.8279(0.01%) 0.8498(0.02%)	220.625(0.03%) 0.03779(0.55%)
	300 2 64.8264(0.01%) 0.8503(0.05%)	220.620(0.03%) 0.03777(0.61%)
	Ref [27] 64.8344	0.8500	220.559	0.03800
	f extra N =3 (α) 64.8269(3.36) f extra N =4 (α) 64.8361(1.97)	0.8517(0.65) 0.8500(1.94)	220.613(2.71) 220.577(1.83)	0.03792(1.88) 0.03773(1.60)
	N x × N y Nu	N u max	y	N u min
	75 2 8.9612(1.54%)	19.5679(11.6%) 0.0473(21.3%)	0.81666(16.6%)
	100 2 8.8980(0.83%)	18.6947(6.61%) 0.0428(9.78%)	0.90300(7.81%)
	150 2 8.8544(0.33%)	18.0240(2.78%) 0.0398(2.05%)	0.96117(1.87%)
	200 2 8.8403(0.17%)	17.7955(1.48%) 0.0391(0.32%)	0.97732(0.22%)
	300 2 8.8310(0.07%)	17.6419(0.60%) 0.0383(1.54%)	0.98440(0.50%)
	Ref [27] 8.8252	17.5360	0.0390	0.97946
	f extra N =3 (α) 8.8245(2.19) f extra N =4 (α) 8.8237(2.13)	17.5162(2.01) 17.4770(1.87)	0.03807(2.42) 0.03847(2.47)	0.9888(2.63) 0.9944(2.31)
	N x × N y u max	y	v max	x
	76 2 64.2827(0.85%) 0.8513(0.16%)	218.983(0.71%) 0.03893(2.46%)
	101 2 64.5255(0.48%) 0.8507(0.08%)	219.734(0.37%) 0.03842(1.12%)
	151 2 64.6988(0.21%) 0.8503(0.04%)	220.206(0.16%) 0.03807(0.18%)
	201 2 64.7601(0.11%) 0.8501(0.02%)	220.372(0.08%) 0.03793(0.18%)
	301 2 64.8010(0.05%) 0.8500(0.00%)	220.482(0.03%) 0.03778(0.57%)
	Ref [27] 64.8344	0.8500	220.559	0.03800
	f extra N =3 (α) 64.8342(2.02) f extra N =4 (α) 64.8393(1.99)	0.8498(1.62) 0.8500(2.21)	220.562(2.14) 220.553(2.26)	0.03764(1.61) 0.03774(1.95)
	N x × N y Nu	N u max	y	N u min
	76 2 8.8362(0.13%)	18.0556(2.96%) 0.0460(17.5%)	0.82666(15.6%)
	101 2 8.8244(0.01%)	17.8409(1.74%) 0.0426(9.37%)	0.90280(7.83%)
	151 2 8.8161(0.10%)	17.6769(0.80%) 0.0403(3.36%)	0.96224(1.76%)
	201 2 8.8139(0.13%)	17.6226(0.49%) 0.0395(1.49%)	0.97792(0.16%)
	301 2 8.8128(0.14%)	17.5853(0.28%) 0.0393(0.91%)	0.98542(0.61%)
	Ref [27] 8.8252	17.5360	0.0390	0.97946
	f extra N =3 (α) 8.8121(2.59) f extra N =4 (α) 8.8114(2.25)	17.5562(2.04) 17.5488(1.91)	0.0391(2.57) 0.0387(2.12)	0.99019(2.54) 0.99742(2.05)
	Table

Table 3 :

 3 Convergence of values with P r = 0.71, Ra = 3.4 × 10 5 and Ma = 10 -1 for BB-BCs. The reference values of Xin and Le Quéré[START_REF] Xin | An extended Chebychev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF] are also presented as the results of RE on three (N = 3) and four (N = 4) meshes. Values in parenthesis are relative differences with reference or convergence order.

	1	θ 1		Nu		∆p 12	τ
	75 × 600 0.052258(7.27%) 0.265875(0.15%) 4.594463(0.33%) 0.000394(78.7%) 3.433715(0.65%)
	100 × 800 0.053549(4.98%) 0.265787(0.12%) 4.587082(0.17%) 0.000709(61.7%) 3.425231(0.40%)
	150 × 1200 0.054752(2.85%) 0.265701(0.08%) 4.582309(0.06%) 0.001065(42.4%) 3.418987(0.22%)
	200 × 1600 0.055291(1.89%) 0.265651(0.06%) 4.580837(0.03%) 0.001251(32.4%) 3.417280(0.17%)
	300 × 2400 0.055794(1.00%) 0.265598(0.04%) 4.579920(0.01%) 0.001443(22.0%) 3.415954(0.13%)
	Ref [43] 0.056356	0.265480		4.579460		0.001850	3.41150
	f extra N =3 (α) 0.056541(1.25) f extra N =4 (α) 0.056539(1.24)	0.265448(0.75) 0.265472(0.85)	4.579335(2.34) 4.579281(2.31)	0.001930(0.82) 0.002032(0.74)	3.415167(2.27) 3.415277(2.14)
	N x × N y u 1		θ 1		Nu	∆p 12
	75 × 600 0.050830(7.29%) 0.039788(6.91%) 0.006336(10.8%) 0.018641(8.53%)
	100 × 800 0.054594(0.43%) 0.042632(0.25%) 0.006860(3.38%) 0.020011(1.81%)
	150 × 1200 0.056746(3.50%) 0.044221(3.47%) 0.007173(1.03%) 0.020817(2.14%)
	200 × 1600 0.057327(4.56%) 0.044617(4.39%) 0.007261(2.27%) 0.021052(3.30%)
	300 × 2400 0.057572(5.00%) 0.044753(4.71%) 0.007307(2.92%) 0.021175(3.90%)
	Ref [43] 0.054828		0.042740		0.007100	0.020380
	f extra N =3 ( α) 0.0577060(2.84) f extra N =4 ( α) 0.0578259(2.68)	0.0448256(3.05) 0.0449381(2.79)	0.0073325(2.64) 0.0073419(2.55)	0.0212455(2.60) 0.0212675(2.56)
	N x × N y u 1	θ 1		Nu		∆p 12	τ
	76 × 601 0.057049(1.23%) 0.265821(0.13%) 4.591869(0.27%) 0.001621(12.3%) 3.441009(0.86%)
	101 × 801 0.056936(1.03%) 0.265701(0.08%) 4.586808(0.16%) 0.001717(7.19%) 3.430037(0.54%)
	151 × 1201 0.056838(0.86%) 0.265621(0.05%) 4.583163(0.08%) 0.001785(3.51%) 3.421374(0.29%)
	201 × 1601 0.056818(0.82%) 0.265586(0.04%) 4.582049(0.06%) 0.001806(2.38%) 3.418490(0.20%)
	301 × 2401 0.056794(0.78%) 0.265557(0.03%) 4.581215(0.04%) 0.001823(1.46%) 3.416565(0.15%)
	Ref [43] 0.056356	0.265480		4.579460		0.001850	3.411500
	f extra N =3 (α) 0.056782(2.26) f extra N =4 (α) 0.056793(2.02)	0.265527(1.64) 0.265539(1.92)	4.580653(2.16) 4.580706(2.10)	0.001834(2.11) 0.001831(2.12)	3.415005(2.03) 3.414373(1.85)
	N x × N y u 1		θ 1		Nu	∆p 12
	76 × 601 0.051456(6.15%) 0.039804(6.87%) 0.006324(10.9%) 0.018610(8.68%)
	101 × 801 0.054653(0.32%) 0.042257(1.13%) 0.006793(4.32%) 0.020167(1.05%)
	151 × 1201 0.056912(3.80%) 0.044062(3.09%) 0.007132(0.45%) 0.020983(2.96%)
	201 × 1601 0.057483(4.84%) 0.044520(4.16%) 0.007232(1.86%) 0.021193(3.99%)
	301 × 2401 0.057689(5.22%) 0.044698(4.58%) 0.007289(2.66%) 0.021279(4.41%)
	Ref [43] 0.054828		0.042740		0.007100	0.020380
	f extra N =3 ( α) 0.057818(2.81) f extra N =4 ( α) 0.058098(2.28)	0.044810(2.74) 0.045040(2.20)	0.007327(2.36) 0.007351(2.12)	0.021321(3.00) 0.021353(2.91)

Table 4 :

 4 Convergence of values with P r = 0.71, Ra = 3.4 × 10 5 and Ma = 10 -1 for ZHCS-BCs. The reference values of Xin and Le Quéré[START_REF] Xin | An extended Chebychev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF] are also presented as the results of RE on three (N = 3) and four (N = 4) meshes. Values in parenthesis are relative differences with reference or convergence order. the difference is smaller and acceptable. It decreases form 12% to 1.5% and a second order accuracy convergence is obtained (α 2.1). Secondly, the oscillatory values do not converge through the reference and present a non negligible relative difference. With the finest meshes, almost a 5% difference is obtained. As shown in table5for BB-BCs, this behaviour is linked to the Mach number (or the time step). The decreasing time step makes the values converge and, with the smallest Mach number Ma = 0.05, their relative difference is about 2%. On the other side, the correct mean values are obtained even for large Mach number. This observation is consistent with the Mach number independence for stationary data.

Table 5 :

 5 Convergence of values with P r = 0.71, Ra = 3.4 × 10 5 and N x × N y = 300 × 2400 for BB-BCs. The reference values of Xin and Le Quéré

	1	θ 1	Nu	∆p 12	τ
	0.2 0.057180(1.46%) 0.265746(0.10%) 4.579715(0.01%) 0.001384(25.1%) 3.426895(0.45%)
	0.1 0.055794(1.00%) 0.265598(0.04%) 4.579920(0.01%) 0.001443(22.0%) 3.415954(0.13%)
	0.05 0.055468(1.58%) 0.265558(0.03%) 4.580093(0.01%) 0.001456(21.3%) 3.413460(0.06%)
	Ref [43] 0.056356	0.265480	4.579460	0.001850	3.411500
	f extra N =3 (α) 0.055368(2.09)	0.265543(1.89)	4.581028(0.24)	0.001460(2.18)	3.412724(2.13)
	Ma u 1	θ 1	Nu	∆p 12	Cpu time (h)
	0.2 0.066978(22.1%) 0.050714(18.6%) 0.007989(12.5%) 0.023792(16.7%) 1700
	0.1 0.057572(5.00%) 0.044753(4.71%) 0.007307(2.92%) 0.021175(3.90%) 3350
	0.05 0.055209(2.04%) 0.043100(0.84%) 0.007114(0.20%) 0.020418(0.19%) 7480
	Ref [43] 0.054828	0.042740	0.007100	0.020380	-
	f extra N =3 ( α) 0.054416(1.99)	0.042466(1.85)	0.007004(1.85)	0.020083(1.79)	-