Hey, what's going on?

Problem Solving Series – Laminar Velocity Profile Between Infinite Parallel Plates

Posted by Syeilendra Pramuditya on October 12, 2012

The Document

The Code

!     Laminar velocity profile between infinite parallel plates
!     By Syeilendra Pramuditya - https://syeilendrapramuditya.wordpress.com
!     October 2012
      program main
      implicit none
      integer j,jmin,jmax,iter,itermax,conv,prnt
      parameter(jmin=0, jmax=10)
      real*8 uold(jmin:jmax),unew(jmin:jmax),L,uinit,dy,errmax,
     &mu,rho,Re,f,mdpdx,alpha,err,ubot,utop,y(jmin:jmax),
     &umath(jmin:jmax)

      !input
      L=1.0d-2
      ubot=0.0d0
      utop=0.0d0
      itermax=1000
      errmax=1.0d-6
      Re=5.0d2

      !props
      mu=1.0d-3
      rho=1.0d3

      !calc internal vars
      uinit=Re*mu/(rho*L) !average x-veloc
      f=24.0d0/Re !friction factor
      mdpdx=(f/L)*0.5d0*rho*uinit*uinit ! -dp/dx
      alpha=0.5d0*mdpdx/mu

      !init all arrays
      do j=jmin,jmax
      uold(j)=0.0d0
      unew(j)=0.0d0
      end do

      !mesh
      dy=L/(jmax-jmin)
      y(jmin)=0.0d0
      do j=jmin+1,jmax
      y(j)=y(j-1)+dy
      end do

      !init guess of x-velocity
      do j=jmin+1,jmax-1
      uold(j)=uinit
      end do

      !init boundary condition
      uold(jmin)=ubot
      uold(jmax)=utop
      unew(jmin)=ubot
      unew(jmax)=utop

      !calculate
      conv=0
      iter=0
      do while((conv.eq.0).and.(iter.le.itermax))
      conv=1
       !calc new u
       do j=jmin+1,jmax-1
       unew(j)=0.5d0*(uold(j+1)+uold(j-1))+alpha*dy*dy
       err=dabs(unew(j)-uold(j))/unew(j)
       if(err.ge.errmax) conv=0 !check convergence
       end do
       !update u
       do j=jmin,jmax
       uold(j)=unew(j)
       end do
      iter=iter+1
      end do

      !exact solution for static walls
      do j=jmin,jmax
      umath(j)=0.5d0*(1.0d0/mu)*mdpdx*(L*y(j)-y(j)*y(j))
      end do

      !result
      prnt=6
      write(prnt,*) 'Re =',Re
      write(prnt,*) 'Solution Converged in',iter,'iterations'
      write(prnt,*) ' '
      write(prnt,*) '      y [m]      u[m/s]  umath[m/s]'
      do j=jmin,jmax
      write(prnt,1000) real(y(j)),real(unew(j)),real(umath(j))
      end do
      write(prnt,*) ' '
      write(prnt,*) 'Result has been saved to fort.10 file'
      prnt=10
      write(prnt,*) 'Re =',Re
      write(prnt,*) 'Solution Converged in',iter,'iterations'
      write(prnt,*) ' '
      write(prnt,*) '      y [m]      u[m/s]  umath[m/s]'
      do j=jmin,jmax
      write(prnt,1000) real(y(j)),real(unew(j)),real(umath(j))
      end do

 1000 format(1P20E12.3)

      stop
      end

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: