Contents
 
 

FORTRAN program
 

      PROGRAM   2D_Burgers_Equation

      Implicit none

      Integer Nbt,Nbx,Nby,n,i,j,m,p

      Real dt,dx,dy,Re,T,Pi,u(0:1100,0:1100),k,v(0:1100,0:1100)

      Character*80 F440

      F440='(210(F10.4,1X))'
      dx=0.001                                  !pas d,espace suivant x
      dy=0.001                                  !pas d'espace suivant y
      dt=0.005                                   !pas de temps
      Re=1000                                  !nombre de Reynolds
      k=5                                            !taux de decroissance du sinus
      p=5                                            !pas de sortie des resultats (spatial)
      write(*,*)'Donner la duree de la simulation: '
      Read(*,*) T  !duree de la simulation
      Nbt=int(T/dt)                         !nombre de pas temporel
      Nbx=int(1/dx)+1                  !nombre de pas spatial en x
      Nby=int(1/dy)+1                  !nombre de pas spatial en y
      Pi=4*atan(1.)                        ! nom bre Pi
      write(*,*)'Nb pas d''espace en x=',Nbx
      write(*,*)'Nb pas d''espace en y=',Nby
      write(*,*)'Nb pas de temps=',Nbt
      write(*,*)'n=0'

      Call Initialiser_vitesses_u_v(u,v,dx,dy,Nbx,Nby,Pi,k)                    !routine d'initialisation (CI) de u et v

      Open(1,file='vit2Du.dat')
      Open(2,file='vit2Dv.dat')

       do j=0,(Nby/p)
                   Write(1,F440) (u(p*i,p*j),i=0,(Nbx/p))                                 !ecriture de la solution initiale pour u
                   Write(2,F440) (v(p*i,p*j),i=0,(Nbx/p))                                 !ecriture de la soution initiale pour v
       end do

      Close(1)
      Close(2)

      write(*,*) 'Tracer les courbes ...'                                                         !trace des courbes (solutions initiales) sous MatLab
      pause
      m=10    !pas de sortie des resultats (temporel)

      do n=1,Nbt
                 Call Calculer_vitesse(u,v,Nbx,Nby,dx,dy,dt,Re,Pi,k,n)       !routine calculant la vitesse u
                 Call Calculer_v(u,v,dx,dy,Nbx,Nby)                                        !routine calculant la vitesse v
                 if (n.eq.m)  then
                            Open(1,file='vit2Du.dat')
                            Open(2,file='vit2Dv.dat')
                            do j=0,(Nby/p)
                                       Write(1,F440) (u(p*i,p*j),i=0,(Nbx/p))              !ecriture de la vitesse u
                                       Write(2,F440) (v(p*i,p*j),i=0,(Nbx/p))              !ecriture de la vitesse v
                            end do
                            Close(1)
                            Close(2)
                            m=m+10
                            write(*,*) 'Tracer les courbes ...'                                   !trace des courbes u et v sous MatLab a l'instant n*dt
                            pause
                 end if
      end do

      END
 

c**************************************************************

      SUBROUTINE Initialiser_vitesses_u_v(u,v,dx,dy,Nbx,Nby,Pi,k)

      Implicit none
      Real dx,dy,k,Pi,u(0:1100,0:1100),v(0:1100,0:1100)
      Integer i,j,Nbx,Nby

      do j=0,Nby
                 do i=0,Nbx
                            u(i,j)= sin(2*Pi*i*dx)*exp(-k*j*dy)
                            v(i,j)=0
                 end do
                u(0,j)=0
                u(Nbx,j)=0
      end do

      return
      End

