Programme mécanique des fluides

**zabbaou** · 21/03/2009, 15h27

Cherche etudiant(e) qui pourrait m'aider à écrire un programme mécanique des fluides en fortran. il s'agit d'une cavité dont en U dont la surface est entrainée (à partir de navier stokes, me faut calculer la vitesse, la pression...)

**FR119492** · 22/03/2009, 18h22

Salut!
S'agit-il d'un exercice, d'une thèse de doctorat, d'un projet industriel ou quoi d'autre?
Jean-Marc Blanc

**zabbaou** · 31/03/2009, 20h47

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
program cavite
implicit none
integer, parameter::nt=100,nx=20,nz=nx,km=1000 !km=iterations maximum!
real,parameter :: eps=1.e-8, Re=100, omg=1.7, h=1./nx, dt=0.25*Re*h**2
real,dimension(0:nx,0:nz+1)::un,u
real,dimension(0:nx+1,0:nz)::wn,w
real,dimension(0:nx+1,0:nz+1)::p,AP,SM,PK ! ce sont des vecteurs
real,dimension(nx,nz) :: psi
integer :: i, j, K, n
real :: t
real :: ududx, wm, wdudz, lapun, um, wdwdz, udwdx, lapwn, dudx, dwdz, sigma! on doit définir toutes les variables
!    open(10,file='psi.out')
! initialisation des champs!
un=0. ; w=0. ; wn=0. ; p=0.
do i=1,nx ; un(i,nz+1)=2.-un(i,nz) ; enddo ; U=UN
!   write(*,*) ' champs U*' ; do j=nz+1,1,-1 ; write(*,fmt='(20F7.3)')(U(i,j),i=0,nx) ; enddo
! valeurs de AP(i,j) 
        AP=4.
        AP(2:nx-1,1)=3; AP(2:nx-1,nz)=3;AP(1,2:nz-1)=3; AP(nx,2:nz-1)=3
        AP(1,1)=2 ; AP(1,nz)=2 ; AP(nx,1)=2 ; AP(nx,nz)=2
DO n=1,nt
t=n*dt
        !resolution [1] calcul de u*!
    DO i=1,nx-1
        Do j=1,nz
            ududx=un(i,j)*(un(i+1,j)-un(i-1,j))/(2*h)
            wm=0.25*(wn(i,j-1)+wn(i+1,j+1)+wn(i,j)+wn(i+1,j))
            wdudz=wm*((un(i,j+1)-un(i,j-1))/(2*h))
            lapun=(un(i+1,j)-2*un(i,j)+un(i-1,j))/(h**2)+(un(i,j+1)-2*un(i,j)+un(i,j-1))/(h**2)
            u(i,j)=un(i,j)-dt*(ududx+wdudz)+(dt/Re)*(lapun)
       enddo
    enddo  
        u(1:nx,0)=-u(1:nx,1)
        u(1:nx,nz+1)=2.-u(1:nx,nz)
        u(0,0:nz)=0
        u(nx,0:nz)=0
  !  write(*,*) ' champs U*' ; do j=nz+1,1,-1 ; write(*,fmt='(20F7.3)')(U(i,j),i=0,nx) ; enddo
        !resolution [2 ] calcul de w*!
    DO i=1,nx
        do j=1,nz-1
            um=0.25*(un(i-1,j)+un(i,j)+un(i-1,j+1)+un(i,j+1))
            wdwdz=wn(i,j)*(wn(i,j+1)-wn(i,j-1))/(2*h)
            udwdx=um*(wn(i+1,j)-wn(i-1,j))/(2*h)
            lapwn=(wn(i+1,j)-2*wn(i,j)+wn(i-1,j))/(h**2)+(wn(i,j+1)-2*wn(i,j)+wn(i,j-1))/(h**2)
            w(i,j)=wn(i,j)-dt*(udwdx+wdwdz)+(dt/Re)*(lapwn)
        enddo
    enddo
        w(0:nx,0)=0
        W(0:nx,nz)=0
        w(0,0:nz)=-w(1,0:nz)
        w(nx+1,0:nz)=-w(nx,0:nz)  
 ! write(*,*) ' champs W*' ; do j=nz,0,-1 ; write(*,fmt='(20F7.3)')(W(i,j),i=1,nx+1) ; enddo        
        !résolution [3]!
    do i=1,nx
        do j=1,nz
         !lappn=(1/h**2)*(pn(i,j-1)+pn(i-1,j)-4*pn(i,j)+pn(i+1,j)+pn(i,j+1))!
            dudx=(u(i,j)-u(i-1,j))/h
            dwdz=(w(i,j)-w(i,j-1))/h
            SM(i,j)=h**2/dt*(dudx+dwdz)
        enddo
    enddo
 sigma=1. ; k=0
