← Back to branch summary

~gabriel1984sibiu/calculix/ccx

« back to all changes in this revision

Viewing changes to CalculiX/ccx_2.11/src/gaspipe_fanno.f

  • Committer: Grevutiu Gabriel
  • Date: 2016-12-30 12:06:41 UTC
  • Revision ID: gabriel1984sibiu@gmail.com-20161230120641-kzmhfy8mn00w3mhg
New upstream version

Show diffs side-by-side

added added

removed removed

Lines of Context:
CalculiX/ccx_2.11/src/gaspipe_fanno.f
 
1
!
 
2
!     CalculiX - A 3-dimensional finite element program
 
3
!     Copyright (C) 1998-2015 Guido Dhondt
 
4
!     
 
5
!     This program is free software; you can redistribute it and/or
 
6
!     modify it under the terms of the GNU General Public License as
 
7
!     published by the Free Software Foundation(version 2);
 
8
!     
 
9
!     
 
10
!     This program is distributed in the hope that it will be useful,
 
11
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
 
12
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 
13
!     GNU General Public License for more details.
 
14
!     
 
15
!     You should have received a copy of the GNU General Public License
 
16
!     along with this program; if not, write to the Free Software
 
17
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
18
!     
 
19
      subroutine gaspipe_fanno(node1,node2,nodem,nelem,lakon,kon,
 
20
     &        ipkon,nactdog,identity,ielprop,prop,iflag,v,xflow,f,
 
21
     &        nodef,idirf,df,cp,r,physcon,dvi,numf,set,
 
22
     &        shcon,nshcon,rhcon,nrhcon,ntmat_,co,vold,mi,ttime,time,
 
23
     &        iaxial)
 
24
!     
 
25
!     pipe with friction losses (Fanno Formulas) GAPF 
 
26
!
 
27
!     author: Yannick Muller
 
28
!     
 
29
      implicit none
 
30
!     
 
31
      logical identity,crit
 
32
      character*8 lakon(*)
 
33
      character*81 set(*)
 
34
!     
 
35
      integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
 
36
     &     ielprop(*),nodef(5),idirf(5),index,iflag,
 
37
     &     inv,ipkon(*),kon(*),icase,kgas,k_oil
 
38
     &     ,nshcon(*),nrhcon(*),ntmat_,i,mi(*),nodea,nodeb,
 
39
     &     nodec,iaxial
 
40
!
 
41
      real*8 prop(*),v(0:mi(2),*),xflow,f,df(5),kappa,R,a,d,l,
 
42
     &     p1,p2,T1,T2,Tt1,Tt2,pt1,pt2,cp,physcon(3),p2p1,km1,dvi,
 
43
     &     kp1,kdkm1,reynolds,pi,e,lambda,lld,kdkp1,T2dTt2,
 
44
     &     T1dTt1,X_t1dTt1,X_t2dTt2,X2_den,X1_den,
 
45
     &     X1,X2,B1,B2,C1,C2,tdkp1,ttime,time,
 
46
     &     pt2zpt1,ks,form_fact,xflow_oil,Tt1dT1,Tt2dT2,
 
47
     &     Pt2zPt1_c,qred_crit,l_neg,Qred,
 
48
     &     expon1,expon2,cte,term1,term2,term3,term4,term5,term6,
 
49
     &     term,phi,M1,M2,qred2,qred1,qred_max1,qred_crit_out,co(3,*),
 
50
     &     shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),vold(0:mi(2),*),
 
51
     &     radius,initial_radius,l_initial
 
52
!
 
53
      if (iflag.eq.0) then
 
54
         identity=.true.
 
55
!     
 
56
         if(nactdog(2,node1).ne.0)then
 
57
            identity=.false.
 
58
         elseif(nactdog(2,node2).ne.0)then
 
59
            identity=.false.
 
60
         elseif(nactdog(1,nodem).ne.0)then
 
61
            identity=.false.
 
62
         endif
 
63
!     
 
64
      elseif (iflag.eq.1)then
 
65
!
 
66
         crit=.false.
 
67
!
 
68
         pi=4.d0*datan(1.d0)
 
69
!     
 
70
         index=ielprop(nelem)
 
71
         kappa=(cp/(cp-R))
 
72
         A=prop(index+1)
 
73
         d=prop(index+2)
 
74
         l=prop(index+3)
 
75
         if(l.lt.0d0) then
 
76
            l_neg=l
 
77
            l=abs(l)
 
78
         else
 
79
            l_neg=l
 
80
         endif
 
81
         ks=prop(index+4)
 
82
         if(lakon(nelem)(2:6).eq.'GAPFA') then
 
83
            icase=0
 
84
         elseif(lakon(nelem)(2:6).eq.'GAPFI') then
 
85
            icase=1
 
86
         endif
 
87
         form_fact=prop(index+5)
 
88
         xflow_oil=prop(index+6)
 
89
         k_oil=nint(prop(index+7))
 
90
!
 
91
         if((lakon(nelem)(2:6).eq.'GAPFF').and.
 
92
     &        (lakon(nelem)(2:7).ne.'GAPFF2')) then
 
93
!            
 
94
            icase=0
 
95
            nodea=nint(prop(index+1))
 
96
            nodeb=nint(prop(index+2))
 
97
c            iaxial=nint(prop(index+3))
 
98
            radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
 
99
     &           co(1,nodea)-vold(1,nodea))**2)
 
100
!
 
101
            initial_radius=dsqrt((co(1,nodeb)-co(1,nodea))**2)
 
102
!
 
103
c            if(iaxial.ne.0) then
 
104
c              A=pi*radius**2/iaxial
 
105
c            else
 
106
              A=pi*radius**2
 
107
c            endif
 
