Variation of stress resultants in concrete structures due to time-dependent creep and shrinkage strains of concrete

Compuf~rs & S~rucrurvr Vol. 18. No. 1. pp 117-126. 1984 Printed in Great Britain.

CM-7949/84 $3.00 t .-xl 0 1984 Pergamon Press LLd

VARIATION OF STRESS RESULTANTS IN CONCRETE STRUCTURES DUE TO TIME-DEPENDENT CREEP AND SHRINKAGE STRAINS OF CONCRETE G.D. STEFANOU University of Patras, Patras, Greece (Received 8 November 1981;receioed for publication 23 March 1982) Abstr&--Ageneral procedure for e~culating the vocation with time of the internal stress resultants and hence the stresses in concrete structures is discussed. In particular, a study is made of the changes in the stress resultants due to time-dependent creep and shrinkage strains of concrete. A general procedure of calculating the variation in the stress resultants due to differential creep strains in concrete structures has been proposed by the author[l]. A similar procedure is followedin this paper to study these variations when creep and shrinkage strains take place simultaneously. The method leads to a system of n-linear differential equations of the form: k=AXItB

the solution of which is performed by a computer using Runge-Kutta numerical procedures. A reinforced concrete portal frame exhibiting creep and shrinkage strains is solved by the proposed method and the results are given in tabular and graphical form. I.INTRODUCTION