c****************************************************************

      SUBROUTINE Calculer_vitesse(u,v,Nbx,Nby,dx,dy,dt,Re,Pi,k,n)

      Implicit none
      Real u(0:1100,0:1100),a(0:1100),b(0:1100),c(0:1100),f(0:1100),
     &     dt,Re,b_prime(0:1100),c_prime(0:1100),f_prime(0:1100),
     &     dy,Pi,k,dx,u1(0:1100,0:1100),v(0:1100,0:1100)
      Integer Nbx,Nby,i,j,n

      write(*,*)'Pas de tps ',n-1,'.5'
      do j=1,Nby-1
                 do i=1,Nbx-1
                              a(i)= -dt*u(i,j)/(4*dx)-dt/(2*Re*dx*dx)
                              b(i)= 1+dt/(Re*dx*dx)
                             c(i)= dt*u(i,j)/(4*dx)-dt/(2*Re*dx*dx)
                              f(i)= u(i,j)-0.5*dt*( v(i,j)*( u(i,j+1)-u(i,j-1) )/(2*dy)
     &                               - ( u(i,j+1)-2*u(i,j)+u(i,j-1) )/(Re*dy*dy) )
                 end do

                 b_prime(1)=b(1)
                 f_prime(1)=f(1)
                 c_prime(1)=c(1)

                 do i=2,Nbx-1
                            c_prime(i)=c(i)
                            b_prime(i)=b(i)-c_prime(i-1)*a(i)/b_prime(i-1)
                            f_prime(i)=f(i)-f_prime(i-1)*a(i)/b_prime(i-1)
                 end do

                 u1(Nbx-1,j)=f_prime(Nbx-1)/b_prime(Nbx-1)

                 do i=Nbx-2,1,-1
                             u1(i,j)=(f_prime(i)-c_prime(i)*u1(i+1,j))/b_prime(i)
                 end do

      end do

      write(*,*)'Pas de tps ',n

      do i=0,Nbx
                 u1(i,0)=sin(2*Pi*i*dx)
                 u1(i,Nby)=sin(2*Pi*i*dx)*exp(-k)
      end do

      do i=1,Nbx-1

                 do j=1,Nby-1
                              a(j)= -dt*v(i,j)/(4*dy)-dt/(2*Re*dy*dy)
                              b(j)= 1+dt/(Re*dy*dy)
                              c(j)= dt*v(i,j)/(4*dy)-dt/(2*Re*dy*dy)
                              f(j)= u1(i,j)-0.5*dt*( u(i,j)*( u1(i+1,j)-u1(i-1,j) )/
     &                               (2*dx)
     &                               - ( u1(i+1,j)-2*u1(i,j)+u1(i-1,j) )/(Re*dx*dx) )
                 end do

                 f(1)=f(1)-a(1)*u1(i,0)
                 f(Nby-1)=f(Nby-1)-c(Nby-1)*u1(i,Nby)
                 b_prime(1)=b(1)
                 f_prime(1)=f(1)
                 c_prime(1)=c(1)

                 do j=2,Nby-1
                             c_prime(j)=c(j)
                             b_prime(j)=b(j)-c_prime(j-1)*a(j)/b_prime(j-1)
                             f_prime(j)=f(j)-f_prime(j-1)*a(j)/b_prime(j-1)
                 end do

                 u(i,Nby-1)=f_prime(Nby-1)/b_prime(Nby-1)

                 do j=Nby-2,1,-1
                            u(i,j)=(f_prime(j)-c_prime(j)*u(i,j+1))/b_prime(j)
                 end do
      end do

      return
      END

c***************************************************************
 

      SUBROUTINE Calculer_v(u,v,dx,dy,Nbx,Nby)

      Implicit none
      Real u(0:1100,0:1100),v(0:1100,0:1100),dx,dy,der_u(0:1100,0:1100)
      Integer i,j,Nbx,Nby
 

      do j=0,Nby
                 der_u(0,j)=(-3*u(0,j)+4*u(1,j)-u(2,j))/(2*dx)
                 do i=1,Nbx-1
                            der_u(i,j)=(u(i+1,j)-u(i-1,j))/(2*dx)
                 end do
                 der_u(Nbx,j)=(3*u(Nbx,j)-4*u(Nbx-1,j)+u(Nbx-2,j))/(2*dx)
      end do

      do i=0,Nbx
                 do j=0,Nby-1
                            v(i,j)=-(der_u(i,j)+der_u(i,j+1))*dy/2
                 end do
                 v(i,Nby)= -(der_u(i,Nby)+der_u(i,Nby-1))*dy/2
      end do

      return
      END

   Contents