108
            d=2*radius
 
109
            l=prop(index+4)
 
110
            if(l.lt.0d0) then
 
111
               l_neg=l
 
112
               l=abs(l)
 
113
            else
 
114
               l_neg=l
 
115
            endif
 
116
            ks=prop(index+5)
 
117
            form_fact=prop(index+6)
 
118
            xflow_oil=prop(index+7)
 
119
            k_oil=nint(prop(index+8))
 
120
!
 
121
         elseif (lakon(nelem)(2:7).eq.'GAPFF2') then
 
122
            write(*,*) nelem,lakon(nelem)(1:6)
 
123
            icase=0
 
124
            nodea=nint(prop(index+1))
 
125
            nodeb=nint(prop(index+2))
 
126
            nodec=nint(prop(index+3))
 
127
c            iaxial=nint(prop(index+4))
 
128
            radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
 
129
     &           co(1,nodea)-vold(1,nodea))**2)
 
130
            initial_radius=dsqrt((co(1,nodeb)-co(1,nodea))**2)
 
131
            d=2*radius
 
132
c           if(iaxial.ne.0) then
 
133
c              A=pi*radius**2/iaxial
 
134
c            else
 
135
               A=pi*radius**2
 
136
c            endif
 
137
            l_initial=dsqrt((co(2,nodec)-co(2,nodeb))**2)
 
138
            l=dsqrt((co(2,nodec)+vold(2,nodec)-
 
139
     &           co(2,nodeb)-vold(2,nodeb))**2)
 
140
            if(l.lt.0d0) then
 
141
               l_neg=l
 
142
               l=abs(l)
 
143
            else
 
144
               l_neg=l
 
145
            endif
 
146
            ks=prop(index+5)
 
147
            form_fact=prop(index+6)           
 
148
            xflow_oil=prop(index+7)
 
149
            k_oil=nint(prop(index+8))
 
150
         endif
 
151
!
 
152
         pt1=v(2,node1)
 
153
         pt2=v(2,node2)
 
154
!
 
155
         if(pt1.ge.pt2) then
 
156
            inv=1
 
157
            Tt1=v(0,node1)+physcon(1)
 
158
            Tt2=v(0,node2)+physcon(1)
 
159
         else
 
160
            inv=-1
 
161
            pt1=v(2,node2)
 
162
            pt2=v(2,node1)
 
163
            Tt1=v(0,node2)+physcon(1)
 
164
            Tt2=v(0,node1)+physcon(1)
 
165
         endif
 
166
!
 
167
         p2p1=pt2/pt1
 
168
         km1=kappa-1.d0
 
169
         kp1=kappa+1.d0
 
170
         kdkm1=kappa/km1
 
171
         tdkp1=2.d0/kp1
 
172
         C2=tdkp1**kdkm1
 
173
!        
 
174
!     incompressible flow
 
175
         xflow=inv*A*dsqrt(d/l*2*Pt1/(R*Tt1)*(pt1-pt2))
 
176
         if(p2p1.gt.C2) then
 
177
           xflow=inv*pt1*a*dsqrt(2.d0*kdkm1*p2p1**(2.d0/kappa)
 
178
     &           *(1.d0-p2p1**(1.d0/kdkm1))/r)/dsqrt(Tt1)
 
179
         else
 
180
            xflow=inv*pt1*a*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
 
181
     &           dsqrt(Tt1)
 
182
         endif
 
183
!
 
184
!     calculation of the dynamic viscosity 
 
185
!     
 
186
         if(dabs(dvi).lt.1E-30) then
 
187
            kgas=0
 
188
            call dynamic_viscosity(kgas,Tt1,dvi)
 
189
         endif  
 
190
!
 
191
         reynolds=dabs(xflow)*d/(dvi*a)
 
192
!
 
193
         call friction_coefficient(l_neg,d,ks,reynolds,form_fact,lambda)
 
194
         xflow=inv*A*dsqrt(d/(lambda*l)*2*Pt1/(R*Tt1)*(pt1-pt2))
 
195
!
 
196
         call pt2zpt1_crit(pt2,pt1,Tt1,Tt2,lambda,kappa,r,l,d,A,iflag,
 
197
     &     inv,pt2zpt1_c,qred_crit,crit,qred_max1,icase)
 
198
!
 
199
         Qred=dabs(xflow)*dsqrt(Tt1)/(A*pt1)
 
200
!
 
201
         if (crit) then
 
202
            xflow=0.5*inv*Qred_crit*Pt1*A/dsqrt(Tt1)
 
203
            if(icase.eq.1) then
 
204
!     
 
205
               call ts_calc(xflow,Tt1,pt1,kappa,r,a,T1,icase)
 
206
               if (inv.eq.1) then 
 
207
                  v(3,node1)=T1
 
208
                  v(3,node2)=T1
 
209
                  if(nactdog(0,node2).eq.1) then
 
210
                     v(0,node2)=T1*(1.d0+km1/(2*kappa))
 
211
                  endif
 
212
               else
 
213
                  v(3,node2)=T1
 
214
                  v(3,node1)=T1
 
215
                  if(nactdog(0,node1).eq.1) then
 
216
                     v(0,node1)=T1*(1.d0+km1/(2*kappa)) 
 
217
                  endif
 
218
               endif
 
219
            endif
 
220
         elseif(Qred.gt.Qred_crit) then
 
221
            xflow=0.5*inv*Qred_crit*pt1*A/dsqrt(Tt1)
 
222
         else
 
223
            xflow=inv*Qred*pt1*A/dsqrt(Tt1)
 
224
         endif
 
225
!
 
226
      elseif (iflag.eq.2)then
 
227
!
 
228
         numf=5
 
