velocity with LBM

[saida29]

hello eveyonerni want to determinate the velocity of a flow in channel. i use for this lattice boltzmann method. the problem is that i get 0 despite that when I view usum which corresponds to the sum function distibution so the speed I find a non-zero value. here is my code. I counted a lot on your support and your understanding. thank you very much.

Code:

rnrn ! declaration de variablesrn rn parameter (n=1000,m=50)rn real f(0:8,0:n,0:m)rn real feq(0:8,0:n,0:m),rho(0:n,0:m)rn real w(0:8), cx(0:8),cy(0:8)rn real u(0:n,0:m), v(0:n,0:m)rn real g(0:8,0:n,0:m), geq(0:8,0:n,0:m),th(0:n,0:m)rn real theta(0:n,0:m),theta_in(0:m)rn real t_in(0:m)rn real Y(0:m)rn integer i,j,k,kkrn character Saida*2910rnrn ! opening filesrnrn Saida=’results/chanflowartth480v2.dat’rn open(12,file=Saida)rn !Saida=’results/chanflowartth200.dat’rn !open(3,file=Saida)rn !Saida=’results/chanflowartth40.dat’rn !open(4,file=Saida)rn !Saida=’results/chanflowartth80.dat’rn !open(5,file=Saida)rn !Saida=’results/chanflowartstrli.dat’rn !open(6,file=Saida)rn Saida=’results/chanflowartxvit200v2.dat’rn open(17,file=Saida)rn !Saida=’results/chanflowartxvit40.dat’rn !open(8,file=saida)rn !Saida=’results/chanflowartxvit120.dat’rn !open(9,file=Saida)rn rn !initialisations de variablesrnrn cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)rn cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)rn w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)rn uo=0.4rn !sumvelo=0.0rn rhof=999.1rn dx=1.0rn dy=dxrn dt=1.0rn tw=0.0rn t_in(m/2)=0.9rn th=0.5rn !H=0.02rn rn rn rn g=0.0rn visco=0.025rn rn pr=1rn alpha=visco/prrn !Re=uo*m/alpharn Re=m*uo/viscorn print *, ‘Re=’, Rern omega=1.0/(3.*visco+0.5)rn omegat=1.0/(3.*alpha+0.5)rn print *,omega,omegatrn mstep=10rn rn do i=0,nrn do j=0,mrn rho(i,j)=rhof rn u(i,j)=0.0rn v(i,j)=0.0rn end dorn end dorn rn ! main looprnrn!/////////////////programme principale//////////////////////////rnrn do kk=1,msteprn call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)rn call streaming(f,n,m)rn call sfbound(f,n,m,uo)rn call rhouv(f,rho,u,v,cx,cy,n,m)rn rn ! collestion for scalarrn call col(u,v,g,geq,th,theta,omegat,w,cx,cy,n,m)rn ! streaming for scalarrn !call streaming(g,n,m)!!!!!!!!!!!! n’existe pasrn call gbound(g,theta_in,Y,w,n,m)rn call tcalcu(g,th,theta,t_in,tw,n,m)rn !call result(u,v,rho,theta,Y,uo,n,m)rn END DOrn call result(u,v,rho,theta,Y,uo,n,m)! end of the main looprnrn !call result(u,v,rho,theta,Y,uo,n,m)rn !do j=0,mrn !print *,theta(12*n/25,j) rn !print *, u(n/5,j),u(n/25,j)rn !end dorn stoprn endrn ! end of the main programrnrn!/////////////////////les subroutines////////////////////////////rnrn !////////////collision////////////////rnrn subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)rn real f(0:8,0:n,0:m)rn real feq(0:8,0:n,0:m),rho(0:n,0:m)rn real w(0:8), cx(0:8),cy(0:8)rn real u(0:n,0:m), v(0:n,0:m)rn DO i=0,nrn DO j=0,mrn t1=(u(i,j)*u(i,j))+(v(i,j)*v(i,j))rn DO k=0,8rn t2=(u(i,j)*cx(k))+(v(i,j)*cy(k))rn feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)rn f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)rn !print*,f(k,i,j)rn END DOrn END DOrn END DOrn rn rn returnrn endrnrn !///////////////////col////////////////rnrn subroutine col(u,v,g,geq,theta,th,omegat,w,cx,cy,n,m)rn real g(0:8,0:n,0:m),geq(0:8,0:n,0:m),theta(0:n,0:m)rn real th(0:n,0:m),t_in(m/2),twrn real w(0:8),cx(0:8),cy(0:8)rn real u(0:n,0:m),v(0:n,0:m)rn rn do i=0,nrn do j=0,mrn rn do k=0,8rn geq(k,i,j)=theta(i,j)*w(k)*(1.0+3.0*(u(i,j)*cx(k)+v(i,j)*cy(k)))rn g(k,i,j)=omegat*geq(k,i,j)+(1.0-omegat)*g(k,i,j)rn rn end dorn end dorn end dorn !print*,thetarn returnrn endrnrn !///////////////////propagation///////////////rnrn subroutine streaming(f,n,m)rn real f(0:8,0:n,0:m)rn ! streamingrn DO j=0,mrn DO i=n,1,-1 !RIGHT TO LEFTrn f(1,i,j)=f(1,i-1,j)rn rn END DOrn DO i=0,n-1 !LEFT TO RIGHTrn f(3,i,j)=f(3,i+1,j)rn END DOrn END DOrn DO j=m,1,-1 !TOP TO BOTTOMrn DO i=0,nrn f(2,i,j)=f(2,i,j-1)rn END DOrn DO i=n,1,-1rn f(5,i,j)=f(5,i-1,j-1)rn END DOrn DO i=0,n-1rn f(6,i,j)=f(6,i+1,j-1)rn END DOrn END DOrn DO j=0,m-1 !BOTTOM TO TOPrn DO i=0,nrn f(4,i,j)=f(4,i,j+1)rn END DOrn DO i=0,n-1rn f(7,i,j)=f(7,i+1,j+1)rn END DOrn DO i=n,1,-1rn f(8,i,j)=f(8,i-1,j+1)rn END DOrn END DOrn returnrn endrnrn !//////////////conditions aux limites////////////////rnrn rn subroutine sfbound(f,n,m,uo)rnrn real f(0:8,0:n,0:m)rnrn do j=0,mrn !bounce back on west boundary!rn rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)rn & +2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))/(1.-uo)rn f(1,0,j)=f(3,0,j)+2.*rhow*uo/3rn f(5,0,j)=f(7,0,j)+rhow*uo/6rn f(8,0,j)=f(6,0,j)+rhow*uo/6rnrn end dornrn !bounce back on south boundary!rn do i=0,nrn f(2,i,0)=f(4,i,0)rn f(5,i,0)=f(7,i,0)rn f(6,i,0)=f(8,i,0)rn end dornrn !bounce back, north boundaryrn do i=0,nrn f(4,i,m)=f(2,i,m)rn f(8,i,m)=f(6,i,m)rn f(7,i,m)=f(5,i,m)rn end dorn rn ! accont for open boundary condition at the outletrn do j=0,mrn f(1,n,j)=2.*f(1,n-1,j)-f(1,n-2,j)rn f(5,n,j)=2.*f(5,n-1,j)-f(5,n-2,j)rn f(8,n,j)=2.*f(8,n-1,j)-f(8,n-2,j)rn end do rnrn returnrn endrnrn !/////////// cond aux limites sur g/////////////rnrn subroutine gbound(g,theta_in,Y,w,n,m)rn real g(0:8,0:n,0:m)rn real w(0:8),theta_in(0:m)rn real Y(0:m)rn ! Boundary conditionsrn ! top boundary condition, theta=0.rn do j=0,mrn g(4,0,j)=-g(2,0,j)rn g(7,0,j)=-g(5,0,j)rn g(8,0,j)=-g(6,0,j)rn g(1,0,j)=-g(3,0,j)rn end dorn ! east boundary condition, outlet cond ouverte.rn do j=0,mrn rn g(8,n,j)=2*g(8,n-1,j)-g(8,n-2,j)rn g(5,n,j)=2*g(5,n-1,j)-g(5,n-2,j)rn g(1,n,j)=2*g(1,n-1,j)-g(1,n-2,j)rn rn end dorn ! inlet west boundary conditions, T=t_in et theta=theta_in dirichletrn do j=0,mrnrn Y(j)=j/float(m)rn theta_in(j)=4*(Y(j)-(Y(j)*Y(j)))rnrn g(8,0,j)=theta_in(j)*(w(8)+w(6))-g(6,0,j)rn g(5,0,j)=theta_in(j)*(w(7)+w(5))-g(7,0,j)rn g(1,0,j)=theta_in(j)*(w(1)+w(3))-g(3,0,j)rn !print*,theta_in(j)rn end dorn !Bottom boundary conditions,Theta=0 rnrn do i=0,nrn g(1,i,0)=-g(3,i,0)rn g(2,i,0)=-g(4,i,0)rn rn rn g(5,i,0)=-g(7,i,0)rn g(6,i,0)=-g(8,i,0)rn rn end dorn returnrn endrnrn!/////////calcul de la temperature///////////rnrn subroutine tcalcu(g,th,theta,t_in,tw,n,m)rn real g(0:8,0:n,0:m),th(0:n,0:m),theta(0:n,0:m)rn real t_in(m/2),twrnrn t_in(m/2)=1.0rn tw=0.0rnrn do j=1,m-1rn do i=1,n-1rn ssumt=0.0rn do k=0,8rn ssumt=ssumt+g(k,i,j)rn end dorn th(i,j)=ssumtrn rn end dorn end dorn rn do j=1,m-1rn do i=1,n-1rn theta(i,j)=(th(i,j)-tw)/(t_in(m/2)-tw)rn end dorn end dorn rn returnrn endrnrn!/////////////calcul de la densite et precisement la vitesse/////////////rnrn subroutine rhouv(f,rho,u,v,cx,cy,n,m)rnrn real f(0:8,0:n,0:m),rho(0:n,0:m)rn real u(0:n,0:m),v(0:n,0:m)rn real cx(0:8),cy(0:8)rnrn do j=0,mrn do i=0,nrn ssum=0.0rn do k=0,8rn ssum=ssum+f(k,i,j)rn end dorn rho(i,j)=ssumrn end dorn end dornrn do j=0,mrn rho(0,j)=(f(0,0,j)+f(2,0,j)+f(4,0,j)rn & +2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))/(1.-uo)rn end dornrn do i=0,nrn do j=0,mrn usum=0.0rn vsum=0.0rn do k=0,8rn usum=usum+f(k,i,j)*cx(k)rn vsum=vsum+f(k,i,j)*cy(k)rn print*, usumrnrn end dorn u(i,j)=usum/rho(i,j)rn v(i,j)=vsum/rho(i,j)rn print*,u(i,j)rn end dorn end dorn rn ! do j=1,mrn ! v(n,j)=0.0rn ! end dorn !do i=1,nrn !u(i,0)=0.0rn !u(i,m)=0.0rn !end dorn rnrn returnrn endrnrn!////////////////calcul de lignes de courant et de resultat/////////////////rnrn subroutine result(u,v,rho,theta,Y,uo,n,m)rn real u(0:n,0:m),v(0:n,0:m),theta(0:n,0:m)rn real rho(0:n,0:m),strf(0:n,0:m)rn real Y(0:m)rn ! streamfunction calculationsrn strf(0,0)=0.rn do i=1,nrn rhoav=0.5*(rho(i-1,0)+rho(i,0))rn strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))rn do j=1,mrn rhom=0.5*(rho(i,j)+rho(i,j-1))rn strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))rn end dorn end dorn rn ! enregistrement de resultatrn rn do j=0,mrn Y(j)=j/float(m)rn rn end dornrn do j=0,mrn rn write(12,*)Y(j),theta(12*n/25,j)rn rn end dorn rn rn rnrn !do j=0,mrn rn !write(3,*)Y(j),theta(n/5,j)rn !end dorn rnrn rnrn !do j=0,mrn !write(4,*)Y(j),theta(n/25,j)rn !end dornrn rn !do j=0,mrn !write(5,*)Y(j),theta(2*n/25,j)rn !end dornrnrn !do j=0,mrn !do i=0,nrn !write(6,*)i,j,strf(i,j)rn !end dorn !end dornrn rnrn do j=0,mrn write(17,*)Y(j),u(100,j)rn !u(n/5,j)rn !print*,u(n/5,j)rn end dorn !pausern rnrn !do j=0,mrn !write(8,*)Y(j),u(n/25,j)rn !end dornrn rnrn !do j=0,mrn !write(9,*)Y(j),u(3*n/25,j)rn !end dornrn rnrn rnrn rn rn returnrn endrn

[mathias]

Dear saida29,rnrnI am n ot really a Fortan expert. To speed up, can you please indicate the line where the velocity is finally computed, shorten the code and comment it.rnrnMathias

[Author]

Posts