CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
liquidpipe.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine liquidpipe (node1, node2, nodem, nelem, lakon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, rho, g, co, dvi, numf, vold, mi, ipkon, kon, set, ttime, time, iaxial)
 

Function/Subroutine Documentation

◆ liquidpipe()

subroutine liquidpipe ( integer  node1,
integer  node2,
integer  nodem,
integer  nelem,
character*8, dimension(*)  lakon,
integer, dimension(0:3,*)  nactdog,
logical  identity,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
integer  iflag,
real*8, dimension(0:mi(2),*)  v,
real*8  xflow,
real*8  f,
integer, dimension(*)  nodef,
integer, dimension(*)  idirf,
real*8, dimension(*)  df,
real*8  rho,
real*8, dimension(3)  g,
real*8, dimension(3,*)  co,
real*8  dvi,
integer  numf,
real*8, dimension(0:mi(2),*)  vold,
integer, dimension(*)  mi,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
character*81, dimension(*)  set,
real*8  ttime,
real*8  time,
integer  iaxial 
)
23 !
24 ! pipe element for incompressible media
25 !
26 ! iflag=0: check whether element equation is needed
27 ! iflag=1: calculate mass flow
28 ! iflag=2: calculate residual and derivative w.r.t.independent
29 ! variables
30 ! iflag=3: output
31 !
32  implicit none
33 !
34  logical identity,flowunknown
35  character*8 lakon(*)
36  character*81 set(*)
37 !
38  integer nelem,nactdog(0:3,*),node1,node2,nodem,iaxial,
39  & ielprop(*),nodef(*),idirf(*),index,iflag,mi(*),
40  & inv,ncoel,ndi,nbe,id,nen,ngv,numf,nodea,nodeb,
41  & ipkon(*),isothermal,kon(*),nelemswirl
42 !
43  real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),a,d,pi,radius,
44  & p1,p2,rho,dvi,friction,reynolds,vold(0:mi(2),*),
45  & g(3),a1,a2,xn,xk,xk1,xk2,zeta,dl,dg,rh,a0,alpha,
46  & coarseness,rd,xks,z1,z2,co(3,*),xcoel(11),yel(11),
47  & yco(11),xdi(10),ydi(10),xbe(7),ybe(7),zbe(7),ratio,
48  & xen(10),yen(10),xgv(8),ygv(8),xkn,xkp,ttime,time,
49  & dh,kappa,r,dkda,form_fact,dzetadalpha,t_chang,
50  & xflow_vol,r1d,r2d,r1,r2,eta, k1, kr, u1,ui, ciu, c1u,
51  & c2u, omega,cinput,un,t
52 !
53  data ncoel /11/
54  data xcoel /0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/
55  data yco /0.5,0.46,0.41,0.36,0.30,0.24,0.18,0.12,0.06,0.02,0./
56  data yel /1.,0.81,0.64,0.49,0.36,0.25,0.16,0.09,0.04,0.01,0./
57 !
58  data ndi /10/
59  data xdi /0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1./
60  data ydi /226.,47.5,17.5,7.8,3.75,1.80,0.8,0.29,0.06,0./
61 !
62  data nbe /7/
63  data xbe /1.,1.5,2.,3.,4.,6.,10./
64  data ybe /0.21,0.12,0.10,0.09,0.09,0.08,0.2/
65  data zbe /0.51,0.32,0.29,0.26,0.26,0.17,0.31/
66 !
67  data nen /10/
68  data xen /0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1./
69  data yen /232.,51.,18.8,9.6,5.26,3.08,1.88,1.17,0.734,0.46/
70 !
71  data ngv /8/
72  data xgv /0.125,0.25,0.375,0.5,0.625,0.75,0.875,1./
73  data ygv /98.,17.,5.52,2.,0.81,0.26,0.15,0.12/
74 !
75  numf=4
76 !
77  pi=4.d0*datan(1.d0)
78  dkda=0.d0
79 !
80  if (iflag.eq.0) then
81  identity=.true.
82 !
83  if(nactdog(2,node1).ne.0)then
84  identity=.false.
85  elseif(nactdog(2,node2).ne.0)then
86  identity=.false.
87  elseif(nactdog(1,nodem).ne.0)then
88  identity=.false.
89  elseif(nactdog(3,nodem).ne.0) then
90  identity=.false.
91  endif
92 !
93  elseif((iflag.eq.1).or.(iflag.eq.2).or.(iflag.eq.3))then
94 !
95  index=ielprop(nelem)
96 !
97  p1=v(2,node1)
98  p2=v(2,node2)
99 !
100  z1=-g(1)*co(1,node1)-g(2)*co(2,node1)-g(3)*co(3,node1)
101  z2=-g(1)*co(1,node2)-g(2)*co(2,node2)-g(3)*co(3,node2)
102 !
103  t=v(0,node1)
104 !
105  if(iflag.eq.1) then
106  inv=0
107  if(nactdog(1,nodem).ne.0) then
108  flowunknown=.true.
109  else
110  flowunknown=.false.
111  xflow=v(1,nodem)*iaxial
112  endif
113  else
114  xflow=v(1,nodem)*iaxial
115  if(xflow.ge.0.d0) then
116  inv=1
117  else
118  inv=-1
119  endif
120  nodef(1)=node1
121  nodef(2)=nodem
122  nodef(3)=node2
123  nodef(4)=nodem
124  idirf(1)=2
125  idirf(2)=1
126  idirf(3)=2
127  idirf(4)=3
128  endif
129 !
130  if((lakon(nelem)(4:5).ne.'BE').and.
131  & (lakon(nelem)(6:7).eq.'MA')) then
132 !
133 ! pipe, Manning (LIPIMA)
134 !
135  if(lakon(nelem)(8:8).eq.'F') then
136  nodea=nint(prop(index+1))
137  nodeb=nint(prop(index+2))
138  xn=prop(index+3)
139  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
140  & co(1,nodea)-vold(1,nodea))**2+
141  & (co(2,nodeb)+vold(2,nodeb)-
142  & co(2,nodea)-vold(2,nodea))**2+
143  & (co(3,nodeb)+vold(3,nodeb)-
144  & co(3,nodea)-vold(3,nodea))**2)
145  a=pi*radius*radius
146  rh=radius/2.d0
147  else
148  a=prop(index+1)
149  rh=prop(index+2)
150  endif
151  xn=prop(index+3)
152  a1=a
153  a2=a
154  dl=dsqrt((co(1,node2)-co(1,node1))**2+
155  & (co(2,node2)-co(2,node1))**2+
156  & (co(3,node2)-co(3,node1))**2)
157  dg=dsqrt(g(1)*g(1)+g(2)*g(2)+g(3)*g(3))
158  if(inv.ne.0) then
159  xk=2.d0*xn*xn*dl*dg/(a*a*rh**(4.d0/3.d0))
160  else
161  xkn=2.d0*xn*xn*dl*dg/(a*a*rh**(4.d0/3.d0))
162  xkp=xkn
163  endif
164  elseif(lakon(nelem)(6:7).eq.'WC') then
165 !
166 ! pipe, White-Colebrook
167 !
168  if(lakon(nelem)(8:8).eq.'F') then
169  nodea=nint(prop(index+1))
170  nodeb=nint(prop(index+2))
171  xn=prop(index+3)
172  radius=dsqrt((co(1,nodeb)+vold(1,nodeb)-
173  & co(1,nodea)-vold(1,nodea))**2+
174  & (co(2,nodeb)+vold(2,nodeb)-
175  & co(2,nodea)-vold(2,nodea))**2+
176  & (co(3,nodeb)+vold(3,nodeb)-
177  & co(3,nodea)-vold(3,nodea))**2)
178  a=pi*radius*radius
179  d=2.d0*radius
180  else
181  a=prop(index+1)
182  d=prop(index+2)
183  endif
184  dl=prop(index+3)
185  if(dl.le.0.d0) then
186  dl=dsqrt((co(1,node2)-co(1,node1))**2+
187  & (co(2,node2)-co(2,node1))**2+
188  & (co(3,node2)-co(3,node1))**2)
189  endif
190  xks=prop(index+4)
191  form_fact=prop(index+5)
192  a1=a
193  a2=a
194  if(iflag.eq.1) then
195 !
196 ! assuming large reynolds number
197 !
198  friction=1.d0/(2.03*dlog10(xks/(d*3.7)))**2
199  else
200 !
201 ! solving the implicit White-Colebrook equation
202 !
203  reynolds=xflow*d/(a*dvi)
204  call friction_coefficient(dl,d,xks,reynolds,form_fact,
205  & friction)
206  endif
207  if(inv.ne.0) then
208  xk=friction*dl/(d*a*a)
209  dkda=-2.5d0*xk/a
210  else
211  xkn=friction*dl/(d*a*a)
212  xkp=xkn
213  endif
214  elseif(lakon(nelem)(6:7).eq.'EL') then
215 !
216 ! pipe, sudden enlargement Berlamont version: fully turbulent
217 ! all section ratios
218 !
219  a1=prop(index+1)
220  a2=prop(index+2)
221  ratio=a1/a2
222  call ident(xcoel,ratio,ncoel,id)
223  if(inv.ge.0) then
224  if(id.eq.0) then
225  zeta=yel(1)
226  elseif(id.eq.ncoel) then
227  zeta=yel(ncoel)
228  else
229  zeta=yel(id)+(yel(id+1)-yel(id))*(ratio-xcoel(id))/
230  & (xcoel(id+1)-xcoel(id))
231  endif
232  if(inv.ne.0) then
233  xk=zeta/(a1*a1)
234  else
235  xkp=zeta/(a1*a1)
236  endif
237  endif
238  if(inv.le.0) then
239  if(id.eq.0) then
240  zeta=yco(1)
241  elseif(id.eq.ncoel) then
242  zeta=yco(ncoel)
243  else
244  zeta=yco(id)+(yco(id+1)-yco(id))*(ratio-xcoel(id))/
245  & (xcoel(id+1)-xcoel(id))
246  endif
247  if(inv.ne.0) then
248  xk=zeta/(a1*a1)
249  else
250  xkn=zeta/(a1*a1)
251  endif
252  endif
253  elseif(lakon(nelem)(4:5).eq.'EL') then
254 !
255 ! pipe, sudden enlargement Idelchik version: reynolds dependent,
256 ! 0.01 <= section ratio <= 0.6
257 !
258  a1=prop(index+1)
259  a2=prop(index+2)
260  dh=prop(index+3)
261  if(dh.eq.0.d0) then
262  dh=dsqrt(4*a1/pi)
263  endif
264  if(inv.eq.0) then
265  reynolds=5000.d0
266  else
267  reynolds=xflow*dh/(dvi*a1)
268  endif
269  if(inv.ge.0) then
270  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
271  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
272  if(inv.ne.0) then
273  xk=zeta/(a1*a1)
274  else
275  xkp=zeta/(a1*a1)
276  endif
277  endif
278  if(inv.le.0) then
279  reynolds=-reynolds
280 !
281 ! setting length and angle for contraction to zero
282 !
283  prop(index+4)=0.d0
284  prop(index+5)=0.d0
285  lakon(nelem)(4:5)='CO'
286  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
287  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
288  lakon(nelem)(4:5)='EL'
289  if(inv.ne.0) then
290  xk=zeta/(a1*a1)
291  else
292  xkn=zeta/(a1*a1)
293  endif
294  endif
295  elseif(lakon(nelem)(6:7).eq.'CO') then
296 !
297 ! pipe, sudden contraction Berlamont version: fully turbulent
298 ! all section ratios
299 !
300  a1=prop(index+1)
301  a2=prop(index+2)
302  ratio=a2/a1
303  call ident(xcoel,ratio,ncoel,id)
304  if(inv.ge.0) then
305  if(id.eq.0) then
306  zeta=yco(1)
307  elseif(id.eq.ncoel) then
308  zeta=yco(ncoel)
309  else
310  zeta=yco(id)+(yco(id+1)-yco(id))*(ratio-xcoel(id))/
311  & (xcoel(id+1)-xcoel(id))
312  endif
313  if(inv.ne.0) then
314  xk=zeta/(a2*a2)
315  else
316  xkp=zeta/(a2*a2)
317  endif
318  endif
319  if(inv.le.0) then
320  if(id.eq.0) then
321  zeta=yel(1)
322  elseif(id.eq.ncoel) then
323  zeta=yel(ncoel)
324  else
325  zeta=yel(id)+(yel(id+1)-yel(id))*(ratio-xcoel(id))/
326  & (xcoel(id+1)-xcoel(id))
327  endif
328  if(inv.ne.0) then
329  xk=zeta/(a2*a2)
330  else
331  xkn=zeta/(a2*a2)
332  endif
333  endif
334  elseif(lakon(nelem)(4:5).eq.'CO') then
335 !
336 ! pipe, sudden contraction Idelchik version: reynolds dependent,
337 ! 0.1 <= section ratio <= 0.6
338 !
339  a1=prop(index+1)
340  a2=prop(index+2)
341  dh=prop(index+3)
342  if(dh.eq.0.d0) then
343  dh=dsqrt(4*a2/pi)
344  endif
345  if(inv.eq.0) then
346  reynolds=5000.d0
347  else
348  reynolds=xflow*dh/(dvi*a2)
349  endif
350  if(inv.ge.0) then
351  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
352  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
353  if(inv.ne.0) then
354  xk=zeta/(a2*a2)
355  else
356  xkp=zeta/(a2*a2)
357  endif
358  endif
359  if(inv.le.0) then
360  reynolds=-reynolds
361  lakon(nelem)(4:5)='EL'
362  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
363  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
364  lakon(nelem)(4:5)='CO'
365  if(inv.ne.0) then
366  xk=zeta/(a2*a2)
367  else
368  xkn=zeta/(a2*a2)
369  endif
370  endif
371  elseif(lakon(nelem)(6:7).eq.'DI') then
372 !
373 ! pipe, diaphragm
374 !
375  a=prop(index+1)
376  a0=prop(index+2)
377  a1=a
378  a2=a
379  ratio=a0/a
380  call ident(xdi,ratio,ndi,id)
381  if(id.eq.0) then
382  zeta=ydi(1)
383  elseif(id.eq.ndi) then
384  zeta=ydi(ndi)
385  else
386  zeta=ydi(id)+(ydi(id+1)-ydi(id))*(ratio-xdi(id))/
387  & (xdi(id+1)-xdi(id))
388  endif
389  if(inv.ne.0) then
390  xk=zeta/(a*a)
391  else
392  xkn=zeta/(a*a)
393  xkp=xkn
394  endif
395  elseif(lakon(nelem)(6:7).eq.'EN') then
396 !
397 ! pipe, entrance (Berlamont data)
398 !
399  a=prop(index+1)
400  a0=prop(index+2)
401  a1=a*1.d10
402  a2=a
403  ratio=a0/a
404  call ident(xen,ratio,nen,id)
405  if(id.eq.0) then
406  zeta=yen(1)
407  elseif(id.eq.nen) then
408  zeta=yen(nen)
409  else
410  zeta=yen(id)+(yen(id+1)-yen(id))*(ratio-xen(id))/
411  & (xen(id+1)-xen(id))
412  endif
413  if(inv.ne.0) then
414  if(inv.gt.0) then
415 ! entrance
416  xk=zeta/(a*a)
417  else
418 ! exit
419  xk=1.d0/(a*a)
420  endif
421  else
422  xkn=1.d0/(a*a)
423  xkp=zeta/(a*a)
424  endif
425  elseif(lakon(nelem)(4:5).eq.'EN') then
426 !
427 ! pipe, entrance (Idelchik)
428 !
429  a1=prop(index+1)
430  a2=prop(index+2)
431  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
432  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
433 !
434 ! check for negative flow: in that case the loss
435 ! coefficient is wrong
436 !
437  if(inv.lt.0) then
438  write(*,*) '*ERROR in liquidpipe: loss coefficients'
439  write(*,*) ' for entrance (Idelchik) do not apply'
440  write(*,*) ' to reversed flow'
441  call exit(201)
442  endif
443 !
444  dh=prop(index+3)
445  if(dh.eq.0.d0) then
446  dh=dsqrt(4*a2/pi)
447  endif
448  if(inv.eq.0) then
449  reynolds=5000.d0
450  else
451  reynolds=dabs(xflow)*dh/(dvi*a2)
452  endif
453 !
454  if(inv.ne.0) then
455  xk=zeta/(a2*a2)
456  else
457  xkn=zeta/(a2*a2)
458  xkp=xkn
459  endif
460  elseif(lakon(nelem)(4:5).eq.'EX') then
461 !
462 ! pipe, exit (Idelchik)
463 !
464  a1=prop(index+1)
465  a2=prop(index+2)
466  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
467  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
468  if(inv.lt.0) then
469  write(*,*) '*ERROR in liquidpipe: loss coefficients'
470  write(*,*) ' for exit (Idelchik) do not apply to'
471  write(*,*) ' reversed flow'
472  call exit(201)
473  endif
474 !
475  dh=prop(index+3)
476  if(dh.eq.0.d0) then
477  dh=dsqrt(4*a1/pi)
478  endif
479  if(inv.eq.0) then
480  reynolds=5000.d0
481  else
482  reynolds=dabs(xflow)*dh/(dvi*a1)
483  endif
484 !
485  if(inv.ne.0) then
486  xk=zeta/(a1*a1)
487  else
488  xkn=zeta/(a1*a1)
489  xkp=xkn
490  endif
491  elseif(lakon(nelem)(4:5).eq.'US') then
492 !
493 ! pipe, user defined loss coefficient
494 !
495  a1=prop(index+1)
496  a2=prop(index+2)
497  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
498  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
499  if(inv.lt.0) then
500  write(*,*) '*ERROR in liquidpipe: loss coefficients'
501  write(*,*) ' for a user element do not apply to'
502  write(*,*) ' reversed flow'
503  call exit(201)
504  endif
505  if(a1.lt.a2) then
506  a=a1
507  a2=a1
508  else
509  a=a2
510  a1=a2
511  endif
512 !
513  dh=prop(index+3)
514  if(dh.eq.0.d0) then
515  dh=dsqrt(4*a/pi)
516  endif
517  if(inv.eq.0) then
518  reynolds=5000.d0
519  else
520  reynolds=dabs(xflow)*dh/(dvi*a)
521  endif
522 !
523  if(inv.ne.0) then
524  xk=zeta/(a*a)
525  else
526  xkn=zeta/(a*a)
527  xkp=xkn
528  endif
529  elseif(lakon(nelem)(6:7).eq.'GV') then
530 !
531 ! pipe, gate valve (Berlamont)
532 !
533  a=prop(index+1)
534  if(nactdog(3,nodem).eq.0) then
535 ! geometry is fixed
536  alpha=prop(index+2)
537  else
538 ! geometry is unknown
539  alpha=v(3,nodem)
540  endif
541  a1=a
542  a2=a
543  dzetadalpha=0.d0
544  call ident(xgv,alpha,ngv,id)
545  if(id.eq.0) then
546  zeta=ygv(1)
547  elseif(id.eq.ngv) then
548  zeta=ygv(ngv)
549  else
550  dzetadalpha=(ygv(id+1)-ygv(id))/(xgv(id+1)-xgv(id))
551  zeta=ygv(id)+dzetadalpha*(alpha-xgv(id))
552  endif
553  if(inv.ne.0) then
554  xk=zeta/(a*a)
555  dkda=dzetadalpha/(a*a)
556  else
557  if(flowunknown) then
558  xkn=zeta/(a*a)
559  xkp=xkn
560  endif
561  endif
562  elseif(lakon(nelem)(6:7).eq.'BE') then
563 !
564 ! pipe, bend; values from Berlamont
565 !
566  a=prop(index+1)
567  rd=prop(index+2)
568  alpha=prop(index+3)
569  coarseness=prop(index+4)
570  a1=a
571  a2=a
572  call ident(xbe,rd,nbe,id)
573  if(id.eq.0) then
574  zeta=ybe(1)+(zbe(1)-ybe(1))*coarseness
575  elseif(id.eq.nbe) then
576  zeta=ybe(nbe)+(zbe(nbe)-ybe(nbe))*coarseness
577  else
578  zeta=(1.d0-coarseness)*
579  & (ybe(id)+(ybe(id+1)-ybe(id))*(rd-xbe(id))/
580  & (xbe(id+1)-xbe(id)))
581  & +coarseness*
582  & (zbe(id)+(zbe(id+1)-zbe(id))*(rd-xbe(id))/
583  & (xbe(id+1)-xbe(id)))
584  endif
585  zeta=zeta*alpha/90.d0
586  if(inv.ne.0) then
587  xk=zeta/(a*a)
588  else
589  xkn=zeta/(a*a)
590  xkp=xkn
591  endif
592  elseif(lakon(nelem)(4:5).eq.'BE') then
593 !
594 ! pipe, bend; values from Idelchik or Miller, OWN
595 !
596  a=prop(index+1)
597  dh=prop(index+3)
598  if(dh.eq.0.d0) then
599  dh=dsqrt(4*a/pi)
600  endif
601  if(inv.eq.0) then
602  reynolds=5000.d0
603  else
604  reynolds=dabs(xflow)*dh/(dvi*a)
605  endif
606  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
607  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
608  if(inv.ne.0) then
609  xk=zeta/(a*a)
610  else
611  xkn=zeta/(a*a)
612  xkp=xkn
613  endif
614  a1=a
615  a2=a
616  elseif(lakon(nelem)(4:5).eq.'LO') then
617 !
618 ! long orifice; values from Idelchik or Lichtarowicz
619 !
620  a1=prop(index+1)
621  dh=prop(index+3)
622  if(inv.eq.0) then
623  reynolds=5000.d0
624  else
625  reynolds=dabs(xflow)*dh/(dvi*a1)
626  endif
627  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
628  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
629  if(inv.ne.0) then
630  xk=zeta/(a1*a1)
631  else
632  xkn=zeta/(a1*a1)
633  xkp=xkn
634  endif
635  a2=a1
636  elseif(lakon(nelem)(4:5).eq.'WA') then
637 !
638 ! wall orifice; values from Idelchik
639 !
640 ! entrance is infinitely large
641 !
642  a1=1.d10*prop(index+1)
643 !
644 ! reduced cross section
645 !
646  a2=prop(index+2)
647  dh=prop(index+3)
648  if(inv.eq.0) then
649  reynolds=5000.d0
650  else
651  reynolds=dabs(xflow)*dh/(dvi*a2)
652  endif
653  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
654  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
655 !
656 ! check for negative flow: in that case the loss
657 ! coefficient is wrong
658 !
659  if(inv.lt.0) then
660  write(*,*) '*ERROR in liquidpipe: loss coefficients'
661  write(*,*) ' for wall orifice do not apply to'
662  write(*,*) ' reversed flow'
663  call exit(201)
664  endif
665  if(inv.ne.0) then
666  xk=zeta/(a2*a2)
667  else
668  xkn=zeta/(a2*a2)
669  xkp=xkn
670  endif
671  elseif(lakon(nelem)(4:5).eq.'BR') then
672 !
673 ! branches (joints and splits); values from Idelchik and GE
674 !
675  if(nelem.eq.nint(prop(index+2))) then
676  a=prop(index+5)
677  else
678  a=prop(index+6)
679  endif
680  a1=a
681  a2=a
682 !
683 ! check for negative flow: in that case the loss
684 ! coefficient is wroing
685 !
686  if(inv.lt.0) then
687  write(*,*) '*ERROR in liquidpipe: loss coefficients'
688  write(*,*) ' for branches do not apply to'
689  write(*,*) ' reversed flow'
690  call exit(201)
691  endif
692  if(inv.ne.0) then
693  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
694  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
695  xk=zeta/(a*a)
696  else
697 !
698 ! here, the flow is unknown. To this end zeta is needed. However,
699 ! zeta depends on the flow: circular argument. Therefore a
700 ! fixed initial value for zeta is taken
701 !
702  zeta=0.5d0
703  xkn=zeta/(a*a)
704  xkp=xkn
705  endif
706 !
707 ! all types of orifices
708 !
709  elseif((lakon(nelem)(4:5).eq.'C1')) then
710  a1=prop(index+1)
711  a2=a1
712  dh=prop(index+2)
713  if(inv.eq.0) then
714  reynolds=5000.d0
715  else
716  reynolds=dabs(xflow)*dh/(dvi*a1)
717  endif
718  zeta=1.d0
719 !
720  a=a1
721  zeta=1/zeta**2
722  if(inv.ne.0) then
723  xk=zeta/(a*a)
724  else
725  xkn=zeta/(a*a)
726  xkp=xkn
727  endif
728 !
729 ! all types of vorticies
730 !
731  elseif((lakon(nelem)(4:4).eq.'V')) then
732 !
733 ! radius downstream
734  r2d=prop(index+1)
735 !
736 ! radius upstream
737  r1d=prop(index+2)
738 !
739 ! pressure correction factor
740  eta=prop(index+3)
741 !
742  if(((xflow.gt.0.d0).and.(r2d.gt.r1d))
743  & .or.((r2.lt.r1).and.(xflow.lt.0d0))) then
744  inv=1.d0
745  p1=v(2,node1)
746  p2=v(2,node2)
747  r1=r1d
748  r2=r2d
749 !
750  elseif(((xflow.gt.0.d0).and.(r2d.lt.r1d))
751  & .or.((r2.gt.r1).and.(xflow.lt.0d0))) then
752  inv=-1.d0
753  r1=r2d
754  r2=r1d
755  p1=v(2,node2)
756  p2=v(2,node1)
757  xflow=-v(1,nodem)*iaxial
758 !
759  nodef(1)=node2
760  nodef(2)=nodem
761  nodef(3)=node1
762 !
763  endif
764 !
765  idirf(1)=2
766  idirf(2)=1
767  idirf(3)=2
768 !
769 ! FREE VORTEX
770 !
771  if((lakon(nelem)(4:5).eq.'VF')) then
772 ! rotation induced loss (correction factor)
773  k1= prop(index+4)
774 !
775 ! tangential velocity of the disk at vortex entry
776  u1=prop(index+5)
777 !
778 ! number of the element generating the upstream swirl
779  nelemswirl=nint(prop(index+6))
780 !
781 ! rotation speed (revolution per minutes)
782  omega=prop(index+7)
783 !
784 ! Temperature change
785  t_chang=prop(index+8)
786 !
787  if(omega.gt.0) then
788 !
789 ! rotation speed is given if the swirl comes from a rotating part
790 ! typically the blade of a coverplate
791 !
792 ! C_u is given by radius r1d (see definition of the flow direction)
793 ! C_u related to radius r2d is a function of r1d
794 !
795  if(inv.gt.0) then
796  c1u=omega*r1
797 !
798 ! flow rotation at outlet
799  c2u=c1u*r1/r2
800 !
801  elseif(inv.lt.0) then
802  c2u=omega*r2
803 !
804  c1u=c2u*r2/r1
805  endif
806 !
807  elseif(nelemswirl.gt.0) then
808  if(lakon(nelemswirl)(2:5).eq.'LPPN') then
809  cinput=prop(ielprop(nelemswirl)+5)
810  elseif(lakon(nelemswirl)(2:5).eq.'LPVF') then
811  cinput=prop(ielprop(nelemswirl)+9)
812  elseif(lakon(nelemswirl)(2:5).eq.'LPFS') then
813  cinput=prop(ielprop(nelemswirl)+7)
814  endif
815 !
816  cinput=u1+k1*(cinput-u1)
817 !
818  if(inv.gt.0) then
819  c1u=cinput
820  c2u=c1u*r1/r2
821  elseif(inv.lt.0) then
822  c2u=cinput
823  c1u=c2u*r2/r1
824  endif
825  endif
826 ! storing the tengential velocity for later use (wirbel cascade)
827  if(inv.gt.0) then
828  prop(index+9)=c2u
829  elseif(inv.lt.0) then
830  prop(index+9)=c1u
831  endif
832 !
833 ! inner rotation
834 !
835  if(r1.lt.r2) then
836  ciu=c1u
837  elseif(r1.ge.r2) then
838  ciu=c2u
839  endif
840 !
841 ! if (iflag.eq.1) then
842  a1=1e-6
843  a2=a1
844  if(inv.ne.0) then
845  xkn=rho/2*ciu**2*(1-(r1/r2)**2)
846  xkp=xkn
847  else
848  xkn=rho/2*ciu**2*(1-(r1/r2)**2)
849  xkp=xkn
850  endif
851  endif
852 !
853 ! FORCED VORTEX
854 !
855  if((lakon(nelem)(4:5).eq.'VS')) then
856 !
857 ! core swirl ratio
858  kr=prop(index+4)
859 !
860 ! rotation speed (revolution per minutes) of the rotating part
861 ! responsible for the swirl
862  omega=prop(index+5)
863 !
864 ! Temperature change
865  t_chang=prop(index+6)
866 !
867  ui=omega*r1
868  c1u=ui*kr
869  c2u=c1u*r2/r1
870 !
871 ! storing the tengential velocity for later use (wirbel cascade)
872  if(inv.gt.0) then
873  prop(index+7)=c2u
874  elseif(inv.lt.0) then
875  prop(index+7)=c1u
876  endif
877 !
878  a1=1e-6
879  a2=a1
880  if(iflag.eq.1)then
881  xflow=0.5d0
882  endif
883 !
884  if(inv.ne.0) then
885  xkn=rho/2*ui**2*((r2/r1)**2-1)
886  xkp=xkn
887  else
888  xkn=rho/2*ui**2*((r2/r1)**2-1)
889  xkp=xkn
890  endif
891  endif
892  endif
893 !
894  if(iflag.eq.1) then
895  if(flowunknown) then
896 !
897  xk1=1.d0/(a1*a1)
898  xk2=1.d0/(a2*a2)
899  xflow=(z1-z2+(p1-p2)/rho)/(xk2-xk1+xkp)
900  if(xflow.lt.0.d0) then
901  xflow=(z1-z2+(p1-p2)/rho)/(xk2-xk1-xkn)
902  if(xflow.lt.0.d0) then
903  write(*,*) '*WARNING in liquidpipe:'
904  write(*,*) ' initial mass flow could'
905  write(*,*) ' not be determined'
906  write(*,*) ' 1.d-10 is taken'
907  xflow=1.d-10
908  else
909  xflow=-rho*dsqrt(2.d0*xflow)
910  endif
911  else
912  xflow=rho*dsqrt(2.d0*xflow)
913  endif
914  else
915 !
916 ! mass flow known, geometry unknown
917 !
918  if(lakon(nelem)(6:7).eq.'GV') then
919  prop(index+2)=0.5d0
920  endif
921  endif
922  elseif(iflag.eq.2) then
923  xk1=1.d0/(a1*a1)
924  xk2=1.d0/(a2*a2)
925 !
926  if(lakon(nelem)(4:4).ne.'V') then
927 !
928  numf=4
929  df(3)=1.d0/rho
930  df(1)=-df(3)
931  df(2)=(xk2-xk1+inv*xk)*xflow/(rho*rho)
932  df(4)=(xflow*xflow*inv*dkda)/(2.d0*rho*rho)
933  f=df(3)*p2+df(1)*p1+df(2)*xflow/2.d0+z2-z1
934 !
935  else if (lakon(nelem)(4:5).eq.'VF') then
936  numf=3
937  if(r2.ge.r1) then
938  f=p1-p2+xkp
939  df(1)=1
940  df(2)=0
941  df(3)=-1
942  elseif(r2.lt.r1) then
943  f=p1-p2-xkp
944  df(1)=1
945  df(2)=0
946  df(3)=-1
947  endif
948  else if (lakon(nelem)(4:5).eq.'VS') then
949  if(((r2.ge.r1).and.(xflow.gt.0d0))
950  & .or.((r2.lt.r1).and.(xflow.lt.0d0)))then
951 !
952  f=p1-p2+xkn
953 ! pressure node1
954  df(1)=1
955 ! massflow nodem
956  df(2)=0
957 ! pressure node2
958  df(3)=-1
959 !
960  elseif(((r2.lt.r1).and.(xflow.gt.0d0))
961  & .or.((r2.gt.r1).and.(xflow.lt.0d0)))then
962 !
963  f=p2-p1+xkn
964 ! pressure node1
965  df(1)=-1
966 ! massflow nodem
967  df(2)=0
968 ! pressure node2
969  df(3)=1
970  endif
971  endif
972 !
973  else if (iflag.eq.3) then
974  xflow_vol=xflow/rho
975  un=dvi/rho
976  if(inv.eq.1) then
977  t=v(0,node1)
978  else
979  t=v(0,node2)
980  endif
981 !
982  write(1,*) ''
983  write(1,55) ' from node',node1,
984  & ' to node', node2,': oil massflow rate = ',xflow,
985  & ' i.e. in volume per time ',xflow_vol
986  55 FORMAT(1x,a,i6,a,i6,a,e11.4,a,e11.4,a)
987  write(1,57)'
988  &Rho= ',rho,', Nu= ',un,', dyn.visc.= ',dvi
989 
990  if(inv.eq.1) then
991  write(1,56)' Inlet node ',node1,': Tt1=',t,
992  & ', Pt1=',p1
993  if(lakon(nelem)(4:5).eq.'EL'.or.
994  & lakon(nelem)(4:5).eq.'CO'.or.
995  & lakon(nelem)(4:5).eq.'EN'.or.
996  & lakon(nelem)(4:5).eq.'EX'.or.
997  & lakon(nelem)(4:5).eq.'US'.or.
998  & lakon(nelem)(4:5).eq.'BE'.or.
999  & lakon(nelem)(4:5).eq.'LO'.or.
1000  & lakon(nelem)(4:5).eq.'WA'.or.
1001  & lakon(nelem)(4:5).eq.'BR')then
1002 
1003  write(1,*)' Element ',nelem,lakon(nelem)
1004  write(1,58)' Re= ',reynolds,' zeta= ',
1005  & zeta
1006 !
1007  elseif((lakon(nelem)(4:5).eq.'C1')) then
1008  write(1,*)' Element ',nelem,lakon(nelem)
1009  write(1,58)' Re= ',reynolds,' cd= ',
1010  & zeta
1011 !
1012  else if(lakon(nelem)(4:5).eq.'FR')then
1013  write(1,*)' Element ',nelem,lakon(nelem)
1014  write(1,59)' Re= ',reynolds,' lambda=
1015  &',friction,' lambda*L/D= ',friction*dl/d
1016 !
1017  else if (lakon(nelem)(4:4).eq.'V')then
1018  write(1,*)' Element ',nelem,lakon(nelem)
1019  write(1,*)' C1u= ',c1u,'m/s ,C2u= '
1020  &,c2u,'m/s',' ,DeltaP= ',xkn
1021  endif
1022 !
1023  write(1,56)' Outlet node ',node2,': Tt2=',t,
1024  & ', Pt2=',p2
1025 !
1026  else if(inv.eq.-1) then
1027 
1028  endif
1029 !
1030  56 FORMAT(1x,a,i6,a,e11.4,a,e11.4,a)
1031  57 FORMAT(1x,a,e11.4,a,e11.4,a,e11.4,a)
1032  58 FORMAT(1x,a,e11.4,a,e11.4)
1033  59 FORMAT(1x,a,e11.4,a,e11.4,a,e11.4)
1034  endif
1035 !
1036  endif
1037 !
1038  xflow=xflow/iaxial
1039  df(2)=df(2)*iaxial
1040 !
1041  return
ident
subroutine ident(x, px, n, id)
Definition: ident.f:26
z1
static double * z1
Definition: filtermain.c:48
zeta_calc
subroutine zeta_calc(nelem, prop, ielprop, lakon, reynolds, zeta, isothermal, kon, ipkon, R, kappa, v, mi, iaxial)
Definition: zeta_calc.f:35
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
r1
static double * r1
Definition: filtermain.c:48
friction_coefficient
subroutine friction_coefficient(l, d, ks, reynolds, form_fact, lambda)
Definition: friction_coefficient.f:26
e11
static double * e11
Definition: radflowload.c:42
Hosted by OpenAircraft.com, (Michigan UAV, LLC)