229
         crit=.false.
 
230
!
 
231
         pi=4.d0*datan(1.d0)
 
232
         e=2.7182818d0
 
233
!
 
234
         kappa=(cp/(cp-R))
 
235
         km1=kappa-1.d0
 
236
         kp1=kappa+1.d0
 
237
         kdkm1=kappa/km1
 
238
         kdkp1=kappa/kp1
 
239
!
 
240
         index=ielprop(nelem)
 
241
         A=prop(index+1)
 
242
         d=prop(index+2)
 
243
!
 
244
         l=prop(index+3)
 
245
         if(l.lt.0d0) then
 
246
            l_neg=l
 
247
            l=abs(l)
 
248
         else
 
249
            l_neg=l
 
250
         endif
 
251
         ks=prop(index+4)
 
252
         if(lakon(nelem)(2:6).eq.'GAPFA') then
 
253
            icase=0
 
254
         elseif(lakon(nelem)(2:6).eq.'GAPFI') then
 
255
            icase=1
 
256
         endif
 
257
         form_fact=prop(index+5)
 
258
         xflow_oil=prop(index+6)
 
259
         k_oil=nint(prop(index+7))
 
260
!
 
261
                  if((lakon(nelem)(2:6).eq.'GAPFF').and.
 
262
     &        (lakon(nelem)(2:7).ne.'GAPFF2')) then
 
263
            icase=0
 
264
            nodea=nint(prop(index+1))
 
265
            nodeb=nint(prop(index+2))
 
266
c            iaxial=nint(prop(index+3))
 
267
            radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
 
268
     &           co(1,nodea)-vold(1,nodea))**2)
 
269
            initial_radius=dsqrt((co(1,nodeb)-co(1,nodea))**2)
 
270
            d=2*radius
 
271
c           if(iaxial.ne.0) then
 
272
c              A=pi*radius**2/iaxial
 
273
c            else
 
274
               A=pi*radius**2
 
275
c            endif
 
276
            l=prop(index+4)
 
277
            if(l.lt.0d0) then
 
278
               l_neg=l
 
279
               l=abs(l)
 
280
            else
 
281
               l_neg=l
 
282
            endif
 
283
            ks=prop(index+5)
 
284
            form_fact=prop(index+6)          
 
285
            xflow_oil=prop(index+7)
 
286
            k_oil=nint(prop(index+8))
 
287
!
 
288
         elseif (lakon(nelem)(2:7).eq.'GAPFF2') then
 
289
            icase=0
 
290
            nodea=nint(prop(index+1))
 
291
            nodeb=nint(prop(index+2))
 
292
            nodec=nint(prop(index+3))
 
293
c            iaxial=nint(prop(index+4))
 
294
            radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
 
295
     &           co(1,nodea)-vold(1,nodea))**2)
 
296
            initial_radius=dsqrt((co(1,nodeb)-co(1,nodea))**2)
 
297
            d=2*radius
 
298
c           if(iaxial.ne.0) then
 
299
c              A=pi*radius**2/iaxial
 
300
c            else
 
301
               A=pi*radius**2
 
302
c            endif
 
303
            l_initial=dsqrt((co(2,nodec)-co(2,nodeb))**2)
 
304
            l=-dsqrt((co(2,nodec)+vold(2,nodec)-
 
305
     &           co(2,nodeb)-vold(2,nodeb))**2)
 
306
            if(l.lt.0d0) then
 
307
               l_neg=l
 
308
               l=abs(l)
 
309
            else
 
310
               l_neg=l
 
311
            endif
 
312
            ks=prop(index+5)
 
313
            form_fact=prop(index+6)         
 
314
            xflow_oil=prop(index+7)
 
315
            k_oil=nint(prop(index+8))
 
316
         endif
 
317
!
 
318
         pt1=v(2,node1)
 
319
         pt2=v(2,node2)
 
320
         xflow=v(1,nodem)*iaxial
 
321
!
 
322
         if((pt1.gt.pt2).or.(xflow.ge.0d0)) then
 
323
            inv=1
 
324
            Tt1=v(0,node1)+physcon(1)
 
325
            call ts_calc(xflow,Tt1,Pt1,kappa,r,a,T1,icase)
 
326
            if(icase.eq.0) then
 
327
               Tt2=Tt1
 
328
               call ts_calc(xflow,Tt2,Pt2,kappa,r,a,T2,icase)
 
329
            else
 
330
               t2=t1
 
331
               Tt2=v(0,node2)+physcon(1)
 
332
            endif
 
333
!     
 
334
            nodef(1)=node1
 
335
            nodef(2)=node1
 
336
            nodef(3)=nodem
 
337
            nodef(4)=node2
 
338
            nodef(5)=node2
 
339
         else
 
340
            inv=-1
 
341
            pt1=v(2,node2)
 
342
            pt2=v(2,node1)
 
343
            xflow=-v(1,nodem)*iaxial
 
344
            Tt1=v(0,node2)+physcon(1)
 
345
            if(icase.eq.0) then
 
346
               Tt2=Tt1
 
347
            else
 
348
               Tt2=v(0,node1)+physcon(1)
 
349
            endif
 
350
!
 
351
            call ts_calc(xflow,Tt1,Pt1,kappa,r,a,T1,icase)
 
352
!
 
353
            call ts_calc(xflow,Tt2,Pt2,kappa,r,a,T2,icase)
 
354
!
 
355
            nodef(1)=node2
 
356
            nodef(2)=node2
 
357
            nodef(3)=nodem
 
358
            nodef(4)=node1
 
359
            nodef(5)=node1
 
360
         endif
 
361
!
 
362
         idirf(1)=2
 
363
         idirf(2)=0
 
364
         idirf(3)=1
 
365
         idirf(4)=2
 
366
         idirf(5)=0
 
367
!     
 
368
         pt2zpt1=pt2/pt1
 
369
!     
 
370
!     calculation of the dynamic viscosity
 
371
!     
 
372
           if(dabs(dvi).lt.1E-30) then
 
373
              kgas=0
 
374
              call dynamic_viscosity(kgas,T1,dvi)
 
375
           endif        
 
376
!           
 
377
           reynolds=dabs(xflow)*d/(dvi*a)
 
378
!
 
379
           if(reynolds.lt.1) then
 
380
              reynolds = 1.d0
 
381
           endif
 
382
!
 
383
!     definition of the friction coefficient for 2 phase flows and pure air
 
384
!     
 
385
!     Friedel's Method
 
386
         if(lakon(nelem)(7:7).eq.'F') then
 
387
!
 
388
            if((k_oil.lt.0).or.(k_oil.gt.12)) then
 
389
               write(*,*) '*ERROR:in gaspipe.f'
 
390
               write(*,*) ' using two phase flow'
 
391
               write(*,*) ' the type of oil is not defined'
 
392
               write(*,*) ' check element ',nelem,' definition'
 
393
               write(*,*) ' Current calculation stops here'
 
394
               call exit(201)
 
395
            elseif(xflow_oil.eq.0.d0) then
 
396
               write(*,*) '*WARNING:in gaspipe.f'
 
397
               write(*,*) ' using two phase flow'
 
398
               write(*,*) ' the oil mass flow rate is NULL'
 
399
               write(*,*) ' check element ',nelem,' definition'
 
400
               write(*,*) ' Only pure air is considered'
 
401
               call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
402
     &              lambda)
 
403
            else
 
404
               call two_phase_flow(Tt1,pt1,T1,Tt2,pt2,T2,xflow,
 
405
     &              xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
 
406
     &              v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
 
407
     &              shcon,rhcon,ntmat_,mi)
 
408
!
 
409
               lambda=lambda*phi
 
410
!
 
411
            endif
 
412
!     
 
413
!     Alber's Method
 
414
!     
 
415
         elseif (lakon(nelem)(7:7).eq.'A') then
 
416
            if((k_oil.lt.0).or.(k_oil.gt.12)) then
 
417
               write(*,*) '*ERROR:in gaspipe_fanno.f'
 
418
               write(*,*) ' using two phase flow'
 
419
               write(*,*) ' the type of oil is not defined'
 
420
               write(*,*) ' check element ',nelem,' definition'
 
421
               write(*,*) ' Current calculation stops here'
 
422
               call exit(201)
 
423
            elseif(xflow_oil.eq.0) then
 
424
               write(*,*) '*WARNING:in gaspipe_fanno.f'
 
425
               write(*,*) ' using two phase flow'
 
426
               write(*,*) ' the oil mass flow rate is NULL'
 
427
               write(*,*) ' check element ',nelem,' definition'
 
428
               write(*,*) ' Only pure air is considered'
 
429
               call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
430
     &              lambda)
 
431
            else
 
432
               call two_phase_flow(Tt1,pt1,T1,Tt2,pt2,T2,xflow,
 
433
     &              xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
 
434
     &              v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
 
435
     &              shcon,rhcon,ntmat_,mi)
 
436
!     
 
437
               call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
438
     &              lambda)
 
439
!
 
440
               lambda=lambda*phi
 
441
!
 
442
            endif
 
443
!     
 
444
!     for pure air
 
445
!     
 
446
         else
 
447
!     
 
448
            phi=1.d0
 
449
            call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
450
     &           lambda)
 
451
         endif
 
452
!     
 
453
         call pt2zpt1_crit(pt2,pt1,Tt1,Tt2,lambda,kappa,r,l,d,A,iflag,
 
454
     &     inv,pt2zpt1_c,qred_crit,crit,qred_max1,icase)
 
455
!
 
456
         Qred1=xflow*dsqrt(Tt1)/(A*Pt1)
 
457
!
 
458
         if(dabs(xflow)*dsqrt(Tt1)/(A*Pt1).gt.qred_max1) then
 
459
            crit=.true.
 
460
         endif
 
461
!
 
462
         Qred2=xflow*dsqrt(Tt2)/(A*Pt2)
 
463
         if(icase.eq.0) then
 
464
            qred_crit_out=dsqrt(kappa/R)*(2/(kappa+1))**(0.5d0*
 
465
     &           (kappa+1)/(kappa-1))
 
466
         else
 
467
            qred_crit_out=R**(-0.5d0)*(2/(kappa+1))**(0.5d0*
 
468
     &           (kappa+1)/(kappa-1))
 
469
         endif
 
470
!     
 
471
!     definition of the coefficients
 
472
!
 
473
         lld=lambda*l/d
 
474
!
 
475
         M2=dsqrt(2/km1*((Tt2/T2)-1))
 
476
         if(icase.eq.0) then
 
477
            if((M2.lt.1)) then
 
478
               crit=.false.
 
479
               if((M2.ge.1.d0).or.(dabs(M2-1).lt.1E-5)) then
 
480
                  pt2=pt1*pt2zpt1_c
 
481
               endif
 
482
            endif
 
483
         elseif (icase.eq.1) then
 
484
            if(M2.lt.1/dsqrt(kappa)) then
 
485
               crit=.false.
 
486
            else
 
487
               crit=.true.
 
488
            endif
 
489
         endif
 
490
!
 
491
!     adiabatic case
 
492
!     
 
493
         if(icase.eq.0) then
 
494
!
 
495
            T2dTt2=T2/Tt2
 
496
            Tt2dT2=1.d0/T2dTt2
 
497
            X_T2dTt2=T2dTt2**(2*kdkm1)
 
498
            T1dTt1=T1/Tt1
 
499
            Tt1dT1=1.d0/T1dTt1
 
500
            X_T1dTt1=T1dTt1**(2*kdkm1)
 
501
!     
 
502
            X2_den=pt2**2*X_T2dTt2
 
503
            X2=t2**2/X2_den
 
504
            X1_den=pt1**2*X_T1dTt1
 
505
            X1=T1**2/X1_den   
 
506
!     
 
507
!            C1=2.d0*cp*A**2*X1_den*(-1.d0+2.d0*kdkm1*T1dTt1)
 
508
!     &           -2.d0*xflow**2*R**2*T1
 
509
            C1=2.d0*cp*A**2*X1_den*(-1.d0+2.d0*kdkm1*(T1dTt1-1))
 
510
     &           -2.d0*xflow**2*R**2*T1
 
511
!     
 
512
!            C2=2.d0*cp*A**2*X2_den*(-1.d0+2.d0*kdkm1*T2dTt2)
 
513
!     &           -2.d0*xflow**2*R**2*T2
 
514
            C2=2.d0*cp*A**2*X2_den*(-1.d0+2.d0*kdkm1*(T2dTt2-1))
 
515
     &           -2.d0*xflow**2*R**2*T2
 
516
!     
 
517
            expon1=(kappa+1)/km1
 
518
            expon2=2*kappa/(km1)
 
519
!     
 
520
            cte=0.5d0*(kappa+1)/kappa
 
521
!
 
522
            term1=pt1**2*T1**expon1*Tt1**(-expon2)*A**2
 
523
!
 
524
            if(.not.crit) then
 
525
               term1=pt1**2*T1**expon1*Tt1**(-expon2)*A**2
 
526
               term2=pt2**2*T2**expon1*Tt2**(-expon2)*A**2        
 
527
!     
 
528
!     simplified version
 
529
               term3=Tt2dT2
 
530
               term4=Tt1dT1
 
531
!     
 
532
               term5=T1**(expon1)*Tt1**(-expon2)*(pt1**2)               
 
533
               term6=T2**(expon1)*Tt2**(-expon2)*(pt2**2)               
 
534
!
 
535
               B1=1/(R*xflow**2)*term1*expon1/T1
 
536
     &              +cte*(-(2/km1)*1/T1)
 
537
!     
 
538
               B2=1/(R*xflow**2)*term2*(-expon1/T2)
 
539
     &              +cte*(2/km1*1/T2)
 
540
!     
 
541
!     residual
 
542
!    
 
543
!     Simplified version
 
544
!     
 
545
               f=1/(R*xflow**2)*(term1-term2)
 
546
     &              +cte*(log(term3)-log(term4)-log(term5)+log(term6))
 
547
     &              -lld 
 
548
     &              +b2/c2*(2*cp*A**2*(Tt2-T2)
 
549
     &              *X2_den-xflow**2*R**2*T2**2)
 
550
     &              +b1/c1*(2*cp*A**2*(Tt1-T1)
 
551
     &              *X1_den-xflow**2*R**2*T1**2)
 
552
!     
 
553
!     pressure node1
 
554
!     
 
555
               df(1)=1/(R*xflow**2)*(term1*2/pt1)
 
556
     &              +cte*(-2/pt1)
 
557
     &              +B1/C1*(4.d0*cp*A**2*(Tt1-T1)*pt1*X_T1dTt1)
 
558
!     
 
559
!     temperature node1
 
560
!     
 
561
               df(2)=1/(R*xflow**2)*term1*(-expon2)/Tt1
 
562
     &              +cte*(expon1*1/Tt1)
 
563
     &              +b1/c1*(2*cp*A**2*X1_den
 
564
     &              *(1.d0-2.d0*kdkm1*(Tt1-T1)/Tt1))
 
565
!     
 
566
!     mass flow
 
567
!     
 
568
               df(3)=-2.d0/(R*(inv*xflow)**3)*(term1-term2)
 
569
     &              +B2/C2*(-2.d0*inv*xflow*R*R*T2**2.d0)
 
570
     &              +B1/C1*(-2.d0*inv*xflow*R*R*T1**2.d0)
 
571
!     
 
572
!     pressure node2
 
573
!     
 
574
               df(4)=1/(R*xflow**2)*(-term2*2/pt2)
 
575
     &              +cte*(2/pt2)
 
576
     &              +B2/C2*(4.d0*cp*A**2*(Tt2-T2)*pt2*X_T2dTt2)
 
577
!     
 
578
!     temperature node2
 
579
!     
 
580
               df(5)=1/(R*xflow**2)*term2*(expon2/Tt2)
 
581
     &              +cte*(-expon1*1/Tt2)
 
582
     &              +b2/c2*(2*cp*A**2*X2_den
 
583
     &              *(1.d0-2.d0*kdkm1*(Tt2-T2)/Tt2))
 
584
!               
 
585
            else
 
586
!
 
587
               term=kappa*term1/(xflow**2*R)
 
588
               B1=expon1*1/T1*(1/kappa*term-1)+cte*1/T1
 
589
!     f=1/kappa*(term1-1)+cte*(log(T1dTt1)-log(2/kp1*term))
 
590
               f=1/kappa*(term-1)+cte*(log(T1dTt1)-log(2/kp1*term))
 
591
     &              -lld
 
592
     &              +b1/c1*(2*cp*A**2*(Tt1-T1)
 
593
     &              *X1_den-xflow**2*R**2*T1**2)
 
594
!     
 
595
!     pressure node1
 
596
!     
 
597
               df(1)=2/pt1*(1/kappa*term-cte)
 
