Preracun fluidovodnih mrez

Vaja pri predmetu Racunalnisko Podprto Konstruiranje

Univerza v Ljubljani
Fakulteta za strojnistvo
Laboratorij za CAD - LECAD
Visokosolski studij
Predavatelj: izr.prof. Dr. Joze Duhovnik, dipl.ing.
Asistent: Mag. Leon Kos, dipl.ing.
 
 
Izdelal: Tomaz Bester
 

Kazalo vsebine:

Abstract

 In this asigmet two Fortran computer progremes have been made. Program Flow calculates fluid flow through pipeline, while program Prikaz presents network in 3D space and exports data in to dxf fiile. Presententation is made wtih PHIGS graphic library. This program allows us to translate, scale and rotete pipeline round x, y and z axis. Program Flow can calculate flows only in loop networks, and can't solwe tree networks.
 

Uvod

Poznamo tri tipe fluidovodnih mrez:
   - mreze v obliki drevesa
   - mreze v obliki zank (prstanov)
   - kombinacije zgornjih dveh tipov

Vaja predstavlja preracun fluidovodne mreze v obliki zanke, ce imamo dolocene karakteristike cevi in dotoke v posameznih vozliscih. Za prstanaste mreze velja, da je vsota vseh tokov, v posameznem vozliscu 0 in vsota vseh padcev tlakov na zakljuceni zanki ravno tako 0. S pomocjo teh enacb lahko izracunamo pretoke po posameznih ceveh.
Izracun je bil izveden s pomocjo racunalniskega programa v Fortranu, program za prikaz pa je poleg le tega za izris uporabljal se PHIGS knjiznico.
 

Mreze v obliki zank

Osnova izracuna fluidovodne mreze je Bernoulijeva enacba za stacionarni tok nestisljivega fluida za sistem s trenjem:

Ker se prerez cevi ne spreminja je hitrost konstantna tako, da lahko napisemo za tlak na izhodu:

Padec tlaka  je funkcija hitrosti fluida, dolzine cevi, kinematicne viskoznosti, premera cevi in hrapavosti cevi. Za izracun padca tlaka si pomagamo z Reynoldsovim stevilom in enacbami za izracun koeficienta trenja pri razlicnih tokovih, lahko pa si pomagamo z diagrami, ki prikazujejo padec tlaka v odvisnosti od masnega pretoka (slike 1-3.):

                                                                          Slika 1.

                                                        Slika 2.

                                                                        Slika 3.

Iz zgornjih slik je razvidno, da je v logaritemskem grafu padec tlaka na tekoci meter cevi linearno odvisen od masnega toka skozi cev, zato lahko enostavno izpeljemo enacbo za padec tlaka v milibar-ih na tekoci meter cevi.

     (3)

        padec tlaka

          masni pretok

A,B         konstantii

Konstante A in B se izracunajo tako, da iz diagrama odcitamo vrednosti za padec tlaka in masni tok za dve tocki in ju vstavimo v gornji enacbi, tako da dobimo sistem dveh enacb z dvema neznankama, iz katerega enaostavno izracunamo neznanki A in B.
Neglede na tip mreze morata biti izpoknjena naslednja dva zakona:

1. Vsota vseh tokov, ki pritekajo in odtekajo  v vozlisce je enaka 0, prav tako pa mora biti vsota vseh tokov, ko pritekajo in odtekajo v sistem 0.

                 (4)

2. Vsota padcev tlakov po vsaki zakljuceni zanki je 0. Staticni tlak ne vpliva na zakljuceno zanko, ker se tlacne razlike v enem obhodu zanke iznicijo saj zacnemo in koncamo v isti tocki.

            (5)

Po gornjih enacbah se lahko zapise sistem enacb za mrezo tokov. Ker padec tlaka ni linearno odvisen od masnega toka sistem poenostavimo tako da enacbo za padec tlaka lineariziramo, tako da jo razvijemo v Taylorjevo vrsto in upostevamo samo linearni clen.

                       (6)

Ker potrebujemo odvisnost padca tlaka od volumskega pretoka zapisemo:

             (7)

        gostota fluida

l            dolzina cevi
 
Ce upostevamo, da je vsota vseh tokov v vozliscu 0 in vsota vseh padcev tlakov na zakljuceni zanki tudi 0, lahko zapisemo naslednje enacbe:

                                                                (8)

                     (9)

        sprememba pretoka v i-ti cevi

           odtok iz i-tega vozlisca

          pretok skozi i-to cev

