salvation has someone who helps me.
I have a program but does not give the same result with book
Mohammad, A. A., Applied Lattice Boltzmann
Method, SURE Print, Dalbrent, Canada, 2007
page 86
the program works so well that flow in a channel without obstacle, but if introduced obstacle stop the program until 1025 iteration stop.
I am looking to run up to 10000 iteration.
!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
uo = 0.1 ! Vitesse d’entrée (m / s) [0.1 - 0.4]
rhoo = 1.0 ! densité de l’air (kg/m³) à 20 °C
! viscd = 0.1 ! viscosité dynamique (Pa)
!visc_lattice = viscd / rhoo ! viscosité cinematique (m²/s)
!visc_lattice = 0.0000154 ! viscosité cinematique (m²/s) à 20 °C
! 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)
!ly = hauteur_lattice* dy ! largeur de canal (m)
hauteur_lattice=80
!Ma = uo / cs ! Nombre de Mach [0.1 - 0.4]
! Re = uo * hauteur_lattice / visc_lattice ! Nombre de Reynolds
Re=400
visc_lattice=(uo*hauteur_lattice)/Re
! Instabilité si 2 <tau <0,506 !!!’
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 *, 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/m, 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
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