Chapter 5 condli

0.1 Fonction

Il est nécessaire de disposer de conditions aux limites dans au moins trois cas principaux :

calcul de termes convectifs (différentiels d’ordre un en espace) au bord : le calcul fait intervenir un flux au bord et demande la donnée au bord de la variable convectée lorsque la frontière est “entrante” au sens des courbes caractéristiques du système (au sens de la vitesse, dans le cas d’une équation isolée portant sur un simple scalaire : interprétation suffisante dans le cadre actuel de Code Saturne¹) ;
calcul de termes de diffusion (différentiels d’ordre deux en espace) : il est alors nécessaire de disposer d’un moyen de connaître la valeur de dérivées spatiales d’ordre un au bord (plus exactement, il est nécessaire de disposer d’information permettant de calculer les termes qui en dépendent, comme les contraintes ou les flux thermiques en paroi) ;
calcul de gradient “cellules” plus standard : il est nécessaire de disposer de la donnée de la variable aux faces de bord (plus généralement, il est nécessaire de disposer d’un moyen de calculer les termes discrets des équations lorsqu’ils font intervenir le gradient dans les cellules de bord, comme par exemple les termes de gradient transposé des équations de Navier-Stokes).

Les considérations présentes concernent uniquement les variables de calcul (vitesse, pression, tenseur de Reynolds, scalaires solution d’une équation de convection-diffusion). Pour ces grandeurs², l’utilisateur doit définir les conditions aux limites en chaque face de bord (usclim ).

Le sous-programme condli permet de traduire les données utilisateur (fournies dans usclim ) au format interne de représentation des conditions aux limites. Des vérifications de complétude et de cohérence sont également menées (dans vericl ). Sont en particulier traitées les conditions aux limites de paroi (clptur ) et de symétrie pour les vitesses et le tenseur de Reynolds (clsyvt ).

Le sous-programme condli fournit en sortie des couples de coefficients A_b et B_b pour chaque variable f et chaque face de bord. Ils sont utilisés pour le calcul des termes discrétisés intervenant dans les équations à résoudre et permettent en particulier de déterminer une valeur de face de bord f_b,int (localisée au “centre” de la face de bord, barycentre de ses sommets) par la relation f_b,int = A_b+B_b f_I′ où f_I′ est la valeur de la variable au point I′, projeté du centre de la cellule jouxtant le bord sur la droite normale à la face de bord et passant par son centre (voir la figure II.0.1).

[height=8cm]./graphics/fluxbord

Figure II.0.1: Cellule de bord.

0.2 Discrétisation

• Notation
On désignera dans la suite par VarScalaire toute variable

- autre que la vitesse, la pression, les grandeurs turbulentes k, ε, R_ij, ϕ, f et ω,
- solution d’une équation de convection-diffusion.

La dénomination VarScalaire pourra en particulier désigner la température, un scalaire passif, une fraction massique ou (sauf mention contraire explicite) la variance des fluctuations d’une autre VarScalaire. Les variables d’état déduites (masse volumique, viscosité...) ne seront pas désignées par VarScalaire.

• Représentation des conditions aux limites standard dans usclim
Des conditions aux limites standardisées peuvent être fournies par l’utilisateur dans usclim . Il est pour cela nécessaire d’affecter un type aux faces de bord des cellules concernées³. Les conditions prévues par défaut sont les suivantes :

• Entrée : correspond à une condition de Dirichlet sur toutes les variables transportées (vitesse, variables turbulentes, VarScalaires...), et à une condition de Neumann homogène (flux nul) sur la pression.
• Sortie :
- - lorsque le flux de masse est effectivement dirigé vers l’extérieur du domaine, ce choix correspond à une condition de Neumann homogène sur toutes les variables transportées et à ∂²P/∂ n∂τ_i=0, pris en compte sous forme de Dirichlet pour la pression (n et (τ_i)_{i ∈ {1,2}} désignent respectivement le vecteur normal de la face de sortie considérée et deux vecteurs normés, orthogonaux entre eux et dans le plan de la face de sortie). Cette condition est appliquée de manière explicite en utilisant le champ de pression et son gradient au pas de temps précédent. En outre, la pression étant définie à une constante près, elle est recalée en un point de sortie pour y conserver la valeur P0 (on évite ainsi toute dérive vers des valeurs très grandes relativement à l’écart maximal de pression sur le domaine)⁴.
- - lorsque le flux de masse est dirigé vers l’intérieur du domaine, situation peu souhaitable a priori⁵, on impose une condition de Dirichlet homogène sur la vitesse (pas sur le flux de masse), à défaut de connaître sa valeur en aval du domaine. La pression est traitée comme dans le cas précédent où le flux de masse est dirigé vers l’extérieur du domaine. Pour les variables autres que la vitesse et la pression, deux cas se présentent :
  - - dans le cas de sorties de type “10” (ISOR10 dans usclim ), on impose une condition de Dirichlet (valeur fournie par l’utilisateur et nulle par défaut) pour représenter la valeur du scalaire introduit dans le domaine par les faces de bord concernées.
  - - dans le cas des sorties de type “9” (ISOR09 dans usclim ), on impose, comme lorsque le flux de masse est sortant, une condition de Neumann homogène (ceci n’est pas une situation souhaitable, puisque l’information portée sur les faces de bord provient alors de l’aval de l’écoulement local).
• Paroi : on se reportera à clptur pour une description du traitement des conditions aux limites de paroi (supposées imperméables au fluide). Brièvement, on peut dire ici qu’une approche par lois de paroi est utilisée afin d’imposer la contrainte tangentielle sur la vitesse. La paroi peut être défilante⁶. Les VarScalaires reçoivent par défaut une condition de Neumann homogène (flux nul). Si l’on souhaite imposer une valeur en paroi pour ces variables (par exemple, dans le cas d’une paroi à température imposée) une loi de similitude est utilisée pour déterminer le flux au bord en tenant compte de la couche limite. Dans le cas des couplages avec , Code Saturnereçoit une température de paroi et fournit un flux thermique. La condition de pression standard est une condition de Neumann homogène⁷
• Symétrie : correspond à des conditions de Neumann homogènes pour les grandeurs scalaires et à des conditions de symétrie classiques pour les vecteurs (vitesse) et les tenseurs (tensions de Reynolds) : voir clsyvt .

• Représentation des conditions aux limites spécifiques dans usclim
On a vu que l’affectation à une face de bord d’un type standard (entrée, sortie, paroi, symétrie) permettait d’appliquer simplement à l’ensemble des variables un assortiment de conditions aux limites cohérentes entre elles pour les types usuels de frontière physique.

Une solution consiste à définir dans usclim , pour chaque face de bord et chaque variable, des conditions aux limites spécifiques⁸ (celles-ci, comme les conditions standards, se ramènent finalement à des conditions de type mixte).

Les deux approches ne sont pas nécessairement incompatibles et peuvent même se révéler complémentaires. En effet, les conditions aux limites standards peuvent être surchargées par l’utilisateur pour une ou plusieurs variables données. Il convient cependant de s’assurer que, d’une façon ou d’une autre, une condition à la limite a été définie pour chaque face de bord et chaque variable.

Des conditions de compatibilité existent également entre les différentes variables (voir vericl ):

- en entrée, paroi, symétrie ou sortie libre, il est important que toutes les composantes de la vitesse aient le même type de condition ;
- lorsque la vitesse reçoit une condition de sortie, il est important que la pression reçoive une condition de type Dirichlet. Pour plus de détails, on se reportera au paragraphe relatif à la condition de sortie pour la pression, page ?? ;
- lorsqu’une des variables de vitesse ou de turbulence reçoit une condition de paroi, il doit en être de même pour toutes ;
- lorsqu’une des composantes R_ij reçoit une condition de symétrie, il doit en être de même pour toutes ;
- lorsqu’une VarScalaire reçoit une condition de paroi, la vitesse doit avoir le même type de condition.

• Représentation interne des conditions aux limites

Objectif

Les conditions fournies par l’utilisateur sont retraduites sous forme de couples de coefficients A_b et B_b pour chaque variable f et chaque face de bord. Ces coefficients sont utilisés pour le calcul des termes discrets intervenant dans les équations à résoudre et permettent en particulier de déterminer une valeur de face de bord f_b,int. Il est important d’insister dès à présent sur le fait que cette valeur est, de manière générale, une simple valeur numérique qui ne reflète pas nécessairement une réalité physique (en particulier aux parois, pour les grandeurs affectées par la couche limite turbulente). On détaille ci-dessous le calcul de A_b, B_b et de f_b,int.

Notations

- On considère l’équation (??) portant sur le scalaire f, dans laquelle ρ représente la masse volumique, u la vitesse, α la conductivité et S les termes sources additionnels. C est défini plus bas.
- Le coefficient α représente la somme des conductivités moléculaire et turbulente (selon les modèles utilisés), soit α=α_m+α_t, avec, pour une modélisation de type viscosité turbulente, α_t=C µ_t/σ_t, où σ_t est le nombre de Prandtl turbulent⁹.
- À titre d’exemple, on précise dans le tableau 5.1 la valeur et l’unité de α pour quelques cas particuliers de f (certains termes de l’équation peuvent alors disparaître : pour la pression, en particulier, l’équation est stationnaire et les termes convectifs sont absents).
- Le coefficient C_p représente la chaleur spécifique, d’unité m² s⁻² K⁻¹=J kg⁻¹ K⁻¹.
- On note λ la conductivité thermique, d’unité kg m s⁻³ K⁻¹=W m⁻¹ K⁻¹.
- Il convient de préciser que C=1 pour toutes les variables hormis pour la température, cas dans lequel on a¹⁰ C=C_p. Dans le code, c’est la valeur de α_m/C que l’utilisateur doit fournir (si la propriété est constante, les valeurs sont affectées dans usini1 à VISCL0 pour la vitesse et à VISLS0 pour les VarScalaires ; si la propriété est variable, ce sont des tableaux équivalents qui doivent être renseignés dans usphyv ).
- Pour la variance des fluctuations d’une VarScalaire, la conductivité α et le coefficient C sont hérités de la VarScalaire associée.

f α

symbole nom unité symbole nom unité

u_i vitesse m s⁻¹ µ ou µ+µ_t viscosité dynamique kg m⁻¹ s⁻¹

p pression kg m⁻¹ s⁻² Δ t pas de temps s

T température K λ conductivité thermique kg m s⁻³ K⁻¹

=W m⁻¹ K⁻¹

H enthalpie m² s⁻² λ/C_p conductivité thermique/C_p kg m⁻¹ s⁻¹

=J kg⁻¹

f VarScalaire unité(f) α conductivité ou diffusivité kg m⁻¹ s⁻¹

Table 5.1: Valeurs et unités de α dans quelques cas particuliers usuels.

Condition de type Dirichlet simple : lorsque la condition est une condition de Dirichlet simple, on obtient naturellement (cas particulier de (??)) : Autres cas : lorsque la condition à la limite porte sur la donnée d’un flux, il s’agit d’un flux diffusif¹¹. On a alors :

Le flux diffusif réel imposé peut être donné

- directement (condition de Neumann), soit φ_réel=φ_imp,ext ou
- déduit implicitement de deux informations imposées : une valeur externe f_imp,ext et un coefficient d’échange h_imp,ext (condition de Dirichlet généralisée).