Za vsako vozlisce lahko napisemo enacbo 8, za vsako zanko pa enacbo 9. S tem dobimo vec enacb kot je potrebno za resitev sistema, zato pri resevanju ne upostevamo ene vozliscne enacbe, saj je dotok v zadnjem vozliscu ze popolnoma dolocen ce poznamo ostale dotoke, prav tako pa ne upostevamo tistih zancnih encb, ki ne doprinesejo nobenega novega padca tlaka v cevi. Sistem enacb lahko zapisemo tudi v matricni obliki tako, da v prvih vrsticah zapisemo zancne enacbe (8), v zadnjih pa vozliscne (9):

                                   (10)

                                    (11)

         za vozliscne enacbe                 (12)
 

        za zancne enacbe                     (13)

Neznanke v sistemu so , s pomocjo katerih izracunamo nove vrednosti za :

                                                        (14)

Postopek se ponavlja dokler ni napaka manjsa od 0.01:

                                                (15)
 

Graficni prikaz mreze

Knjiznica omogoca risanje in prikaz na izhodnih enotah. Navodilo za risanje s pomocjo PHIGS knjiznice se nahaja na spletni strani http://www.lecad.uni-lj.si/documents/phigsd/lib-users.html
Da lahko s pomocjo PHIGS knjiznice narisemo mrezo in jo poleg tega tudi poljubno transliramo, skaliramo in rotiramo, moramo najprej izracunati transformirane kordinate vozlsc, ki jih dobomo tako da stare koordinate mnozimo z transformacijsko matriko:

Translacija

                                                         (16)

Skalacija

                                                              (17)

Rotacija okoli x osi

               (18)

Rotacija okoli y osi

                 (19)

Rotacija okoli z osi

                  (20)

Ce opravimo vec transformacij lahko dolocimo skupno transformacijsko matriko tako, da mnozimo transformacijaske matrike za posamezne transformacije.
 


Program

Za izracun pretokov po posameznih ceveh mreze in graficni prikaz rezultatov sta bila narejena dva programa v Fortranu. Prvi je flow.f, ki iz datoteke Podatki.txt prebere podatke o mrezi, izracuna pretoke po mrezah in zapise rezultate v datoteko results.bin.
Datoteka Podatki.txt:
 
junctions32 # i, x, y, z, Qi
1   5.0   2.0   0.0  -0.10
2   5.0   6.0   0.0   0.0
3   4.1   6.0   0.0   0.0
4   4.1   6.0   1.0   0.0
5   4.0   6.0   1.0   0.0
6   4.0   6.0   0.0   0.0
7   0.0   6.0   0.0   0.0
8   0.0   3.1   0.0   0.0
9   0.0   3.1   1.0   0.0
10  0.0   3.0   1.0   0.0
11  0.0   3.0   0.0   0.0
12  0.0   0.1   0.0   0.0
13  0.0   0.1   1.0   0.0
14  0.0   0.0   1.0   0.0
15  0.0   0.0   0.0   0.0
16  5.0   0.0   0.0   0.0
17  5.0   1.9   0.0   0.10
18  6.0   0.0   0.0   0.0
19  7.9   0.0   0.0   0.0
20  7.9   0.0   1.0   0.0
21  8.0   0.0   1.0   0.0
22  8.0   0.0   0.0   0.0
23  9.0   0.0   0.0   0.0
24  9.0   1.9   0.0   0.0
25  9.0   1.9   1.0   0.0
26  9.0   2.0   1.0   0.0
27  9.0   2.0   0.0   0.0
28  9.0   4.9   0.0   0.0
29  9.0   4.9   1.0   0.0
30  9.0   5.0   1.0   0.0
31  9.0   5.0   0.0   0.0
32  6.0   5.0   0.0   0.0
pipes32 # i, zaY, kon, l(m), a, b
1   1   2    4.0  0.5663  110.0
2   2   3    1.9  0.5663  110.0
3   3   4    1.0  0.5663  110.0
4   4   5    0.1  0.5663  110.0
5   5   6    1.0  0.5663  110.0
6   6   7    4.0  0.5663  110.0
7   7   8    2.9  0.5663  110.0
8   8   9    1.0  0.5663  110.0
9   9   10   0.1  0.5663  110.0
10  10  11   1.0  0.5663  110.0
11  11  12   2.9  0.5663  110.0
12  12  13   1.0  0.5663  110.0
13  13  14   0.1  0.5663  110.0
14  14  15   1.0  0.5663  110.0
15  15  16   5.0  0.5663  110.0
16  16  17   1.9  0.5663  110.0
17  18  17   2.1  0.5663  110.0
18  19  18   2.0  0.5663  110.0
19  20  19   1.0  0.5663  110.0
20  21  20   0.1  0.5663  110.0
21  22  21   1.0  0.5663  110.0
22  23  22   1.0  0.5663  110.0
23  24  23   1.9  0.5663  110.0
24  25  24   1.0  0.5663  110.0
25  26  25   0.1  0.5663  110.0
26  27  26   1.0  0.5663  110.0
27  28  27   2.9  0.5663  110.0
28  29  28   1.0  0.5663  110.0
29  30  29   0.1  0.5663  110.0
30  31  30   1.0  0.5663  110.0
31  32  31   3.0  0.5663  110.0
32   1  32   3.2  0.5663  110.0
loops 1 # çt. cevi {smer cevi ...}
32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 -17 -18 -19 -20 -21 -22 -23 -24 -25 -26 -27 -28 -29 -30 -31 -32
constants 1000 1.01e-6  # gostota, viskoznost
solve

