Le code et sa construction

Cette seconde partie détaille les méthodes numériques de notre code 3D dans un repère cylindrique.

Les schémas numériques et la discrétisation de l'équation de la chaleur sont exposés, les conditions de stabilité rappelées. On a démontré également sa convergence en maillage et enfin, on a rappelé les inputs et outputs pour les utilisateurs à l'aide d'un diagramme.

 

Equation de la chaleur et Discrétisation

Notre étude étant un problème de conduction instationnaire, le modèle régissant la physique, ici, n'est autre que l'équation de la chaleur :

En faisant l'hypothèse d'une conductivité uniforme $\lambda$ dans le piston:

$ \rho C_p\Large{ \frac{\partial T}{\partial t}} $=$    \nabla .( \lambda \nabla T )  = \lambda \Delta  T$ 

Soit en cylindrique :

$ \Large{ \frac{\partial T}{\partial t}} $=$\Large{ \frac{\lambda}{\rho C_p}(\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial T}{\partial r})+\frac{1}{r²}\frac{\partial^2 T}{\partial \theta²}+\frac{\partial^2 T}{\partial z²})}$=$\Large{ \frac{\lambda}{\rho C_p}(\frac{\partial^2 T}{\partial r²}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{1}{r²}\frac{\partial^2 T}{\partial \theta²}+\frac{\partial^2 T}{\partial z²})}$

En optant pour un schéma explicite, centré, on obtient :

Les dérivées selon $e_r$ : $ \Large{\Delta_r = \frac{T^n_{i+1,j,k}-T^n_{i-1,j,k}}{2rdr}+\frac{T^n_{i+1,j,k}-2T^n_{i,j,k}+T^n_{i-1,j,k}}{dr²}}$

La dérivée selon  $e_\theta$ : $ \Large{\Delta_\theta =  \frac{T^n_{i,j+1,k}-2T^n_{i,j,k}+T^n_{i,j-1,k}}{r² d\theta²}} $

La dérivée selon $e_z$ : $ \Large{\Delta_z =  \frac{T^n_{i,j,k+1}-2T^n_{i,j,k}+T^n_{i,j,k-1}}{ dz²}} $

$ \Large{T^{n+1}_{i,j,k}=T^{n}_{i,j,k} +\frac{\lambda}{\rho C_p}(\Delta_r+\Delta_\theta+\Delta_z) }$

Maillage 3D cylindrique

Pour un maillage 3D cylindrique, se pose comme problème la définition et le calcul de l'axe central, un axe singulier dans notre maillage.

Ce point, s'il est maillé, est connecté à une multitude d'autres points (proportionnel au pas d'espace $d\theta$ selon $e_\theta$).

En volumes finis, cela impliquerait des flux infinis car la surface d'échange est nulle.
De même en différences finies, cela implique de définir un flux au centre, que nous ne pouvons définir. L'objectif du projet étant de faire un calcul pour une physique dissymétrique, nous ne pouvons pas utiliser la condition de symétrie pour l'axe central.

La seule solution proposé par [7] est de ne pas définir l'axe central dans notre maillage mais de créer un cylindre creux de diamètre $dr$.

 

Ce maillage, pour être stable, doit respecter un nombre de Fourier $\leq$ 0.5.

=> Fo = $dt\Large ( \frac{D}{dr²}+\frac{D}{dz²}+\frac{D}{(r_{min}*dth)²})$ $\leq$ 0.5.  avec D: coefficient de diffusion thermique

 

 

 

Consistence en maillage du code 3D cylindrique

Pour vérifier la consistance de notre schéma numérique sur un maillage 3D, nous avons utilisé une fonction gaussienne pour l'initialisation de notre problème. 

Elle est définie de la façon suivante sous matlab :

T0=400; % Kelvin

for k=1:Nz
    for j=1:Nth+1
        for i=1:Nr
            r=(2*i-1)*dr/2;
            T(i,j,k)=(T0-293.15)*exp(-r*r*4*4/D/D)+293.15;
        end
    end
end

En fixant des transferts nuls aux parois, la gaussienne doit se diffuser au cours du temps, et la température tendra vers la température moyenne. Cette température moyenne doit se conserver au cours du temps. 

Nous nous sommes placés dans les configurations suivantes, en faisant varier uniquement le nombre de points selon $e_r$:

  • Nr=[25,50,500]
  • Nth=2
  • Nz=4
  • Fo=0.45
  • tf=100*tcycle

 

La température moyenne  $\Large \frac{\int_{S(z=h)} T dS}{S}$ doit rester constante au cours du temps. En fonction du maillage et de son raffinement, cette contrainte physique est plus ou moins respecter. Un maillage avec 25 points selon $ e_r$ est moins précis qu'avec 500 points, mais en fonction de l'erreur que nous pouvons tolérer, cela peut convenir.

 

 

 

 

 

 

Il n'existe pas de solution analytique, nous considérons donc que 500 points selon $ e_r$ est la solution quasi-exacte de notre problème. Nous pouvons observer que pour un faible nombre de points selon $e_r$, les résultats sont assez proches de la solution exacte.

 

 

 

 

 

 

Quelque soit le nombre de points (>25) selon $e_r$, notre code converge bien vers la solution stationnaire qui est équivalente à la température moyenne de notre problème de diffusion thermique pur.

 

 

 

 

 

 

 

Fonctionnement du code

L'utilisateur doit donc fournir au code un certain nombre de données afin de l'executer et voir l'influence de certains paramètres sur la température de surface du piston). La pression et température des gaz nous ont été fourni par Continental en fonction de l'angle vilebrequin (CA) sur un cycle (4 Temps) pour des régimes moteur et charge différents. En outre, l'utilisateur peut changer manuellement le coefficient de Swirl, le SOI (qui correspond à l'angle vilebrequin à partir duquel on commence l'injection).

Note: ce schéma n'est pas exhaustif , d'autres paramètres seront demandés par le code.