Selon le type de condition (Dirichlet ou Neumann) et en prenant pour hypothèse la conservation du flux dans la direction normale au bord, on peut alors écrire (voir figure II.0.1) :

Le rapport entre le coefficient h_b et le coefficient h_int rend compte de l’importance de la traversée de la zone proche du bord et revêt une importance particulière dans le cas des parois le long desquelles se développe une couche limite (dont les propriétés sont alors prises en compte par h_b : se reporter à clptur ). Dans le cadre plus simple considéré ici, on se limitera au cas h_b=h_int et f_b,ext=f_b,int=f_b. La relation (??) s’écrit alors :

En réarrangeant, on obtient la valeur de bord :

Conclusion : on notera donc les conditions aux limites de manière générale sous la forme : avec A_b et B_b définis selon le type des conditions :

Remarques

- La valeur f_I′ est calculée en utilisant le gradient cellule de f, soit : f_I′=f_I+II′grad f_I.
- Il reste à préciser la valeur de h_int. Il s’agit d’une valeur numérique, n’ayant a priori aucun rapport avec un coefficient d’échange physique, et dépendante du mode de calcul du flux diffusif dans la première maille de bord. Ainsi h_int=α/I′F (l’unité s’en déduit naturellement).
- On rappelle que dans le code, c’est la valeur de α_m/C que l’utilisateur doit fournir. Si la propriété est constante, les valeurs sont affectées dans usini1 à VISCL0 pour la vitesse (viscosité dynamique moléculaire µ en kg m⁻¹ s⁻¹) et à VISLS0 pour les VarScalaires (par exemple, pour la température et l’enthalpie, λ/C_p en kg m⁻¹ s⁻¹). Si la propriété est variable en espace ou en temps, ce sont des tableaux équivalents qui doivent être renseignés dans usphyv . En outre, la variance des fluctuations d’une VarScalaire hérite automatiquement la valeur de α_m/C de la VarScalaire associée (Code Saturne 1.1 et suivantes).
- On rappelle également, car ce peut être source d’erreur, que dans le code, on a :
- - pour la température α_m=λ et C=C_p
- - pour l’enthalpie α_m=λ/C_p et C=1

Exemples de cas particuliers

- Dans le cas d’une condition de Dirichlet, l’utilisateur est donc conduit à fournir deux données : f_imp,ext et h_imp,ext. Pour obtenir une condition de Dirichlet simple (sans coefficient d’échange) il suffit d’imposer h_imp,ext=+∞. C’est le cas d’utilisation le plus courant (en pratique, h_imp,ext=10³⁰ ).
- Dans le cas d’une condition de Neumann, l’utilisateur fournit une seule valeur φ_imp,ext (nulle pour les conditions de Neumann homogènes).

Valeur et unité des données à fournir
Le tableau 5.2 permet d’identifier les unités de h (h_imp,ext ou h_int) et φ_imp,ext pour quelques variables classiques. La variable d représente une longueur (égale, pour h_int, à I′F). Il est également important de noter qu’un flux “sortant” doit être positif.

f h

symbole nom unité homogène à unité

u_i vitesse m s⁻¹ (µ+µ_t)/d kg m⁻² s⁻¹

p pression kg m⁻¹ s⁻² (Δ t)/d s m⁻¹

T température K (λ+C_pµ_t/σ_t)/d kg s⁻³ K⁻¹

=W m⁻² K⁻¹

H enthalpie m² s⁻²=J kg⁻¹ (λ/C_p+µ_t/σ_t)/d kg m⁻² s⁻¹

f scalaire unité(f) α/d kg m⁻² s⁻¹

f φ_imp,ext

symbole nom unité homogène à unité

u_i vitesse m s⁻¹ (µ+µ_t) grad (u_i).n kg m⁻¹ s⁻²

p pression kg m⁻¹ s⁻² (Δ t) grad (p).n kg m⁻² s⁻¹

T température K (λ+C_pµ_t/σ_t) grad (T).n kg s⁻³=W m⁻²

H enthalpie m² s⁻²=J kg⁻¹ (λ/C_p+µ_t/σ_t) grad (H).n kg s⁻³=W m⁻²

f scalaire unité(f) α grad (f).n kg m⁻² s⁻¹ unité(f)

Table 5.2: Valeur et unité de h et φ_imp,ext dans quelques cas particuliers usuels.

• Condition de sortie pour la pression
On précise ici la condition de sortie appliquée à la pression dans le cas des sorties standards. Il est nécessaire d’imposer une condition de type Dirichlet (accompagnée d’une condition de type Neumann homogène sur les composantes de la vitesse). On la calcule à partir des valeurs de la variable au pas de temps précédent.

- En raisonnant sur une configuration simple (de type canal, avec une sortie plane, perpendiculaire à l’écoulement), on peut faire l’hypothèse que la forme des profils de pression pris sur les surfaces parallèles à la sortie est inchangée aux alentours de celle-ci (hypothèse d’un écoulement établi, loin de toute perturbation). Dans cette situation, on peut écrire ∂²P/∂n∂τ_i=0 (n est le vecteur normal à la sortie, τ_i représente une base du plan de sortie).
- Si, de plus, on peut supposer que le gradient de pression pris dans la direction perpendiculaire aux faces de sortie est uniforme au voisinage de celle-ci, le profil à imposer en sortie (valeurs p_b) se déduit du profil pris sur un plan amont (valeurs p_amont) en ajoutant simplement la constante R=d grad (p).n (où d est la distance entre le plan amont et la sortie), soit p_b=p_amont+R (le fait que R soit identique pour toutes les faces de sortie est important afin de pouvoir l’éliminer dans l’équation (??) ci-dessous).
- Avec l’hypothèse supplémentaire que les points I′ relatifs aux faces de sortie sont sur un plan parallèle à la sortie, on peut utiliser les valeurs en ces points (p_I′) pour valeurs amont soit p_amont=p_I′=p_I+II′.grad p.
- Par ailleurs, la pression étant définie à une constante près (écoulement incompressible) on peut fixer sa valeur arbitrairement en un point A (centre d’une face de sortie choisie arbitrairement¹²) à p₀ (valeur fixée par l’utilisateur, égale à P0 et nulle par défaut), et donc décaler le profil imposé en sortie en ajoutant :
R₀=p₀−(p_amont,A+R)=p₀−(p_I′,A+R).
- On obtient donc finalement :

On constate donc que la condition de pression en sortie est une condition de Dirichlet dont les valeurs sont égales aux valeurs de la pression (prises au pas de temps précédent) sur le plan amont des points I′ et recalées pour obtenir P0 en un point de sortie arbitraire.

0.3 Mise en œuvre

• Mode de prescription des conditions aux limites dans usclim
L’utilisateur fournit pour chaque face un moyen de déterminer les conditions aux limites de toutes les variables de calcul. Comme on l’a vu, plusieurs méthodes sont envisageables.

La plus simple consiste à renseigner, pour chaque face, un code désignant une condition à la limite standard (dans le tableau ITYPFB de dimension NFABOR) et les informations complémentaires éventuelles. Les valeurs de ITYPFB suivantes sont prédéfinies¹³ : IENTRE (entrée), ISOR09 (sortie de type 9), ISOR10 (sortie de type 10), ISYMET (symétrie), IPAROI (paroi). Dans le cas des entrées (et des sorties de type 10), il est nécessaire de fournir une valeur de Dirichlet.

L’utilisateur peut également renseigner un entier (tableau ICODCL de dimension NFABOR×NVAR) désignant le type de condition à appliquer à une variable donnée : les valeurs 1 (pour Dirichlet) et 3 (pour Neumann) sont souvent utilisées, plus rarement 5 (pour les conditions de paroi) et exceptionnellement 4 (pour la symétrie d’une vitesse ou du tenseur de Reynolds), 9 ou 10 (pour les sorties de type 9 ou 10). Lorsque ICODCL est renseigné pour une variable donnée, il prend alors le pas sur ITYPFB sur la face considérée. L’utilisateur doit dans ce cas également fournir, suivant les cas, zéro, un ou deux réels (dans le tableau RCODCL de dimension NFABOR×NVAR× 3). Pour une face et une variable f données, les 3 valeurs de RCODCL désignent respectivement les grandeurs f_imp,ext, h_imp,ext et φ_imp,ext.

On indique dans le tableau 5.3 le mode de traitement des variables pour les différents types de conditions aux limites standards. La liste des valeurs complémentaires à fournir est également précisée. Le tableau 5.4 propose le même type d’information pour les conditions aux limites spécifiques plus complexes, définies à partir de la valeur de ICODCL. Enfin, le tableau 5.5 synthétise les valeurs admissibles de ICODCL pour chaque variable de calcul.

ITYPFB type variables traitées type de condition données complémentaires

IENTRE entrée variables transportées Dirichlet RCODCL(.,.,1) valeur d’entrée nulle par défaut

pression Neumann homogène

ISYMET symétrie vitesse symétrie vecteur

tensions de Reynolds symétrie tenseur

autres variables Neumann homogène

IPAROI paroi pression Neumann homogène^a

vitesses lois de paroi (Dirichlet)^b RCODCL(.,.,1) vitesse de défilement nulle par défaut

grandeurs turbulentes lois de paroi^c

autres variables

transportées Neumann homogène^d

ISOR09 sortie pression ∂²P/∂ n∂τ_i=0

vitesses Neumann homogène

(flux de masse sortant)

Dirichlet homogène

(flux de masse entrant)

autres variables

transportées Neumann homogène

ISOR10 sortie pression ∂²P/∂ n∂τ_i=0

vitesses Neumann homogène

(flux de masse sortant)

Dirichlet homogène

(flux de masse entrant)

autres variables Neumann homogène

transportées (flux de masse sortant)

Dirichlet^e RCODCL(.,.,1) valeur d’entrée nulle par défaut

(flux de masse entrant)

^aLe gradient peut être extrapolé en paroi selon les valeurs de EXTRAP : voir gradrc
^bvoir clptur
^cvoir également clptur
^dLa condition de Neumann homogène est obtenue lorsque l’utilisateur ne précise rien d’autre que le type de frontière IPAROI ; d’autres conditions sont possibles : le tableau 5.4 les précise.
^eLa valeur est imposée par l’utilisateur. Il n’y a pas d’implicitation temporelle.
Table 5.3: Conditions aux limites standards.

ICODCL variables type de condition données complémentaires nécessaires et suffisantes

1 toutes Dirichlet RCODCL(.,.,1) f_imp,ext valeur nulle par défaut

RCODCL(.,.,2) h_imp,ext coef. d’échange =RINFIN=10³⁰ par défaut

3 toutes Neumann RCODCL(.,.,3) φ_imp,ext valeur nulle par défaut

4 vitesse Symétrie vecteur

R_ij Symétrie tenseur

5 vitesse Loi de paroi RCODCL(.,.,1) vitesse de défilement nulle par défaut

k, R_ij, ε, ϕ, f, ω Loi de paroi

VarScalaire Loi de paroi RCODCL(.,.,1) f_imp,ext valeur en paroi nulle par défaut

(sauf variance) RCODCL(.,.,2) h_imp,ext coef. d’échange =RINFIN=10³⁰ par défaut

9 vitesse Neumann homogène

(flux de masse sortant)

Dirichlet homogène

(flux de masse entrant)

10 k, R_ij, ε, ϕ, f, ω Dirichlet RCODCL(.,.,1) f_imp,ext valeur nulle par défaut

VarScalaire (flux de masse entrant) RCODCL(.,.,2) h_imp,ext coef. d’échange =RINFIN=10³⁰ par défaut

Neumann RCODCL(.,.,3) φ_imp,ext valeur nulle par défaut

(flux de masse sortant)

Table 5.4: Conditions aux limites spécifiques.

Variable Valeurs de ICODCL admissibles

Vitesse U 1 3 4 5 9

Pression p 1 3

Variable scalaire de turbulence k, ε, ϕ, f, ω 1 3 5 10

Tenseur de Reynolds R_ij 1 3 4 5 10

VarScalaire (hormis variances) 1 3 5 10

Variance des fluctuations d’une VarScalaire 1 3 10

Table 5.5: Valeurs admissibles de ICODCL pour chaque variable.

Il est important de connaître également les compatibilités à assurer entre les valeurs de ICODCL associées aux différentes variables (Cf. Points à traiter).

- Si ICODCL vaut 4, 5 ou 9 pour une composante de la vitesse, la même valeur de ICODCL doit être associée à toutes les composantes (symétrie, paroi ou sortie libre).
- Si ICODCL vaut 9 pour une composante de la vitesse, la valeur de ICODCL associée à la pression doit être 1 (sortie libre).
- Si ICODCL vaut 5 pour une composante de la vitesse, pour k, ε, ϕ, f, ω ou pour une des composantes du tenseur de Reynolds R_ij, la valeur de ICODCL associée à toutes ces variables doit être 5 (paroi).
- Si ICODCL vaut 4 pour une composante de la vitesse ou pour une des composantes du tenseur de Reynolds R_ij, la valeur de ICODCL associée à toutes ces variables doit être 4 (symétrie).
- Si ICODCL vaut 5 pour une variable VarScalaire, la valeur de ICODCL associée à toutes les composantes de la vitesse doit être 5 (paroi).

• Conversion des données utilisateur
Une première étape (TYPECL) permet essentiellement de convertir les données des utilisateurs. Plus précisément, les actions suivantes sont réalisées :

- les faces de bord sont triées selon le type de condition à la limite qui leur a éventuellement été affecté (valeurs de ITYPFB) ;
- les types ITYPFB (si l’utilisateur en a défini) sont convertis en quadruplets définis pour chaque face IFAC et chaque variable IVAR (hormis lorsque l’utilisateur a affecté une valeur à ICODCL) :
(ICODCL, RCODCL(IFAC,IVAR,1), RCODCL(IFAC,IVAR,2), RCODCL(IFAC,IVAR,3))
- pour la pression, la condition de sortie se traduit par une condition de type Dirichlet : on impose aux faces de bord la valeur de la pression p au point I′, calculée à partir des données du pas de temps précédent (et le calcul du gradient de pression est donc réalisé au préalable). Cette valeur p_I′ est conservée dans le tableau COEFU(.,1), de dimension NFABOR, puis affectée à RCODCL(.,IPR,1) accompagnée d’un recalage (identique pour toutes les faces) destiné à assurer une valeur de Dirichlet de pression fixe P0 sur la première face de sortie rencontrée (Cf. Points à traiter).
- L’indicateur IDIRCL (tableau de dimension NVAR) est annulé pour repérer les variables pour lesquelles le système diffusif à résoudre n’est pas nécessairement inversible (le problème sera traité dans matrix ). Il s’agit en pratique, dans le cadre actuel, de la variable pression lorsqu’il n’y a aucune face de sortie (cavité fermée par exemple). Plus précisément, IDIRCL est annulé pour les variables qui n’ont aucune condition de Dirichlet¹⁴ et pour lesquelles l’indicateur ISTAT est positionné à zéro. Cet indicateur (de valeur 0 ou 1) est, dans le bilan explicite, le coefficient multiplicatif du terme de dérivée temporelle de la variable résolue ; il est nul pour la pression.
- la valeur du flux de masse est imprimée pour les ensembles de faces de bord de même type ITYPFB (il existe un type par défaut appliqué aux faces de bord pour lesquelles l’utilisateur n’a pas renseigné ITYPFB).

• Vérifications
Une vérification des conditions aux limites est ensuite menée (vericl ) (l’utilisateur est prévenu et le programme interrompu en cas de problème).

- Le sous-programme s’assure que toutes les faces de bord ont reçu une condition à la limite pour toutes les variables.
- L’admissibilité des conditions aux limites (selon la nature de la variable à laquelle elles sont appliquées) est vérifiée.
- La cohérence des conditions aux limites entre les différentes variables est vérifiée (voir les tableaux 5.4 et 5.5).

• Calculs préliminaires
Différentes grandeurs sont ensuite construites pour préparer le traitement ultérieur des conditions aux limites.

- quand la connaissance de la distance à la paroi est nécessaire (R_ij−ε avec termes d’écho de paroi, modèle LES de Smagorinsky avec amortissement de van Driest, modèle k−ω SST) et que la méthode de calcul “ancienne” a été choisie (|ICDPAR|=2, non compatible avec le parallélisme et la périodicité), pour chaque phase IPHAS, on détermine, pour chaque cellule IEL, le numéro de la face de paroi la plus proche (cette information est stockée dans IA(IIFAPA(IPHAS)-1+IEL)) ; cette distance est inutile pour le traitement des conditions aux limites, mais condli est le sous-programme dans lequel la réalisation de ce calcul est la plus simple.
- en cas de couplage avec , on détermine, pour le scalaire couplé f (habituellement la température), les valeurs f_I′ dans les cellules de bord au moyen de la formule f_I′=f_I+II′grad f_I (ou au moyen de f_I′=f_I si on traite le premier pas de temps¹⁵ ou que la reconstruction a été desactivée par ITBRRB=0¹⁶). Ces valeurs sont stockées dans TBORD (tableau de dimension NFABOR) qui est utilisé en sortie de condli pour transmission d’information à .
- si des faces portent des conditions de symétrie ou de paroi, on construit les composantes de la vitesse u_j,I′ dans les cellules de bord (valeurs stockées dans le tableau local COEFU de dimension NFABOR× 3) au moyen de la formule u_j,I′=u_j,I+II′grad u_j_I (ou au moyen de u_j,I′=u_j,I au premier pas de temps) ; en outre, si le modèle de turbulence est le R_ij−ε, on construit les composantes du tenseur de Reynolds dans les cellules de bord (valeurs stockées dans le tableau local RIJIPB de dimension NFABOR× 6) de la même manière. Les deux tableaux COEFU et RIJIPB sont utilisés dans clptur et clsyvt .

• Conditions aux limites de paroi
On détermine alors (clptur ) les conditions aux limites de paroi pour la vitesse et les grandeurs turbulentes. Les VarScalaires (excepté les variances) peuvent également recevoir un traitement particulier prenant en compte la couche limite si la condition qui leur est appliquée est de type Dirichlet (ICODCL=5). On traite donc en particulier dans clptur les conditions portant sur la température lorsque la paroi est à température imposée. Par contre, les conditions de flux sont traitées plus tard dans condli .

Le sous-programme clptur remplit les tableaux COEFA et COEFB pour les faces de paroi et les variables traitées. Ces tableaux représentent les coefficients A_b et B_b. Plus précisément, pour une face IFAC et une variable IVAR donnée, les valeurs renseignées sont COEFA(IFAC,ICLVAF) et COEFB(IFAC,ICLVAF) avec ICLVAF=ICLRTP(IVAR,ICOEFF). Deux autres valeurs sont également renseignées :
COEFA(IFAC,ICLVAR) et COEFB(IFAC,ICLVAR) avec ICLVAR=ICLRTP(IVAR,ICOEF).
Elles ne diffèrent des précédentes qu’en k−ε et représentent des coefficients A_b et B_b particuliers pour le calcul des gradients intervenant dans le terme de production turbulente (il y a donc bien deux types de conditions aux limites pour la vitesse dans ce cas, et chacun est utilisé lors du traitement d’un terme particulier : on se reportera à clptur ).

• Conditions aux limites de symétrie
On détermine également (clsyvt ) les conditions aux limites de symétrie pour les vecteurs (vitesse) et les tenseurs (de Reynolds). Elles nécessitent des rotations pour s’exprimer dans le repère lié au bord et le traitement lourd qui en résulte a donc été encapsulé. Le sous-programme complète les tableaux COEFA et COEFB aux faces de symétries pour la vitesse et les tensions de Reynolds (contrairement aux faces de paroi en k−ε, les valeurs COEFA(IFAC,ICLVAF) et COEFB(IFAC,ICLVAF) sont égales à COEFA(IFAC,ICLVAR) et COEFB(IFAC,ICLVAR)).

• Autres conditions aux limites pour la vitesse
On détermine, pour la vitesse, les conditions aux limites en sortie. On traite en outre les conditions restantes de type Dirichlet et Neumann. Ici également, on fournit une valeur de COEFA et de COEFB pour chaque face concernée (les valeurs COEFA(IFAC,ICLVAF) et COEFB(IFAC,ICLVAF) sont égales à COEFA(IFAC,ICLVAR) et COEFB(IFAC,ICLVAR)). Notons que, si le flux de masse est entrant à travers une face de sortie, on impose un Dirichlet homogène sur la vitesse (COEFA(IFAC,ICLVAF)=0, COEFB(IFAC,ICLVAF)=0) au lieu de la condition de Neumann homogène (COEFA(IFAC,ICLVAF)=0, COEFB(IFAC,ICLVAF)=1), lorsque le flux de masse est sortant.

• Autres conditions aux limites
On détermine enfin, pour la pression, les grandeurs turbulentes et les autres scalaires, les conditions aux limites de type Dirichlet et Neumann restantes. Ici également, on fournit donc une valeur de COEFA et de COEFB pour chaque face concernée (les valeurs COEFA(IFAC,ICLVAF) et COEFB(IFAC,ICLVAF) sont égales à COEFA(IFAC,ICLVAR) et COEFB(IFAC,ICLVAR)).

0.4 Points à traiter

• Représentation des conditions par une valeur de face
Bien que la méthode utilisée permette une simplicité et une homogénéité de traitement de toutes les conditions aux limites, elle est relativement restrictive au sens où une seule valeur ne suffit pas toujours pour représenter les conditions à appliquer lors du calcul de termes différents.

Ainsi, en k−ε a-t-il été nécessaire, lors du calcul des conditions aux limites de paroi, de disposer de deux couples (A_b, B_b) afin de prendre en compte les conditions à appliquer pour le calcul de la contrainte tangentielle et celles à utiliser lors du calcul du terme de production (et un troisième jeu de coefficients serait nécessaire pour permettre le traitement des gradients intervenant dans les termes de gradient transposé, dans vissec ).

Peut-être pourrait-il être utile de mettre en place une méthode permettant d’utiliser (au moins en certains points stratégiques du code) directement des forces, des contraintes ou des flux, sans passer nécessairement par le calcul d’une valeur de face.