DO WHILE ( sigma > eps .and. K < KM)
  PK=p ; k=k+1 ; sigma=0.
  do i=1,nx
   do j=1,nz    
   P(i,j)=(1.-omg)*P(i,j)+omg*(P(i,j-1)+P(i-1,j)+P(i,j+1)+P(i+1,j)-SM(i,j))/AP(i,j)      
   sigma = sigma + ((P(i,j)-Pk(i,j))**2)
   enddo    
  enddo
sigma=sqrt(sigma)
 ! print *,' k sigma ',n,k,sigma
enddo
  print *,t,k,sigma
! write(*,*) ' champs P' ; do j=nz,1,-1 ; write(*,fmt='(20F7.3)')(p(i,j),i=1,nx) ; enddo 
!conditions limites sur p 
 !       p(1:nx,0)=p(1:nx,1)
 !       p(1:nx,nz+1)=p(1:nx,nz)
 !       p(0,1:nz)=p(1,1:nz)
 !       p(nx+1,1:nz)=p(nx,1:nz)
 
            !résolution 4 pour obtenir Vn+1!
    Do i=1,nx-1
        DO j=1,nz
            u(i,j)=u(i,j)-dt*(p(i+1,j)-p(i,j))/h ! la variable à gauche écrase celle de droite!
        enddo 
    enddo
            !résolution 5!
    do i=1,nx
        do j=1,nz-1
            w(i,j)=w(i,j)-dt*(p(i,j+1)-p(i,j))/h
        enddo
    enddo
   ! reinit
   UN=U ; WN=W
enddo
!
 call gnu(psi,x,z,h,nx,nz)
!
end program
 
!    subroutine gnup(psi,x,z,h,nx,nz)
!        implicit none
!            integer :: nx,nz,i,j
!            real :: h
!            real, dimension(nx,nz) :: psi
!            real, dimension(nx) :: x
!            real, dimension(nz) :: z
!            do j=1,nz
!                do i=1,nx
!                    write(20,fmt='(2F8.3,2E14.5)') x(i),z(j),psi(i,j)
!                enddo 
!                    write(20,*) 
!            enddo 
!    end subroutine gnup
!
!    subroutine calpsi(psi,u,w,h,nx,nz)
!    integer :: nx,nz,i,j
!    real :: h
!    real, dimension(nx,nz) :: psi
!    real,dimension(0:nx,0:nz+1):: u
!    real,dimension(0:nx+1,0:nz)::w
!    !calcul pour psi qui est la même grille que p (décalé-décalé)!
!    psi(1,1)=0
!    do i=2,nx
!    psi(i,1)=psi(i-1,1)-w(i,1)*h!on travail sur la première ligne
!    enddo
!    do j=2,nz
!     psi(1,j)=psi(1,j-1)+u(1,j)*h!on travail sur la première colonne!
!    enddo
!    do i=2,nx
!        do j=2,nz
!            psi(i,j)=psi(i,j-1)+u(i,j)*h          
!        enddo    
!    enddo
!    end subroutine    
 
!

**kango** · 31/03/2009, 22h31

bonsoir,

qu'est ce que tu entends par "ne tourne pas" ?

- ça compile pas
- ça compile mais tu as une erreur à l'exécution
- pas d'erreur à l'exécution mais le résultat est incorrect

**zabbaou** · 01/04/2009, 10h29

ça ne complile pas

**kango** · 01/04/2009, 10h35

et quel est le message d'erreur que tu obtiens à la compilation ?

**FR119492** · 02/04/2009, 17h21

Salut!

D'abord, toutes mes excuses pour ne pas t'avoir répondu plus tôt: je réfléchissais à ton problème.

Comme tu commences à t'en rendre compte, l'équation de Navier-Stokes est une chose assez monstrueuse qui constituerait à elle-même une thèse de doctorat en mathématiques appliquées.

Lorsque je travaillais dans l'industrie de l'aluminium, pour tenter d'optimiser la situation magnétohydrodynamique des fours d'électrolyse, j'ai été confronté à un problème analogue: la direction avait entendu de l'équation de Navier-Stokes et m'a sommé de la résoudre. Après mûre réflexion, je me suis rendu compte que

c'était impossible dans les délais impartis;
ça ne servait à rien.

J'ai eu alors recours à l'astuce suivante: j'ai remplacé l'équation de Navier-Stokes par celle de Poisson, dont l'intégration ne présentait aucune difficulté. De cette manière, j'ai calculé les tourbillons dans les fours existants, en ajustant la viscosité de manière à ce que les résultats coïncident plus ou moins avec les valeurs mesurées. Ensuite, j'ai fait de même avec les données d'autres fours. De cette manière, j'ai pu optimiser la disposition des conducteurs autour des fours. Gain de l'opération: 0,5 kWh par kg d'Al par des modifications mineures sur des fours existants et plus du double sur de nouveaux prototypes.

Comme quoi, une supercherie bien réfléchie, ça peut être payant!

Si tu penses qu'une astuce analogue peut convenir dans ton cas, ça vaut la peine d'essayer. Sinon, j'ai quelques idées sur ce qu'il faut faire ou ne pas faire avec Navier-Stokes, et on pourra y revenir par la suite.

Jean-Marc Blanc

**Leflix** · 08/04/2009, 15h00

Envoyé par FR119492

Salut!
Après mûre réflexion, je me suis rendu compte que

c'était impossible dans les délais impartis;
ça ne servait à rien.

Jean-Marc Blanc

Desolé de te contredire Jean-Marc mais contrairement a ce que tu dis resoudre les equations de Navier-Stokes ne sert pas a rien.
Si les avions vols c'est grace aux equations de Navier-Stokes....(entre autre). Toutes les souffleries qu'on appelle "numeriques" sont basees sur la resolution numerique de ces equations (ce que tu cherches a faire zabbaou).
Jean-Marc resoudre l'equation de Poisson ne t'apportera absolument rien, car ca n'a rien a voir ! J'irais pas faire cuire un poulet dans un des fours que tu as optimise !!

Navier-Stokes te permet de calculer les composantes de la vitesse du fluide puis le champ de pression. L'equation de Poisson te permet au mieux de calculer la fonction de courant (si tu disposes encore du rotationelle de la vitesse.. et pour calculer le rotationnelle il te faut la vitesse...) Bon bref je ne vais pas te faire un cours de meca flu ici....

par contre la ou je suis d'accord avec toi , c'est que c'est extremement delicat a faire...
et des gens ont developpe de tels programmes dans le cadre d'une these de 3 ans !!
Je ne veux pas t'effrayer zabbaou mais ce que l'ont t'a demande comme travail c'est au minimum 3 mois de boulot pour un master1 en mecanique...et encore tu as interet a etre bon en mecanique des fluides et methodes numeriques pour savoir de quoi l'on parle...

**FR119492** · 08/04/2009, 15h29

Salut!

J'irais pas faire cuire un poulet dans un des fours que tu as optimise

Navier-Stokes ou pas, tu risquerais surtout d'être intoxiqué par les composants fluorés les plus divers.

