Discussion :

[ruyan]

Bonjour à tous!

Je rencontre des difficultés en programmation fortran 77 et j'espère que vous pourrez m'aiguiller...

Je voudrais simplement modéliser un problème de conduction instationnaire dont voici les principales données:

Soit un puits contenant de l'eau.
Le puits fait 500m de profondeur et 18cm de large.
Le puits est entouré sur toute sa profondeur par le terrain environnant depuis la surface.
Je considère donc le diamètre du terrain autour de 3m de diamètre comme condition limite.
Je me place en coordonnées cylindrique et je veux connaitre la température
J'ai un terrain de 500m de profondeur avec un gradient géothermique de 3°C tous les 100m.


Mon but est de calculer la température de l'eau dans le puits.


Pour ce faire j'ai dans un premier temps considéré le problème en 1D dans le sens de la profondeur (Z) en faisant un bilan sur une maille quelconque, la maille en tête de puits et la maille de fond et en considérant que la chaleur apportée par le terrain est constante.

Maintenant je veux pouvoir modéliser le flux de chaleur apporté par le terrain en 2D (profondeur+radial)

Comment faire,notamment pour programmer les conditions limites?


Merci à tous!

[afrancisco]

salut,

Tu as un problème axisymétrique stationnaire, donc tu fais un maillage rectangulaire, de largeur 3m et de hauteur 503m.
- Dedans tu as continuité de la température et un gradient de température imposé
- sur les noeuds de l'axe d'axisymétrie le flux est constant = 0 (dérivée thermique nulle)
- ailleurs sur les bords, tu as des températures imposées

Ces conditions sont imposées :
- sur le système linéaire constitué, soit par élimination, soit par pénalisation
- ou bien, explicitement sur une résolution itérative type gauss-seidel

voilà, à vue d'oeuil, c'est de ce style

PS: préfère fortran 90/95 !

[ruyan]

Bonjour Afranscisco!

Merci pour ta réponse rapide mais au fait je dois rajouter une difficulté au problème:

En fait, j'aimerais faire évaporer l'eau de ce puits en créant une dépression en tête du puits.
En créant une dépression,l'eau du puits va entrer en ébullition et s'évaporer...
Ce phénomène est dit endothermique, c'est-à dire qu'il nécessite l'apport de chaleur pour se réaliser.
De ce fait,la seule source de chaleur disponible est la chaleur amenée par le terrain et l'eau des mailles sous jacente.
En cédant de la chaleur, le terrain va se refroidir ainsi que l'eau en surface.


Dans mon modèle,j'impose un débit d'évaporation (exemple 0.05 kg/s) que j'ai choisis de telle sorte à éviter que le terrain et l'eau ne gèle.

je pense que ça change un peu les conditions aux limites et un autre de mes problème érside également dans le maniement des tableaux avec 'DIMENSION'...

Je ne sais pas si tu as tout saisi mais si toi ou quelqu'un d'autre pouvais m'aider ce serait cool!!

[afrancisco]

salut,

pour l'évaporation, c'est simplement un flux thermique sortant par le haut imposé.
L'aspect instationnaire peut venir de la perte d'eau dans le puits, en fonction du temps, compensée par de l'air, mais ça devient un peu vicieux ....

Pour les tableaux, je n'ai pas fait de f77 depuis très très longtemps ; en f90/95 les tableaux s'allouent/desallouent dynamiquement avec l'attribut allocatable.
Mais dans ton problème, il n'y a pas lieu de changer les bornes ; et en plus tu peux ne stocker que les termes non nuls (5 diagonales) ce qui allège incroyablement ton système.
Par contre, il faut être à l'aise avec les indices, lors de la résolution...

bon courage

[gregory2526]

salut.
je n'ai pas beaucoup de temps pour te donner une réponse précise mais si tu tire de la chaleur du terrain, il faudrait que tu pense à éloigner ta condition limite verticale, c'est à dire mettre plus de 3m de terrain autour du puit.
ça te permettrai de conserver une CL de dirichlet. parcequ'à mon avis, 3m de terrain ne sont pas suffisants.
Sinon, peut-être qu'une CL de type mixte pourrait arranger le problème. Tu considère qu'en frontière, le flux est proportionnel à la température. En calculant le flux de température dû à l'absorption de chaleur par le puit, tu peux en déduire la nouvelle température de ta condition limite. Il te suffit de prendre un coefficient du genre conductivité thermique (ou qq chose comme ça) comme coefficient de proportionnalité...

Je ne sais pas si ça peut t'aider (la reflexion à été rapide, j'ai peut-être dis des bêtises ou mal compris ton problème)

[ruyan]

Merci pour vos réponses.

J'aurais cependant une autre question:

Le fait que l'eau du puits s'évapore fait que l'interface eau/air ambiant diminue au fur et à mesure du processus.
Comme la température du puits est plus élevée en profondeur,cela veut dire que plus l'interface descend,plus l'évaporation va être intense.

Je ne sais malheureusement pas comment faire pour prendre en compte ce phénomène (recalculer le maillage à chaque itération?)

[gregory2526]

Peut-être que pour faire simple, tu pourais rajouter une variable dans ton problème qui serait la quantité d'eau dans une cellule (1 pour une cellule remplie d'eau et 0 pour de l'air). L'interface entre l'eau et l'air pourrait être représenté par les cellules dont la fraction d'eau serait comprise entre 0 et 1. Et tu pourrait dire que après calcul de l'évaporation, si cette variable est supérieure à 0.5, l'interface est située à l'interface entre la cellule et la cellule au dessus, et si elle est plus petite que 0.5, l'interface serait située sur la face entre cette cellule et la cellule en dessous.
Je ne sait pas si cette approximation peut te convenir. ça a l'interet de t'éviter de recalculer un maillage.
Après faut trouver un moyen pour que cela soit pris en compte seulement dans le puit, et pas dans le sol.
Voila.
C'est ce à quoi je pense. Mais si le calcul du maillage n'est pas trop long ou trop compliqué, tu peut tenter... ce sera probablement plus juste pour la position de l'interface. Par contre, il te faudra interpoller tes valeurs dans le nouveau maillage en fonction de l'ancien, et ça risque d'introduire des erreurs ..

Bon courage

[genteur slayer]

ta seule variable est-elle la température???
je veux dire par là que tu n'as pas d'écoulement de fluide (air ou eau)?
tu veux connaitre la tempértature où? seulement dans l'eau oubien partout!!! sinon pour savoir où est l'eau c simple elle est dans le domaine définit par z<h



parce que dans ce cas là, tes seuls inconnues sont T et h où T est la température en tous pointsd et h la hauteur d'eau dans ton puit.

on a dh/dt=div(q) où q est ton flux 's'évaporant' que tu calcul de façon annexe (dsl je suis pas callé en transition de phase) q est lié à T qui lui est calculer avec une équation du style div(grad(T))=0 (equation de poisson) et les conditions aux bords qui vont bien (surment dépendant de h)ensuite, tu te choisit un schémas numérique par exemple element/volumue finis ou même machins finis (hybride entre les deux précedent) et tu fais résoudre!!!

CONSEIL (à tous):avant de programmer, il faut poser les équations et les conditions limites, si déjà là, tout est bien posé, la programmation se fait les doights dans le nez!!!

P.S. il faudrais un module LaTeX dans ce forum genre de balise [LaTeX][\LaTeX] affin de pouvoir posé des équations propres....

[Butterfly83]

Essaies unmaillage radial isosurface c est le mieu pour ton problème.
Rentre tes condition initiales dans un vecteur entrées et fait tourner...

[gregory2526]

Un maillage radial isosurface ?
c quoi ?
ça ressemble a un nom de type de maillage 3D de logiciel commercial ou qq chose comme nça non?

[Butterfly83]

A priori les 2d, si j ai bien compris sont l axe de profondeur et la vue de face.
Un maillage radial isosurface c est tout simple, tu regardes ton disque de face et tu fait un maillage en cercle concentrique, le premier reflexe serrait de faire des cercles espacés de 50cm par exemple, le maillage isosurface dit qu en fait toutes tes maille aurront la meme surface donc de plus en plus rapproché au fur et à mesure que l on s éloigne du centre du cylindre.

Voilà, et a priori les calcul des méthodes numériques sur ce genre de maillage est plus agréable...

[genteur slayer]

le maillage isosurface radial est inconnu au bataillon!!! c quoi exactement???????

de plus je suis pas sûr qu'un cas 3D soit profondément nécéssaire, tu as une axisymétrie il me semble non? donc du 2D carré grossier doit suffire à moins que ton puis ne soit pas vertical ou qu'il chauffe plus d'un coté que de l'autre.... en fait cela dépend surtout de tes condition aux bords (normal tu essaye de résoudre un problème aux limites :p )

cela dit, étant matheux, si tu peux me donner du grain à moudre, des fois, on se rend compte que des tas de choses se simplifie et que le problème de départ se réduit à une équation de domaine et une sur l'interface eau/air.

il me semble vu ton problème que tu résout l'equation de poisson dans l'eau avec des condition de bord de type neuman explicit avec une onde de choc sur l'interface, le choc étant quatifié par le flux de vapeur à travers l'interface et faisant bouger cette même interface....

[gregory2526]

J'ai fait une recherche sur google pour maillage isosurface radial : rien (ni en anglais d'ailleurs).
Et puis en 2D, le puit devient rectangle, et non disque...

Je suis d'accord avec genteur slayer, du carré 2D s'est bien suffisant. Surtout si le maillage est fait à la main (c'est bien plus simple à programmer)

Par contre j'ai une question.. Qu'est ce que tu entend pas Neuman explicite ?
je connait Neuman homogène (flux nul) ou non, mais pas l'explicite.

[genteur slayer]

c'est juste une manière de coder la condition de neumann, c'est à dire que tu connais le flux de chaleur à priori et non que tu le calcul à posteriori...
en gros cela veut dire qu'il est connu

moi je bosse sur l'erosion d'un sol par un écoulement, et cela ressemble pas mal à de la transition de phase.... et j'ai quelque mal à faire bouger ma ligne d'interface sans être obliger de remaillé à chaque pas de temps, mais il existe pas mal de méthode pour le faire normalement.... tu peux aussi chercher par là.

[gregory2526]

je m'excuse d'avance pour cette réponse qui ne fait pas avancer le problème de ruyan, mais ça m'intrigue cette appellation de Neuman explicite.

Ce que je connait c'est la condition limite de Neuman (homogène si le flux est nul, non homogène si le flux est non nul) dans lequel la valeur du flux est donnée au départ.

Est-ce que dans ce cas Neuman explicite serait pour la Cl de type Neuman que j'ai écrite ci-dessus et un "Neuman implicite" serait une sorte de CL de type mixte, dans laquelle la valeur du flux serait calculée en utilisant la valeurs de la grandeur dans le domaine ? par exemple, on aurait Ft=a*t, où Ft serait le flux de température, t la température calculée dans le domaine et a un coefficient de proportionnalité ?

[genteur slayer]

oui c'est à peu près cela... disons que l'explicit est connu et donc bel et bien imposée!!! tandis que l'implicite est une inconnue qu'il faut calculer (avec une méthode x ou y, mais c souvent très délicat car, resoudre un pb aux limite sans connaitre les limites...)
une des methodes possible est la belle Newton-Raphson: on suppose le flux d'une certaine valeur, on calcul et on fait un calcul d'érreur qui permet d'adapter ce flux et recommencer le calcul (en une dizaine d'itération on retrouve notre flux (quand ça marche bien))

[ruyan]

Bonjour à tous!

Merci pour vos réponse qui me donnent des pistes à explorer.

Voilà quelques infos supplémentaires pour genter slayer:
en effet je cherche à calculer le profil de température le long du puits en fonction du temps ainsi que la hauteur de l'eau dans le temps.

Mes notions de programmation ne sont que rudimentaires et je ne connais que les différences finis...

Comme l'a dit genter slayer je me trouve avec un problème axisymétrique et je cherche à calculer la température T(r,z,t) où z:profondeur,r:coordonéée radiale et t:temps.

Est-ce que quelqu'un pourrais me montrer un exemple de bout de code pour résoudre ce problème?(surtout au niveau de l'interface qui me laisse perplexe...)

Merci d'avance pour vos réponse.

[genteur slayer]

et bien comme qui dirait si tu veux vraiment modélisé l'interface.... c un problème qui n'est que partiellement résolu... cela dit dans ton cas il doit y avoir des astuces:
-on n'a pas besoins vraiment d'avoir une hauteur d'eau dépendante de r, on peut donc dire que ton domaine est de la forme [0,R]x[0,H(t)] donc ta hauteur d'eau est donné par H(t); dH/dt est donné par ta loi d'évaporation, et dépend de la température.

-Pour la température, à un tenps donné t; on résout \Laplacien T=0 avec comme conditions aux bords que le flux est nul pour z=0 et r=0, le flux est égal à g sur r=R et inconnu sur z=h(t) donc on écris la formulation intégrale du problème on fait l'intégration par partie, on ajoute les conditions de bord... en diférences finie, je crois que cela se fait.

-une fois que tu as le champ de T° à l'instant t, tu peux définir la quantité d'eau qui s'évapore, et mettre à jour H

-passe à l'instant suppérieur

{c un schémas d'euler explicit tout simple}

[ruyan]

Merci pour tes précieux conseils...

Cependant pourquoi est-ce que le flux est nul en z=0,r=0?

[gregory2526]

Ton problème est axisymétrique. Il y a donc une symétrie en r=0 (axe de symétrie). Dans ce cas, le flux d'une grandeur est nulle. Il y a probablement une démonstration mathématique pour cela mais je ne suis pas capable de te la donner.
tu peux te représenter ça par le fait que si tu a un plan de symétrie, et que le flux est non nul, tu aura une certaine quantité de la grandeur qui traversera ce plan. Pour conserver la symétrie, il faudra que la même quantité de cette grandeur le traverse dans la direction opposée.

Je ne suis pas sur que mon explication est très claire...