• Difficultés liées à la donnée d’un flux (entrée des données peu intuitive)
L’utilisateur est invité à fournir, dans le cas d’une condition de Neumann, une valeur de flux en W m⁻² pour la température ou l’enthalpie. Par souci de cohérence, lorsqu’il souhaite fournir une condition de Neumann sur une autre variable, il est amené à fournir une grandeur similaire (en fait homogène à α grad f.n, si f est la variable), ce qui est peu naturel pour la vitesse, la pression ou un scalaire, dans le cas où l’on souhaite imposer la valeur du gradient normal à la face et non pas la valeur de (µ+µ_t)grad u_j.n, de Δ t grad (p).n ou de α/C grad (f).n.

• Condition de sortie en pression
La condition de pression en sortie se traduit par p_f=p_I′+R1 et le profil obtenu correspond au profil amont pris aux points I′ et recalé pour obtenir p₀ en un point A arbitraire. Ce type de condition est appliqué sans précautions, mais n’est pas toujours justifié (une condition de Dirichlet basée sur la valeur calculée directement aux faces de bord pourrait être plus adaptée). Les hypothèses sont en particulier mises en défaut dans les cas suivants :

- la sortie est proche d’une zone où l’écoulement n’est pas établi en espace (ou varie en temps) ;
- la sortie n’est pas une surface perpendiculaire à l’écoulement ;
- le gradient de pression dans la direction normale à la sortie n’est pas le même pour toutes les faces de sortie (dans le cas de sortie multiples, par exemple, le gradient n’est probablement pas le même au travers de toutes les sorties) ;
- les points I′ ne sont pas sur une surface parallèle à la sortie (cas des maillage irréguliers par exemple).

On pourrait également tester la méthode d’Orlansky (équation de convection sur la pression).

Par ailleurs, en l’absence de condition de sortie, il pourrait peut-être se révéler utile de fixer une valeur de référence sur une cellule donnée ou de ramener la moyenne de la pression à une valeur de référence (avec le décalage du spectre, on assure l’inversibilité de la matrice à chaque pas de temps, mais il faudrait vérifier si la pression n’est pas susceptible de dériver au cours du calcul).

• Termes non pris en compte
Les conditions aux limites actuelles semblent causer des difficultés lors du traitement du terme de gradient transposé de la vitesse dans l’équation de Navier-Stokes (terme traité de manière explicite en temps). Il est possible de “débrancher” ce terme en positionnant le mot clé IVISSE à 0. Sa valeur par défaut est 1 (les termes en gradient transposé sont pris en compte).
On remarquera que l’indicateur IVISSE active non seulement les termes de gradient transposé en ( ^tgrad v), mais également le terme en − 2/3 grad ( µ_tot div v) avec :

1: sauf en module compressible, cf. cfxtcl
2: Les autres grandeurs (propriétés physiques par exemple) font l’objet d’un traitement différent qui ne sera pas détaillé ici (par exemple, pour la masse volumique, l’utilisateur définit directement les valeurs aux bord, information qui est conservée telle quelle ; on pourra se reportera à usphyv ou phyvar ).
3: L’affectation d’un type se fait en renseignant le tableau ITYPFB.
4: Lorsqu’il n’y a pas de sortie, le spectre des valeurs propres de la matrice est décalé d’une valeur constante afin de rendre le système inversible : voir matrix .
5: Un message indique à l’utilisateur combien de faces de sortie voient un flux de masse entrer dans le domaine.
6: On doit alors fournir les composantes de la vitesse de la paroi.
7: On pourra se reporter à gradrc pour la condition conduisant à l’extrapolation de la pression au bord, condition pilotable par l’utilisation de la variable EXTRAP.
8: Les conditions aux limites spécifiques sont codées en renseignant directement les tableaux ICODCL et RCODCL pour chaque face de bord et chaque variable : des exemples sont fournis dans usclim .
9: Le nombre de Prandtl turbulent est sans dimension et, dans certains cas usuels, pris égal à 0,7.
10: Plus exactement, on a C=C_p pour toutes les VarScalaires f que l’on souhaite traiter comme la température pour les conditions aux limites. Ces VarScalaires sont repérables par l’utilisateur au moyen de l’indicateur ISCSTH=1. Par défaut cet indicateur est positionné à la valeur 0 pour toutes les VarScalaires (qui sont alors traitées comme des scalaires passifs avec C=1) hormis pour la variable thermique éventuelle (ISCALT^ième VarScalaire), pour laquelle on a ISCSTH=1 : on suppose par défaut que la variable thermique est la température et non l’enthalpie. Si l’on souhaite résoudre en enthalpie, il faut positionner ISCSTH à la valeur 2 pour la variable thermique. Pour le compressible, la variable thermique est l’énergie, identifiée par ISCSTH=3. On se reportera à cfxtcl pour le traitement des conditions aux limites.
11: En effet, le flux total sortant du domaine est donné par la somme du flux convectif (si la variable est effectivement convectée) et du flux diffusif. Néanmoins, pour les parois étanches et les symétries, le flux de masse est nul et la condition se réduit à une contrainte sur le flux diffusif. De plus, pour les sorties (flux de masse sortant), la condition à la limite ne porte que sur le flux diffusif (souvent une condition de Neumann homogène), le flux convectif dépendant des conditions amont (il n’a donc pas besoin de condition à la limite). Enfin, aux entrées, c’est le plus souvent une condition de Dirichlet simple qui est appliquée et le flux diffusif s’en déduit.
12: première face de sortie rencontrée en parcourant les faces de bord dans l’ordre naturel induit par la numérotation interne au code
13: L’utilisateur peut définir d’autres types (i.e. affecter à ITYPFB d’autres valeurs entières), mais elles ne recouvrent pas de conditions aux limites par défaut. Les faces de bord ainsi repérées sont cependant traitées comme un ensemble particulier lors de l’impression d’informations de type flux de masse par exemple.
14: Le test contient une référence à ICODCL = 2 qui est un héritage d’une version dans laquelle les conditions de Dirichlet avec coefficient d’échange étaient distinguées des conditions de Dirichlet sans coefficient d’échange.
15: On ne dispose pas encore de conditions aux limites permettant le calcul du gradient lors du passage dans condli au premier pas de temps.
16: ce qui est le cas par défaut