Na prvi vrstici datoke je zapisano stevili vozlisc, v naslednjih vsticah pa koordinate vozlisc in odtok iz vozlisc. Za vozlisci je zapisano stevilo cevi, v naslednjih vrsticah pa zacetno in koncno vozlisce cevi, dolzina cevi in koeficienta A in B. Za cevmi je zapisano stevilo zank, v naslednjih vrsticah pa najprej stevilo cevi v zanki nato pa zaporedne stevilke cevi v zanki, z ustreznim predznakom. V zadnji vrstici je zapisana gostota fluida.
Program flow glede na podatke prebrane iz vhodne datoteke naredi matrike iz enacbe 10, ki jo nato resi s pomocjo podprogramov ludcmp in lubksb. Iterativni postopek se ponavlja dokler napaka izracunana po enacbi 15 ni manjsa od 0.01. Na koncu program zapise podatke o koordinatah vozlisc ter odtokih iz njih in podatke o zacetnih in koncnih vozliscih cevi ter pretokih po ceveh.
Za graficni prikaz mreze je bil narejen program Prikaz. Program prebere podatke o mrezi iz datoteke Results.bin in nam s pomocjo graficne knjiznice PHIGS izrise mrezo z oddtoki iz vozlisc, ter pertoki po ceveh. Pred izrisom mreze imamo moznost izbirati med naslednjimi geometrijskimi trensformacijami:
- Translacija
- Skalacija
- Rotacija okoli x osi
- Rotacija okoli y osi
- Rotacija okoli z osi
Geometrijske transformacije so izvedene po enecbah (16 - 20). Poleg geometrijskih transformacij nam program Prikaz omogoca tudi izvoz geometrijsjkih podatkov v DXF format, saj po zelji kreira datoteko results.dxf.
Primer:
Mreza na slliki 4. ima 32 vozlisc in cevi ter eno samo zanko, vhodni podatki za to mrezo so podani v datoteki podatki.txt


                                                                          Slika 4.

Za primer na sliki lahko pretoke izracunane s programom flow preverimo, tako da po enacbii 3 izracunamo padce tlakov po obeh polovicah zanke. Ti padci tlakov morajo biti enakii.

l1 =28.9m                            dolzina prve polovice zanke
l2 =22.4m                            dolzina druge polovice zanke
Q1 =0,4641m3 /s                pretok po prvi polovici zanke
Q2 =0,5359m3 /s                pretok po drugi polovici zanke

 

V prvem primeru sta koeficienta A in B enaka za vse cevi, ce pa program uporabimo za iizracun pretokov v cevovodu, kjer so poleg cevi tudi elementi z lokalnimi izgobami, moralmo v cevovodu lokalne izgube ponazoriti s cevjo dolzine 1m in ustreznima koeficientoma A in B. Podatki iz datoteke podatki2.txt je nam prikazujejo tak primer.

