R Fourneaux

M Vermalle

F Albarel

I Morange

T Graillon

V Amodru

T Cuny

H Dufour

T Brue

MD, Ph. Frederic Castinetti
email: frederic.castinetti@ap-hm.fr

ACROMEGALY IN REMISSION : A VIEW FROM THE PARTNER 1

Keywords: Acromegaly, Surgery, Quality of life, Pituitary, Relative, Patient Related Outcome

Purpose A relative can be an asset in dealing with chronic illnesses such as acromegaly where quality of life is altered even after remission. However, it has been shown that quality of life of caregivers can also be impacted. Our main objective was to explore the consequences of acromegaly in remission in the patient-relative dyad in a matter of quality of life and self-esteem.

Methods

In this observational study, to better characterize the perception of the disease by the partner, patient's body image and self-esteem were evaluated from the patient's point of view (n=27) and from the relative's using the same questionnaires with modified instructions. The patient and the partner were also asked to fulfil quality of life, anxiety/depression and coping strategies questionnaires.

Results

The relative had an overall accurate estimation of the patient's body image using Stunkard figurines. However, there were wide variations between the patient's and the relative's answers for self-esteem and body perception. The relative's quality of life was not altered and was significantly higher in the social domain than for the patient.

Conclusion Despite inter-individual changes between the patient's and the relative's view, our results

show that the relative should be educated in all the steps of the management of acromegaly, to help him/her better understand the disease and support the patient.

INTRODUCTION

In chronic diseases, health-related quality of life has become an increasingly important aspect of patient care [START_REF] Webb | Patientcentred outcomes with pituitary and parasellar disease[END_REF]. For instance, patients with active as well as controlled acromegaly report an impaired quality of life, linked to the duration of the disease, the necessity of a prolonged medical treatment, a history of radiotherapy or persistent symptoms such as joints complaints [START_REF] Andela | Quality of life (QoL) impairments in patients with a pituitary adenoma: a systematic review of QoL studies[END_REF][START_REF] Biermasz | Morbidity after long-term remission for acromegaly: persisting joint-related complaints cause reduced quality of life[END_REF][START_REF] Andela | More concerns and stronger beliefs about the necessity of medication in patients with acromegaly are associated with negative illness perceptions and impairment in quality of life[END_REF]. The partner (family member or close friend) can also be impacted by the chronic disease of the patient [START_REF] Weitzner | The impact of pituitary disease on the family caregiver and overall family functioning[END_REF]. Indeed, it has been thoroughly described with dementia and cancer patients [START_REF] Farina | Factors associated with the quality of life of family carers of people with dementia: A systematic review[END_REF][START_REF] Stenberg | Review of the literature on the effects of caring for a patient with cancer[END_REF][START_REF] Tan | Burden, emotional distress and quality of life among informal caregivers of lung cancer patients: An exploratory study[END_REF]. Furthermore, when evaluated with various generic scales such as HAD and SF-36, chronic kidney disease patients and relatives showed the same level of anxiety and depression [START_REF] Pereira B Dos | Beyond quality of life: a cross sectional study on the mental health of patients with chronic kidney disease undergoing dialysis and their caregivers[END_REF].

To the best of our knowledge, the patient-relative dyad in pituitary disease has only been studied in two qualitative studies [START_REF] Dunning | Pituitary disease -perspectives of patients and partners[END_REF] [START_REF] Andela | The partner's perspective of the impact of pituitary disease: Looking beyond the patient[END_REF]. For instance, Andela et al. conducted focused groups interviews to explore the partners' perspective. They found that the partners had negative beliefs about medication and felt that they had to make adaptations of their own behaviour for their ill partner [START_REF] Andela | The partner's perspective of the impact of pituitary disease: Looking beyond the patient[END_REF]. In both studies, partners asked for more information about the disease and additional guidance.

However, consequences on the dyad in acromegaly has never been specifically studied from a quantitative point of view.

The main objective of this study was to determine the partner's perception of acromegaly in remission via an original approach in which we asked him/her to imagine what the patient had answered in terms of body image and self-esteem. Secondary objectives were to determine the quality of life of the partner, as well as anxiety, and depression (comparing it with the patient's characteristics), and the coping strategies of the dyad.

PATIENTS AND METHODS

A prospective monocentric, non-interventional study was conducted in the Tertiary reference Center of La Conception Hospital, Assistance Publique Hôpitaux de Marseille, France.

Patients and partner/family member

Patients aged 18-80 years, in remission and/or controlled by a pharmacological treatment for 1-10 years, were recruited while doing their follow-up at our department between September 2019 and June 2020. Remission was defined as follows: Normal age and sex matched -insulin-like growth factor (IGF-1) levels, random growth hormone (GH) < 1 ng/mL and/or a nadir GH level after OGTT < 0,4 ng/mL. Patients were considered controlled if they had normal IGF-1 levels and random GH < 1 ng/mL on medical treatment. Patients with known cognitive deficiency could not be included. The partner and/or family member was chosen by the patient.

Data collection and conduct of the study

After being selected by the endocrinologist in charge of each patient in the Department, patients were contacted by mail. Informed written consent was obtained from all patients to participate in the study, which was approved by the ethics committee of Aix Marseille University. The following data were collected for each patient: socio-demographic data (including age, gender, education level, profession, marital status), medical data (medical history, comorbidities), acromegaly data (therapeutic course, current treatment, hormonal deficits).

Patients and relatives were asked to fill in the following questionnaires:

-Quality of life assessment using : First, the French version of the generic scale WhoQoL-BREF developed in 1998 [START_REF]Development of the World Health Organization WHOQOL-BREF quality of life assessment. The WHOQOL Group[END_REF] and translated in 2010 by Baumann et al. [START_REF] Baumann | The WHOQOL-BREF questionnaire: French adult population norms for the physical health, psychological health and social relationship dimensions[END_REF]. It allows the evaluation of quality of life in several areas like "physical health", "psychological health", "social relations" and "environment". Secondly, the only disease-specific scale in acromegaly called AcroQoL [START_REF] Webb | Acromegaly Quality of Life Questionnaire (ACROQOL) a new health-related quality of life questionnaire for patients with acromegaly: development and psychometric properties[END_REF] for patients only.

-Anxiety and depression evaluation with the Hospital Anxiety and Depression scale (HAD) developed by Zigmond et al. in 1983 [15] and validated in French by Lepine et al. in 1985 [16]. This scale has the particularity of studying both the anxiety and the depression components. It includes 14 items rated from 0 to 3. Seven questions refer to the anxiety component and seven others to the depressive component, giving two scores.

-Self-esteem was evaluated with Rosenberg's questionnaire, developed in 1965 [17] and translated in 1992 [START_REF] Vallieres | TRADUCTION ET VALIDATION CANADIENNE-FRANÇAISE DE L'ÉCHELLE DE L'ESTIME DE SOI DE ROSENBERG*[END_REF].

-Body image perception using the Image Body Questionnaire (IBQ), a French questionnaire developed by Bruchon-Schweitzer in 1990 [START_REF] Koleck | The body-image questionnaire: an extension[END_REF] and figurines by Stunkard, Sorensen and Schlusinger (patients were asked to choose the figurine that they believed to be closest to their representation of themselves) [START_REF] Stunkard | Use of the Danish Adoption Register for the study of obesity and thinness[END_REF].

