一维, 二维热传导代码大全

sillybear

http://www.sali.freeservers.com/engineering/cfd/cfd_codes.html

sillybear

The model equation to be solved is of the form:
d2T/dx2 = - S(x)  
***************************************************************************
*  THIS PROGRAM SOLVES 1-D POISSON EQUATION WITHIN THE DOMAIN 0<x<1 FOR   *
*  DIRICHLET BOUNDARY CONDITION USING FINITE DIFFERENCE METHOD (UNIFORM   *
*       GRID) AND DETERMINE GLOBAL ERROR BY COMPARING THE NUMERICAL       *
*                       AND ANALYTICAL SOLUTION                           *
***************************************************************************
      program pg1
      integer nd,i,n,nm1,nm2
      parameter (nd=9000)
      doubleprecision S,L(nd),D(nd),U(nd),C(nd),T(nd)
      doubleprecision T1,Tn,x,dx,q1,qn
*Defining Source Function S on the Right-Hand-Side of the Poisson equation.
*Setting S = 0 will change the Poisson equation into Laplace equation.
      S(x) = x
*
      write(6,*) '            ***************************************'
      write(6,*) 'Enter the Number of GRID POINTS:'
      read(*,*) n
      nm1 = n - 1
      nm2 = n - 2
      write(6,*) 'Enter the Value of Dependent Variable at the LEFT Boun
     $dary:'
      read(*,*) T(1)
      T1 = T(1)
      write(6,*) 'Enter the Value of Dependent Variable at the RIGHT Bou
     $ndary:'
      read(*,*) T(n)
      Tn = T(n)
      write(6,*) '             ****************************************'
      dx = (1.0d0 - 0d0)/nm1
*
*Construction of Elements in the Tridiagonal Matrix
      do 10 i = 2,nm1
         L(i) = 1.d0
         D(i) = -2.d0
         U(i) = 1.d0
10    continue
*
*Construction of RHS Vector
      x = 0d0
      do 20 i = 2,nm1
         x = x + dx
         C(i) = -S(x)*dx**2
20    continue
      C(2)   = C(2) - L(2)*T(1)
      C(nm1) = C(nm1) - U(nm1)*T(n)
*
*Solution of the Equation Using TDMA
      CALL TDMA(nm1,L,D,U,C,T)
      T(1) = T1
      T(n) = Tn
*
      write(6,*) ''
      write(6,*) 'The Numerical Solution is:'
      write(6,*)''
      write(6,30) (T(i),i=1,n)
30    format('',7(2x,e10.4))
*
*Calculation of Analytical Solution and % Global Error
*Note: If analytical solution is known, the comment can be removed from
*      the CALL statement below.
C     CALL ANALYT(n,dx,T1,Tn,T)
*
*Calculation of Heat Fluxes
      q1 = -(-T(3) + 4*T(2) - 3*T(1)) / (2*dx)
      qn = -(3*T(n) - 4*T(nm1) + T(nm2)) / (2*dx)
      write(6,*) ''
      write(6,*) 'The Heat Flux at LEFT Boundary is:',q1
      write(6,*) 'The Heat Flux at RIGHT Boundary is:',qn
      stop
      end

***SUBROUTINE: ANALYTICAL SOLUTION***************************************
*
      subroutine ANALYT(n,dx,T1,Tn,T)
      integer nd,i,n
      parameter (nd=9000)
      doubleprecision TA,T(nd),TANL(nd),E(nd)
      doubleprecision x,dx,T1,Tn,maxerr
      parameter (eps=1.0d-8)
*Defining the Analytical Solution (solution depends on boundary conditions)
      TA(x,T1,Tn) = T1 + (Tn - T1 + 1.0d0/6.0d0)*x - x**3/6.0d0
*      
      x = 0d0
      do 10 i = 1,n
         x = x + dx
         TANL(i) = TA(x-dx,T1,Tn)
         E(i) = 0d0
         if (TANL(i).gt.eps) then
            E(i) = abs(TANL(i) - T(i))*100.0d0 / TANL(i)
         endif
10    continue
*
      write(6,*)''
      write(6,*) 'The Analytical Solution is:'
      write(6,*)''
      write(6,20) (TANL(i),i=1,n)
20    format('',7(2x,e10.4))
c     write(6,*)''
c     write(6,*) 'The Global Error:'
c     write(6,20) (E(i),i=1,n)
*
*Calculation of Maximum Global Error
      maxerr = 0d0
      do 30 i = 1,n
         if (E(i).gt.maxerr) then
            maxerr = E(i)
         endif
30    continue
      write(6,*)''
      write(6,*) 'The Maximum % Global Error is  ',maxerr
      return
      end
***SUBROUTINE: TDMA********************************************************
*
      subroutine TDMA(nm1,L,D,U,C,X)
      integer nd,i,nm1
      parameter (nd=9000)
      doubleprecision L(nd),D(nd),U(nd),C(nd),X(nd),P(nd),Q(nd),temp
*
      L(2) = 0d0
      U(nm1) = 0d0
*
*Forward Elimination
      do 10 i = 2,nm1
         temp = D(i) + L(i)*P(i-1)
         P(i) = -U(i) / temp
         Q(i) = (C(i) - L(i)*Q(i-1)) / temp
10    continue
*
*Back Substitution
      do 20 i = nm1,2,-1
         X(i) = P(i)*X(i+1) + Q(i)
20    continue
      return
      end
***************************************************************************