junctions32 # i, x, y, z, Qi
1   5.0   2.0   0.0  -0.10
2   5.0   6.0   0.0   0.0
3   4.1   6.0   0.0   0.0
4   4.1   6.0   1.0   0.0
5   4.0   6.0   1.0   0.0
6   4.0   6.0   0.0   0.0
7   0.0   6.0   0.0   0.0
8   0.0   3.1   0.0   0.0
9   0.0   3.1   1.0   0.0
10  0.0   3.0   1.0   0.0
11  0.0   3.0   0.0   0.0
12  0.0   0.1   0.0   0.0
13  0.0   0.1   1.0   0.0
14  0.0   0.0   1.0   0.0
15  0.0   0.0   0.0   0.0
16  5.0   0.0   0.0   0.0
17  5.0   1.9   0.0   0.10
18  6.0   0.0   0.0   0.0
19  7.9   0.0   0.0   0.0
20  7.9   0.0   1.0   0.0
21  8.0   0.0   1.0   0.0
22  8.0   0.0   0.0   0.0
23  9.0   0.0   0.0   0.0
24  9.0   1.9   0.0   0.0
25  9.0   1.9   1.0   0.0
26  9.0   2.0   1.0   0.0
27  9.0   2.0   0.0   0.0
28  9.0   4.9   0.0   0.0
29  9.0   4.9   1.0   0.0
30  9.0   5.0   1.0   0.0
31  9.0   5.0   0.0   0.0
32  6.0   5.0   0.0   0.0
pipes32 # i, zaY, kon, l(m), a, b
1   1   2    4.0  0.5663  110.0
2   2   3    1.9  0.5663  110.0
3   3   4    1.0  0.5663  110.0
4   4   5    0.1  0.5663  110.0
5   5   6    1.0  0.5626  0.0269
6   6   7    4.0  0.5663  110.0
7   7   8    2.9  0.5663  110.0
8   8   9    1.0  0.5663  110.0
9   9   10   0.1  0.5663  110.0
10  10  11   1.0  0.5626  0.0269
11  11  12   2.9  0.5663  110.0
12  12  13   1.0  0.5663  110.0
13  13  14   0.1  0.5626  0.0296
14  14  15   1.0  0.5663  110.0
15  15  16   5.0  0.5663  110.0
16  16  17   1.9  0.5663  110.0
17  18  17   2.1  0.5663  110.0
18  19  18   2.0  0.5663  110.0
19  20  19   1.0  0.5663  110.0
20  21  20   0.1  0.5663  110.0
21  22  21   1.0  0.5626  0.0296
22  23  22   1.0  0.5663  110.0
23  24  23   1.9  0.5663  110.0
24  25  24   1.0  0.5663  110.0
25  26  25   0.1  0.5663  110.0
26  27  26   1.0  0.5626  0.0296
27  28  27   2.9  0.5663  110.0
28  29  28   1.0  0.5663  110.0
29  30  29   0.1  0.5663  110.0
30  31  30   1.0  0.5626  0.0296
31  32  31   3.0  0.5663  110.0
32   1  32   3.2  0.5663  110.0
loops 1 # çt. cevi {smer cevi ...}
32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 -17 -18 -19 -20 -21 -22 -23 -24 -25 -26 -27 -28 -29 -30 -31 -32
constants 1000 1.01e-6  # gostota, viskoznost
solve

Pretoki izracunani s programom Flow za podatke iz datoteke podatki2.txt

Q1 =0,4728m3 /s                pretok po prvi polovici zanke
Q2 =0,5272m3 /s                pretok po drugi polovici zanke

Program Flow

      PROGRAM FLOW
 
      parameter (NP=100)
      character*80 line
      real a(NP,NP),dQ(NP),d
      integer indx(NP)
      real lambda,lambda1,PI
      real Q(NP),l(NP),visk,gost,x(NP),y(NP),z(NP)
      real Qd(NP),C(NP),aa(NP),bb(NP)
      integer i,j,pipes,loops,junctions,stcevi(NP),zac(NP),kon(NP)
      integer sm(NP,NP),vrst
      PI=3.141592654
      open(1,file="podatki.txt")
      open(2,file="results.bin",status="unknown",
     *access="direct",form="UNFORMATTED",
     *recl=58)
      open(3,file="results.txt",status="unknown",
     *form="FORMATTED")
10    read(1,"(a)",end=50) line
      if (line(1:9).eq."junctions")then
        read(line(10:12),"(I10)",err=103)junctions
        Qd0=0
          do 30 i=1,junctions
          read(1,*) j, x(j), y(j), z(j), Q(j)
          Qd0=Qd0+abs(Q(j))
 
30    continue
      endif
 
      if (line(1:5).eq."pipes")then
        read(line(6:8),"(I10)",err=102)pipes
         write(*,*) "Pipes=",pipes
        do 20 i=1,pipes
        read(1,*) j, zac(j), kon(j), l(j), aa(j), bb(j)
          if(zac(j).lt.junctions) sm(zac(j),j)=1
          if(kon(j).lt.junctions) sm(kon(j),j)=-1
20      continue
      endif
 
      if (line(1:5).eq."loops")then
        read (line(6:8),"(I10)",err=104)loops
         vrst=junctions-1
           do 40 ii=1,loops
             read (1,*) j, (stcevi(m), m=1,j)
               do 40 i=1,j
                 if (stcevi(i).gt.0)then
                   sm(vrst+ii,stcevi(i))=1
                 else
                   sm(vrst+ii,-stcevi(i))=-1
                 endif
40         continue
      endif
 
      if (line(1:9).eq."constants")then
        read(line(11:15),"(f4)",err=105)gost
        read(line(16:22),"(e10.5)",err=106)visk
      endif
      goto10
 
50    close(1)

      Qd0=Qd0/pipes
      print*,"Qd0=",Qd0
        do 90 i=1,pipes
90      Qd(i)=Qd0
c      Print*,"Smerna matrika SM"
c      do 60 i=1,pipes
c        do 70 j=1,pipes
c70        print "(I3,$)",sm(i,j)
c60    print *, " "
c      pause
200   print*,"MATRIKA A"
      do 80 i=1,junctions-1
      do 80 j=1,pipes
80    a(i,j)=sm(i,j)
      do 120 i=junctions,pipes
        do 120 j=1,pipes
      a(i,j)=(1/aa(j))*(((abs(Qd(j))*gost)/bb(j))**(1/aa(j)-1))*
     a(1/bb(j))*sm(i,j)*gost*l(j)
      if (Qd(j).lt.0) a(i,j)=-a(i,j)
120   continue
      do 150 i=1,junctions-1
        dQ(i)=-Q(i)
        do 150 j=1,pipes
          dQ(i)=dQ(i)-sm(i,j)*Qd(j)
150   continue
      do 175 i=junctions,pipes
        dQ(i)=0
          do 170 j=1,pipes
          if (Qd(j).ge.0) then
            dQ(i)=dQ(i)-(((abs(Qd(j))*gost)/bb(j))**(1/aa(j)))*sm(i,j)
     s*l(j)
          else
            dQ(i)=dQ(i)+(((abs(Qd(j))*gost)/bb(j))**(1/aa(j)))*sm(i,j)
     s*l(j)
      endif
170   continue
175   continue
      do 180 i=1,pipes
180     write(3,*),"dQ(",i,")=",dQ(i)
      n=pipes
      call ludcmp(a,n,NP,indx,d)
      call lubksb(a,n,NP,indx,dQ)
      do 13 i=1,n
13      if (abs(dQ(i)/Qd(i)).gt.0.01) goto 210
      goto 300
210   do 212 i=1,n
        Qd(i)=Qd(i)+dQ(i)
212   continue
      goto 200

102   stop
103   stop
104   stop
105   stop
106   stop
300   print*," "
        do i=1,n
          print*,"Qd(i)=",Qd(i)
        enddo
      write(2,rec=1)junctions
        do i=1,junctions
          write(2,rec=i+1)i,x(i),y(i),z(i),Q(i)
        enddo
      write(2,rec=junctions+2)pipes
        do j=1,pipes,10
          write (2,rec=junctions+j+2)j,zac(j),kon(j),Qd(j)
        enddo
      read(2,rec=1)junctions
      read(2,rec=junctions+2)pipes
        do i=1,pipes
          read(2,rec=junctions+i+2)j,zac(j),kon(j),Qd(j)
        enddo
        close (unit=3)
        close (unit=2)
      end
 
 
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      PARAMETER (NMAX=100,TINY=1.0E-20)
      DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11      CONTINUE
        IF (AAMAX.EQ.0) PAUSE "SINGULAR MATRIX"
        VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,N
        IF (J.GT.1)THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF(I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
            A(I,J)=SUM
            ENDIF
14        CONTINUE
        ENDIF
        AAMAX=0
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
          A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX)THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
16      CONTINUE
      IF (J.NE.IMAX)THEN
        DO 17 K=1,N
          DUM=A(IMAX,K)
          A(IMAX,K)=A(J,K)
          A(J,K)=DUM
17      CONTINUE
        D=-D
        VV(IMAX)=VV(J)
      ENDIF
      INDX(J)=IMAX
C        IF (J.NE.N)THEN
          IF (A(J,J).EQ.0)A(J,J)=TINY
          IF (J.NE.N)THEN
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
18        CONTINUE
        ENDIF
19    CONTINUE
C      IF(A(N,N).EQ.0)A(N,N)=TINY
      RETURN
      END
 
 
 
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
      DIMENSION A(np,np),INDX(N),B(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0) THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
11        CONTINUE
        ELSE IF (SUM.NE.0) THEN
          II=I
        ENDIF
        B(I)=SUM
12    CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF (I.LT.N) THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
13        CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END

Program Prikaz

      Program prikaz
      include 'phigsdef.f'
      parameter (NP=100)
      integer zac(NP),kon(NP),i,j,junctions,pipes,o
      real u(3) /0.0,0.0,100.0/
      real v(3) /100.0,0.0,0.0/
      real xxx(NP),yyy(NP),zzz(NP)
      real xx(NP),yy(NP),zz(NP),Q(NP),Qd(NP),Tx,Ty,Tz
      real x(2,NP),y(2,NP),z(2,NP),PI,Sx,Sy,Sz,Rx,Ry,Rz
      real  ViewMappingMatrix(3,3),xp(2),yp(2),xt(2,NP)
      real yt(2,NP),zt(2,NP),xpr,ypr,xxt(NP),yyt(NP),zzt(NP)
      real  WindowLimits(4) /-300.0, 400.0, -300.0, 400.0/
      real  ViewportLimits(4) /0.0, 1.0, 0.0, 1.0/
      real  ClipLimits(4)     /0.0, 1.0, 0.0, 1.0/
      integer wkid,ErrorReturn
      character*5 pretok
      character*2 voz
      character*5 dotok
      character*2 cev
      PI=3.1415927
      open(1,file="results.bin",status="unknown",
     *access="direct",form="UNFORMATTED",
     *recl=58)
      open(2,file="results.dxf",status="unknown",
     *access="sequential",form="FORMATTED")
      read(1,rec=1)junctions
      print*,"junctions=",junctions
      do i=1,junctions
        read(1,rec=i+1)j,xx(i),yy(i),zz(i),Q(j)
c        print*,j,  xx(i),  yy(i),  zz(i),  Q(i)
      enddo
      read(1,rec=junctions+2)pipes
      do i=1,pipes
        read(1,rec=junctions+i+2)j,zac(j),kon(j),Qd(j)
c        print*,j,zac(j),kon(j),Qd(j)
      enddo
      do i=1,pipes
        x(1,i)=xx(zac(i))
        y(1,i)=yy(zac(i))
        z(1,i)=zz(zac(i))
        x(2,i)=xx(kon(i))
        y(2,i)=yy(kon(i))
        z(2,i)=zz(kon(i))
        xt(1,i)=x(1,i)
        yt(1,i)=y(1,i)
        zt(1,i)=z(1,i)
        xt(2,i)=x(2,i)
        yt(2,i)=y(2,i)
        zt(2,i)=z(2,i)
        xxt(i)=xx(i)
        yyt(i)=yy(i)
        zzt(i)=zz(i)
        xxx(i)=xx(i)
        yyy(i)=yy(i)
        zzz(i)=zz(i)
      enddo
      sx=1
      sy=1
      sz=1
      rx=0
      ry=0
      rz=0
      tx=0
      ty=0
      tz=0
10    continue
      print*,"          *******************************************"
      print*,"          *                                         *"
      print*,"          *          1. TRANSLACIJA                 *"
      print*,"          *                                         *"
      print*,"          *          2. SKALACIJA                   *"
      print*,"          *                                         *"
      print*,"          *          3. ROTACIJA V X OSI            *"
      print*,"          *                                         *"
      print*,"          *          4. ROTACIJA V Y OSI            *"
      print*,"          *                                         *"
      print*,"          *          5. ROTACIJA V Z OSI            *"
      print*,"          *                                         *"
      print*,"          *          6. IZRIS                       *"
      print*,"          *                                         *"
      print*,"          *          7. IZVOZ V DXF DATOTEKO        *"
      print*,"          *                                         *"
      print*,"          *          8. KONEC                       *"
      print*,"          *                                         *"
      print*,"          *******************************************"
      read(*,*),o

      if (o.eq.1)then
        print*,"Translacija po x osi"
        read(*,*),Tx
        print*,"Translacija po y osi"
        read(*,*),Ty
        print*,"Translacija po z osi"
        read(*,*),Tz
          do i=1,pipes
            xt(1,i)=x(1,i)+Tx
            yt(1,i)=y(1,i)+Ty
            zt(1,i)=z(1,i)+Tz
            xt(2,i)=x(2,i)+Tx
            yt(2,i)=y(2,i)+Ty
            zt(2,i)=z(2,i)+Tz
          enddo
          do i=1,junctions
            xxt(i)=xx(i)+Tx
            yyt(i)=yy(i)+Ty
            zzt(i)=zz(i)+Tz
          enddo
      endif
      if (o.eq.2)then
        print*,"Skalacija po x osi"
        read(*,*),Sx
        print*,"Skalacija po y osi"
        read(*,*),Sy
        print*,"Skalacija po z osi"
        read(*,*),Sz
          do i=1,pipes
            xt(1,i)=x(1,i)*Sx
            yt(1,i)=y(1,i)*Sy
            zt(1,i)=z(1,i)*Sz
            xt(2,i)=x(2,i)*Sx
            yt(2,i)=y(2,i)*Sy
            zt(2,i)=z(2,i)*sz
          enddo
          do i=1,junctions
            xxt(i)=xx(i)*Sx
            yyt(i)=yy(i)*Sy
            zzt(i)=zz(i)*sz
          enddo
      endif
      if (o.eq.3)then
        print*,"Rotacija v x osi"
        read(*,*),Rx
        Rx=Rx*PI/180
          do i=1,pipes
            xt(1,i)=x(1,i)
            yt(1,i)=y(1,i)*cos(Rx)-z(1,i)*sin(Rx)
            zt(1,i)=y(1,i)*sin(Rx)+z(1,i)*cos(RX)
            xt(2,i)=x(2,i)
            yt(2,i)=y(2,i)*cos(Rx)-z(2,i)*sin(Rx)
            zt(2,i)=y(2,i)*sin(Rx)+z(2,i)*cos(RX)
          enddo
          do i=1,junctions
            xxt(i)=xx(i)
            yyt(i)=yy(i)*cos(Rx)-zz(i)*sin(Rx)
            zzt(i)=yy(i)*sin(Rx)+zz(i)*cos(RX)
          enddo
      endif
      if (o.eq.4)then
        print*,"Rotacija v y osi"
        read(*,*),Ry
        Ry=Ry*PI/180
          do i=1,pipes
            xt(1,i)=x(1,i)*cos(Ry)-z(1,i)*sin(Ry)
            yt(1,i)=y(1,i)
            zt(1,i)=x(1,i)*sin(Ry)+z(1,i)*cos(Ry)
            xt(2,i)=x(2,i)*cos(Ry)-z(2,i)*sin(Ry)
            yt(2,i)=y(2,i)
            zt(2,i)=x(2,i)*sin(Ry)+z(2,i)*cos(Ry)
          enddo
          do i=1,junctions
            xxt(i)=xx(i)*cos(Ry)-zz(i)*sin(Ry)
            yyt(i)=yy(i)
            zzt(i)=xx(i)*sin(Ry)+zz(i)*cos(Ry)
          enddo
      endif
      if (o.eq.5)then
        print*,"Rotacija v z osi"
        read(*,*),Rz
        Rz=Rz*PI/180
          do i=1,pipes
            xt(1,i)=x(1,i)*cos(Rz)-y(1,i)*sin(Rz)
            yt(1,i)=x(1,i)*sin(Rz)+y(1,i)*cos(Rz)
            zt(1,i)=z(1,i)
            xt(2,i)=x(2,i)*cos(Rz)-y(2,i)*sin(Rz)
            yt(2,i)=x(2,i)*sin(Rz)+y(2,i)*cos(Rz)
            zt(2,i)=z(2,i)
          enddo
          do i=1,junctions
            xxt(i)=xx(i)*cos(Rz)-yy(i)*sin(Rz)
            yyt(i)=xx(i)*sin(Rz)+yy(i)*cos(Rz)
            zzt(i)=zz(i)
          enddo
      endif
      if (o.eq.6)then
      call popph(1,0)
      wkid=1
      call popwk(wkid,"",WK151024)
      call pevmm(WindowLimits, ViewportLimits,
     *ErrorReturn, ViewMappingMatrix)
      do 50 i=1,3
        do 60 j=1,3
          write(*,*) ViewMappingMatrix(i,j)
60      continue
50    continue
      call pswkw(wkid, 0.0, 1.0, 0.0, 1.0)
      call pswkv(wkid, 0.0, 0.15, 0.0, 0.15)
      call psplci(4)

      do i=1,pipes
        xp(1)=xt(1,i)
        yp(1)=yt(1,i)
        xp(2)=xt(2,i)
        yp(2)=yt(2,i)

      call ppl(2,xp,yp)
      enddo
      call pstxci(2)
      call pschh (12)
      call pstxal(paleft,0)
      do i=1,pipes
        xpr=(xt(1,i)+xt(2,i))/2
        ypr=(yt(1,i)+yt(2,i))/2
        write(pretok,"(f5.3)")Qd(i)
        call ptx(xpr,ypr,pretok)
        xpr=(xt(1,i)+xt(2,i))/2
        ypr=(yt(1,i)+yt(2,i))/2+18
        write(cev,"(i2)")i
        call ptx(xpr,ypr,cev)
      enddo
      call pstxci(4)
      do i=1,junctions
        xpr=(xxt(i))
        ypr=(yyt(i)+18)
        write(voz,"(I2)")i
        call ptx(xpr,ypr,voz)
      enddo
      call pstxci(6)
      do i=1,junctions
        xpr=(xxt(i))
        ypr=(yyt(i))
        write(dotok,"(f5.3)")Q(i)
        call ptx(xpr,ypr,dotok)
      enddo
      call psplci(3)
      call ppl(3,u,v)
      call pstxci(3)
      call ptx(100.0,0.0,"X")
      call ptx(0.0,100.0,"Y")
      pause
      call pclwk(wkid)
      call pclph()
      endif