Of note, we studied the relative's perception of the patient's body image using the same questionnaires (Rosenberg's questionnaire, IBQ and figurines by Stunkard, Sorensen and Schlusinger) with modified instructions, and without knowing the patient's answers. This technique has already been used in the literature to study the perception of the patient's quality of life by the relative in chronic diseases [START_REF] Gundy | The influence of proxy perspective on patient-proxy agreement in the evaluation of health-related quality of life: an empirical study[END_REF], particularly in cancer and neurological diseases like strokes where it can be useful since the patient's communication is impaired [START_REF] Sneeuw | Assessing quality of life after stroke. The value and limitations of proxy ratings[END_REF] [START_REF] Sneeuw | The use of significant others as proxy raters of the quality of life of patients with brain cancer[END_REF].

For the patient-relative dyad:

We also explored the coping strategies, defined as cognitive and behavioural efforts to deal with a stressful situation, using the Brief Cope situational scale. The version we used was translated and validated in French by Muller et al. in 2003 [24]. It includes 14 different coping dimensions (active coping, planning, using instrumental support, using emotional support, venting, behavioural disengagement, self-distraction, self-blame, positive reframing, humour, denial, acceptance, religion and substance use). We used a French version regrouping these dimensions in 4 coping strategies :

Seeking social support (including emotional support, instrumental support, venting and religion), positive thinking (including humour, positive reframing and acceptance), avoidance (including behavioural disengagement, self-distraction, substance use, denial and self-blame), problem solving (including active coping and planning) [START_REF] Baumstarck | Assessment of coping: a new french four-factor structure of the brief COPE inventory[END_REF].

Statistical analysis

We first performed a descriptive analysis of the population. Quantitative variables were expressed as median with interquartile ranges or means with standard deviation; qualitative variables were expressed as proportions and percentages. We considered each Qol score as a variable of interest: the four Qol scores of the WHOQoL-BREF (physical, psychological, social and environmental scores) and the total score of the AcroQoL. The same analysis was performed for the partner. Body image questionnaire, Stunkard figurines score and Rosenberg's questionnaire evaluations made by the partner were correlated with the evaluations of the patients. Data analyses were performed using Prism version 8.4.3 (GraphPad software, LLC). p-values < 0.05 were considered significant and all statistical tests were two-tailed.

RESULTS

Thirty-five patients were contacted; twenty-seven patients (15 women and 12 men) agreed to participate, and were enrolled in this study. The epidemiological and socio-cognitive characteristics of the patients are reported in table 1. Most patients had familial support, and were educated to a middle to high level. Three patients (n=3/27, 11.1%) had a disabled worker status. Regarding the partner and/or family member, 77.7% of the relatives were partners, 22.2% were children.

We asked the relative to determine the self-esteem score that the patient would obtain by filling in the same questionnaire as he/she was the patient. Eight relatives considered the patient to have a very low or low self-esteem (vs. 6 when completed by the patient), while 13 (vs. 12 when fulfilled by the patient) considered the patient to have a high or very high self-esteem (p=0.348). We asked every relative to depict the Stunkard body that would be given by the patient: the median estimated Stunkard body score was 5 [START_REF] Webb | Patientcentred outcomes with pituitary and parasellar disease[END_REF][START_REF] Andela | Quality of life (QoL) impairments in patients with a pituitary adenoma: a systematic review of QoL studies[END_REF][START_REF] Biermasz | Morbidity after long-term remission for acromegaly: persisting joint-related complaints cause reduced quality of life[END_REF][START_REF] Andela | More concerns and stronger beliefs about the necessity of medication in patients with acromegaly are associated with negative illness perceptions and impairment in quality of life[END_REF][START_REF] Weitzner | The impact of pituitary disease on the family caregiver and overall family functioning[END_REF][START_REF] Farina | Factors associated with the quality of life of family carers of people with dementia: A systematic review[END_REF][START_REF] Stenberg | Review of the literature on the effects of caring for a patient with cancer[END_REF][START_REF] Tan | Burden, emotional distress and quality of life among informal caregivers of lung cancer patients: An exploratory study[END_REF], which was not significantly different compared to the patient's answers (p= 0.171 in comparison with the patient). The same non-significant difference was observed for the IBQ (p=0.549). However, as shown in Figure 1, though the overall differences in these parameters were not significant between the patient's and the relative's view, there was a wide variation when taken from an individual viewpoint: except for Stunkard scores (r=0.792, p<0.001), both Rosenberg and IBQ scores were indeed not correlated between the patient and the relative (r=0.203 and r=0.080, respectively).

Regarding the socio-cognitive characteristics of the partner, the quality of life median scores of the WhoQoL-BREF were 69 [31-100] for physical, 69 for psychological, 75 [44-88] for social relationships and 75 [44-88] for environmental. When comparing with patients' Who-QoL (data not shown), there was no significant difference for all but one parameter, social relationships, for which there was a median gain of 19 points in the relative questionnaire (p=0.0435). Two family members (7.4%) had a HAD-anxiety score above 8 (likely anxious), while 1 (3.7%) had a score above 11 (certain anxiety). According to the HAD depression scale, 1 family member (3.7%) had a depression score above 8 (likely depressive), and 1 (3.7%), a score above 11 (certain depressive). Regarding these scores, we found no significant difference between the patient and the family member (p=0.247).

Finally, the adaptation strategy preferentially adopted by the dyad was positive thinking. Problem solving came in second position, and seeking social support in third.

DISCUSSION

Illness perceptions can be different between a patient suffering from a chronic disease and his/her relative. In acromegaly, we showed that the patient's relative had a global accurate estimation of the patient's body image using the Stunkard scale. However, even though there was no statistical difference between the patient's IBQ and Rosenberg scores and the partner's, there were major differences when taken individually as shown in Figure 1. To our knowledge, this facet of the dyad has never been evaluated in our field. Using the same method for advanced cancer patients, caregivers seem to give an adequate evaluation of the patient's body dissatisfaction and associate it to the weight loss due to their disease [START_REF] Rhondali | Association between body image dissatisfaction and weight loss among patients with advanced cancer and their caregivers: a preliminary report[END_REF]. Body image is only mentioned by the patient in the few studies available about the patient-relative dyad in pituitary disease. Therefore, we speculate that family members of acromegaly patients are not fully aware of the consequences of body changes on their relative's self-esteem. As a matter of fact, in the focus groups conducted by Andela et al., partners expressed viewing their relative differently but for other reasons such as changes in their relationship (partner becoming the counsellor or carer) or differences in coping mechanisms. Both patients and partners reported difficulties in communicating about the disease with doctors and their social network, resulting in a decrease of social interactions [START_REF] Andela | The partner's perspective of the impact of pituitary disease: Looking beyond the patient[END_REF]. Accordingly, seeking social support only came in third place on our results of the Brief cope. This may be explained by the fact that acromegaly is a rare disease, little known by the medical community and general population.

Dunning et al. explored experiences of patients with pituitary disease and their partners by monitoring chat room discussions and they observed that patients discussed more openly about physical changes, appearance and feelings among their peers than with their family or friends [START_REF] Dunning | Pituitary disease -perspectives of patients and partners[END_REF].

Moreover, in an online survey conducted on the Carenity website, patients expressed the need to share experiences with other patients in discussion groups and that acromegaly be better known from their relatives [START_REF] Albarel | Daily life, needs and expectations of patients with acromegaly in France: An on-line survey[END_REF]. Taken together, this suggests that we should use body image as a starting point to discuss the particularity of acromegaly either in consultation or in education programs [START_REF] Albarel | Evaluation of an individualized education program in pituitary diseases: a pilot study[END_REF].

Relatives of acromegaly patients do not seem to have an altered quality of life. We also found a difference of 19 points between the patient's and the relative's social WhoQoL score which contrasts with previous findings [START_REF] Weitzner | The impact of pituitary disease on the family caregiver and overall family functioning[END_REF]. However, it should be interpreted with caution considering the high variability of answers and the size of our cohort. Indeed, acromegaly is a rare disease and finding a family member willing to participate reduces the number of participants. A review of the literature published by Stenberg et al. identified difficulties experienced by family caregivers of cancer patients like fatigue, sleep disturbance, anxiety or mood disturbances highlighting the burden of caregiving responsibilities [START_REF] Stenberg | Review of the literature on the effects of caring for a patient with cancer[END_REF]. Using specific scales like the Caregiver Burden Scale (CBS) and the Caregivers Quality of Life Index-Cancer (CQOLC), another study found that caregiver burden increased with anxiety, depression and poorer quality of life and that these results were similar to those reported in other chronic conditions like chronic heart failure [START_REF] Tan | Burden, emotional distress and quality of life among informal caregivers of lung cancer patients: An exploratory study[END_REF]. Furthermore, some findings in chronic kidney disease show that patients and caregivers have similar scores in the mental and social domains of quality of life score [START_REF] Pereira B Dos | Beyond quality of life: a cross sectional study on the mental health of patients with chronic kidney disease undergoing dialysis and their caregivers[END_REF]. We cannot fully compare our results as we are the first to analyse the relative's quality of life in pituitary disease. Nonetheless, we can hypothesize that results might be different around the diagnostic of acromegaly when there are more outpatient visits, exams, and hospitalisations and therefore a higher burden for the partner.

Our study has a main limitation, which is the fact that we were not allowed to analyse the medical history of the partners/family members as they were not inpatients of our institution. We thus cannot exclude that comorbidities might have biased the results of the 2 nd part of our study, ie. their sociocognitive characteristics. However, the way we evaluated our main criterion, was based on a direct comparison of the perception of the disease by the partner and the patient, was not biased by this medical history. Moreover, as a pilot study, we only included a small number of patients, and our results will have to be confirmed by studies on a larger number of patients. Finally, we decided to include patients controlled by medical treatments or cured by surgery: even if we acknowledge that quality of life is probably different between these 2 groups of patients, this should not modify our main criterion.

To conclude, our study is the first to analyse the consequences of acromegaly of the patient-relative dyad in a matter of quality of life and self-esteem. Partners appeared to have an accurate estimation of the patient's body image. However, they were not fully aware of the consequences of body changes on their relative's self-esteem. It is thus proposed that discussing body image may be used as a starting point in dealing with the patient's perception of the disease. Both patients and relatives reported a decrease in social interactions although seeking social support was not a preferential coping mechanism. Education programs should therefore consider a specific approach not only centred on the patient, but also on the familial environment. In a disease such as acromegaly in which physical and psychological sequelae can likely remain for years after remission, endocrinologists will probably increase the overall acceptance of the disease by involving the relatives. Overestimation by the relative: 50% for the Rosenberg, 45.8% for the IBQ and 19.2% for Stunkard's figurines. Underestimation: 42.3%, 54.1% and 38.4% respectively. Accurate estimation: 7.7%, 0% and 42.3% respectively.

LEGEND TO TABLES AND FIGURES

Table 1 :Table 2 :Figure 1 :

 121 Figure 1: Correlation scores between the patient and the relative for Rosenberg, Body image and

Availability of data: Data collected for the study, including individual participant data and a data dictionary defining each field in the set, will be made available upon reasonable request to the corresponding author.

Declarations

Funding: This work was financed by an unrestricted educational grant received from IPSEN.

Competing interests:

The authors have nothing to disclose.

Ethics approval: Informed written consent was obtained from all patients to participate in the study, which was approved by the ethics committee of Aix Marseille University.


Christian Ghiaus
email: christian.ghiaus@insa-lyon.fr

Causality issue in the heat balance method for calculating the design heating and cooling load

Keywords: thermal network, state space model, proper and improper transfer function, linear algebra

Causality issue in the heat balance method for calculating the design heating and cooling load

Christian Ghiaus

Introduction

Estimating the heating or cooling load means to find out what is the power needed by the building in order to maintain the indoor temperature at a desired value when the heat gains and losses (due to weather, ventilation rates, occupancy, usage of electrical appliances, etc.) vary. Thermal load calculation is one of the most important tasks of building design. It is used to improve the architectural design by estimating the energy performance of buildings, to size the heating, cooling and ventilation systems, to conceive the control strategies for lowering the peak demand and the energy consumption, to analyze the life cycle cost. Thermal simulation of buildings by employing the heating and cooling load is used to check if a building design complies with standards and regulations.

It is generally accepted that the heat balance method is fundamental for the thermal load calculation. Other procedures are seen as variants or simplifications of this method. Therefore, heat balance method is used by simulation software such as TRNSYS (TRNSYS 2007), EnergyPlus (EnergyPlus 2012), ESP-r [START_REF] Clarke | Energy Simulation in Building Design 2 nd edition[END_REF], CLIMA [START_REF] Skeiker | Advanced software tool for the dynamic analysis of heat transfer in buildings: applications to Syria[END_REF], etc. It is the recommended method of ASHRAE [START_REF] Bellenger | Nonresidential cooling and heating load calculation procedures[END_REF]) and the basis of the CEN standard for calculation of the design heat load [START_REF] Cen | Heating systems in buildings -Method for calculation of the design heat load EN 12831[END_REF] as well as of other standards related to thermal performances of buildings (e.g. CEN 2007a[START_REF] Cen | Thermal performance of buildings -Sensible room cooling load calculation -General criteria and validation procedures EN 15255[END_REF].

The principle of the method is to do the heat balance on the air volume of the room. This implies calculating the conductive, convective and radiative heat balance for each surface of the room and the convective heat balance for the room air (e.g. Wang The result of the heat balance is an input-output model. A dynamical input-output system is causal when its output depends only on past inputs and noncausal if its output depends also on future inputs. All physical systems need to be causal [START_REF] Romagnoli | Introduction to Process control 2nd edition[END_REF]). An acausal system is not physically realizable because it requires prediction (Morari and Zafiriou 1989).

When used for calculating the indoor temperature, the heat balance method is a direct method (i.e. it respects the physical causality): the output, i.e. the indoor temperature, depends on more inputs, including the thermal load. However, when it is used to estimate the load, the heat balance method, in its actual form, results in an acausal system: the output of the physical system, i.e. the indoor temperature, becomes an input of the algorithm for calculating the thermal load. The practical consequence of the acausality is that, if the indoor room temperature has a step variation, the calculated load tends to infinity when the time step for calculation tends to zero. Mathematically, if the thermal model is linear (or piece-wise linearized for one or more time steps during the simulation), it is characterized by an improper transfer function which grows unbounded as the frequency approaches infinity. Practically, when the indoor temperature varies, the value of the calculated thermal load will depend on the value of the time step used in simulation. The physical explanation is rather intuitive: an infinite power is needed to change the indoor temperature instantaneously. Therefore, the sizing of the HVAC system, which is determined by the calculated load, will depend on the time step used in simulation. Practitioners noticed this problem and introduced load limitations in the heat balance method (e.g. CEN 2007b).

The way in which the heat balance method for load calculation is usually presented makes it difficult to notice that the model used for estimating the load results in an improper transfer function that does not respect the causality. However, if the whole system of equations is written in state space form and then transformed into a set of transfer functions, it becomes evident that some transfer functions are improper, which reflects their acausality.

Principle of the heat balance method

The most important assumption of the heat balance (HB) method is that the temperature of the air in the thermal zone is homogeneous. The other important assumptions are that the surfaces (walls, windows, floors, ceilings, etc.) have uniform temperature, uniform long-wave and shortwave irradiation and diffuse radiation [START_REF] Bellenger | Nonresidential cooling and heating load calculation procedures[END_REF]). The heat conduction is usually treated in one dimension but two or three dimensions may be handled. Based on these assumptions, heat balance is done for each wall for outside face, nodes in the solid walls and inside face. Finally, the convective heat balance is done for the room air (Figure 1): [START_REF] Bellenger | Nonresidential cooling and heating load calculation procedures[END_REF] where:

is the room air temperature considered homogeneous [°C or K];

-rate of increase of heat stored in the room air volume of mass and If the thermal capacity of the air in the zone is neglected, then equation (1) becomes:

(2)

In the heat balance method, the load, , is determined from the inside air heat balance. If the heat capacity of the room air is not neglected, then

If the heat capacity of the air is neglected, then (4) In equations ( 3) and (4) the room air temperature, , is considered known because it is the desired temperature of the indoor air.

HVAC g v ci a a a Q Q q q dt d c m ! ! + + + = q a q dt d c m a a a q a m a c å - = i a ci i i ci h S q ) ( q q i S ci q i h ) ( a o a v v T c m q q - = ! v m ! o T g Q ! HVAC Q ! 0 = + + + HVAC g v ci Q Q q q ! ! HVAC Q ! g v ci a a a HVAC Q q q dt d c m Q ! ! - - - = q g v ci HVAC Q q q Q ! ! - - - = page 4/22
Finding indoor air temperature from equations (1) or ( 2) is a simulation (or direct) problem.

Finding the load from equation ( 1) or ( 2) is a control (or inverse) problem. If the model for heat transfer is linear (or piecewise linear), than equations (3) or (4) result in an improper transfer function. In order to demonstrate this, we will first obtain a state space formulation of the thermal model. Then, from the state space form we will obtain the transfer functions and we will demonstrate that equations (3) or (4) end in an improper transfer function. A simple numerical example, which does not reduce the generality of the problem, will be used to illustrate these steps.

Differential algebraic equation representation of the thermal model

The thermal model of the room is described by a system of differential algebraic equations. In order to obtain this system, an algorithm is proposed which is inspired by [START_REF] Strang | Introduction to Applied Mathematics[END_REF][START_REF] Strang | Introduction to Linear Algebra[END_REF],

) (t a q a q HVAC Q !
who introduced it for systems in equilibrium, and by the Modified Nodal Analysis method, frequently used in electrical circuits modeling and analysis. As compared to electrical circuits, the thermal circuits do not have an element equivalent to inductance. Thermal circuits may be arranged in branches, containing thermal resistances and temperature sources, and nodes, to which thermal capacities and flow sources are connected.

! ! = # ! $ = % "#$ -% " + ( ! ) " % ̇" = $ ! -$ !%$ + + "

Figure 2 Typical branch and node in thermal networks

The temperature difference over a resistance is [START_REF] Clarke | Energy Simulation in Building Design 2 nd edition[END_REF] where and are the temperatures in the nodes connected by the branch , and is the temperature source on the branch . Writing equation [START_REF] Clarke | Energy Simulation in Building Design 2 nd edition[END_REF] for all branches, we obtain: [START_REF] Degiovanni | Transmission de l'énergie thermique -Conduction[END_REF] where:

is the vector of temperature drops over the thermal resistances, is the vector of temperature in the nodes, is the vector of temperature sources on the branches.

The matrix is a difference operator for the temperatures while the matrix is the incidence matrix of the network representing the thermal circuit. The rows of the matrix correspond to the branches containing the resistances (i.e. the edges of the graph) and the columns correspond to the nodes representing the temperatures (i.e. the vertices of the graph). The elements of the incidence matrix have the value:

(7) k R l k l k b e q q - + = -1 1 - l q l q k k b k b Aθ e + - = [ ] T m k e e e ... ... 1 = e [ ] T n l q q q ... ... 1 = θ [ ] T m k b b b ... ... 1 = b A - A A k R l q A ï î ï í ì - + = l k l k l k kl R R R a q q q node from exists through passing flow if , 1 node in enters through passing flow if , 1 node to connected not if , 0 page 6/22
The heat flux in each branch is [START_REF] Energyplus | EnergyPlus Engineering Reference[END_REF] Expressed for all fluxes in matrix form, equation ( 8) gives [START_REF] Ghiaus | Calculation of optimal thermal load of intermittently heated buildings[END_REF] where:

is the vector of heat rates in the branches and is a diagonal matrix of thermal conductivities.

The balance of heat rates in a node states that the variation in time of the energy accumulated in the thermal capacity, , is equal to the algebraic sum of heat rates entering the node, and the heat flow sources connected to the node:

Writing the balance equation for all nodes in matrix form gives: [START_REF] Li | Applicability of calculation methods for conduction transfer function of building constructions[END_REF] where:

is a diagonal matrix of thermal capacities;

, the transpose of the incidence matrix, is a matrix operator which makes the algebraic sum of heat transfer rates ;

is the vector of heat rate sources in each temperature node.

By substituting

from equation ( 6) into equation ( 9) and then from equation ( 9) into equation [START_REF] Li | Applicability of calculation methods for conduction transfer function of building constructions[END_REF], we obtain:

(12) or ( 13
)
where we note and .

k k k e R q 1 - = e G q = [ ] T m k q q q ... ... 1 = q ú ú ú û ù ê ê ê ë é = - - 1 1 1 0 0 0 0 0 0 m R R ! G l q l l C q ! å l l q l f l l l l l f q C + = å q ! f q A θ C + = T ! ú ú ú û ù ê ê ê ë é = n C C 0 0 0 0 0 0 1 ! C T A q [ ] T n l f f f ... ... 1 = f b Aθ e + - = Ge q = f Gb A GAθ A θ C + + - = T T ! f b K Kθ θ C + + = b ! GA A K T - º G A K T b º page 7/22
Since some elements of the diagonal matrix are zero, the relations ( 12) and (13) represent a system of differential algebraic equations (DAE): the equations corresponding to zero elements on the diagonal of matrix are algebraic while the equations corresponding to non-zero elements on the diagonal of matrix are differential.

State space representation of the thermal model

Equation ( 13) can be written in state space form by eliminating the temperatures corresponding to the nodes without thermal capacity, i.e. the nodes for which . Equation ( 13) can be written in bloc-matrix form as: [START_REF] Morari | Robust Process Control[END_REF] where:

and correspond to the nodes without thermal capacity; and correspond to the nodes with thermal capacity; is the bloc of the partitioned matrix for which the elements are non-zero; , , , are blocs of the partitioned matrix obtained according to the partitioning of the matrix ; , are blocs of the partitioned matrix obtained according to the partitioning of the matrix ;

, are identity matrices.

The vector needs to be eliminated from equation [START_REF] Morari | Robust Process Control[END_REF]. By multiplying the first row of equation ( 14) by , we obtain: [START_REF] Romagnoli | Introduction to Process control 2nd edition[END_REF] By replacing the second row with the sum of the two rows, equation [START_REF] Romagnoli | Introduction to Process control 2nd edition[END_REF] becomes:

All coefficients of in equation ( 16) are zero which implies that this equation does not depend on temperatures of nodes without thermal capacity, [START_REF] Strang | Introduction to Applied Mathematics[END_REF] Since the matrix in equation ( 17) is invertible, the state space representation of the modeled thermal network is: [START_REF] Strang | Introduction to Linear Algebra[END_REF] or (19) where the state matrix is [START_REF] Wang | A novel and simple building load calculation model for building and system dynamical simulation[END_REF] the input matrix is [START_REF] Xu | A simplified dynamic model for existing buildings using CTF and thermal network models[END_REF] and the input vector is:

C C C l 0 = ll C ú û ù ê ë é ú û ù ê ë é + ú û ù ê ë é + ú û ù ê ë é ú û ù ê ë é = ú û ù ê ë é ú û ù ê ë é C b b C C C f f I 0 0 I b K K θ θ K K K K θ θ C 0 0 0 0 22 11 2 1 0 22 21 12 11 0 ! ! 0 θ 0 f C θ C f C C C 11 K 12 K 21 K 22 K K C 1 b K 2 b K b K C 11 I 22 I 0 θ 1 11 21 - - K K ú û ù ê ë é ú û ù ê ë é- + ú û ù ê ë é- + ú û ù ê ë é ú û ù ê ë é - - = ú û ù ê ë é ú û ù ê ë é - - - C b b C C C f f I 0 0 K K b K K K K θ θ K K K K K K θ θ C
! ! [ ] [ ] [ ] ú û ù ê ë é - + + - + ú û ù ê ë é + - = ú û ù ê ë é - - - C b b C C C f f I K K b K K K K θ θ K K K K 0 θ θ C 0 0 22
(22)
Considering that the output is , the state space model is:

(23)
[ ] ú û ù ê ë é - + + - + + - = - - - C b b C C C f f I K K b K K K K θ K K K K θ C 0 22 1 11 21 2 1 1 11 21 22 12 1 11 21 ) ( ) ( ! C C [ ] ú ú ú û ù ê ê ê ë é - + - + + - = - - - - - C b b C C C C f f b I K K K K K K C θ K K K K C θ 0 22 1 11 21 2 1 1 11 21 1 22 12 1 11 21 1 ) ( ! u B θ A θ S C S C + = ! ) ( 22 12 1 11 21 1 K K K K C A + - = - - C S [ ] 22 1 11 21 2 1 1 11 21 1 I K K K K K K C B - - - - + - = b b C S [ ] T C f f b u 0 = a q u D θ C u B θ A θ S C S a S C S C + = + = q ! a q a q S C a q C θ S D a q a q a q 0 θ 0 θ ( ) [ ] ÷ ÷ ÷ ø ö ç ç ç è ae ú ú ú û ù ê ê ê ë é + - = + + - = - - C b b f f b 0 I K θ K K f I b K θ K K θ 0 11 1 C 12 1 11 0 11 1 C 12 1 11 0 u D θ C θ S C S + = 0 12 1 11 K K C - - = S page 9/22 , & = --$$ #$ [-'$ 0 $$ 1] (27) and (28)
Then can be extracted from the vector of temperatures .

Transfer function representation of the thermal model

The relation between the inputs, , and the output, , can be expressed as a set of transfer functions or a transfer matrix. Applying Laplace transform for zero initial conditions, the set of equations (23) become:

(29

)
where is the complex variable. From the first equation of (29), we obtain . Eliminating the state vector from (29), we obtain

is the transfer matrix in which each element is a transfer function relating the output to a particular input from the input vector . Conduction transfer functions (CFT) are widely used in thermal load calculation. The coefficients of the transfer functions can be calculated precisely [START_REF] Li | Applicability of calculation methods for conduction transfer function of building constructions[END_REF].

Transfer function model in the case of non-negligible thermal capacity of the indoor air

We will consider a thermal network (Figure 3 a) that is an epitome of the heat balance method (Figure 1). Since the aim is to present the principle, the model from Figure 3 keeps the essential aspects of the heat balance method at a minimum number of temperature nodes. Normally, the wall would consist of several layers, each layer discretized in order to obtain temperatures at several nodes inside it. Generally, different surfaces would be considered for the room walls, windows, doors, floors, ceilings, etc. Thermal bridges can be also included in the model by treating heat conduction in two and three dimensions. But the simplicity of the model shown in Figure 3 does not reduce its generality for the problem presented hereafter.

Any thermal network can be easily put to comply with the typical branch and node shown in Figure 2. For example, the thermal network of Figure 3b is the equivalent of the network

[ ] T C f f b u 0 = a q 0 θ u a q u D θ C u B θ A θ S C S a S C S C s + = + = q w s j s + = u B θ A I S C S s = - ) ( C θ u D B A I C ] ) ( [ 1 S S S S a s + - = - q S S S S s D B A I C H + - = -1
) ( a q u represented in Figure 3a. For a numerical example, let us consider that the elements of the circuit shown in Figure 3 have the values given in Table 1.

For the thermal network presented in Figure 3, the temperature differences for each thermal resistance are: q q q q q q q q 5 4 3

= e T o o T T ] 0 0 0 [ = b T w a si so ] [ q q q q = θ ci w w co v w a si so R R R R R 2 1 0 1 1 0 1 0 1 0 1 0 0 1 0 0 0 1 0 1 0 0 ú ú ú ú ú ú û ù ê ê ê ê ê ê ë é - - - = q q q q A b Aθ e + - = page 11/22
The heat transfer rates in the network branches can be expressed as e R q q e R q q e R q q e R q q e R q q ci ci w w the set of equations ( 36) or (37), representing the heat flow in each branch, can be written as (38)

q q q ci w w co v T q q q q q ] [
The heat balance in each temperature node gives:

(39)

It is worth noting that the 3 rd equation, , where and from the system of equations ( 39) is, in fact, the same as equation ( 1) used in the heat balance method. In matrix form, the system of equations ( 39 ;

q = T w a si so ] [ q q q q = θ ï ï î ï ï í ì - = + + + = + - = + - = 4 3 5 1 5 4 3 2 0 0 q q C Q Q q q C Q q q Q q q w w HVAC g a a i o q q ! ! ! ! ! ! HVAC g a a Q Q q q C ! ! ! + + + = 5 1 q v q q º 1 ci q q º 5 ú ú ú ú û ù ê ê ê ê ë é + + ú ú ú ú ú ú û ù ê ê ê ê ê ê ë é ú ú ú ú û ù ê ê ê ê ë é - - - = ú ú ú ú û ù ê ê ê ê ë é ú ú ú ú û ù ê ê ê ê ë é 0 0 1 1 0 0 1 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 4 3 2 1 HVAC g i o w a si so w a Q Q Q Q q q q q q C C ! ! ! ! ! ! q q q q T HVAC g i o Q Q Q Q ] 0 [ ! ! ! ! + = f ú ú ú ú û ù ê ê ê ê ë é = w a C C 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C ú ú ú ú û ù ê ê ê ê ë é - - - = 0 1 1 0 0 1 0 0 0 1 1 1 0 0 0 0 0 1 1 0 T A page 14/22
f q A θ C + = T ! Ge q = b Aθ e + - = f Gb A GAθ A θ C + + - = T T ! f b K Kθ θ C + + = b ! ú ú ú ú ú ú û ù ê ê ê ê
G A K T b ú û ù ê ë é - - = 125 9 . 2 0 0 0 0 0 9 . 2 250 0 1 b K ú û ù ê ë é - = 0 9 . 2 9 . 2 0 0 125 0 0 0 3 . 38 2 b K [ ] ú ú ú û ù ê ê ê ë é - + - + + - = - - - - - C b b C C C C f f b I K K K K K K C θ K K K K C θ 0 22
C C ú û ù ê ë é = 1 0 0 1 22 I u B θ A θ S C S C + = ! ú û ù ê ë é × - × × × - = + - = - - - - - - 6
I K K K K K K C B - - - - + - = b b C S page 16/22
The input vector is (52) where , and . However, we can define the input vector by taking only the elements 1, 2, 6, 7, and 8 from the input vector given by ( 52),

and, therefore, only the columns 1, 2, 6, 7, and 8 from the matrix ,

the state space model corresponding to the thermal circuit from Figure 3 with the parameters given in Table 1 is

(56)
where , is given by (53), is given by (50), is given by (54), and is given by (55).

The transfer matrix corresponding to the state space model (56) is obtained by substituting and in equation (31):

T b T 0 f T C f [ ] T HVAC g i o o o T C Q Q Q Q T T ] 0 0 0 0 [ 0 ! ! ! ! + = = f f b u T o o T T ] 0 0 0 [ = b T i Q Q ] [ 0 0 ! ! = f T HVAC g C Q Q ] 0 [ ! ! + = f u [ ] T HVAC g i o ow ov Q Q Q Q T T ! ! ! ! + = u S B ú û ù ê ë é × × × × × × = - - - - - - 0 10 7 .
0 1 = S C C S a S C S C θ C u B θ A θ = + = q ! [ ] T w a C q q = θ u S A S B S C S A S B S C page 17/22
By using the transfer functions from the transfer matrix, the output is or (57) All transfer functions , …, are strictly proper (they respect the physical causality).

It is worth noting that the heat balance on air volume, expressed by equation (1) as the 3 rd equation of the system of equations (39), was used to obtain the expression (57). However, this is not evident at first sight.

Transfer function model in the case of negligible thermal capacity of the indoor air

In the case of negligible capacity in the node , i.e. , the partitioned matrix from equation (46) becomes and the partitioned matrix from equation (47) becomes

T s s s s s s s s s s s s s H H H H H ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú û ù ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë é + × + × × + × + × + × × + × + × + × × + × + × × + × + × × + × = = - - - - -
Q Q Q Q T T H H H H H ! ! ! ! ] [ 5 4 3 2 1 u H q ) ( 5 4 3 2 1 HVAC g i o o o a Q Q H Q H Q H T H T H ! ! ! ! + + + + + = q 1 H 5 H a q 0 = a C K ú û
Q Q Q Q T T H H H H H ! ! ! ! ] [ 5 4 3 2 1 u H q ) ( 5 4 3 2 1 HVAC g i o o o a Q Q H Q H Q H T H T H ! ! ! ! + + + + + = q page 19/22
As in the case of non-negligible thermal capacity of the room air, the heat balance equation ( 2), which is one of the equations from the system of equations (39), was used to obtain the output . However, the relation between equation (2) and equation ( 62) is not self-evident.

Transfer function for the thermal load

In order to obtain the thermal load, the current method uses equation ( 3) or (4). By using the transfer functions, this is equivalent to express from equation ( 57) or ( 62):

(63)

The transfer function , corresponding to input , is improper. The transfer functions corresponding to the influence of the outdoor temperature, , through ventilation and corresponding to the influence of the heat flux on the indoor walls, , are proper. They correspond to inputs connected directly to the node , i.e. not through a node containing a thermal capacity. The functions corresponding to the influence of the outdoor temperature, , through walls and corresponding to influence of the heat flux on the outdoor wall, , are strictly proper. They correspond to inputs connected to the node through nodes containing at least a thermal capacity.

It results that this method yields an improper function for input . This is true not only for the two epitomes of non-negligible and negligible indoor air capacity presented above, but for any model of this type. For any thermal network constructed according to the heat balance method (Figure 1):

• the transfer function relating the load to the temperature in node will be always improper,

• the transfer functions relating the inputs to the output through branches which do not contain at least a thermal capacity will be strictly proper and • the transfer functions connecting the inputs to the output by branches containing at least a thermal capacity will be proper.

Consequences of the improper transfer function

The models of multiple input -multiple output systems (MIMO) describe a causal relation from input to the output if and only if the transfer matrix is proper, i.e. if in any entry of the transfer matrix the degree of the denominator is larger than or equal to the degree of the numerator (Mareels and Polderman 1996). A system with an improper transfer matrix is not physically realizable because it requires prediction [START_REF] Morari | Robust Process Control[END_REF]. Furthermore, an improper transfer function amplifies the high frequencies.

q HVAC Q ! g i o o o a HVAC Q Q H H Q H H T H H T H H H Q ! ! ! ! - - - - - = - - - - - 4 1 5 3 1 5 2 1 5 1 1 5 1 5 q 1 5 - H a q 1 1 5 H H - o T 4 1 5 H H - i Q ! a q 2 1 5 H H - o T 3 1 5 H H - o Q ! a q a q a q a q a q
Heat balance method for load calculation results in an improper transfer function. The practical consequence is that this method cannot be used to calculate the load, , when the space temperature, , varies. For example, if the air space temperature has a step change (which is frequently the case in simulations), the high frequencies present in the signal of will be amplified. If the simulation time step tends to zero, the calculated load will tend to infinity, which technologically is obviously not achievable. One solution is to limit the value of the load [START_REF] Cen | Thermal performance of buildings -Sensible room cooling load calculation -General criteria and validation procedures EN 15255[END_REF], which implies that the simulation will not be used to calculate the maximum load; however, it can be useful for energy consumption estimation. Another solution is to introduce a term in the transfer function of in order to make it proper. This is achievable by considering the thermal inertia of the HVAC system or by introducing a feed-back control loop, in which case the maximum load depends on the controller type and settings (Ghiaus and Hazyuk 2010).

Conclusions

Heat balance method used for simulation of the time variation of the air space temperature as a function of inputs has strong physical basis. The most restrictive assumption is that the air in the thermal zone has a uniform temperature. The heat balance results in a set of differential algebraic equations which, after piece-wise linearization, can be written in state space form and then transformed in a set of transfer functions. These models respect the causality: the outputs depend on inputs through the state variables. Consequently, the transfer functions are proper.

Thermal load calculation is an inverse problem. One input of the physical system, the thermal load, becomes an output of the model and one output of the physical system, the indoor temperature, becomes an input of the model used in calculation. The method used to solve this problem does not respect the causality and results in an improper transfer function. Therefore, the high frequencies present in the signal of the indoor temperature will be amplified. The practical consequence is that, if the indoor temperature varies, the calculated load will depend on the simulation time step: smaller the time step, larger the calculated load. If the room temperature has a step variation, the calculated load will tend to infinity when the simulation step tends to zero. This effect is less evident if the simulation time step is one hour (the usual value in building simulation software) but becomes evident if the time step is reduced.

The problem of non causality can be solved by introducing a filter which will make the transfer function proper. This can be achieved by transforming the load calculation into a control problem. Then, the load will be the output of the control algorithm.

The issue of causality is not self-evident for complex systems. Therefore, the claim of simulation systems which are equation-based (i.e. in which the equations describe equalities, not assignments) that there is a benefit in the ability to choose arbitrarily the inputs and the outputs of the model (Modelica 2000[START_REF] Dymola | Dymola. Dynamic Modeling Laboratory[END_REF]) is questionable if this choice results in improper transfer functions.

  and Chen 2000, Bellenger et al. 2001, Hernandez et al. 2003, Xu and Wang 2008, Li et al. 2009).

  rate due to mass flow of outdoor air at temperature introduced into the space by infiltration and ventilation [W]; -convective heat flow gained from the internal loads [W]; -heat flow to/from the HVAC system [W].

Figure 1

 1 Figure 1 Heat balance method (adapted after[START_REF] Bellenger | Nonresidential cooling and heating load calculation procedures[END_REF]

Table 1 .

 1 Parameter of the thermal network shown in Figure3

Figure 3

 3 Figure 3 Epitome thermal network for the heat balance method: a) usual representation; b) representation in the form of typical branches (according to Figure 2)

  the matrix operator for differences of flows in temperature nodes is the transpose of the incidence matrix, , equation (40) becomes

Ge

  in equation (31) the values of the matrices and given above, we obtain the matrix transfer:As for the case of non-negligible thermal capacity of the indoor air, all transfer functions , …, are proper or strictly proper. By using the transfer functions, the output is or

  a

Acknowledgements

This work was conducted in the framework of the International Energy Agency Annex 53 Total energy use in buildings: analysis and evaluation methods and in the research project Fault detection and diagnosis of energy systems (AIDE-3D) financed by the French National Research Agency (ANR). Special thanks are given to Ion Hazyuk (INSA Lyon, France), Eric Blanco (Ecole Centrale de Lyon, France), Ad van der Aa (Cauherg-Huygen Raadgevende Ingenieurs BV, The Netherlands) and Ian Hensen (Eindhoven University of Technology, The Netherlands) for their interest, discussions, and support.


Florian Miralles
email: florian.miralles@umontpellier.fr

Bastien Sauvage
email: bastien.sauvage@inria.fr

Stephen Wornom
email: stephen.wornom@inria.fr

Bruno Koobus
email: bruno.koobus@umontpellier.fr

Alain Dervieux
email: alain.dervieux@inria.fr

Modélisation hybride RANS/DVMS autour d'obstacles statique et tournants

Keywords: modèles hybrides de turbulence, modèle variationnel multi-échelle pour LES, méthode chimera, écoulement autour d'un cylindre, machine tournante, grilles non structurées hybrid turbulence models, variational multiscale model for LES, chimera method, cylinder flow, rotating machine, unstructured grids

. Ce travail a été soutenu par le projet ANR NORMA dans le cadre de la subvention ANR-19-CE40-0020-01.

the space of large scales V h the space of small scales Inria 1 Introduction

In this paper, hybrid models are evaluated for the simulation of massively separated flows around fixed and moving geometries, with the objective of applying them to aeroacoustic problems involving complex industrial flows at high Reynolds numbers. For this purpose, the turbulence models need to behave preferentially like large eddy simulation models (LES) in order to introduce little dissipation and thus better capture the small scales and the fluctuations of the unsteady flows considered. However, these flows involve very thin boundary layers, which current computers and softwares cannot compute in LES mode, only in RANS mode. To accommodate both needs, LES and RANS, the hybrid turbulence approach, combining LES and RANS in a somewhat zonal manner, is considered by many research teams.

In this work, beside the classical Delayed Detached Eddy Simulation model (DDES, [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF]), we study an approach which combines zonally this last model and the dynamic variationnal multiscale model (DVMS, [START_REF] Koobus | A variational multiscale method for the large eddy simulation of compressible turbulent flows on unstructured meshes-application to vortex shedding[END_REF]). In this hybrid DDES/DVMS approach [START_REF] Moussaed | Simulation of the flow past a circular cylinder in the supercritical regime by blending rans and variationalmultiscale les models[END_REF], the DVMS model is preferentially activated in the wake regions where this latter model introduces less dissipation than the LES component of DDES. A hybrid RANS/DVMS model [START_REF] Itam | Combining a ddes model with a dynamic variational multiscale formulation[END_REF], where the RANS component is Goldberg's k -ε model [START_REF] Goldberg | A wall-distance-free k -ε model with enhanced near-wall treatment[END_REF], is also applied in this study. A smooth blending function, which is based on the value of a blending parameter, is employed in these hybrid strategies, for switching from RANS to DVMS or DDES to DVMS. In [START_REF] Itam | Combining a ddes model with a dynamic variational multiscale formulation[END_REF], these models are applied to flows around a cylinder at moderate Reynolds number, starting from Reynolds number 20000. In many industrial cases, and particularly cases with rotating geometries like propellers, the Reynolds number is much higher.

In the present paper, these hybrid models are applied to the flow around a circular cylinder in the supercritical regime, namely at Reynolds numbers 1 million and 2 million. At Reynolds number 1 million, the supercritical regime shows turbulent flow separation on both the upper and lower surface of the cylinder with a laminar-turbulent transition in the boundary layers located between the stagnation point and the separation one. At Reynolds number 2 million, the high transition regime shows a fully turbulent boundary layer on one side of the cylinder, and a partially laminar, partially turbulent boundary layer on the other side. This benchmark, which contains many characteristics encountered in industrial flows, is challenging due to the complex physics of the flow and the considered high Reynolds numbers, and remains a difficult calculation to perform. The second application concerns the flow at Reynolds number 1.8 million around a cross shaped mixing device rotating inside a cylinder.

This work is part of a cooperation program with Keldysh Institute of Applied Mathematics of Moscow focussing on the "Efficient simulation of noise of rotating machines", see e.g. [START_REF] Garanzha | Moving deforming mesh generation based on quasi-isometric functional[END_REF][START_REF] Bakhvalov | On using artificial viscosity in edge-based schemes on unstructured grids[END_REF][START_REF] Duben | Ebr schemes with curvilinear reconstructions of variables in the near-wall region[END_REF].

The remainder of this paper is organized as follows. Section 2 presents the hybrid turbulence models used in this work. Section 3 describes the numerical scheme. In Section 4, applications are presented. The obtained results are analyzed and compared with those obtained in other numerical studies and with experimental data. Finally, concluding remarks are drawn in Section 4.

Turbulence modeling

RANS

First, we want to specify that RANS stands for unsteady RANS throughout this document. When comparing with other published results, we use the notation that appear in the published Miralles articles notable URANS which is shorthand for Unsteady RANS, which the same as RANS in this document. Steady RANS will be denoted as SRANS and Steady LES as SLES.

In this work, and as far as the closure of the RANS equations is concerned, we use the low Reynolds k -ε model proposed in Goldberg et al. [START_REF] Goldberg | A wall-distance-free k -ε model with enhanced near-wall treatment[END_REF]. This model was designed to improve the predictions of the standard k -ε one for adverse pressure gradient flows, including separated flows. The equations of this model are recalled in this document, in dimensionless frame :

µ t = ρc µ f µ (k, ε) k 2 ε , (1)
where the damping function is defined as follow :

f µ (k, ε) = 1 -e aµRt(k,ε) 1 -e - √ Rt(k,ε) max 1, c τ R t (k, ε) , R t (k, ε) = Re k 2 ε . (2)
function R µ represent the turbulent Reynolds number. Transports equation on k and ε can be written as :

∂ρk ∂t + ∇ • (ρ uk) = τ : ∇ u + ∇ • 1 Re + µ t σ k ∇k -ρε, (3)
∂ρε ∂t + ∇ • (ρ uε) = ∇ • 1 Re + µ t σ ε ∇ε + c (1) ε ε k τ : ∇ u -c (2) ε ρ ε 2 k + C (2) (k, ε) C (1) (k, ε). (4)
Since to control production and dissipation of dissipation rate, functions C (1) and C (2) are introduced. These functions are defined such that :

C (1) (k, ε) = k ε max 1, c τ √ R t , C (2) (k, ε) = 3 10 ρ max √ k, 1 Re ε 1/4
εC (1) max (∇k • ∇ k ε , 0).

The total system, including the fluids variables and turbulent variables, can be written as :

   ∂W ∂t + ∇ • F (W ) = τ RAN S (W ), +Boundary conditions, Q(0, x) = Q 0 (x), ∀x ∈ Ω f (5)
Where τ RAN S means the RANS closure terms, defined such that :

τ RAN S (W ) =     ρ 0 , ρu 0 , ρE 0 , ρk τ : ∇ u -ρ , ρ c (1) ε ε k τ : ∇ u -c (2) ε ρ ε 2 k + C (2) (k, ε) C (1) (k, ε)     (6)
In practice, for compressible flow, we use the constants summarize in the table below :

c (1) ε c (2) ε σ k σ ε 1.44 1.92 1.0 1.3 Inria 2.

DDES model

The DDES approach used in this work is based on the Goldberg's k -ε model [START_REF] Goldberg | A wall-distance-free k -ε model with enhanced near-wall treatment[END_REF]. In other words, this RANS model is introduced in the DDES formulation [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF] l DDES where

l DDES = k 3/2 ε -f d × max(0, k 3/2 ε -C DDES ∆) with f d = 1 -tanh((8r d ) 3 )
and

r d = ν t + ν max( √ u i,j u i,j , 10 -10 )K 2 d 2 w .
K denotes the von Karman constant (K = 0.41), d w the wall-normal distance, u i,j the x jderivative of the ith-component of the velocity u, and the model constant C DDES is set to the standard value 0.65. We define here the DDES k -ε Goldberg source term :

τ DDES (W ) =     ρ 0 , ρu 0 , ρE 0 , ρk τ : ∇ u -D DDES k , ρ c (1)
ε ε k τ : ∇ u -c (2) ε ρ ε 2 k + C (2) (k, ε) C (1) (k, ε)     (7)
Symbol ∆ holds for the filter width, in the type of LES which we use, the filter width is fixed to the local mesh size. The definition of the local mesh size is not an easy question when (highly) anisotropic meshes are used. Two main options are (1) the third root of the grid element volume and (2) the largest local edge. However, our experience is that, when combined with devices like VMS and Dynamic-LES, the choice in defining the local mesh size is not so sensitive. Our filter width is defined as the third root of volume of the grid element T, then :

∆ = ∆ T = T dx 1/3 . ( 8
)
Miralles

High Reynolds computation

In order to catch the turbulent boundary layer, we use a wall law modeling, designed by Reichardt [START_REF] Reichardt | Vollständige darstellung der turbulenten geschwindigkeitsver-teilung in glatten leitungen[END_REF]. By this way, the no slip boundary condition on the wall become :

U + = 1 κ ln 1 + κy + + 7.8 1 -exp -y + 11 - -y + 11 exp -y + 3 (9)
Also, for turbulent variables, we use the work of Jaeger [START_REF] Jaeger | An extended k-epsilon finite element model[END_REF], with the following condition on the wall :

k = u 2 f C µ y + 10 2 , ε = Re u 4 f 10κ y + 10 2 + 0.2 C µ 1 - y + 10 2 .
The Reichardt wall law give us a smooth transition from linear to logarithmic boundary layer.

LES-like model : DVMS

In this section, we present the DVMS model which is preferred to the classical LES approach in our hybrid strategies because of some specific properties, as will be specified hereafter, that allow this model to have a better behavior in some regions of the flow, such as shear layers and wakes.

The Variational Multiscale (VMS) model for the large eddy simulation of turbulent flows has been introduced in [START_REF] Hughes | Large eddy simulation and the variational multiscale method[END_REF] in combination with spectral methods. In [START_REF] Koobus | A variational multiscale method for the large eddy simulation of compressible turbulent flows on unstructured meshes-application to vortex shedding[END_REF], an extension to unstructured finite volumes is defined. This method is adapted in the present work. Let us explain this VMS approach in a simplified context. Assume the mesh is made of two embedded meshes. On the fine mesh we have a P 1 -continuous finite-element approximation space V h with the usual basis functions Φ i vanishing on all vertices but vertex i. Let be V 2h its embedded coarse subspace V 2h . Let V h be the complementary space : V h = V 2h ⊕ V h . The space of small scales V h is spanned by only the fine basis functions Φ i related to vertices which are not vertices of V 2h . Let us denote the compressible Navier-Stokes equations by :

∂W ∂t + ∇ • F (W ) = 0 where W = (ρ, ρu, E). The VMS discretization writes for W h = W i Φ i : ∂W h ∂t , Φ i + (∇ • F (W h ), Φ i ) = -τ LES (W h ), Φ i ∀a i ∈ Ω h , Q h (0, a i ) = Q 0 (a i ), ∀a i ∈ Ω h (10)
Here the convective and viscous flux do not contains any closure equations. Moreover, the closure terms is defined as the following :

τ LES (W h ), φ i := 0, M S (W h , Φ i ), M H (W h , Φ i ) T (11)
and

M S (W h , Φ i ) := T ∈Ω h T µ sgs t P ∇Φ i dx (12)
Inria here P = 2S -1 3 T r(S)Id . The basis function are written as Φ k = Φ k -Φ k where :

Φ k = V ol(C k ) j∈I k V ol(C j ) j∈I k Φ j (13)
I k means the index set of cells wich are contains in macro-cell C k obtained by agglomeration process, see also [START_REF] Lallemand | Unstructured multigridding by volume agglomeration : Current status[END_REF]. For the subgrid terms in relation with heat transfer, we write :

M H (W h , Φ i ) = T ∈Ω h T C p µ sgs t P r t ∇T • ∇Φ i dx ( 14
)
For a test function related to a vertex of V 2h , the RHS vanishes, which limits the action of the LES term to small scales. In practice, embedding two unstructured meshes V h and V 2h is a constraint that we want to avoid. The coarse level is then built from the agglomeration of vertices/cells as sketched in Fig. 1. It remains to define the model term τ LES (W h ). This term represents the subgrid-scale (SGS) stress term, acting only on small scales W h , and computed from the small scale component of the flow field by applying either a Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF] or a WALE SGS model [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF], the constants of these models being evaluated by the Germano-Lilly dynamic procedure [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF]. The resulting model is denoted DVMS in this paper. It has been checked [START_REF] Moussaed | Impact of dynamic subgrid-scale modeling in variational multiscale large-eddy simulation of bluff body flows[END_REF] that combining VMS and dynamic procedure effectively brings improved predictions.