598
     &              +B1/C1*(4.d0*cp*A**2*(Tt1-T1)*pt1*X_T1dTt1)
 
599
!     
 
600
!     temperature node1
 
601
!     
 
602
               df(2)=expon2*1/Tt1*(-1/kappa*term+1)-cte*1/Tt1
 
603
     &              +b1/c1*(2*cp*A**2*X1_den
 
604
     &              *(1.d0-2.d0*kdkm1*(Tt1-T1)/Tt1))
 
605
!     
 
606
!     mass flow
 
607
!     
 
608
               df(3)=2.d0/(inv*xflow)*(-term/kappa+cte)
 
609
     &              +B1/C1*(-2.d0*inv*xflow*R*R*T1**2.d0)
 
610
!     
 
611
!     pressure node2
 
612
!     
 
613
               df(4)=0.d0
 
614
!     
 
615
!     temperature node2
 
616
!     
 
617
               df(5)=0.d0
 
618
!
 
619
            endif
 
620
!     
 
621
!     isothermal icase
 
622
!     
 
623
         elseif(icase.eq.1) then
 
624
               T2dTt2=T2/Tt2
 
625
               Tt2dT2=1.d0/T2dTt2
 
626
               X_T2dTt2=T2dTt2**(2*kdkm1)
 
627
               T1dTt1=T1/Tt1
 
628
               Tt1dT1=1.d0/T1dTt1
 
629
               X_T1dTt1=T1dTt1**(2*kdkm1)
 
630
!     
 
631
               X2_den=pt2**2*X_T2dTt2
 
632
               X2=t2**2/X2_den
 
633
               X1_den=pt1**2*X_T1dTt1
 
634
               X1=T1**2/X1_den   
 
635
!     
 
636
               C1=2.d0*cp*A**2*X1_den*(1.d0-2.d0*kdkm1*(Tt1dT1-1.d0))
 
637
     &              +2.d0*xflow**2*R**2*T1
 
638
!     
 
639
               C2=2.d0*cp*A**2*X2_den*(1.d0-2.d0*kdkm1*(Tt2dT2-1.d0))
 
640
     &              +2.d0*xflow**2*R**2*T2
 
641
!     
 
642
               expon1=(kappa+1)/km1
 
643
               expon2=2*kappa/(kappa-1)
 
644
!     
 
645
               cte=0.5d0*(kappa+1)/kappa
 
646
!     
 
647
               term1=pt1**2*T1**expon1*Tt1**(-expon2)*A**2
 
648
               term2=pt2**2*T2**expon1*Tt2**(-expon2)*A**2
 
649
!     
 
650
               term5=T1**(expon1)*Tt1**(-expon2)*(pt1**2*A**2)
 
651
               term6=T2**(expon1)*Tt2**(-expon2)*(pt2**2*A**2)
 
652
!
 
653
            if(.not.crit) then
 
654
               B1=1/(R*xflow**2)*term1*expon1/T1
 
655
     &              -expon1/T1   
 
656
!     
 
657
               B2=1/(R*xflow**2)*term2*(-expon1/T2)
 
658
     &              +expon1/T2  
 
659
!     
 
660
!     Simplified version
 
661
!     
 
662
               f=1/(R*xflow**2)*(term1-term2)
 
663
     &              +(-log(term5)+log(term6))            
 
664
     &              -lld 
 
665
     &              +b2/c2*(2*cp*A**2*(Tt2-T2)
 
666
     &              *X2_den-xflow**2*R**2*T2**2)
 
667
     &              +b1/c1*(2*cp*A**2*(Tt1-T1)
 
668
     &              *X1_den-xflow**2*R**2*T1**2)
 
669
!     
 
670
!     pressure node1
 
671
!     
 
672
               df(1)=1/(R*xflow**2)*(term1*2/pt1)
 
673
     &              +(-(2/pt1))
 
674
     &              +B1/C1*(4.d0*cp*A**2*(Tt1-T1)*pt1*X_T1dTt1)
 
675
!     
 
676
!     temperature node1
 
677
!     
 
678
               df(2)=1/(R*xflow**2)*term1*(-expon2)/Tt1
 
679
     &              +(expon2/Tt1)            
 
680
     &              +b1/c1*(2*cp*A**2*X1_den
 
681
     &              *(1.d0-2.d0*kdkm1*(Tt1-T1)/Tt1))
 
682
!     
 
683
!     mass flow
 
684
!     
 
685
               df(3)=-2.d0/(R*xflow**3)*(term1-term2)
 
686
     &              +B2/C2*(-2.d0*inv*xflow*R*R*T2**2.d0)
 
687
     &              +B1/C1*(-2.d0*inv*xflow*R*R*T1**2.d0)
 
688
!     
 
689
!     pressure node2
 
690
!     
 
691
               df(4)=1/(R*xflow**2)*(-term2*2/pt2)
 
692
     &              +(2/pt2)
 
693
     &              +B2/C2*(4.d0*cp*A**2*(Tt2-T2)*pt2*X_T2dTt2)
 
694
!     
 
695
!     
 
696
!     temperature node2
 
697
!     
 
698
               df(5)=1/(R*xflow**2)*term2*(expon2/Tt2)
 
699
     &              +(-expon2/Tt2)            
 
700
     &              +b2/c2*(2*cp*A**2*X2_den
 
701
     &              *(1.d0-2.d0*kdkm1*(Tt2-T2)/Tt2))
 
702
               
 
703
            else
 
704
               term=term1/(xflow**2*R)
 
705
               B1=expon1/T1*(term-1)
 
706
!     alternate critical equation
 
707
!  
 
708
               f=term-1-log(term)            
 
709
     &              -lld 
 