f			α
symbole	nom	unité	symbole	nom	unité
u_i	vitesse	m s⁻¹	µ ou µ+µ_t	viscosité dynamique	kg m⁻¹ s⁻¹
p	pression	kg m⁻¹ s⁻²	Δ t	pas de temps	s
T	température	K	λ	conductivité thermique	kg m s⁻³ K⁻¹
					=W m⁻¹ K⁻¹
H	enthalpie	m² s⁻²	λ/C_p	conductivité thermique/C_p	kg m⁻¹ s⁻¹
		=J kg⁻¹
f	VarScalaire	unité(f)	α	conductivité ou diffusivité	kg m⁻¹ s⁻¹

f			h
symbole	nom	unité	homogène à	unité
u_i	vitesse	m s⁻¹	(µ+µ_t)/d	kg m⁻² s⁻¹
p	pression	kg m⁻¹ s⁻²	(Δ t)/d	s m⁻¹
T	température	K	(λ+C_pµ_t/σ_t)/d	kg s⁻³ K⁻¹
				=W m⁻² K⁻¹
H	enthalpie	m² s⁻²=J kg⁻¹	(λ/C_p+µ_t/σ_t)/d	kg m⁻² s⁻¹
f	scalaire	unité(f)	α/d	kg m⁻² s⁻¹

f			φ_imp,ext
symbole	nom	unité	homogène à	unité
u_i	vitesse	m s⁻¹	(µ+µ_t) grad (u_i).n	kg m⁻¹ s⁻²
p	pression	kg m⁻¹ s⁻²	(Δ t) grad (p).n	kg m⁻² s⁻¹
T	température	K	(λ+C_pµ_t/σ_t) grad (T).n	kg s⁻³=W m⁻²
H	enthalpie	m² s⁻²=J kg⁻¹	(λ/C_p+µ_t/σ_t) grad (H).n	kg s⁻³=W m⁻²
f	scalaire	unité(f)	α grad (f).n	kg m⁻² s⁻¹ unité(f)

`ITYPFB`	type	variables traitées	type de condition	données complémentaires
`IENTRE`	entrée	variables transportées	Dirichlet	RCODCL(.,.,1)	valeur d’entrée	nulle par défaut
		pression	Neumann homogène
`ISYMET`	symétrie	vitesse	symétrie vecteur
		tensions de Reynolds	symétrie tenseur
		autres variables	Neumann homogène
`IPAROI`	paroi	pression	Neumann homogène^a
		vitesses	lois de paroi (Dirichlet)^b	RCODCL(.,.,1)	vitesse de défilement	nulle par défaut
		grandeurs turbulentes	lois de paroi^c
		autres variables
		transportées	Neumann homogène^d
`ISOR09`	sortie	pression	∂²P/∂ n∂τ_i=0
		vitesses	Neumann homogène
			(flux de masse sortant)
			Dirichlet homogène
			(flux de masse entrant)
		autres variables
		transportées	Neumann homogène
`ISOR10`	sortie	pression	∂²P/∂ n∂τ_i=0
		vitesses	Neumann homogène
			(flux de masse sortant)
			Dirichlet homogène
			(flux de masse entrant)
		autres variables	Neumann homogène
		transportées	(flux de masse sortant)
			Dirichlet^e	RCODCL(.,.,1)	valeur d’entrée	nulle par défaut
			(flux de masse entrant)

`ICODCL`	variables	type de condition	données complémentaires nécessaires et suffisantes
1	toutes	Dirichlet	`RCODCL(.,.,1)`	f_imp,ext	valeur	nulle par défaut
			`RCODCL(.,.,2)`	h_imp,ext	coef. d’échange	=`RINFIN`=10³⁰ par défaut
3	toutes	Neumann	`RCODCL(.,.,3)`	φ_imp,ext	valeur	nulle par défaut
4	vitesse	Symétrie vecteur
	R_ij	Symétrie tenseur
5	vitesse	Loi de paroi	`RCODCL(.,.,1)`		vitesse de défilement	nulle par défaut
	k, R_ij, ε, ϕ, f, ω	Loi de paroi
	VarScalaire	Loi de paroi	`RCODCL(.,.,1)`	f_imp,ext	valeur en paroi	nulle par défaut
	(sauf variance)		`RCODCL(.,.,2)`	h_imp,ext	coef. d’échange	=`RINFIN`=10³⁰ par défaut
9	vitesse	Neumann homogène
		(flux de masse sortant)
		Dirichlet homogène
		(flux de masse entrant)
10	k, R_ij, ε, ϕ, f, ω	Dirichlet	`RCODCL(.,.,1)`	f_imp,ext	valeur	nulle par défaut
	VarScalaire	(flux de masse entrant)	`RCODCL(.,.,2)`	h_imp,ext	coef. d’échange	=`RINFIN`=10³⁰ par défaut
		Neumann	`RCODCL(.,.,3)`	φ_imp,ext	valeur	nulle par défaut
		(flux de masse sortant)

Variable		Valeurs de `ICODCL` admissibles
Vitesse	U	1	3	4	5	9
Pression	p	1	3
Variable scalaire de turbulence	k, ε, ϕ, f, ω	1	3		5		10
Tenseur de Reynolds	R_ij	1	3	4	5		10
VarScalaire (hormis variances)		1	3		5		10
Variance des fluctuations d’une VarScalaire		1	3				10