Si les avions vols c'est grace aux equations de Navier-Stokes

Parce que tous les avions d'un même modèle ont des caractéristiques plus ou moins analogues. En revanche, des fours de construction identique avaient des comportements très différents. C'est pourquoi une approximation plus que grossière était suffisante.
Jean-Marc Blanc

**El Zeru** · 13/05/2009, 22h51

Il ne faut pas non plus exagérer sur le temps qu'il faut pour écrire un code Navier-Stokes. Tout dépend de ce que l'on veut.
On peut faire un code 2D incompressible monophasique sur maillage cartésien très rapidement si l'on a pas à réinventer la poudre (les discrétisations)

Si vous voulez du 3D multiphasique avec tension de surface, maillage non structuré, objets complexes mobiles et modèles de turbulence à gogo,oui là ça prendra du temps...

**genteur slayer** · 15/05/2009, 16h16

je tiend à dire que perso, je développe un code navier-stokes incompressible turbulent diphasique avec surfaces mobile par érosion suportant des conditions limites aussi diverses que variées (adherence, glissement, débit imposé, pression imposée, symétrie, périodicité, symétrie central... etc...) le tout en 2D, 3Daxi et 3D, et bin les gars, j'en chie!!!! surtout qu'il faut des schéma de résolution d'ordre 5 minimum pour avoir une érosion acceptable... en clair je fais fumer toutes les machines de la boite qui n'ont rien à faire (ah oui, petit p.s. tout cela est parallèle, évidement) le seul truc, j'ai un maillage cartésien fixe, sauf qu'il est décalé... ce serai trop beau!!!

pour son problème au petit, je suis pas sûr que navier-stokes soit obligatoire, une bète equation de stokes devrai suffir, voir peut-être plus simple encore...

**El Zeru** · 15/05/2009, 21h59

Tu fais de la décomposition de domaine ? Si oui, tu adaptes comment tes schémas aux limites des sous domaines ? Tu utilises quelle méthode pour prendre en compte tes obstacles ? (IBM ?)

**genteur slayer** · 19/05/2009, 11h16

pour les obstacles, c'est méthode de pénalisation (c.f. Angot 1999) en gros tu rajoute un terme source Mu/K*u où Mu c'est ta visco, K la perméabilité du poreux(l'obstacle est vu comme un poreu) et u ton champ de vitessedu coup tu met ce truc à 0 lorque tu est dans la partie fluide, et non nul dans la partie obstacle tu repère tes milieu avec un méthyode Level Set.

pour les sous domaines bin c simple tu va chercher l'info qu'il te manque dans le processus d'à coté il suffit de faire une "condition limite" "intérieur" (qui considère qu'elle est au milieu du domaine)

**El Zeru** · 20/05/2009, 16h46

Haaaa mais la pénalisation Darcy, c'est dépassé Monsieur ! Jette donc un oeil la dessus :
http://www.trefle.u-bordeaux1.fr/per...thou_FVCA5.pdf

C'est un extension à l'ordre 2. Bien plus précise !

**Leflix** · 22/05/2009, 23h18

Envoyé par El Zeru

Haaaa mais la pénalisation Darcy, c'est dépassé Monsieur ! Jette donc un oeil la dessus :
http://www.trefle.u-bordeaux1.fr/per...thou_FVCA5.pdf

C'est un extension à l'ordre 2. Bien plus précise !

Mouaih....