710
     &              +b1/c1*(2*cp*A**2*(Tt1-T1)
 
711
     &              *X1_den-xflow**2*R**2*T1**2)
 
712
!     
 
713
!     pressure node1
 
714
!     
 
715
               df(1)=2/pt1*(term-1)
 
716
     &              +B1/C1*(4.d0*cp*A**2*(Tt1-T1)*pt1*X_T1dTt1)
 
717
!     
 
718
!     temperature node1
 
719
!     
 
720
               df(2)=expon2/Tt1*(-term+1)            
 
721
     &              +b1/c1*(2*cp*A**2*X1_den
 
722
     &              *(1.d0-2.d0*kdkm1*(Tt1-T1)/Tt1))
 
723
!     
 
724
!     mass flow
 
725
!     
 
726
               df(3)=2/xflow*(-term+1)
 
727
     &              +B1/C1*(-2.d0*inv*xflow*R*R*T1**2.d0)
 
728
!     
 
729
!     pressure node2
 
730
!     
 
731
               df(4)=0.d0
 
732
!     
 
733
!     temperature node2
 
734
!     
 
735
               df(5)=0.d0
 
736
!     
 
737
            endif
 
738
         endif
 
739
!
 
740
!     output
 
741
!
 
742
      elseif(iflag.eq.3) then
 
743
!
 
744
         pi=4.d0*datan(1.d0)
 
745
         e=2.7182818d0
 
746
!     
 
747
         kappa=(cp/(cp-R))
 
748
         km1=kappa-1.d0
 
749
         kp1=kappa+1.d0
 
750
         kdkm1=kappa/km1
 
751
         kdkp1=kappa/kp1
 
752
!         
 
753
         index=ielprop(nelem)
 
754
         A=prop(index+1)
 
755
         d=prop(index+2)    
 
756
         l=prop(index+3)
 
757
!
 
758
         lambda=0.5
 
759
!
 
760
         if(l.lt.0d0) then
 
761
            l_neg=l
 
762
            l=abs(l)
 
763
         else
 
764
            l_neg=l
 
765
         endif
 
766
         ks=prop(index+4)
 
767
         if(lakon(nelem)(2:6).eq.'GAPFA') then
 
768
            icase=0
 
769
         elseif(lakon(nelem)(2:6).eq.'GAPFI') then
 
770
            icase=1
 
771
         endif
 
772
         form_fact=prop(index+5)
 
773
         xflow_oil=prop(index+6)
 
774
         k_oil=nint(prop(index+7))
 
775
!
 
776
         pt1=v(2,node1)
 
777
         pt2=v(2,node2)
 
778
!     
 
779
         if(xflow.ge.0d0) then
 
780
            inv=1
 
781
            xflow=v(1,nodem)*iaxial
 
782
            Tt1=v(0,node1)+physcon(1)
 
783
            call ts_calc(xflow,Tt1,Pt1,kappa,r,a,T1,icase)
 
784
            if(icase.eq.0) then
 
785
               Tt2=Tt1
 
786
               call ts_calc(xflow,Tt2,Pt2,kappa,r,a,T2,icase)
 
787
            else
 
788
               T2=T1
 
789
               Tt2=v(0,node2)+physcon(1)
 
790
               call tt_calc(xflow,Tt2,Pt2,kappa,r,a,T2,icase,iflag)
 
791
            endif
 
792
!      
 
793
!            call ts_calc(xflow,Tt1,Pt1,kappa,r,a,T1,icase)
 
794
!!
 
795
!            call ts_calc(xflow,Tt2,Pt2,kappa,r,a,T2,icase)
 
796
!     
 
797
         else
 
798
            inv=-1
 
799
            pt1=v(2,node2)
 
800
            pt2=v(2,node1)
 
801
            xflow=v(1,nodem)*iaxial
 
802
!           
 
803
            Tt1=v(0,node2)+physcon(1)
 
804
            if(icase.eq.0) then
 
805
               Tt2=Tt1
 
806
            else
 
807
               Tt2=v(0,node1)+physcon(1)
 
808
            endif
 
809
!               
 
810
            call ts_calc(xflow,Tt1,Pt1,kappa,r,a,T1,icase)
 
811
            call ts_calc(xflow,Tt2,Pt2,kappa,r,a,T2,icase)
 
812
!     
 
813
         endif
 
814
!     
 
815
         pt2zpt1=pt2/pt1
 
816
!     
 
817
!     calculation of the dynamic viscosity 
 
818
!     
 
819
         if(dabs(dvi).lt.1E-30) then
 
820
            kgas=0
 
821
            call dynamic_viscosity(kgas,T1,dvi)
 
822
         endif
 
823
! 
 
824
         reynolds=dabs(xflow)*d/(dvi*a)
 
825
!
 
826
        if(reynolds.lt.1.d0) then
 
827
            reynolds= 1.d0
 
828
         endif
 
829
!     
 
830
!     definition of the friction coefficient for 2 phase flows and pure air
 
831
!     
 
832
!     Friedel's Method
 
833
         if(lakon(nelem)(7:7).eq.'F') then
 
834
!     
 
835
            if((k_oil.lt.0).or.(k_oil.gt.12)) then
 
836
               write(*,*) '*ERROR:in gaspipe.f'
 
837
               write(*,*) ' using two phase flow'
 
838
               write(*,*) ' the type of oil is not defined'
 
839
               write(*,*) ' check element ',nelem,' definition'
 
840
               write(*,*) ' Current calculation stops here'
 
841
               call exit(201)
 
842
            elseif(xflow_oil.eq.0) then
 
843
               write(*,*) '*WARNING:in gaspipe.f'
 
844
               write(*,*) ' using two phase flow'
 
845
               write(*,*) ' the oil mass flow rate is NULL'
 
