- #1

- 1

- 0

Please help me find where is the mistake. I copied the program from an old book.

DIMENSION S(50),V(50) ! Element stress and velocity

DIMENSION U(51) ! Nodal displacement

! Double arrays for tabular or graph representation of stress discontinuity

DIMENSION X(102),SP(102)

INTEGER J,JP

WRITE(*,*) 'Number of elements N<=50' ! Imput data

READ(*,*) N

WRITE(*,*) 'Number of time steps NT'

READ(*,*) NT

NP=2*N

DO 10 J=1,N ! Initial state of the system at T=O.

S(J)=0.

V(J)=0.

10 U(J)=0.

DX=1./N ! Length of element

DT=DX ! Time increment

DO 20 I=1,NT ! Beginning of time loop

T=I*DT ! Current time

S1=-1. ! Boundary conditions at free end, J=1

V1=V(1)-S1+S(1)

U(1)=U(1)+V1*DT ! Displacement at J=1

DO 30 J=2,N ! Beginning of element loop

V2=0.5*(V(J)+V(J-1)+S(J)-S(J-1)) ! V and S at right border

S2=S(J-1)+V2-V(J-1) ! of the element J-1

U(J)=U(J)+V2*DT ! Displacement at node J

V(J-1)=V1+V2-V(J-1) ! New values of V and S

S(J-1)=S1+S2-S(J-1) ! for element J-1

V1=V2 ! V and S at left border

S1=S2 ! of element J

30 CONTINUE ! End of element loop

V2=0. ! Boundary conditions at fixed

S2=S(N)+V2-V(N) ! end, J=N+1

V(N)=V1+V2-V(N) ! New values of V and S

S(N)=S1+S2-S(N) ! for element J=N

20 CONTINUE ! End of time loop

! To output of stress discontinuity along the rod length

X(1)=0.

X(2)=0.

SP(1)=S(1)

DO 40 JP=3,NP+2

X(JP)=X(JP-2)+DX

SP(JP-1)=S(JP/2)

40 CONTINUE

OPEN(5,FILE='STRESS') ! Element stress vs. rod length at I=NT

WRITE(5,*)'Element stress vs. rod length'

WRITE(5,*)' X S(X)'

WRITE(5,50)(X(JP+1),SP(JP),JP=1,NP)

50 FORMAT(2X,2(E12.5,2X))

OPEN(7,FILE='RESULTS') ! Control results for NT=2*N

WRITE(7,*) 'Control results for NT=2*N'

WRITE(7,*) 'Input data: N=',N,' ','NT=',NT

WRITE(7,*) 'Benchmark: T=2, U(1 )=2, S(N)=-2, S2=-2'

WRITE(7,*) 'Calculated: T=',T,' ','U(1)=',U(1)

WRITE(7,*) 'S(N)=',S(N),' ','S2=',S2

STOP

END ! End of Program

Thank you very much.