The general problem of calculating the intern& distribution of stress in concrete structures due to timedependent creep or shrinkage strains can be treated in the following manner. It is assumed that the frame exhibits time-dependent stresses when its members obey to a sum of different creep functions (o(t). Sufficient number of cuts are made in the frame[21. The arbitrary stress resultants at the cuts of a frame with n-indeterminacies are given by eqn (1). [x]‘=[x,x~...xK...x”].

(1)

At time (t) time function cp(t) assumes the values: qPl(O, cpz(%..,~&~,..*, qA(f), and at time (t + dt) the change in these functions is: dq,, dq2,. . . , d4o,,.. . , dq,. Due to load [P] the relative displacements at the cuts in the direction of arbitrary constant Xk will be given in terms of time functions 9eft). At time interval dt the relative displacements at the cuts and hence the changes of these displacements result from the summation of relative displacements of all remaining members. The displacements (S& at the cut due to the external load [PI in the direction of arbitrary stress resultant X, must be considered. These displacements are governed by time function cpz(t). The change of the relative displacement d URtxti in the direction of Xk is given by the following expression:

the distribution of stresses in concrete structures due to time-dependent strains is given in Ref. [1] and a brief outline is given in Appendix 1. 2.FOR~LATIONOFS~KA~EPROBL~

To introduce the effect of shrinkage strains in the calculation of time-varying stress resultants we assume that shrinkage strains are developing in a similar manner as creep strains, i.e. “creep analogy” holds also for shrinkage 12, 31, At time interval dt, the relative displacement URcxti between the faces of the arbitrary cut in the direction of stress resultant X, due to shrinkage strains of all members of the structure is governed by a time function 9&t), will be given by eqn (3).

where dpE is the change of time function v*(t) in the interval dt, and 9e- is the value of time function 9&t) for time t = a. It follows that the total relative displacements Z U, of the faces of one cut in the direction of the arbitrary stress resultant Xk due to shrinkage strains of all members of the structure in the interval dt, will be given by the sum

In matrix form eqn (4) becomes: = tsKo,,sKo,...sKo,,...6K”.~l

A full discription of the general problem of calculating 117

G.

118

D. STEFANOU

Considering now the total relative displacements URtXtinin the directions of all n-stress resultants due to shrinkage strains of all members of the structure in the interval dt we shall obtain the n-relative displacements in all arbitrary cuts.

(c) Due to creep strains. Ww1[&,21.

. . b&,1.

. .

r [xl 101. . . [O]. . . [O]WI [xl

u

R(X&

= -

&“,l &o, .

RIO,,.. ho,*

s20.1 820.2 . .

szo.*.

. f

. . S2O.h

=

. . . LOI . . . [Ol

.

X

6KO.C.

. . . . . .

1

~,,,.A

. .

*

(8)

. . . WI

LOI WI . 1 .I01 . . . WI_

6.0.2,

In matrix form eqn (5) becomes:

uR(Xkh

.

[Ol WI . * * 1-w

SKO.1 SK,,, .

S”,,

[At.Jl*

.2..

-[AK0,3]*[zz..

. s]’

(5a)

_

4. FORMULATIONOF THE TIME-DEPENDENTCREEP ANDSRRINK-

(d) Due to shrinkage strains. They are given by eqns (5) or (Sa). The resultant relative displacements are obtained by superposition, assuming continuity of the displacements is ensured at the cuts. This will give rise to the following compatibility system. [AKJ

[d-V

+

[[AK~.II

[AK~.z]

. . . [AK~.cI

. . . [A~~.*11

AGE PROBLEM

The total solution for the time interval dt is obtained by considering the effects on the stress resultants of all relative displacements due to: (a) The external loads [PI; (b) The changes [dX] in the stress resultants; (c) The creep strains; and (d) The shrinkage strains. This is possible by assuming that compatibility of discontinuity is maintained and be equating the sum of all four relative displacements to zero. These displacements are given by the following expressions: (a) Due to loads [P] QIO.I bo.2 ho,, . . . [AKOLI

.

’

+ [fjKO.l

(9) aKO.2 . . . 8K0.5.

. . 6

~20.2 . . .

~KO.1 SK.2 . . . .

-0 0

= 0

(b) Due to loads [dXj:

s*z.. s ... . . . . . . .

&I

I

_o_

Ip..

..I

.

.-.

[AKJ

. . .

Im

.

.

6x1 SK7

[AKJ

.

..*...

.

+ [AKA

.

..a..

.

.

.

.

.

Dividing both sides of eqn (9) by time function we obtain the final system of n-linear differential equations eqn (10). A step-by-step solution of eqn (10) with predetermined steps will give, at any desired time interval, the changes in the internal stress resultants of the structural system and hence the distribution of stresses in its members.

.

SK,.

. . . . . 6

!1 WI ml

WI

PI

[xl

PI

:

:

:

(01 ro1 G-l

*

- [AKo.sI

(10)

Variation of stress resultants in concrete structures

Equation (10) can also be written in the

form:

[A~pl[fi + bw1

119 APPENDIX 1

Numerical example A reinforced concrete one-bay two-storey portal shown in Fig.

1 is analysed. Frame data: uniformly distributed load, p, = p2 = 8.0tlm. Dimensions: length I = lO.Om, height h = 5.0m, beam section

[b] rbl..[xl

11’I b,

I+1 cp1-

b x d = 40 x 100cm, column section b X d 40 x 70 cm. Concrete strength at 28-days B300= 300kg/cm*. Values of coefficients: z = 1, cpO = 2, d, = 40cm, T = 20°C a, = I5 days, K, & K, from Figs. 2 and 3 of Ref. [3, 41. Age of concrete of 1st storey t = t, = 15days. At time t > 1, the time function for concrete is pi(r > t,) as shown in Fig. 4. At time t = t2= 30 days the age of concrete of 1st storey is 30 days and that of the 2nd is 15days. At time t > t2the time function for concrete is cpz(t > fJ. The time functions cp,(t> I,) and cp2(f> tz)

1

5. DISCUSSIONAND CONCLUSIONS

An attempt has been made in this paper to investigate the variation of stress resultants in concrete structures as a result of time dependent creep and shrinkage strains of concrete. Although several simplifying assumptions had to be made regarding creep and shrinkage behaviour of structural concrete, it has been possible to predict the variations in the stress resultants due to variations of creep and shrinkage strains with time. It can be seen from the output of the attached programme that the initial values of the stress resultants obtained from the solution of the linear system of differential equations, eqn (lo), ‘uctuate in a way which allows the structure to be in equilibrium at all time intervals. It i also shown that due to variations in the time dependent deformational characteristics of concrete redistribution of stress takes place in the structure. Furthermore, an asymptotic change in the stress resultants is observed leading to constant values. Experimental evidence showed that constant values could be reached at an approximate time of 4 years after loading.

Portal

Fig. 1.

185 170 155

15-

140 125 I IO Y,

'O-

124

Acknowledgemenrs-The author wishes to express his appreciation to Dipl. Ing., Ch. Dimas, and E. Karaoglanis for their assistance in the preparation of this work.

frame

2

10

Days, 1. Fig. 2. Theoretical age of concrete, (time I in days, coeff. K,, DIN4227).

REFERENCES

1. G. D. Stefanou, A general method of calculating stress redistribution due to differential creep in concrete structures. Comput. 2.

3. 4. 5. 6.

Structures,

198 1.

P. B. Morice, Linear Structural Analysis. Thames & Hudson, London (1%2). F. Leonhhardt and E. Monnig, Grundbagen zur Eemessung in Stohlbetonban, Sweite Au&-e. Springer-Verlag, Berlin (1972). H. Rusch, Stahlbeton-Spar&ton. Werkstofeingenschajten und Beemessungoerfahren. Verner-Verlag. GmbH., Dusseldorf (1973). DIN 4227/1972.For prestressed concrete. DIN 1045/1972.For plain and reinforced concrete.

Days, t.- a. Fig. 3. Theoretical duration (t - a) in days. Course of delayed elastic deformation coeff. K..

G.

D. STEFANOU

P t/m

P t/m

Porta 1 frame Time,

AX

months

Fig. 5.

Fig. 4. Time functions @, and @z. are calculated according to DIN 4227 and DIN 1045N of Refs. I5, 61. They are given by eqn (Al) and are shown in Fig. 4. ~(0 = c~n*(&~ - &,,I + 0.4*JL - a,

(Al)

follows: X, = R,’ , X2 =

M,,X3 = R,?, X4= Mz2.

(A3)

To calculate time function cps(1),(I= 1,2), eqns (AIHA3) and Figs. 2 and 3 are used. Results of function cp&t)are shown in Fig. 4. Values of cpOare taken from Table 1. The calculation of normalised time function F,(t) follows and the results are shown in Fig. 6. Time function F,(t), is. obtained if function dp* is divided by drp,, Ft = dqddrp, = pdcp,, (5 = 2,3,. . A),as shown in Fig. 7. Dividing the variation of stress resultant dJ$ by dq, we define time function cp](t) as a dependence function of stress resultants Xi, (i = 1, 2,. . . , n), X, = dXJdq,, (k = 1, 2,. . . , n). The analysis of the frame assumes the following stages of calculations.

To calculate the changes in the stress resultants due to the time-dependent strains, which are governed by time functions p,(t), and &t), the frame is cut as shown in Fig. 5, and the coefficients

(a) Stage 1. 15
The effective depth is calculated from eqn (A2).

W’) Coefficient C of eqn (A2) is taken from Table 1 and the effective age is calculated from eqn (A3) rw=~~t(Ttlo)

sKO,fr&p,PSKp,(K~= 1,2,3,4 and 5= L2) are calculated. The stress resultants are replaced by the arbitrary constants as

(A41

Table 1. Basic flow and basic shinkage dependent on the position of the construction member Position of Average relative oonstr. member ggy?.

Conslstencs

Field 82

Basic flow

BaSiC

@O 1

2

3

Coefficient C according

Shrinkage

to paragraph 5

E 510 4

5

0, 8

+10.10-5

30

90

1. 3

-10.10+

5,o

!%x~ly

70

2, 0

-25.10-~

1,5

indry atmorphera i.e. in 0y interiora

40

3. 0

-40.10’5

1.0

1

In the water

2

in verg dqmp atmosphere immediately over water

3 4

121

Variation of stress resultants in concrete structures

01 ’

I

0512

I

I

I

I

I

I

3

4

5

6

f

6

months

Tnne,

(a)

"05

I

2

3

4

5

Time,

6

7

8

months

(b)

(e)

Fig. 6. Values of normalised time functions F!(t), Fz(t)

Fig. 7. Time functions

(b) Stage 2. t = 30 days For this interval there exist stress resultants X, and X, in the second storey of the frame and variations by AX, and AX, in stress resultants X,” and X,’ respectively. These are calculated from compatibility system (AS).

- 2193 18274.8 --2193 7309.9

1180.2 - 6578.9 877.2 -2193

- 6578.9 t 58479.5 -8771.9 18274.8

9(t)

and normalised time functions F(t).

Superposition of eqns (A4) and (AS) will eventually give the final values of the initial conditions for the stress resultants.

1

877.2 -8771.9 2057.4 -2193

The mean value of time function F(r) throughout time intervals dt = 15 days is assumed to be constant. Variations of the stress resultants X,, X,, Xs, X,, of the analysed frame for a period of 8 months are shown in tabular form on the attached output of BASE I computer programme which was written to solve the linear differential system of eqn (10) using Runge-Kutta numerical procedures and graphically in Fig. 8.

100

*..-

_ .

..a

**.................

-..

\

1

1

\

Xiii‘)

. __*._........

a.TTFlr.

35

\

*

70

L

105

1

140

TN

1

17s

DAY'S

I

I

'I

240 X,0

\ ‘.

‘\

‘.._

\

‘\ ‘1. ‘--

-._____._

-_ --__

---_ -

-.-.-.

X4 (TM)

-----_____ _.-.-.

___._____.__ ---r:LY’_____

-.

-.--‘-‘-

---------___,_____

x, (J’)

I t

Fig. 8. Graphical representation of the changes in the stress resultants of the frame due to creep and shrinkage.

G. D. STEFANOU

122

APPENDIX ****9*+,3+**+,+++**+***~*~~**~*~*~~~*~*~~***~*~~**~~~~~~~** + CHANGES OF STRESS RESULTANTS DUE TO CREEP AND SHRINKAGE * ~++*9+***++t9++**+++*++9++++++9++9++3+++**~~~~~~*~**~***~*~ INITIAL VALUES ________-_-_---_ -4.150

-27.600

+1.750

-71.480

MATRIX Cl -_____-_ -7309,.9 +2193.0 -18274.8 +2193.0 -97820.3 +6578.9 -18274.8 +6578.9 -58479.8 +2193.0 -877.2 -8771.9

+2193.0 -877.2 +8771.9 -2057.4

MATRIX C2 --______ +o.o +o.o +o.o +o.o

+o.o

+o.o

+D.O +7309.9 -2193.0

+o.o +2193.0 +1180.2

MATRIX C3 ____-_-+7309.9 -21.9 +18274.8 -2193.0 +1180.2 -6578.9 +18274.8 -6578.9 +51169.6 +877.2 -6578.9 -2193.0

-2193.0 +877.2 -6578.9 +877.2

+o.o +o.o +o.o +o.o

MATRIX C4 _--____-219298.2 +97820.3 -657894.7 +87719.3 MATRIX CS -_-___--+o.o +o.o -219298.2 +97820.3 MATRIX C6 ____---__ +0.0 +219298.2 +o.o -97820.3 +657894.7 -219298.3 -87719.3 -97820.3

123

Variation of stress resultants in concrete structures X2L(TM) Xl(T) STEP F2 ______-_____________~~~~~~~-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~---~---_-____

X3(T)

X4( TM)

-15.6249 -27.5179 .42880 3.29400 -27.1618 -27.4339 .46160 3.29400 -38.7346 -27.3484 .49440 3.29400 -50.3166 -27.2613 .5272O 3.29400 -61.8804 -27.1729 .56000 3.29400 ~~~~~~~----------_~~~_~~~~-~-------~~~~~~~~~~~-----------------------_

+4.3803 +7.1165 +9.9520 +12.8797 +15.8923

-82.8899 -93.7841 -104.1258 -113.8795 -123.0109

-67.1219 -27.1200 .59040 1.68700 -72.2846 -27.0670 .62080 1.68700 -77.3665 -27.0137 -65120 1.68700 -82.3657 -26.9603 .68160 1.68700 -87.2805 -26.9067 .71200 1.68700 __-______________-__~~~~~~~~~-~----~~~~~~~~~~~~~----------------------

+17.2522 +18.6175 +19.9871 +2 1.3509 +22.7351

-125.2507 -127.3004 -129.1586 -130.8239 -132.2952

-89.3173 -26.8809 .72840 1.36300 -91.3314 -26.8551 .74480 1.36300 -93.3226 -26.8292 .76120 1.36300 -95.2910 -26.8034 .77760 1.36300 -97.2364 -26 -7776 -79400 1.36300 ~~~-~~-~~--~-----~~~~~~~~~~~~~~~~~-~~~~~~~----~~~-~~----______________

+23.3132 +23.8905 +24.4669 +25.0422 +25.6164

-132.4446 -132.5549 -132.6263 -132.6586 -132.6521

-98.9190 -26.7534 .81000 1.2SOOO -100.5808 -26.7293 .82600 1.25000 -102.2219 -26.7052 .84200 1.25000 -103.8422 -26.6810 .85800 1.25000 -105.4418 -26.6569 .87400 1.25000 ______________----__~~~~~-~~~~~~~~~~--~~~--~~~~~~~~~~----_--___--___-_

+26.1198 +26.6218 +27.1223 +27.6213 +28.1187

-132.4267 -132.1699 -131.8817 -131.5622 -131.2115

-106.8336 -26.6354 .88840 1.23000 -108.2088 -26.6139 .90280 1.23000 -109.5673 -26.5924 -91720 1.23000 -110.9093 -26.5710 -93160 1.23000 -112.2347 -26.5495 .94600 1.23000 ~~~~~~------------~~~~~~-~~~~~~~~~~~~~~~~~~--------------_---~~~~-____

+28.5563 +28.9924 +29.4271 +29.8604 +30.2921

-130.8419 -130.4479 -130.0296 -129.5871 -129.1205

+30.6051 +30.9173 +3 1.2286 +31 .S389 +31.8483

-128.7253 -128.3174 -127.8966 -127.4630 -127.0166

-118.4187 -26.4408 1.02080 1.10000 -119.9125 -26.4116 1.04160 1.10000 -121.3750 -26.3823 1.06240 1.10000 -122.8063 -26.3532 1.08320 1.10000 -124.2046 -26.3241 1.10400 1.10000 -----------------~---__--------____---____-~~~~~~~~-~~~~~~~--_________

+32.3795 +32.9071 +33 . 4310 +33.9511 +34.4675

-125.9425 -124.8292 -123.6769 -122.4861 -121.2569

-125.1855 -26.3022 1.12000 1.05000 -126.1469 -26.2804 1.13600 1.05000 -127.0908 1.15200 -26.2585 1.05000 -128.0173 -26.2368 1.16800 1.05000 -128.9264 -26.2151 1.18400 1.05000 -------__-------__------~~~~~~~~---~~~~~~~~~~___----_---_------~~~----

+34.8388 +35.2078 +3!i. 5744 +35.9387 +36.3006

-120.2218 -119.1665 -118.0911 -116.9959 -115.8808

-26.2075 -26.1999 -26.1923 -26.1848 -26.1772

+36.4267 +36.5525 +36.6781 +36 . 8033 +36.9282

-115.4859 -115.0886 -114.6889 -114.2868 -113.8824

-26.1504 -26.1238 -26.0972 -26.0707 -26.0443

+37.3605 +37.7891 +3s. 1137 +39.6344 +39.0510

-112.3895 -110.8680 -109.3181 -107.7401 -106.1343

1.20000 1.20000 1.20000 1.20000 1.20000

1.05000 1.05000 1.05000 1.05000 1.05000

.95680 .96760 .97840 .98920 1.00000

1.18960 1.19520 I.20080 1.20640 1.21200

-113.1847 -114.1256 -115.0574 -115.9800 -116.8934

-129.2405 -129.5525 -129.8624 -130.1702 -130.4759

1.03000 1.23200 -131 .S136 -132.5252 1.25200 1.03000 -133.5106 1.27200 1.03000 -134.4701 1.29200 1.03000 1.31200 1.03000 -135.4036 -----------------~~~---~___---~~__-~----_--------------_____-------___

-26.5336 -26.5177 -26.5019 -26.4860 -26.4702

124

G.D.STEFANOU 10

!

Q++*+**9*9+*+**+++***++**++++*+*9+9*+*+*9*+*9+*++**994**,*91++~4+

20 30 40 SO 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610

!

+

VARIATION

OF

STRESS

RESULTANTS

DUE

TO

CREEP

AND

SHRINKAGE

! *1o**999*+99**9**+9+*~~~*~**~~~~~**~*~~~~~~~~*~~**~~~***~~~~****~ OPTION BASE 1 REAL C1(4,4).C2(4,4),C3(4,4),C4(4),CS(4),C6(4,2),C7~2) REAL DlC4) ,D2(4) ,03(4) ,D4(4,4),D(4,4)rW(20) REAL A(4,4),B(4),Yn(4),DYn(4),D12(4),C43(4),F(lO),H~ll) REAL Plot(4~100)~Xronos(iOO~~Xr~ll) ! I

*++*99

INPUT

OF

DATA

AND

MATRIX

FORMULATION

+

aOr.41

ire***+

INPUT “ENTER PRINTER CODE Irr(O:Thermikos,l6:Othonh)“,Irr PRINTER IS Irr PRINT “+++9+3++++3++++++++~e****~~*****~e~e**~~***~*~e~*~~*~~*~~~~” PRINT “+ CHANGES OF STRESS RESULTANTS DUE TO CREEP AND SHRINKAGE PRINT “+++9+++++39*e+e++++e*~*******~*~~e~~*~***~*~~**~~*~~e~*” READ Yn(*) ! INITIAL VALUES DATA -4.15, -27.6r1.75, -71.48 PRINT ” INITIAL VALUES “I LIN(l),"-----------------",LINcl) PRINT USING 220:Yn(+) IMAGE 4cSDDD.DDD, 2X) // READ F(+) ! NORMALIZED FUNCTION DATA 3.294,1.687,1.363r1.250~1.230,1.200,1.1,1.05,1.05,1.03 READ H(+l ! TIME FUNCTION DATA .396,.56~.712,.794~.874~.946,1,1.104,1.184,1.212,1.312 READ Xr (*) DATA 30,45~60,75,90>105~120,150,180,210~240 READ Cl(+) DATA -7309.9,2193.,-18274.8,2193. DATA 2193.~ -97820.3,6578.9,-877.2 DATA -18274.8,6578.9,-58473.818771.9 -877.2,-8771.9,-2057.4 DATA 2193.J PRINT USING 350:Cl(*) IMAGE ~~MATRIX ~~**/~~--------**/~(~(sDDDDD.D,~)/)/ READ CZ<+) DATA O,O,O,D DATA O,O,O,O DATA 0,0,7309.9,2193. DATA O,O,-2193.11180.2 PRINT USING 42O:C2(*) “/4C4(SDDDDD.D,X)/)/ IMAGE “MATRIX CZ”/“-------READ C31*) DATA 7309.9,-21.93, 18274.8,-2193 DATA -2193,1180.2,-6578.9~877.2 DATA 18274.8,~6578.9,51169.6,-6578.9 DATA -2193,877.2,-6578.9~877.2 PRINT USING 49O:C3(*) IMAGE “MATRIX C3”/“---------“/4(4(SDDDDD.D, X)/1/ READ C4(*) DATA -219298.2,97820.3,-657894.7a87719.3 PRINT USING 53O:C4(*) IMAGE “MATRIX C4”/“---------“/4((SDDDDD.D,X)f)/ READ CS(+) DATA O,O,-219298.2897820.3 PRINT USING 570:CS(*) IMAGE “MATRIX C5”/“---------“/~((SDDDDDD.D.E)/)/ READ C6(*) DATA +219298.2,0,-97820.3,O DATA +657894.7,-219298.3,-87719.3,-97820.3 PRINT USING 62O:CC(Jo

+”

Variation to20

IMAGE

“pj&)TRI

.‘.7’0

PRINT

I’

,:I

_.

x

(-6”

,“_______--

“/4(2(SDDDDDD.D,X)/)/ Xl(T)

STEP

F2

125

of stress resultants in concrete structures

X2Z[TN)

x

X3CT‘,

tm)” 640 __._@I 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 180 890 900 910 920 930 940 950 960 970 980 990 1000 1010 1020 1030 1040 1050 1060 1070 1080 1090 1100 1110 1120 1130 1140 1150 1160 _-_(I

I

!

++**++

INTEGRATION

UF

THE

SYSTEM

++a+***

Fl=O w=H1’1) Xrono=XrCl) FOR I=1 TO 10 Hstep=(H(I+l)-H(I))/5 Xronstep=~Xr~I+l~-Xro)/5 H=H+Hstep ,,33+,step+.33333 Hlt=Hst,eP*.16666 HS=HsteP*.S H125=HstePi+.125 GOSAE Su br FOR J=l TO 4 WCJ+4)=Yn(J) Y~(J)=WCJ+~)+H~~*DY~(J) WCJ+S)=OynCJ) NEXT J GOSUB Subr FOR J=l TO 4 Yn~J)=W~J+4)+Hl6*(W(J+s)+DYno) NEXT J FOR J=l TO 4 Yn(J)=W(J+4)+H125+(Wo+3+DYn[J)) WCJ+~~)=DY~(J) NEXT J GOSUB Subr FOR J=l TO 4 YnCJ)=WCJ+4l+H5*(W(J+8)-3*W~J+l2~+4aDyn[J~~ W(J+12)=Dyn(J) NEXT J GOSUB Su br FOR J=l TO 4 YnCJ)=W(J+4)+H16*~W~J+8)+Duno+49W(J+12~> NEXT J

I ! ++**++ ,

PRINT

FINAL

RESULTS

*+++*1(

PRINT USING lOLOiF~I)~H,Yn~+) IMAGE D.SDa3X,D.50,3X,4(SDDDDD.DDDD.2X) Pl=P1+1 Xrono=Xrono+XronstePp H=H+Hstep Xronos(Pl)=Xrono FOR J=l TO 4 PlotCJ,Pll=Yn(J) NEXT J IF H<=H(I+l) THEN 750 H=H-Hstep PRINT ~-----------------------~~~-~~~~-~~---__-__-__-____________

--_-_---

G. D.

126 1170 1180 1190 1200 1210 1220 1230 1240 1250 1260 1270 1280 1290 1300 1310 1320 1330 1340 1350 1360 1370 1380 1390 1400 1410 1420 1430 1440 1450 1460 1470 1480 1490 1500 1510 1520 1530

NEXT !

I

9+?+++9

PLOT

RESULTS

INPUT “ENTER PLOTTER PLOTTER IS Iplot$ GRAPHICS SCALE 30,240,-400,100 AXES 30,100,30,0 FOR J=l TO 4 FOR Plo=l TO Pl PLOT Xronos(Plo),Plot(J,Plo) NEXT Plo PENUP LINE TYPE J NEXT J IF IPlotS=“GRAPHICS” PAUSE END !

STEFANOU

**+I++*

COMPUTE

Subr: C7( 1)=1.825 F2=F( I) C7 ( 2) =l . S446+F2 MAT Dl =Cb+C7 NAT D2=(F2)*CS MAT D12=Dl+D2 MAT D3=D12+C4 MAT D4( F2) aC2 MAT D4=D4+C3 MAT D=INV(Cl) MAT A=D+D4 MAT B=D+D3 MAT G=A+Yn MAT Dvn=G+B RETURN

9+99+9 CODE

IPlotS~GRAPHICS:Othonh~9872A:Plotter~”~I~lot~

THEN

DUMP

DIFFERENTIAL

GRAPHICS

EQUATIONS

OF

THE

SYSTEM

WW-+**