CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   troble in code for wave equation by lassonen method (https://www.cfd-online.com/Forums/main/128100-troble-code-wave-equation-lassonen-method.html)

 nargess_md January 3, 2014 08:32

troble in code for wave equation by lassonen method

Hi every body
I have written a code in FORTRAN 95 to compute the solution of wave equation dT/dt=a* d^2T/dx^2 by laasonen method, but there is a problem:
the problem is a wall 1ft thick and infinite in other directions. it's initially in constant temperature 100 F and suddenly temperature of two sides is increased and maintained at 300 F. I wrote the following code to find the solution using the Laasonen method.. but there is a problem that my answers contains numbers more than 300 F near the wall. can sb help me about that:

program lassonen
implicit none
integer,parameter::mmax=100,mmin=-100
integer::i,j,n,m,k
real,parameter::alpha=0.1
real::diff,a,b,c,delta_t,delta_x,total_t,x,P,L
real,dimension(mmin:mmax)::G,H,D
real,dimension(mmin:mmax)::T,U

open(unit=3,file="1.3.txt")
!------------------------------------------------------------
!-----------------------------getting the input--------------
!------------------------------------------------------------
write(*,*)'enter delta_x:'
write(*,*)'enter delta_t:'
write(*,*)'enter total length:'
write(*,*)'enter total time:'

diff=(alpha*delta_t/(delta_x**2))!diffusion number
!-------------------------------------------------------------
!------------------------defining matrix k elements-----------
!-------------------------------------------------------------

a=diff
b=-(1+2*diff)
c=diff

P=x/delta_x
L=(total_t/delta_t)+1

m=P+1

!------------------------------------------------------------
!-----------------------initial condition--------------------
!------------------------------------------------------------
do k=2,m-1
T(k)=100
end do
!------------------------------------------------------------
!-----------------------boundary condition-------------------
!------------------------------------------------------------
!do n=1,L
T(1)=300
T(m)=300
!end do

!-------------------------------------------------------------
!-------------------------D deffinition-----------------------
!-------------------------------------------------------------
do n=2,L

D(2)=-T(2)-a*T(1)
do i=3,m-2
D(i)=-T(i)
end do
D(m-1)=-T(m-1)-c*T(m)

!-------------------------------------------------------------
!------------------------H & G deffinition--------------------

H(1)=0
G(1)=300
do i=2,m-1
H(i)=c/(b-a*H(i-1))
G(i)=(D(i)-a*G(i-1))/(b-a*H(i-1))
!write(*,*)i, H(i), G(i)
end do

do j=m-1,2,-1
T(j)=-H(j)*T(j+1)+G(j)
write(*,*) j, T(j)
end do
-

do i=1,m
write(3,*) i, T(i)
end do

end do

end program lassonen

I really gave up.
I have no idea to fix this.