A key property of the VMS formulation is that the modeling of the dissipative effects of the unresolved structures is only applied on the small resolved scales, as sketched in Fig. 2. This property is not satisfied by LES models which damp also the other resolved scales. Important consequences are that a VMS model introduces less dissipation than its LES counterpart (based on the same SGS model) and that the backscatter transfer of energy from smallest scales to large scales is not damped by the model. These VMS models then generally allow a better behavior near walls, in shear layers and in the presence of large coherent structures.

SGS model

To obtain the SGS viscosity needed to close the above VMS formulation, we have to define the SGS viscosity coefficient µ sgs t . A first option is the widely used Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF] is first considered. In the adopted VMS formulation this writes :

µ sgs t = ρ(C s ∆) 2 | ˇ S|, with | ˇ S| = ˇ S : ˇ S, (15)
where C s is the Smagorinsky coefficient. A typical value for the Smagorinsky coefficient for shear flows is C s = 0.1, which is used herein.

Miralles

A second SGS model is the Wall-Adapting Local Eddy -Viscosity (WALE) SGS model proposed by Nicoud and Ducros [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]. The eddy-viscosity term in the VMS formulation is defined as follows :

µ sgs t = ρ(Cw∆) 2 τ W : τ W 3/2 Š : Š 5/2 + (τ W : τ W ) 5/2 . ( 16
)
with the tensor invariant :