ordre 2 c'est vite dit !!!
pour votre info les gars les deux methodes pré-cités en matière de gestion d'obtsacles proviennent de la meme culture...(ces gens en question ont travaille ensemble.... lol labo Master qui est ensuite devenu Trefle..)
De toute facon (sauf dans le cas d'obstacles mobiles) il n'y a rien de plus precis et de mieux que le BFC (Body Fitted Coordinate)

Mais bon la on s'eloigne de la preoccupation initiale du jeune...

**El Zeru** · 23/05/2009, 14h40

Ha oui, tu parles surement d'Angot et Caltagirone, mes directeurs de thèses

Donc oui, je connais bien tout ça. Nous avons formulé la méthode à l'ordre deux comme une extension de la méthode à l'ordre un, en particulier la VPM de la thermique où tu ajoutes un terme B(u-uD), car la formulation Darcy, bien que plus physique pour Navier, ne permet pas d'imposer naturellement une solution autre que 0. Je précise qu'on est en lagrangien augmenté pour le couplage vitesse pression, en projection c'est une autre histoire.
Pour l'ordre deux, pour le moment nous l'avons tout le temps obtenu, (cf le lien d'avant FVCA5) sans faire toutefois des millions de validations. Nous avons aussi la preuve théorique dans les cartons. Effectivement on ne peut pas rivaliser avec le body fitted pour la précision. Après on peut se poser la question de la précision et de la complexité de la discrétisation des opérateurs dans ce cadre là. Je travaille aussi sur le Code_Saturne d'EDF qui est en VF non-structuré et c'est un code plus complexe que Thétis, le code du labo en structuré.
Dans le cas du Couette cylindrique, on est 100 fois plus précis sur maillage polaire que sur Cartésien. Mais dans ce cas là, la discrétisation sur maillage polaire est à la fois simple et précise. A voir sur maillage plus tordu.

Si vous allez au SMAI, on pourra en discuter

a+

**Leflix** · 23/05/2009, 20h04

Envoyé par El Zeru

Ha oui, tu parles surement d'Angot et Caltagirone, mes directeurs de thèses

yes

Donc oui, je connais bien tout ça. Nous avons formulé la méthode à l'ordre deux comme une extension de la méthode à l'ordre un, en particulier la VPM de la thermique où tu ajoutes un terme B(u-uD), car la formulation Darcy, bien que plus physique pour Navier, ne permet pas d'imposer naturellement une solution autre que 0.

oui c'est normal !!

Je précise qu'on est en lagrangien augmenté pour le couplage vitesse pression, en projection c'est une autre histoire.

Les autres methodes d'IBM marchent aussi en methode de Projection , mais peut etre moins bien que le lagrangien augmente....justement a cause du terme supplementaire de Darcy dans l'equation de NS.

Pour l'ordre deux, pour le moment nous l'avons tout le temps obtenu, (cf le lien d'avant FVCA5) sans faire toutefois des millions de validations. Nous avons aussi la preuve théorique dans les cartons. Effectivement on ne peut pas rivaliser avec le body fitted pour la précision. Après on peut se poser la question de la précision et de la complexité de la discrétisation des opérateurs dans ce cadre là. Je travaille aussi sur le Code_Saturne d'EDF qui est en VF non-structuré et c'est un code plus complexe que Thétis, le code du labo en structuré.

oui c sur que c'est plus que c'est plus complexe de discretiser les equations sur un maillage non-orthogonal. Mais ca reste toutefois tres abordable en utilisant le principe de la "deferred correction"

Dans le cas du Couette cylindrique, on est 100 fois plus précis sur maillage polaire que sur Cartésien. Mais dans ce cas là, la discrétisation sur maillage polaire est à la fois simple et précise. A voir sur maillage plus tordu.

Si vous allez au SMAI, on pourra en discuter

a+

Oui c normal que sur un maillage polaire pour un ecoulement entre deux cylindres ca marche mieux que sur un maillage cartesien avec une IBM. Ceci etant avec la discretisation utilisant la "deferred correction" c'est tranquille.
Voir le bouquin de Ferziger et Peric "computational fluid dynamic". C'est la bible ! Peric ancien prof de l'universite de Hambourg est aujourd'hui le directeur du developpement de Star-CD !!

Le code saturne utilise cette discretisation "a la Peric" (deferred correction)

Bon mais peut etre faudrait il ouvrir une nouvelle discussion intitulee " Comment vaincre Navier-Stokes ?" car je crains que l'on s'eloigne du sujet original....

Programme mécanique des fluides

Fortran

Discussions similaires

Partager

Partager