!programme PRINCIPAL
parameter (n=1000, m=40) ! dimensions du canal
REAL*8 f(0:8,0:n,0:m), feq(0:8,0:n,0:m)! fonction de distribution des particules et la fonction d’équilibre
REAL8 cx(0:8), cy(0:8), w(0:8) ! ! les vitesses et les dimensions de reseau
REAL8 u(0:n,0:m), v(0:n,0:m), rho(0:n,0:m) ! velocidade-x, velocidade-y e densidade do fluido
REAL*8 Ma, Lx, hauteur_lattice,tau ! Número de Mach, comprimento e largura do canal
cx(
= (/0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0/) !composante x de la vitesse du réseau de D2Q9
cy( = (/0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0/) !composante y de la vitesse du réseau de D2Q9
w( = (/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./) ! facteurs de pondération: 16/36, 4/36, 1/36
!visc_physique = viscd / rhoo ! viscosité cinematique (m²/s)
!visc_lattice = 0.0000154
! dx = 1 ! distance entre les sites de réseau dans la direction x (m)
! dy = dx ! distance entre les sites de réseau dans la direction y (m)
! dt = 1 ! pas de temps (s)
!cs = dx / dt ! vitesse du son dans le réseau (m / s)
!Lx = n* dx ! longueur de canal (m)
!Ma = uo / cs ! Nombre de Mach [0.1 - 0.4]
!ly = hauteur_lattice* dy ! largeur de canal (m)
uo = 0.1 ! Vitesse d’entrée (m / s) [0.1 - 0.4]
rhoo = 1.0 ! densité de l’air (kg/m³) à 20 °C
! hauteur_lattice de canal (LU)
hauteur_lattice=80
! Nombre de Reynolds ! Re = uo * hauteur_lattice / visc_lattice
Re=400
visc_lattice=(uo*hauteur_lattice)/Re
! Instabilité si 2 <tau <0,506 !!!’
! stabilité tau [0.506 - 2]
omega=1./(3.*visc_lattice+0.5)
tau= 1./omega
!tau =(3*visc_lattice)+0.5 ! Relaxation time
print *, "tau = ", tau ! sortie: vitesse d’entrée
!sumvelo = 0.0 ! réinitialisation (initialisation) somme
!omega = 1.0/tau ! facteur de relaxation: flux
!mstep = 8000 ! nombre de pas de temps
mstep = 10000 ! nombre de pas de temps
print , "**** unités SI *****"
print *, "Lx = ", Lx !sortie: longueur de canal
print *, "hauteur_lattice = ", hauteur_lattice ! sortie: largeur de canal
print *, "rhoo = ", rhoo ! sortie: densité
print *, "mu = ", viscd ! sortie: viscosité dynamique
print *, "visc_lattice = ", visc_lattice ! saída: viscosidade cinemática
print *, "Uo = ", uo ! sortie: vitesse d’entrée
print *, "Cs = ", cs ! sortie: vitesse du son dans le réseau
print *, "Re = ", Re ! saída: número de Reynolds
print *, "Ma = ", Ma ! saída: número de Mach
print *, "dx = ", dx ! saída: distância entre sítios da rede na direção x
print *, "dy = ", dy ! saída: distância entre sítios da rede na direção y
print *, "dt = ", dt ! saída: passo de tempo
print *, "omega = ", omega ! saída: fator de relaxação
print *, “Appuyez sur une touche pour continuer”
read *
do j = 0, m ! conditions initiales: densité, vx, vy
do i = 0, n
rho(i,j) = rhoo
u(i,j) = 0.0
v(i,j) = 0.0
enddo
enddo
do j = 1, m-1 ! conditions aux limites: entrée (pas glisser pour j = 0 & j = m)
u(0,j) = uo
v(0,j) = 0.0
enddo
!!!
!!!
! main loop
do kk = 1, mstep ! boucle principale
call collision(u,v,f,feq,rho,omega,w,cx,cy,n,m) ! l'étape de collision
call streaming(f,n,m) ! l'étape de propagation
call boundcond(f,n,m,uo) ! conditions aux limites
call calcrhouv(f,rho,u,v,cx,cy,n,m) ! l'obtention la vitesse du fluide
m2=m/2 ! milieu d’enregistrement (la direction y)
print *,kk,v(0,m2),rho(0,m2),u(n,m2),v(n,m2),rho(n,m2) ! sortie (écran): Les valeurs de point median
enddo
call output(u,v,rho,uo,n,m)! sortie des valeurs (de fichiers)
end
! kk,
!!!
!! sous-programme à l’étape de collision
!!!
subroutine collision(u,v,f,feq,rho,omega,w,cx,cy,n,m)! l’étape de collision
REAL8 f(0:8,0:n,0:m), feq(0:8,0:n,0:m) ! fonction de distribution de particules: Réseau D2Q9
REAL8 w(0:8), cx(0:8), cy(0:8) ! des poids, la vitesse (x et y) du réseau D2Q9
REAL*8 rho(0:n,0:m), u(0:n,0:m), v(0:n,0:m) ! la densité, la vitesse (x et y) macroscopique
do i = 0, n
do j = 0, m
t1 = u(i,j)*u(i,j) + v(i,j)*v(i,j) ! vecteur de vitesse: Module
do k = 0, 8 ! vitesse de reseau D2Q9
t2 = u(i,j)*cx(k) + v(i,j)*cy(k) ! produit scalaire: v.c
feq(k,i,j) = rho(i,j) * w(k) * (1.0 + 3.0*t2 + 4.50*t2*t2-1.50*t1)! la fonction d'équilibre
f(k,i,j) = omega*feq(k,i,j) + (1.- omega)*f(k,i,j) ! fonction de distribution
enddo
enddo
enddo
return
end
!!!
! sous-programme à l’étape de propagation
!!!
subroutine streaming(f,n,m) ! etapa de propagação
REAL*8 f(0:8,0:n,0:m) ! fonction de distribution de particules:
! Para fins de varredura --> direções 1 & 3:isoladas, direções 5, 6, 7, 8, além de 2 & 4 (obrigatórias): conjunta
! Regras de ouro:função progressiva <-> varredura regressiva & função regressiva <-> varredura progressiva
do j = 0, m ! varredura progressiva (irrelevante, poderia ser regressiva)
do i = n, 1, -1 ! varredura regressiva obrigatória
f(1,i,j) = f(1,i-1,j) ! f1: função progressiva horizontal (não há “componente” vertical)
enddo
do i = 0, n-1 ! varredura progressiva obrigatória
f(3,i,j) = f(3,i+1,j) ! f3: função regressiva horizontal (não há “componente” vertical)
enddo
enddo
do j = m, 1, -1 ! varredura regressiva obrigatória
do i = 0, n ! varredura progressiva (irrelevante, poderia ser regressiva)
f(2,i,j) = f(2,i,j-1) ! f2: função progressiva vertical (não há “componente” horizontal)
enddo
do i = n, 1, -1 ! varredura regressiva obrigatória
f(5,i,j) = f(5,i-1,j-1) ! f5: função progressiva horizontal,progressiva vertical
enddo
do i = 0, n-1 ! varredura progressiva obrigatória
f(6,i,j) = f(6,i+1,j-1) ! f6: função regressiva horizontal,progressiva vertical
enddo
enddo
do j = 0, m-1 ! varredura progressiva obrigatória
do i = 0, n ! varredura progressiva (irrelevante, poderia ser regressiva)
f(4,i,j) = f(4,i,j+1) ! f4: função regressiva vertical (não há “componente” horizontal)
enddo
do i = 0, n-1 ! varredura progressiva obrigatória
f(7,i,j) = f(7,i+1,j+1) ! f7: função regressiva horizontal,regressiva vertical
enddo
do i = n, 1, -1 ! varredura regressiva obrigatória
f(8,i,j) = f(8,i-1,j+1) ! f8: função progressiva horizontal,regressiva vertical
enddo
enddo
return
end
!!!
! Sous-programme pour des conditions aux limites
!!!
subroutine boundcond(f,n,m,uo) ! conditions aux limites: fonction de distribution des particules
REAL*8 f(0:8,0:n,0:m) ! fonction de distribution de particules:
! frontière ouest: entrée
do j = 0, m
!-----------------------------
rhow =(f(0,0,j) + f(2,0,j) + f(4,0,j) + 2.(f(3,0,j)+f(6,0,j)+f(7,0,j)))/(1.- uo) !—>
f(1,0,j) = f(3,0,j) + 2.rhowuo/3. !—>
f(5,0,j) = f(7,0,j) + rhowuo/6. !—>
f(8,0,j) = f(6,0,j) + rhow*uo/6. !-----------------------------
enddo
!!!
!!!chang!!!
! frontière sud: rebondir contre le mur(bounce back)
do i = 0, n
f(2,i,0) = f(4,i,0)
f(5,i,0) = f(7,i,0)
f(6,i,0) = f(8,i,0)
enddo
! frontière nord: rebondir contre le mur(bounce back)
do i = 0, n
f(4,i,m) = f(2,i,m)
f(8,i,m) = f(6,i,m)
f(7,i,m) = f(5,i,m)
enddo
! frontière orientale: Sortie
do j = 1, m
f(1,n,j) = 2.*f(1,n-1,j) - f(1,n-2,j)
f(5,n,j) = 2.*f(5,n-1,j) - f(5,n-2,j)
f(8,n,j) = 2.*f(8,n-1,j) - f(8,n-2,j)
enddo
!!!
do i=0,40
f(2,i,20)=f(4,i,20)
f(5,i,20)=f(7,i,20)
f(6,i,20)=f(8,i,20)
end do
do j=0,20
f(1,40,j)=f(3,40,j)
f(5,40,j)=f(7,40,j)
f(8,40,j)=f(6,40,j)
end do
!!!
return
end
!!!
!
!Sous-programme pour des quantités macroscopiques
!
!!!
subroutine calcrhouv(f,rho,u,v,cx,cy,n,m) ! évaluation de variables macroscopiques
REAL8 f(0:8,0:n,0:m) ! fonction de distribution de particules:
REAL8 rho(0:n,0:m), u(0:n,0:m), v(0:n,0:m) ! densité du fluide,les vitesses (composantes x et y)
REAL*8 cx(0:8), cy(0:8) ! velocidades da rede D2Q9
do j = 0, m ! varredura vertical
do i = 0, n ! varredura horizontal
ssum = 0.0 ! inicialização do somatório para avaliar a densidade
do k = 0, 8 ! varredura: direções da rede D2Q9
ssum = ssum + f(k,i,j) ! somme p / évaluer la densité: SSUM = F0 + F1 + ... + f8
enddo
rho(i,j) = ssum ! densité du fluide à la position (i, j)
enddo
enddo
do i = 0, n ! varredura horizontal
do j = 0, m ! varredura vertical
usum = 0.0 ! inicialização do somatório para avaliar a velocidade-x
vsum = 0.0 ! inicialização do somatório para avaliar a velocidade-y
do k = 0, 8 ! varredura: direções da rede D2Q9
usum = usum + f(k,i,j)cx(k) ! somatório vx: usum = f0c0,x + f1c1,x + … + f8c8,x
vsum = vsum + f(k,i,j)cy(k) ! somatório vy: vsum = f0c0,y + f1c1,y + … + f8c8,y
enddo
u(i,j) = usum / rho(i,j) ! velocidade-x na posição (i,j)
v(i,j) = vsum / rho(i,j) ! velocidade-y na posição(i,j)
enddo
enddo
do j = 0, m
v(n,j) = 0.0 ! correction p / profil vitesse développé: la direction y
enddo
!!!
do j=0,20
do i=0,40
u(i,j)=0.0
v(i,j)=0.0
end do
end do
!!!
return
end
! sous-routine pour sorties (enregistrement) des résultats
subroutine output(u,v,rho,uo,n,m) ! sortie des résultats: le fichier journal
REAL8 u(0:n,0:m), v(0:n,0:m), rho(0:n,0:m) ! vitesse en x, vitesse en y, la densité
REAL8 stf(0:n,0:m) ! fonction actuelle
!n5= n/20
!n10 = n/10
!n50 = n/2
!n95 = 95*n/100
m2 = m/2
open(2,file=‘result.dat’)
write(2,)“VARIABLES =X, Y, U, V”
write(2,)“ZONE “,“I=”,n+1,“J=”,m+1,”,”,“F=BLOCK”
n100 = n/10
n200 = n/5
n500 = n/2
open(3,file=‘u_x different_100_200_500.dat’)
do j = 0,m
write(3,50) j, u(n/10,j)/uo, u(n/5,j)/uo, u(n/2,j)/uo
enddo
do j=0,m
write(2,*)(i,i=0,n)
end do
do j=0,m
write(2,*)(j,i=0,n)
end do
do j=0,m
write(2,*)(u(i,j),i=0,n)
end do
do j=0,m
write(2,*)(v(i,j),i=0,n)
end do
!open (8, file=‘excel_si.dat’) ! enregistrement de profils de vitesse dans différentes positions en x de CANAL
!do j = 0, m
! write(8,50) j, u(0,j), u(n5,j), u(n10,j), u(n50,j), u(n95,j)
!enddo
!write(10,*) ‘ux velocity at the vertical center geometry’
!do j = 0, m
! write(10,*) j, u(n/2,j)/uo
!end do
open (9, file=‘vmax.dat’) ! enregistrement de l’évolution de la vitesse maximale (vitesse à l’axe central)
do i = 0, n
write(9,50) i, u(i,m2)
enddo
50 format (I4, 6(2X,F11.8))
return
end