τ W i,j = 1 2 ǧ2 i,j + ǧ2 j,i - 1 2 δ i,j ǧ2 k,k , ǧ2 i,j := k ∂ ǔi ∂x k ∂ ǔi ∂x k . ( 17
)
In this model Cw is set as 0.5, as indicated [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF].

DVMS

Moreover, in this work, the dynamic procedure, which provides a tuning of the SGS dissipation in space and time, is combined with the VMS approach, which limits its effects to the smallest resolved scales, so that the resulting DVMS model enjoys synergistic effects. In their original formulations, Cs and Cw appearing in the expression of the viscosity of the Smagorinsky and WALE SGS model (Eqs. 15 and 16 respectively) are set to a constant over the entire flow field and in time. In the dynamic model [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF], this constant is replaced by a dimensionless parameter C(x, t) that is allowed to be a function of space and time. The dynamic approach also provides a systematic way for adjusting this parameter in space and time by using information from the resolved scales. After the introduction of the grid filter, denoted by overline and tilde, tilde being Favre averaging, f = ρf ρ , a second filter is considered, having a larger width than the grid one, which is called the test-filter and denoted by a hat. The test-filter is applied to the grid-filtered Navier-Stokes equations, and then, the subtest-scale stress is defined as follows :

T test = ρu ⊗ u - ρu ⊗ ρu ρ , and
τ sgs = ρu ⊗ u -ρ u ⊗ u (18)
Then we consider the resolved turbulent tensor :