846
               write(*,*) ' check element ',nelem,' definition'
 
847
               write(*,*) ' Only pure air is considered'
 
848
               call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
849
     &              lambda)
 
850
            else
 
851
               call two_phase_flow(Tt1,pt1,T1,Tt2,pt2,T2,xflow,
 
852
     &              xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
 
853
     &              v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
 
854
     &              shcon,rhcon,ntmat_,mi)
 
855
!
 
856
               call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
857
     &              lambda)
 
858
!
 
859
            endif
 
860
!     
 
861
         elseif (lakon(nelem)(7:7).eq.'A') then
 
862
            if((k_oil.lt.0).or.(k_oil.gt.12)) then
 
863
               write(*,*) '*ERROR:in gaspipe.f'
 
864
               write(*,*) ' using two phase flow'
 
865
               write(*,*) ' the type of oil is not defined'
 
866
               write(*,*) ' check element ',nelem,' definition'
 
867
               write(*,*) ' Current calculation stops here'
 
868
               call exit(201)
 
869
            elseif(xflow_oil.eq.0) then
 
870
               write(*,*) '*WARNING:in gaspipe.f'
 
871
               write(*,*) ' using two phase flow'
 
872
               write(*,*) ' the oil mass flow rate is NULL'
 
873
               write(*,*) ' check element ',nelem,' definition'
 
874
               write(*,*) ' Only pure air is considered'
 
875
               call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
876
     &              lambda)
 
877
            else
 
878
               call two_phase_flow(Tt1,pt1,T1,Tt2,pt2,T2,xflow,
 
879
     &              xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
 
880
     &              v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
 
881
     &              shcon,rhcon,ntmat_,mi)
 
882
!     
 
883
               call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
884
     &              lambda)
 
885
!
 
886
            endif
 
887
!     
 
888
!     for pure air
 
889
!     
 
890
         else
 
891
            phi=1.d0
 
892
            call friction_coefficient(l_neg,d,ks,reynolds,form_fact,
 
893
     &           lambda)
 
894
         endif
 
895
!     
 
896
         call pt2zpt1_crit(pt2,pt1,Tt1,Tt2,lambda,kappa,r,l,d,A,iflag,
 
897
     &     inv,pt2zpt1_c,qred_crit,crit,qred_max1,icase)
 
898
!     
 
899
!     definition of the coefficients 
 
900
!     
 
901
         M1=dsqrt(2/km1*((Tt1/T1)-1))
 
902
         M2=dsqrt(2/km1*((Tt2/T2)-1))
 
903
!     
 
904
            write(1,*) ''
 
905
            write(1,55) ' from node ',node1,
 
906
     &           ' to node ', node2,':   air massflow rate = ',xflow,
 
907
     &           ' , oil massflow rate = ',xflow_oil
 
908
!     
 
909
            if(inv.eq.1) then
 
910
               write(1,53)'       Inlet node ',node1,' :    Tt1 = ',Tt1,
 
911
     &              ' , Ts1 = ',T1,'  , Pt1 = ',Pt1,
 
912
     &              ' ,M1 = ',M1
 
913
               write(1,*)'             Element ',nelem,lakon(nelem)
 
914
               write(1,57)'             eta = ',dvi,' kg/(m*s) , Re = '
 
915
     &              ,reynolds
 
916
               write(1,58)'             PHI = ',phi,' , LAMBDA = ',
 
917
     &              lambda,
 
918
     &              ', LAMBDA*l/d = ',lambda*l/d,' , ZETA_PHI = ',
 
919
     &              phi*lambda*l/d
 
920
               write(1,53)'       Outlet node ',node2,' :    Tt2 = ',
 
921
     &              Tt2,
 
922
     &              '  , Ts2 = ',T2,'  , Pt2 = ',Pt2,
 
923
     &              ' , M2 = ',M2
 
924
!    
 
925
            else if(inv.eq.-1) then
 
926
               write(1,53)'       Inlet node ',node2,':    Tt1= ',Tt1,
 
927
     &              ' , Ts1= ',T1,' , Pt1= ',Pt1,
 
928
     &              ' , M1= ',M1
 
929
               write(1,*)'             Element ',nelem,lakon(nelem)
 
930
               write(1,57)'             eta = ',dvi,' kg/(m*s) , Re = '
 
931
     &              ,reynolds
 
932
               write(1,58)'             PHI = ',phi,' , LAMBDA = ',
 
933
     &              lambda,
 
934
     &              ', LAMBDA*l/d = ',lambda*l/d,' , ZETA_PHI = ',
 
935
     &              phi*lambda*l/d
 
936
               write(1,53)'       Outlet node ',node1,' :    Tt2 = ',
 
937
     &              Tt2,
 
938
     &              '  , Ts2 = ',T2,'  , Pt2 =',Pt2,
 
939
     &              ' , M2 = ',M2
 
940
            endif
 
941
      endif
 
942
!     
 
943
 55   format(1X,a,i6,a,i6,a,e11.4,a,a,e11.4,a)
 
944
 53   format(1X,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
 
945
 57   format(1X,a,e11.4,a,e11.4)
 
946
 58   format(1X,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
 
947
 77   format(3x,a,e11.1,a,e11.4,a,e11.4,a,e11.4,a,e11.4,a,e11.4
 
948
     &     ,a,e11.4)     
 
949
 78   format(a,i4,a,a,a)
 
950
 79   format(3X,a,i4,a,i4,a,i4)
 
951
 80   format(3x,a,e11.4,a,e11.4,a,e11.4)
 
952
!     
 
953
      xflow=xflow/iaxial
 
954
      df(3)=df(3)*iaxial
 
955
!     
 
956
      return
 
957
      end