1     format (I1)
2     format (I2)
3     format (a3)
4     format (a4)
5     format (a6)
6     format (a7)
7     format (a8)
8     format (f8.3)
9     format (a5)
11    format (a1)
12    format (a2)
      if (o.eq.7)then
        write(2,1)0
        write(2,6)"SECTION"
        write(2,1)2
        write(2,7)"ENTITIES"
        do i=1,pipes
          write(2,1)  0
          write(2,4)"LINE"
          write(2,1)  8
          write(2,1)0
          write(2,2) 10
          write(2,8)xxx(zac(i))
          write(2,2) 20
          write(2,8)yyy(zac(i))
          write(2,2) 30
          write(2,8)zzz(zac(i))
          write(2,2) 11
          write(2,8)xxx(kon(i))
          write(2,2) 21
          write(2,8)yyy(kon(i))
          write(2,2) 31
          write(2,8)zzz(kon(i))
        enddo
        do i=1,junctions
          write(2,1)0
          write(2,4)"TEXT"
          write(2,1)  8
          write(2,1)0
          write(2,2) 10
          write(2,8)xxx(i)
          write(2,2) 20
          write(2,8)yyy(i)+1
          write(2,2) 30
          write(2,8)zzz(i)
          write(2,2)40
          write(2,1)5
          write(2,1)1
          if (i.lt.10) then
            write(voz,"(i1)")i
            write(2,11)voz
          elseif (i.lt.100.and.i.gt.9) then
            write(voz,"(i2)")i
            write(2,12)voz
          endif
          write(2,1)0
          write(2,4)"TEXT"
          write(2,1)  8
          write(2,1)0
          write(2,2) 10
          write(2,8)xxx(i)
          write(2,2) 20
          write(2,8)yyy(i)-5
          write(2,2) 30
          write(2,8)zzz(i)
          write(2,2)40
          write(2,1)5
          write(2,1)1
          if ((q(i).lt.10).and.(q(i).ge.0))then
            write(dotok,"(f4.3)")q(i)
            write(2,4)dotok
          elseif ((q(i).lt.0).or.(q(i).ge.10)) then
            write(dotok,"(f5.3)")q(i)
            write(2,9)dotok
          endif
        enddo
        do i=1,pipes
          write(2,1)0
          write(2,4)"TEXT"
          write(2,1)  8
          write(2,1)0
          write(2,2) 10
          write(2,8)(xxx(zac(i))+xxx(kon(i)))/2
          write(2,2) 20
          write(2,8)(yyy(zac(i))+yyy(kon(i)))/2-6
          write(2,2) 30
          write(2,8)(zzz(zac(i))+zzz(kon(i)))/2
          write(2,2)40
          write(2,1)5
          write(2,1)1
          if (qd(i).lt.10.and.qd(i).ge.0)then
            write(pretok,"(f4.3)")qd(i)
            write(2,4)pretok
          elseif (qd(i).lt.0.or.qd(i).ge.10) then
            write(pretok,"(f5.3)")qd(i)
            write(2,9)pretok
          endif
          write(2,1)0
          write(2,4)"TEXT"
          write(2,1)  8
          write(2,1)0
          write(2,2) 10
          write(2,8)(xxx(zac(i))+xxx(kon(i)))/2
          write(2,2) 20
          write(2,8)(yyy(zac(i))+yyy(kon(i)))/2+1
          write(2,2) 30
          write(2,8)(zzz(zac(i))+zzz(kon(i)))/2
          write(2,2)40
          write(2,1)5
          write(2,1)1
          if (i.lt.10) then
            write(cev,"(i1)")i
            write(2,11)cev
          elseif (i.lt.100.and.i.gt.9) then
            write(cev,"(i2)")i
            write(2,12)cev
          endif
        enddo
        write(2,1)  0
        write(2,5)"ENDSEC"
        write(2,1)  0
        write(2,3)"EOF"
      endif
      if (o.eq.8) stop
      if (o.lt.1.or.o.gt.8)then
      print*,"NapaYen vnos. Poskusi znova."
      endif
      do i=1,pipes
        x(1,i)=xt(1,i)
        y(1,i)=yt(1,i)
        z(1,i)=zt(1,i)
        x(2,i)=xt(2,i)
        y(2,i)=yt(2,i)
        z(2,i)=zt(2,i)
      enddo
      do i=1,junctions
        xx(i)=xxt(i)
        yy(i)=yyt(i)
        zz(i)=zzt(i)
      enddo
      goto 10
      close (1)
      CLOSE (2)
      end
 
 

 

NASLOV:
Tomaz Bester
Zoisova 44
4000 Kranj
Slovenija
+386 64 225-847