L = ρ u ⊗ u - ρu ⊗ ρu ρ (19)
This tensor can be interpreted as a intermediate grid filter tensor. We related these quantity by the following equation :

L = T -τ , (20)
Inria By Smagorinsky hypothesis [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF], the deviatoric part of subgrid scale tensor is proportional to stress tensor P = 2 S -2 3 T r(S)Id :

dev τ = τ - 1 3 T r(τ )Id = -ρ(C s ∆) 2 g( u) P . ( 21
)
dev T = T - 1 3 T r(T )Id = -ρ(C s ∆) 2 g( u) P (22)
and where g( u) denotes the contribution to the SGS viscosity depending on the gradient velocity that appears in 15, explicitely g( u) = | S| for the Smagorinsky model, and in 16,

explicitely g( u) = (τ W :τ W ) 3/2 ( Š: Š) 5/2 +(τ W :τ W ) 5/2
, for the WALE model. Then, from the resolved field, we can explicitely compute L. And, moreover L = devL , we have :

devL = L - 1 3 T r(L)Id, (23)
= T -τ - 1 3 T r(T -τ )Id, (24)
= T -

1 3 T r(T )Id -τ - 1 3 T r(τ )Id , (25)
= -(C ∆) 2 ρ g( u) P + (C∆) 2 ρg( u) P (26)
= (C∆) 2   ρg( u) P - ∆ ∆ 2 ρ g( u) P   B . (27)
Thenceforth, the probleme become to find C∆ solution of equation below :

L -(C∆) 2 B = 0 ⇔ min L -(C∆) 2 B 2 . ( 28
)
By convexity and a first optimalality order, the least square problem provide a uniqueness solution :

(C∆) 2 = B : L B : B . (29)
In practice this filter is choosen such that ∆ ∆ > 1. In our case ∆ is defined by :

∆ = 1 N ai Ta i dx 1/3 , ∆ = Ta i dx 1/3
.

Thus ∆ ∆ 2 = (N ai ) 2/3 greater than one.
Note that all quantities in the right-hand side of Eq. 29 are known from the LES computation. Note also that we preferred to dynamically compute (C∆) 2 , instead of C as done in the original dynamic procedure, in order to partially overcome difficulties in the definition of the filter width for inhomogeneous and unstructured grids. Finally, as done also in [START_REF] Holmen | Sensitivity of the scale partition for variational multiscale large-eddy simulation of channel flow[END_REF], the classical dynamic procedure previously briefly outlined, which involves all the resolved scales, is used herein. Once (C∆) 2 is dynamically computed, it is injected in Eq. 15 or 16 to obtain the SGS viscosity used in the VMS approach.

A possible drawback of the dynamic procedure based on the Germano-identity [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF] when applied to a SGS model already having a correct near-wall behavior, as the WALE one, is the introduction of a sensitivity to the additional filtering procedure. A simple way to avoid this inconvenient is to have a sensor able to detect the presence of the wall, without a priori knowledge of the geometry, so that the dynamic SGS model adapts to the classical constant of the model, which is equal to 0.5 in the near wall region for the WALE model, and compute the constant dynamically otherwise. We adopt the sensor proposed in [START_REF] Baya Toda | Is the dynamic procedure appropriate for all sgs models[END_REF], having the following expression :

f SV S = τ W : τ W 3/2 (τ W : τ W ) 3/2 + ˇ S : ˇ S 3 (30)
This parameter has the properties to behave like y + 3 near a solid wall, to be equal to 0 for pure shear flows and to 1 for pure rotating flows. It should be noticed that the implementation of the dynamic SGS models in our software has been optimized so that the additional cost of the resulting dynamic LES and VMS models, in the case of an implicit time-marching scheme, which is our default option, is less than 1% compared to their non-dynamic counterparts.

Hybrid RANS/DVMS

The central idea of the hybrid RANS/DVMS model [START_REF] Moussaed | Simulation of the flow past a circular cylinder in the supercritical regime by blending rans and variationalmultiscale les models[END_REF] applied in this study is to combine the mean flow field obtained by the RANS component with the application of the DVMS model wherever the grid resolution is adequate. In this hybrid model, the k -ε model of Goldberg (subsection 2.1) is used as the RANS component. First, let us write the semi-discretization of the RANS equations :

∂ W h ∂t , Φ i + (∇ • F ( W h ), Φ i ) = -τ RAN S ( W h ), Φ i .
A natural hybridation writes :

∂W h ∂t , Φ i + (∇ • F (W h ), Φ i ) = -θ τ RAN S ( W h ), Φ i -(1 -θ) τ LES (W h ), Φ i
where W h denotes now the hybrid variables and θ is the blending function which varies between 0 and 1 and is defined by :

θ = 1 -f d (1 -tanh(ξ 2 )) with ξ = ∆ l RAN S or ξ = µ SGS µ RAN S ,
the shielding function f d being defined as in Subsection 2.2. The blending function θ allows for a progressive switch from RANS to LES (DVMS in our case) where the grid resolution becomes fine enough to resolve a significant part of the local turbulence scales or fluctuations, i.e. computational regions suitable for LES-like simulations. Additionally, this blending function, thanks to the shielding function f d , prevents the activation of the LES Inria mode in the boundary layer. By way of example, the isocontours of the blending function θ for the flow past a circular cylinder are shown in Fig. 3. We can see in particular that in the boundary layer the RANS model is activated, while in the wake the LES approach is recovered. However, the near-wake is RANS whereas it should be LES. This indicates that the blending function needs futher improvement or that the near-wake mesh is too coarse for LES computation. ) : for θ = 1 (in red) the RANS model is activated, wherever 0 < θ < 1 additional resolved fluctuations are computed, and in the limit θ → 0 (in blue) the LES approach (DVMS in this work) is recovered.

Hybrid DDES/DVMS

The key idea of the proposed hybrid DDES/DVMS model [START_REF] Itam | Combining a ddes model with a dynamic variational multiscale formulation[END_REF] is to use the DVMS approach instead of the LES component of DDES in locations where this component is expected to be activated, especially in wake regions where the DVMS approach allows more accurate prediction.

Assuming that the semi-discretization of the DDES equations writes :

∂W h ∂t , Φ i + (∇ • F (W h ), Φ i ) = -τ DDES (W h ), Φ i ,
the hybrid equations are then defined by :

∂W h ∂t , Φ i + (∇ • F (W h ), Φ i ) = -θ τ DDES (W h ), Φ i -(1 -θ) τ LES (W h ), Φ i
where W h denotes the hybrid variables and θ is the blending function defined as in the previous subsection.

3 Numerical discretization

Numerical scheme

The governing equations are discretized in space using a mixed finite-volume/finite-element method applied to unstructured tetrahedral grids. The adopted scheme is vertex-centered, i.e.

Miralles

all degrees of freedom are located at the vertices. The diffusive terms are discretized using P1 Galerkin finite-elements on the tetrahedra, whereas finite-volumes are used for the convective terms. The numerical approximation of the convective fluxes at the interface of neighboring cells is based on the Roe scheme [START_REF] Roe | Approximate Riemann solvers parameters vectors and difference schemes[END_REF] with low-Mach preconditioning [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF]. To obtain second-order accuracy in space, the Monotone Upwind Scheme for Conservation Laws reconstruction method (MUSCL) [START_REF] Van Leer | Towards the ultimate conservative scheme. iv : a new approach to numerical convection[END_REF] is used, in which the Roe flux is expressed as a function of reconstructed values of W at each side of the interface between two cells. A particular attention has been paid to the dissipative properties of the resulting scheme, since this is a key point for its successful use in LES, and therefore in simulations performed with a hybrid turbulence model. The numerical dissipation provided by the scheme used in the present work is made of sixth-order space derivatives [START_REF] Camarri | A low diffusion muscl scheme for les on unstructured grids[END_REF] and thus is concentrated on a narrow-band of the highest resolved frequencies. This is expected to limit undesirable damping of the large scales by numerical dissipation. Moreover, a parameter γ directly controls the amount of introduced numerical viscosity and can be explicitly tuned in order to reduce it to the minimal amount needed to stabilize the simulation. Time advancing is carried out through an implicit linearized method, based on a second-order accurate backward difference scheme and on a first-order approximation of the Jacobian matrix [START_REF] Martin | A second-order defect correction scheme for unsteady problems[END_REF]. The resulting numerical discretization is second-order accurate both in time and space. It should be noted that the spatial discretization used in this work leads to a superconvergent approximation, i.e. the accuracy can be well above second-order for some Cartesian meshes. One can also add that the objective is not a higher-order convergence but a strong reduction of dissipation and a certain reduction of the dispersion in the general case of a non-Cartesian (but not too irregular) mesh.

Inria

Mesh adaptation for rotating machines

Our numerical model has been extended to rotor-stator simulation with a Chimera technique. The Chimera method aims at solving partial differential equations by decomposition into subdomains with overlap in order to avoid having to use a global mesh. This method allows the communication between the computational subdomains thanks to the overlapping of the subdomains. In our case we consider a decomposition in two domains as shown in Fig. 4,a A mesh adaptation loop based on a metric depending of the Mach number has been combined with the Chimera algorithm. The Transient Fixed Point mesh adaptation algorithm [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems : Application to CFD simulations[END_REF][START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for timedependent problems[END_REF][START_REF] Guégan | An L ∞ -L p mesh adaptive method for computing unsteady bi-fluid flows[END_REF] is applied to each domain from a metric field evaluated on the whole domain.

Applications

Flow past a cylinder

The predictions of the flow around a circular cylinder are presented. Two Reynolds numbers, 10 6 and 2×10 6 , based on the cylinder diameter, D, and on the freestream velocity, are considered. Only a few numerical investigations have been performed for Reynolds numbers higher than 10 6 . This interval is inside the supercritical regime which appears at Reynolds number higher than 2 × 10 5 and for which the separation becomes turbulent [START_REF] Sumer | Hydrodynamics Around Cylindrical Structures[END_REF]. The computational domain is such that -15 ≤ x, y ≤ 15, and -1 ≤ z ≤ 1 where x, y and z denote the streamwise, transverse and spanwise directions respectively, the cylinder axis being located at x = y = z = 0. The mesh contains 4.8 millions nodes. In order to control the computational costs, a wall law is applied in the close vicinity of the wall. For accuracy purpose, the Reichardt analytical law [START_REF] Hinze | Turbulence[END_REF], which gives a smooth matching between linear, buffer and logarithmic regions, is chosen. As the y + normalized distance is generally subject to large variations in complex flows, the wall law is combined with low Reynolds modeling which locally damps the fully-turbulent model in regions in which the wall law does not cover the buffer zone. The behavior of the hybrid models presented in Section 2 are first investigated in terms of flow bulk coefficients. From Table ??, it can be noted that the predictions of these coefficients are globally in good agreement with the experimental data and the numerical results in the literature. The lift time fluctuations are nevertheless better predicted by the RANS/DVMS model compared to the other two hybrid models. This good behavior is confirmed by the correct prediction of the distribution of the mean pressure coefficient, see Fig. 6. On the other hand, it can be observed from Table ?? and Fig. 6 that Smagorinsky and WALE SGS models give globally comparable results. It is also worth noting that the RANS model predicts rather well the bulk coefficients. Fig. 7 shows instantaneous isocontours of the vorticity magnitude for each of the turbulence models used for this benchmark. As one might expect, the unsteady RANS model is more dissipative and captures less flow detail than the hybrid models, with in particular a much more damped and regular wake. The main outputs obtained by the RANS, DDES, DDES/DVMS and RANS/DVMS models are summed up in Table 2. Regarding the two hybrid models using the DVMS approach, only the results obtained with the Smagorinsky SGS model are shown, the WALE SGS model giving very similar results. It can be noticed that all models overpredict the separation angle, and that the fluctuations of the lift coefficient are overpredicted by the RANS and DDES/DVMS model. The RANS/DVMS approach gives on the whole the most satisfactory results. In Fig. 10, which shows the distribution of the mean pressure coefficient, the numerical results obtained with both hybrid models RANS/DVMS and DDES/DVMS are in very good agreement with the experimental results, in particular RANS/DVMS.

Flow around a rotating cross

We present a preliminary computation obtained with the combination of the mesh adaptation, the Chimera approximation and the DDES model. They are used for simulating the mixing process obtained with the rotation of a cross in a cylindric box. The Reynolds number is 1.8 million, based on the thickness of the blades. The number of vertices of the initial mesh is 10.466. After two turns, and five mesh adaptations, the number of vertices of the adapted mesh is 1.05 million. In Fig. 12a, we present a view of this final adapted mesh. Mesh size normal to the cylinder is yet not small enough to allow the creation of strong non-2D features. This is verified from the examination of the Q-factor, Fig. 12b. Fig. 13 shows details of mesh and velocity magnitude. It seems from this first calculation that despite the hign Reynold number, the flow separates at blade angle is a rather stable mode.

Conclusion

Several hybrid strategies, based on the DDES, RANS and DVMS models, are first evaluated for the simulation of the supercritical flow around a circular cylinder. This benchmark is characterized by turbulent boundary-layer separation. The cylinder flow is first computed at Reynolds 10 6 . The predictions of the main flow parameters are in overall good agreement with experimental data and the results of simulations in the literature, especially the RANS/DVMS model. A second series of computations of the flow past a circular cylinder is carried out at Reynolds 2 × 10 6 .

The predictions obtained by the hybrid models used to perform this benchmark are a little less accurate compared to the previous Reynolds number, more specifically for the mean separation angle which is overestimated. Again, the RANS/DVMS model is the hybrid model among those applied which gives the best results overall. It is also to be noted that the pressure distribution is well predicted by the hybrid models for both Reynolds numbers, especially by the RANS/DVMS model. Then, the DDES formulation is applied to the flow around a rotating cross inside a cylinder with a Reynolds number of 1.8 × 10 6 . Preliminary results using mesh adaptation with meshes of 1 million vertices show that separations are more stable than in the flows around the cylinder. RANS/DVMS computations are planned. The new tool is also being applied to other rotating devices like propellers and will produce mesh adapted results for these geometries in a near future.

  by replacing the D RAN S k = ρε dissipation term in the k equation by D DDES k = ρ k 3/2

Figure 1 -

 1 Figure 1 -Building the VMS coarse level

Figure 2 -

 2 Figure 2 -VMS principle.

Figure 3 -

 3 Figure 3 -Isocontours of the blending function for the circular cylinder benchmark (Reynolds 10 6 ) : for θ = 1 (in red) the RANS model is activated, wherever 0 < θ < 1 additional resolved fluctuations are computed, and in the limit θ → 0 (in blue) the LES approach (DVMS in this work) is recovered.

  fixed domain in red and a rotating domain in blue. For computations we start by locating the boundary nodes of domain 1 in domain 2, and reciprocally we locate the boundary nodes of domain 2 in domain 1, then the aerodynamic values of each boundary node are determined by interpolation and each domain performs its calculation with the new interpolated values.

Figure 4 -

 4 Figure 4 -Definition of rotor and stator from an initial mesh : in red, the stator, in blue the rotor, and in gray, the overlap region.

Figure 5 -

 5 Figure 5 -Computational domain

Figure 6 -

 6 Figure 6 -Distribution over the cylinder surface of the mean pressure coefficient. Comparaison between experimental data and numerical results at Reynods 10 6 .

Figure 7 -Figure 8 -

 78 Figure 7 -Circular cylinder, Reynolds 10 6 : instantaneous isocontours of the vorticity magnitude. From top to bottom : RANS, DDES, DDES/DVMS and RANS/DVMS.

Figure 9 -

 9 Figure 9 -Instantaneous vorticity magnitude. From top to bottom : DDES, DDES/DVMS and RANS/DVMS.

Figure 10 -

 10 Figure 10 -Distribution over the cylinder surface of the mean pressure coefficient. Comparaison between experimental data and numerical results at Reynods 2 × 10 6 .

Figure 11 -

 11 Figure 11 -Lift and drag coefficient fluctuation. From top to bottom : DDES/DVMS and RANS/DVMS.

Inria

Figure 12 -

 12 Figure12-DDES computation of a rotating mixing device with a mesh adaptative Chimera approach (horizontal cut). View of the mesh and of Q-factor (colored with velocity magnitude) after the vertical upper blade.

Figure 13 -

 13 Figure13-DDES computation of a rotating mixing device with a mesh adaptative Chimera approach (horizontal cut). Partial view of the mesh and velocity magnitude.

Table des matières

 des RANS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 DDES model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 High Reynolds computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 LES-like model : DVMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 SGS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 DVMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Hybrid RANS/DVMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Hybrid DDES/DVMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	NOMENCLATURE
	1 Introduction W	[-]	flow variables
	2 Turbulence modeling ρ [kg/m 3 ]	density
	2.1 5 Conclusion	u E k ε W h W h F τ RAN S [-] [m/s] [J/m 3 ] [m 2 /s 2 ] [m 2 /s 3 ] [-] [-] [-] [-] τ LES [-] τ DDES [-] ν [m 2 /s] ν t [m 2 /s] C d [-] C l [-] C p b [-] θ [deg.] ∆ [m]	velocity vector total energy per unit volume turbulence kinetic enery dissipation rate of k discrete flow variables small resolved scales of W h Reynolds average convective and diffusive fluxes RANS closure term LES closure term DDES closure term kinematic viscosity of the fluid turbulent kinematic viscosity mean drag coefficient r.m.s. of the lift coefficient mean base pressure coefficient mean separation angle LES filter width
		l RAN S	[m]	RANS turbulence length scale
		µ SGS	[kg/m/s] LES eddy viscosity
		µ RAN S [kg/m/s] RANS eddy viscosity
		τ	[-]	Reynolds stress tensor
		Sets
		Ω f		Fluid domain
		Ω h		Collection of Tetrahedron
		J k		Index set wich are contains in macro cell C k
		C k		Agglomeration cell
		V 2h
	RR n°9464

3 Numerical discretization 3.1 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Mesh adaptation for rotating machines . . . . . . . . . . . . . . . . . . . . . . . . 4 Applications 4.1 Flow past a cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Flow around a rotating cross . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 -

 1 Bulk coefficient of the flow around a circular cylinder at Reynolds number 1M, C d holds for the mean drag coefficient, C l is the root mean square of lift time fluctuation, C p b is the pressure coefficient at cylinder basis, L r is the mean recirculation lenght, θ is the mean separation angle.

		Mesh size y + w	y + m	C d	C l	-C pb	Lr	θ	St
	Present simulation
	RANS k -ε	4.8M	20 100	0.23	0.06	0.23	0.69 132 0.46
	DDES k -ε Goldberg	4.8M	20 100	0.20	0.04	0.22	0.87 138 0.13
	DDES/ DVMS
	k -ε / cubic Smagorinsky	4.8M	20 100	0.20	0.02	0.22	0.82 135 0.42
	k -ε / cubic WALE	4.8M	1	100	0.20	0.02	0.26	0.80 132 0.58
	RANS / DVMS
	k -ε / cubic Smagorinsky	4.8M	20 100	0.24	0.05	0.22	0.62 133 0.42
	k -ε / cubic Smagorinsky	4.8M	1	100	0.25	0.09	0.25	0.64 132 0.46
	k -ε / cubic WALE	4.8M	1	100	0.26	0.11	0.22	0.65 134 0.42
	Other simulations
	SRANS Catalano [5]	2.3M	-	-	0.39	-	0.33	-	-	-
	SLES Catalano [5]	2.3M	-	-	0.31	-	0.32	-	-	-
	URANS Catalano [5]	2.3M	-	-	0.40	-	0.41	1.37	-	0.31
	LES Kim [18]	6.8M	-	-	0.27	0.12	0.28	-	108	-
	Measurements
	Exp. [30]				0.24	-	0.33	-	-	-
	Exp. [29]				0.22	-	-	-	-	-
	Exp. [35]				0.25	-	0.32	-	-	-
	Exp. [9]				-	-	-	-	130	-
	Exp. [38]				0.2-0.4 0.1-0.15 0.2-0.34	-

Table 2 -

 2 Bulk coefficient of the flow around a circular cylinder at Reynolds number 2 × 10 6 . The subscript S holds for Smagorinsky SGS model.

RR n°9464

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

This work was supported by the ANR NORMA project, grant ANR-19-CE40-0020-01 of the French National Research Agency. The authors gratefully acknowledge GENCI for granted access to HPC resources through CINES (grants 2020-A0092A05067 and 2021-A0102A06386) and IDRIS (grant 2021-A0112A05067).

This work was supported by the ANR project NORMA under grant ANR-19-CE40-0020-01.