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

Go to the source code of this file.

Functions/Subroutines

subroutine restrictor (node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, cp, r, physcon, dvi, numf, set, shcon, nshcon, rhcon, nrhcon, ntmat_, mi, ttime, time, iaxial, co, vold)
 

Function/Subroutine Documentation

◆ restrictor()

subroutine restrictor ( integer  node1,
integer  node2,
integer  nodem,
integer  nelem,
character*8, dimension(*)  lakon,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
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  cp,
real*8  r,
real*8, dimension(*)  physcon,
real*8  dvi,
integer  numf,
character*81, dimension(*)  set,
real*8, dimension(0:3,ntmat_,*)  shcon,
integer, dimension(*)  nshcon,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer  ntmat_,
integer, dimension(*)  mi,
real*8  ttime,
real*8  time,
integer  iaxial,
real*8, dimension(3,*)  co,
real*8, dimension(0:mi(2),*)  vold 
)
24 !
25 ! pressure loss element with partial total head loss
26 !
27 ! author: Yannick Muller
28 !
29  implicit none
30 !
31  logical identity,crit,isothermal
32  character*8 lakon(*)
33  character*81 set(*)
34 !
35  integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
36  & ielprop(*),nodef(*),idirf(*),index,iflag,iaxial,
37  & inv,ipkon(*),kon(*),kgas,icase,k_oil,nshcon(*),
38  & nrhcon(*),ntmat_,mi(*)
39 !
40  real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,d,
41  & tt1,tt2,pt1,pt2,cp,physcon(*),km1,dvi,co(3,*),
42  & kp1,kdkm1,reynolds,kdkp1,ttime,time,
43  & pt2pt1,pt1pt2,pt1pt2_crit,qred_crit,qred1,qred2,zeta,
44  & a1,a2,root, expon1,expon2,expon3,fact1,fact2,sqrt,pi,
45  & pt2_lim,m2,m1,xflow_oil,t1,t2,phi,vold(0:mi(2),*),
46  & shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),zeta_phi,aeff,
47  & c2,tdkp1
48 !
49  phi=0.d0
50  index=ielprop(nelem)
51 !
52  if(iflag.eq.0) then
53  identity=.true.
54 !
55  if(nactdog(2,node1).ne.0)then
56  identity=.false.
57  elseif(nactdog(2,node2).ne.0)then
58  identity=.false.
59  elseif(nactdog(1,nodem).ne.0)then
60  identity=.false.
61  endif
62 !
63  elseif(iflag.eq.1)then
64 !
65 ! complementing the properties of restrictor elements
66 !
67  if(lakon(nelem)(2:5).eq.'REEX')then
68  prop(index+2)=100000*prop(index+1)
69  elseif(lakon(nelem)(2:7).eq.'REWAOR')then
70  prop(index+1)=100000*prop(index+2)
71  elseif(lakon(nelem)(2:5).eq.'REEN')then
72  prop(index+1)=100000*prop(index+2)
73  prop(index+4)=0.5d0
74  endif
75 !
76  isothermal=.false.
77  kappa=(cp/(cp-r))
78  kp1=kappa+1d0
79  km1=kappa-1d0
80 !
81 ! defining surfaces for branch elements
82 !
83  if(lakon(nelem)(2:6).eq.'REBRJ') then
84  if(nelem.eq.nint(prop(index+2))) then
85  a1=prop(index+5)
86  a2=a1
87  elseif(nelem.eq.nint(prop(index+3)))then
88  a1=prop(index+6)
89  a2=a1
90  endif
91  zeta=1.d0
92  elseif(lakon(nelem)(2:6).eq.'REBRS') then
93  if(nelem.eq.nint(prop(index+2))) then
94  a1=prop(index+5)
95  a2=a1
96  elseif(nelem.eq.nint(prop(index+3)))then
97  a1=prop(index+6)
98  a2=a1
99  endif
100  zeta=1.d0
101  elseif(lakon(nelem)(2:5).eq.'REUS' ) then
102 !
103 ! for other Restrictor elements
104 !
105  a1=prop(index+1)
106  a2=prop(index+2)
107  zeta=prop(index+4)
108 !
109 ! change in flow direction
110 !
111  if(v(2,node1).ge.v(2,node2))then
112  zeta=prop(index+4)
113  else
114  zeta=prop(index+7)
115  if(zeta.le.0.d0)then
116  zeta=prop(index+4)
117  endif
118  endif
119 !
120  if(a1.gt.a2) then
121  a1=a2
122  endif
123  else
124  a1=prop(index+1)
125  a2=prop(index+2)
126  zeta=1.d0
127  endif
128 !
129  pt1=v(2,node1)
130  pt2=v(2,node2)
131 !
132  if(pt1.ge.pt2) then
133  inv=1
134  tt1=v(0,node1)-physcon(1)
135  tt2=v(0,node2)-physcon(1)
136  else
137  inv=-1
138  pt1=v(2,node2)
139  pt2=v(2,node1)
140  tt1=v(0,node2)-physcon(1)
141  tt2=v(0,node1)-physcon(1)
142  endif
143 !
144  pt1pt2=pt1/pt2
145  pt2pt1=1/pt1pt2
146  km1=kappa-1.d0
147  kp1=kappa+1.d0
148  kdkm1=kappa/km1
149 !
150  if(.not.isothermal) then
151  pt1pt2_crit=(0.5d0*kp1)**(zeta*kdkm1)
152  else
153  pt1pt2_crit=0.5d0*(3*kappa-1)**(zeta*kdkm1)
154  endif
155 !
156  if(pt1pt2.gt.pt1pt2_crit) then
157  crit=.true.
158  pt1pt2=pt1pt2_crit
159  endif
160 !
161  if(a1.le.a2) then
162 !
163  qred1=dsqrt(kappa/r)*pt1pt2**(-0.5d0*kp1/(kappa*zeta))
164  & *dsqrt(2.d0/km1*(pt1pt2**(km1/(kappa*zeta))-1d0))
165 !
166  qred2=pt1pt2*a1/a2*qred1
167 !
168  if(.not.isothermal) then
169  qred_crit=dsqrt(kappa/r)*(1.d0+0.5d0*km1)
170  & **(-0.5d0*kp1/km1)
171  else
172  qred_crit=dsqrt(1/r)*(1+0.5*km1/kappa)
173  & **(-0.5d0*kp1/km1)
174  endif
175 !
176  if(qred2.lt.qred_crit) then
177  if((qred1.gt.qred_crit).or.(pt1pt2.gt.pt1pt2_crit)) then
178  xflow=inv*a1*pt1*qred_crit/dsqrt(tt1)
179  else
180  xflow=inv*a1*pt1*qred1/dsqrt(tt1)
181  endif
182  else
183  call pt2_lim_calc(pt1,a2,a1,kappa,zeta,pt2_lim)
184 !
185  xflow=inv*a2*pt2_lim*qred_crit/dsqrt(tt2)
186 !
187  endif
188 !
189  else
190  qred_crit=dsqrt(kappa/r)*(1.d0+0.5d0*km1)
191  & **(-0.5d0*kp1/km1)
192  call pt2_lim_calc(pt1,a2,a1,kappa,zeta,pt2_lim)
193 !
194  xflow=inv*a2*pt2_lim*qred_crit/dsqrt(tt2)
195  endif
196 !
197  pt2pt1=pt2/pt1
198  km1=kappa-1.d0
199  kp1=kappa+1.d0
200  kdkm1=kappa/km1
201  tdkp1=2.d0/kp1
202  c2=tdkp1**kdkm1
203  if(a1.gt.a2) then
204  aeff=a2
205  else
206  aeff=a1
207  endif
208  if(pt2pt1.gt.c2) then
209  xflow=inv*pt1*aeff*dsqrt(2.d0*kdkm1*pt2pt1**(2.d0/kappa)
210  & *(1.d0-pt2pt1**(1.d0/kdkm1))/r)/dsqrt(tt1)
211  else
212  xflow=inv*pt1*aeff*dsqrt(kappa/r)*tdkp1**(kp1/(2.d0*km1))/
213  & dsqrt(tt1)
214  endif
215  if(lakon(nelem)(2:5).ne.'RECO') then
216  xflow=0.75*xflow
217  else
218  xflow=xflow
219  endif
220 !
221  elseif(iflag.eq.2)then
222 !
223 ! complementing the properties of restrictor elements
224 !
225  if(lakon(nelem)(2:5).eq.'REEX')then
226  prop(index+2)=100000*prop(index+1)
227  elseif(lakon(nelem)(2:7).eq.'REWAOR')then
228  prop(index+1)=100000*prop(index+2)
229  elseif(lakon(nelem)(2:5).eq.'REEN')then
230  prop(index+1)=100000*prop(index+2)
231  prop(index+4)=0.5d0
232  endif
233 !
234  numf=4
235  isothermal=.false.
236  pi=4.d0*datan(1.d0)
237  kappa=(cp/(cp-r))
238  km1=kappa-1.d0
239  kp1=kappa+1.d0
240  kdkm1=kappa/km1
241  kdkp1=kappa/kp1
242 !
243  pt1=v(2,node1)
244  pt2=v(2,node2)
245 !
246  if(pt1.ge.pt2) then
247  inv=1
248  else
249  inv=-1
250  endif
251 !
252 ! defining surfaces and oil properties for branches elements
253 !
254  if(lakon(nelem)(2:6).eq.'REBRJ') then
255  if(nelem.eq.nint(prop(index+2))) then
256  a1=prop(index+5)
257  a2=a1
258  xflow_oil=prop(index+9)
259  k_oil=nint(prop(index+11))
260  elseif(nelem.eq.nint(prop(index+3)))then
261  a1=prop(index+6)
262  a2=a1
263  xflow_oil=prop(index+10)
264  k_oil=nint(prop(index+11))
265  endif
266  elseif(lakon(nelem)(2:6).eq.'REBRS') then
267  if(nelem.eq.nint(prop(index+2))) then
268  a1=prop(index+5)
269  a2=a1
270  if(lakon(nelem)(2:8).eq.'REBRSI1') then
271  xflow_oil=prop(index+11)
272  k_oil=nint(prop(index+13))
273  else
274  xflow_oil=prop(index+9)
275  k_oil=nint(prop(index+11))
276  endif
277  elseif(nelem.eq.nint(prop(index+3)))then
278  a1=prop(index+6)
279  a2=a1
280  if(lakon(nelem)(2:8).eq.'REBRSI1') then
281  xflow_oil=prop(index+12)
282  k_oil=nint(prop(index+13))
283  else
284  xflow_oil=prop(index+10)
285  k_oil=nint(prop(index+11))
286  endif
287  endif
288  else
289 !
290 ! for other Restrictor elements
291 !
292  if(inv.gt.0.d0) then
293  a1=prop(index+1)
294  a2=prop(index+2)
295  else
296  a1=prop(index+2)
297  a2=prop(index+1)
298  endif
299 !
300  if(lakon(nelem)(2:5).eq.'REEL') then
301  xflow_oil=prop(index+4)
302  k_oil=nint(prop(index+5))
303  elseif((lakon(nelem)(2:7).eq.'RELOID').or.
304  & (lakon(nelem)(2:5).eq.'REUS').or.
305  & (lakon(nelem)(2:5).eq.'REEN').or.
306  & (lakon(nelem)(2:5).eq.'REEX').or.
307  & (lakon(nelem)(2:7).eq.'REWAOR').or.
308  & (lakon(nelem)(2:7).eq.'RELOLI')) then
309  xflow_oil=prop(index+5)
310  k_oil=nint(prop(index+6))
311  elseif((lakon(nelem)(2:5).eq.'RECO').or.
312  & (lakon(nelem)(2:7).eq.'REBEMA').or.
313  & (lakon(nelem)(2:7).eq.'REBEMI').or.
314  & (lakon(nelem)(2:8).eq.'REBEIDC')) then
315  xflow_oil=prop(index+6)
316  k_oil=nint(prop(index+7))
317  elseif(lakon(nelem)(2:8).eq.'REBEIDR') then
318  xflow_oil=prop(index+8)
319  k_oil=nint(prop(index+9))
320  endif
321  endif
322 !
323  if(pt1.gt.pt2) then
324  inv=1
325  xflow=v(1,nodem)*iaxial
326  tt1=v(0,node1)-physcon(1)
327  tt2=v(0,node2)-physcon(1)
328 !
329  icase=0
330  call ts_calc(xflow,tt1,pt1,kappa,r,a1,t1,icase)
331  call ts_calc(xflow,tt2,pt2,kappa,r,a2,t2,icase)
332 !
333  nodef(1)=node1
334  nodef(2)=node1
335  nodef(3)=nodem
336  nodef(4)=node2
337 !
338  elseif(pt1.eq.pt2) then
339  inv=1
340  xflow=v(1,nodem)*iaxial
341  tt1=v(0,node1)-physcon(1)
342  tt2=v(0,node2)-physcon(1)
343 !
344  pt2=pt2-0.01*pt2
345  icase=0
346  call ts_calc(xflow,tt1,pt1,kappa,r,a1,t1,icase)
347  call ts_calc(xflow,tt2,pt2,kappa,r,a2,t2,icase)
348 !
349  nodef(1)=node1
350  nodef(2)=node1
351  nodef(3)=nodem
352  nodef(4)=node2
353 !
354  else
355  inv=-1
356  pt1=v(2,node2)
357  pt2=v(2,node1)
358  xflow=-v(1,nodem)*iaxial
359  tt1=v(0,node2)-physcon(1)
360  tt2=v(0,node1)-physcon(1)
361  icase=0
362  call ts_calc(xflow,tt1,pt1,kappa,r,a1,t1,icase)
363  call ts_calc(xflow,tt2,pt2,kappa,r,a2,t2,icase)
364  nodef(1)=node2
365  nodef(2)=node2
366  nodef(3)=nodem
367  nodef(4)=node1
368  endif
369 !
370  idirf(1)=2
371  idirf(2)=0
372  idirf(3)=1
373  idirf(4)=2
374 !
375 ! calculation of the dynamic viscosity
376 !
377  if(lakon(nelem)(2:3).eq.'RE') then
378  icase=0
379  endif
380 !
381  if(dabs(dvi).lt.1e-30) then
382  write(*,*) '*ERROR in restrictor: '
383  write(*,*) ' no dynamic viscosity defined'
384  write(*,*) ' dvi= ',dvi
385  call exit(201)
386  endif
387 !
388 ! Reynolds number calculation
389 !
390  if(lakon(nelem)(2:5).eq.'REBR') then
391  d=dsqrt(4d0*a1/pi)
392  reynolds=dabs(xflow)*d/(dvi*a1)
393  else
394  d=prop(index+3)
395  if(a1.le.a2) then
396  reynolds=dabs(xflow)*d/(dvi*a1)
397  else
398  reynolds=dabs(xflow)*d/(dvi*a2)
399  endif
400  endif
401 
402  if(xflow_oil.lt.1.d-10) then
403  xflow_oil=0.d0
404  endif
405  if(lakon(nelem)(2:7).eq.'REBEMI') then
406 !
407 ! BEND MILLER with oil
408 !
409  if(xflow_oil.ne.0.d0) then
410 !
411  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
412  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
413  & v,dvi,cp,r,k_oil,phi,zeta,nshcon,nrhcon,
414  & shcon,rhcon,ntmat_,mi,iaxial)
415 !
416  zeta=phi*zeta
417  else
418 
419  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
420  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
421  phi=1.d0
422  endif
423  elseif(lakon(nelem)(2:7).eq.'RELOID') then
424 !
425 ! long orifice idelchick with oil
426 !
427  if(xflow_oil.ne.0.d0) then
428 !
429  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
430  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
431  & v,dvi,cp,r,k_oil,phi,zeta,nshcon,nrhcon,
432  & shcon,rhcon,ntmat_,mi,iaxial)
433  zeta=phi*zeta
434 
435  else
436 
437  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
438  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
439  phi=1.d0
440  endif
441 !
442  else
443 !
444 ! every other zeta elements with/without oil
445 !
446  if(xflow_oil.ne.0.d0) then
447  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
448  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
449  & v,dvi,cp,r,k_oil,phi,zeta,nshcon,nrhcon,
450  & shcon,rhcon,ntmat_,mi,iaxial)
451  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
452  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
453  zeta=phi*zeta
454  else
455  phi=1.d0
456  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
457  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
458  zeta=phi*zeta
459  endif
460  endif
461 !
462  if(dabs(zeta).lt.1.d-10) zeta=0.d0
463 !
464  if(zeta.lt.0.d0) then
465  pt1=v(2,node1)
466  pt2=v(2,node2)
467  xflow=v(1,nodem)*iaxial
468  tt2=v(0,node2)
469  tt1=v(0,node1)
470  call ts_calc(xflow,tt1,pt1,kappa,r,a1,t1,icase)
471  call ts_calc(xflow,tt2,pt2,kappa,r,a2,t2,icase)
472 !
473  nodef(1)=node1
474  nodef(2)=node1
475  nodef(3)=nodem
476  nodef(4)=node2
477 !
478  endif
479 !
480  if(.not.isothermal) then
481  pt1pt2_crit=(0.5d0*kp1)**(zeta*kdkm1)
482  else
483  pt1pt2_crit=0.5d0*(3*kappa-1)**(zeta*kdkm1)
484  endif
485  pt1pt2=pt1/pt2
486 !
487 ! Mach number caclulation
488 !
489  m1=dsqrt(2d0/km1*(tt1/t1-1d0))
490  if((1.d0-m1).le.1.d-6) then
491  if(zeta.gt.0.d0) then
492  call limit_case_calc(a2,pt1,tt2,xflow,zeta,r,kappa,
493  & pt2_lim,m2)
494 !
495  endif
496  else
497  m2=dsqrt(2d0/km1*(tt2/t2-1d0))
498  endif
499 !
500 ! Section A1 smaller than or equal to section A2
501 ! or for all BRANCHES ELEMENTS
502 !
503  if(a1.le.a2) then
504 !
505 ! definition of the reduced mass flows
506 !
507  if(zeta.gt.0.d0) then
508 !
509  qred1=dsqrt(kappa/r)*pt1pt2**(-0.5d0*kp1/(kappa*zeta))
510  & *dsqrt(2.d0/km1*(pt1pt2**(km1/(kappa*zeta))-1d0))
511 !
512  elseif(zeta.lt.0.d0) then
513 !
514  qred1=dabs(xflow)*dsqrt(tt1)/(pt1*a1)
515 !
516  endif
517 !
518  qred2=pt1pt2*a1/a2*qred1
519 !
520  if(.not.isothermal) then
521  qred_crit=dsqrt(kappa/r)*(1.d0+0.5d0*km1)
522  & **(-0.5d0*kp1/km1)
523  else
524  qred_crit=dsqrt(1/r)*(1+0.5*km1/kappa)
525  & **(-0.5d0*kp1/km1)
526  endif
527 !
528 ! icase zeta greater than zero
529 !
530  if(zeta.gt.0.d0) then
531 !
532 ! definition of the coefficients
533 !
534  if((pt1pt2.ge.0.9999999999d0).and.
535  & (pt1pt2.le.1.000000000111111d0))then
536  pt1=1.0001d0*pt2
537  pt1pt2=pt1/pt2
538  endif
539 !
540  sqrt=dsqrt(r*tt1/kappa)
541  expon1=-0.5d0*kp1/(zeta*kappa)
542  fact1=pt1pt2**expon1
543  expon2=km1/(zeta*kappa)
544  fact2=pt1pt2**expon2
545  expon3=1d0/(zeta*kappa)
546  root=2d0/km1*(fact2-1d0)
547 !
548  if(qred2.lt.qred_crit) then
549 !
550  if((qred1.lt.qred_crit)
551  & .and.(pt1pt2.lt.pt1pt2_crit))then
552 !
553 ! section 1 is not critical
554 !
555 ! residual
556 !
557  f=xflow*sqrt/(a1*pt1)-fact1*dsqrt(root)
558 !
559 ! pressure node1
560 !
561  df(1)=-xflow*sqrt/(a1*pt1**2)+
562  & fact1/pt1*dsqrt(root)
563  & *(-expon1-expon3*fact2/root)
564 !
565 ! temperature node1
566 !
567  df(2)=0.5d0*xflow*dsqrt(r/(kappa*tt1))/(a1*pt1)
568 !
569 ! mass flow
570 !
571  df(3)=inv*sqrt/(a1*pt1)
572 !
573 ! pressure node2
574 !
575  df(4)=fact1/pt2*dsqrt(root)*
576  & (expon1+expon3*fact2/root)
577 !
578  else
579 !
580 ! section1 is critical
581 !
582  f=xflow*sqrt/(pt1*a1)-dsqrt(r/kappa)*qred_crit
583 !
584 ! pressure node1
585 !
586  df(1)=-xflow*sqrt/(a1*pt1**2)
587 !
588 ! temperature node1
589 !
590  df(2)=0.5d0*xflow*dsqrt(r/kappa)
591  & /(pt1*a1*dsqrt(tt1))
592 !
593 ! mass flow
594 !
595  df(3)=inv*sqrt/(a1*pt1)
596 !
597 ! pressure node2
598 !
599  df(4)=0.d0
600 !
601  endif
602 !
603  else
604 !
605 ! section A2 critical
606 !
607  call pt2_lim_calc(pt1,a2,a1,kappa,zeta,pt2_lim)
608  pt1pt2=pt1/pt2_lim
609 !
610  fact1=pt1pt2**expon1
611 !
612  fact2=pt1pt2**expon2
613 !
614  root=2d0/km1*(fact2-1d0)
615 !
616  f=xflow*sqrt/(a1*pt1)-fact1*dsqrt(root)
617 !
618 ! pressure node1
619 !
620  df(1)=-xflow*sqrt/(a1*pt1**2)+
621  & fact1/pt1*dsqrt(root)
622  & *(-expon1-expon3*fact2/root)
623 !
624 ! temperature node1
625 !
626  df(2)=0.5d0*xflow*dsqrt(r/(kappa*tt1))/(a1*pt1)
627 !
628 ! mass flow
629 !
630  df(3)=inv*sqrt/(a1*pt1)
631 !
632 ! pressure node2
633 !
634  df(4)=0
635 !
636  endif
637 !
638 ! icase zeta less than zero
639 !
640  elseif(zeta.lt.0.d0) then
641 !
642  expon1=-kp1/(zeta*kappa)
643  fact1=pt1pt2**expon1
644  expon2=km1/(zeta*kappa)
645  fact2=pt1pt2**expon2
646  expon3=1d0/(zeta*kappa)
647  root=2d0/km1*(fact2-1d0)
648 !
649  if(qred1.lt.qred_crit) then
650 !
651 ! section 1 is not critical
652 !
653 ! residual
654 !
655  f=xflow**2*r*tt1/(a1**2*pt1**2*kappa)
656  & -fact1*root
657 !
658 ! pressure node1
659 !
660  df(1)=-2*xflow**2*r*tt1/(a1**2*pt1**3*kappa)
661  & -1/pt1*fact1*(expon1*root
662  & +2/(zeta*kappa)*fact2)
663 !
664 ! temperature node1
665 !
666  df(2)=xflow**2*r/(a1**2*pt1**2*kappa)
667 !
668 ! mass flow
669 !
670  df(3)=2*xflow*r*tt1/(a1**2*pt1**2*kappa)
671 !
672 ! pressure node2
673 !
674  df(4)=-(1/pt2*fact1)
675  & *(-expon1*root-2/(zeta*kappa)*fact2)
676 !
677  else
678 !
679 ! section1 is critical
680 !
681  f=xflow**2*r*tt1/(a1**2*pt1**2*kappa)
682  & -r/kappa*qred_crit**2
683 !
684 ! pressure node1
685 !
686  df(1)=-2*xflow**2*r*tt1/(a1**2*pt1**3*kappa)
687 !
688 ! temperature node1
689 !
690  df(2)=xflow**2*r/(a1**2*pt1**2*kappa)
691 !
692 ! mass flow
693 !
694  df(3)=2*xflow*r*tt1/(a1**2*pt1**2*kappa)
695 !
696 ! pressure node2
697 !
698  df(4)=0.d0
699 !
700  endif
701 !
702  elseif(zeta.eq.0.d0) then
703 !
704 ! zeta = 0
705 !
706  f=pt1/pt2-1.d0
707 !
708  df(1)=1/pt2
709  df(2)=0.d0
710  df(3)=0.d0
711  df(4)=-pt1/pt2**2
712 !
713  endif
714 !
715  else
716 !
717 ! A1 greater than A2
718 !
719  qred2=dabs(xflow)*dsqrt(tt2)/(a2*pt2)
720 !
721  qred1=1/pt1pt2*a2/a1*qred2
722 !
723  qred_crit=dsqrt(kappa/r)*(1.d0+0.5d0*km1)
724  & **(-0.5d0*kp1/km1)
725 !
726 ! definition of the coefficients
727 !
728  if(zeta.gt.0.d0) then
729 !
730  sqrt=dsqrt(r*tt1/kappa)
731 !
732  if((pt1pt2.ge.0.9999999999d0).and.
733  & (pt1pt2.le.1.000000000111111d0))then
734  pt1=1.0001d0*pt2
735  pt1pt2=pt1/pt2
736  endif
737 !
738  expon1=-0.5d0*kp1/(zeta*kappa)
739  fact1=pt1pt2**expon1
740  expon2=km1/(zeta*kappa)
741  fact2=pt1pt2**expon2
742  expon3=1d0/(zeta*kappa)
743  root=2d0/km1*(fact2-1d0)
744 !
745  if(pt1pt2.ge.pt1pt2_crit) then
746  pt1pt2=pt1pt2_crit
747  pt2=pt1/pt1pt2_crit
748  endif
749 !
750  if((qred2.lt.qred_crit)
751  & .and.(pt1/pt2.lt.pt1pt2_crit)) then
752 !
753 ! section 2 is not critical
754 !
755 ! residual
756 !
757  f=xflow*sqrt/(a2*pt2)-fact1*dsqrt(root)
758 !
759 ! pressure node1
760 !
761  df(1)=-fact1/pt1*dsqrt(root)
762  & *(expon1+0.5*dsqrt(2/km1)*expon2*fact2/root)
763 !
764 ! temperature node1
765 !
766  df(2)=0.5d0*xflow*sqrt/(a2*pt2*tt1)
767 !
768 ! mass flow
769 !
770  df(3)=inv*sqrt/(a2*pt2)
771 !
772 ! pressure node2
773 !
774  df(4)=-xflow*sqrt/(a2*pt2**2)
775  & -fact1/pt2*dsqrt(root)*
776  & (-expon1-0.5*dsqrt(2/km1)*expon2*fact2/root)
777 !
778  else
779  write(*,*)
780  & '*WARNING in restrictor: A1 greater than A2'
781  write(*,*) ' critical flow in element',nelem
782 !
783 ! section2 is critical
784 !
785  pt2=pt1/pt1pt2_crit
786 !
787  f=xflow*dsqrt(tt1)/(pt2*a2)-qred_crit
788 !
789 ! pressure node1
790 !
791  df(1)=0
792 !
793 ! temperature node1
794 !
795  df(2)=0.5d0*xflow/(a2*pt2*dsqrt(tt2))
796 !
797 ! mass flow
798 !
799  df(3)=inv*dsqrt(tt1)/(a2*pt2)
800 !
801 ! pressure node2
802 !
803  df(4)=-xflow*dsqrt(tt1)/(a2*pt2**2)
804 !
805  endif
806 !
807  elseif(zeta.eq.0.d0) then
808 !
809  qred1=dabs(xflow)*dsqrt(tt1*kappa/r)/(a1*pt1)
810  qred2=dabs(xflow)*dsqrt(tt2*kappa/r)/(a2*pt2)
811  qred_crit=dsqrt(kappa/r)*(1.d0+0.5d0*km1)
812  & **(-0.5d0*kp1/km1)
813 !
814  f=pt1/pt2-1.d0
815 !
816  df(1)=1/pt2
817 !
818  df(2)=0.d0
819 !
820  df(3)=0.d0
821 !
822  df(4)=-pt1/pt2**2
823 !
824  endif
825  endif
826 !
827  elseif((iflag.eq.3).or.(iflag.eq.4)) then
828 !
829 ! complementing the properties of restrictor elements
830 !
831  if(lakon(nelem)(2:5).eq.'REEX')then
832  prop(index+2)=100000*prop(index+1)
833  elseif(lakon(nelem)(2:7).eq.'REWAOR')then
834  prop(index+1)=100000*prop(index+2)
835  elseif(lakon(nelem)(2:5).eq.'REEN')then
836  prop(index+1)=100000*prop(index+2)
837  prop(index+4)=0.5d0
838  endif
839 !
840  isothermal=.false.
841  pi=4.d0*datan(1.d0)
842  kappa=(cp/(cp-r))
843  km1=kappa-1.d0
844  kp1=kappa+1.d0
845  kdkm1=kappa/km1
846  kdkp1=kappa/kp1
847 !
848  pt1=v(2,node1)
849  pt2=v(2,node2)
850  if(pt1.ge.pt2) then
851  inv=1
852  else
853  inv=-1
854  endif
855 !
856 ! defining surfaces for branches elements
857 !
858  if(lakon(nelem)(2:6).eq.'REBRJ') then
859  if(nelem.eq.nint(prop(index+2))) then
860  a1=prop(index+5)
861  a2=a1
862  xflow_oil=prop(index+9)
863  k_oil=nint(prop(index+11))
864  elseif(nelem.eq.nint(prop(index+3)))then
865  a1=prop(index+6)
866  a2=a1
867  xflow_oil=prop(index+10)
868  k_oil=nint(prop(index+11))
869  endif
870  elseif(lakon(nelem)(2:6).eq.'REBRS') then
871  if(nelem.eq.nint(prop(index+2))) then
872  a1=prop(index+5)
873  a2=a1
874  if(lakon(nelem)(2:8).eq.'REBRSI1') then
875  xflow_oil=prop(index+11)
876  k_oil=nint(prop(index+13))
877  else
878  xflow_oil=prop(index+9)
879  k_oil=nint(prop(index+11))
880  endif
881  elseif(nelem.eq.nint(prop(index+3)))then
882  a1=prop(index+6)
883  a2=a1
884  if(lakon(nelem)(2:8).eq.'REBRSI1') then
885  xflow_oil=prop(index+12)
886  k_oil=nint(prop(index+13))
887  else
888  xflow_oil=prop(index+10)
889  k_oil=nint(prop(index+11))
890  endif
891  endif
892  else
893 !
894 ! for other Restrictor elements
895 !
896  a1=prop(index+1)
897  a2=prop(index+2)
898  if(lakon(nelem)(2:5).eq.'REEL') then
899  xflow_oil=prop(index+4)
900  k_oil=nint(prop(index+5))
901  elseif((lakon(nelem)(2:7).eq.'RELOID').or.
902  & (lakon(nelem)(2:5).eq.'REUS').or.
903  & (lakon(nelem)(2:5).eq.'REEN').or.
904  & (lakon(nelem)(2:5).eq.'REEX').or.
905  & (lakon(nelem)(2:7).eq.'REWAOR').or.
906  & (lakon(nelem)(2:7).eq.'RELOLI')) then
907  xflow_oil=prop(index+5)
908  k_oil=nint(prop(index+6))
909  elseif((lakon(nelem)(2:5).eq.'RECO').or.
910  & (lakon(nelem)(2:7).eq.'REBEMA').or.
911  & (lakon(nelem)(2:7).eq.'REBEMI').or.
912  & (lakon(nelem)(2:8).eq.'REBEIDC')) then
913  xflow_oil=prop(index+6)
914  k_oil=nint(prop(index+7))
915  elseif(lakon(nelem)(2:8).eq.'REBEIDR') then
916  xflow_oil=prop(index+8)
917  k_oil=nint(prop(index+9))
918  endif
919  endif
920 !
921  if(pt1.ge.pt2) then
922  inv=1
923  xflow=v(1,nodem)*iaxial
924  tt1=v(0,node1)-physcon(1)
925  tt2=v(0,node2)-physcon(1)
926  icase=0
927  call ts_calc(xflow,tt1,pt1,kappa,r,a1,t1,icase)
928  call ts_calc(xflow,tt2,pt2,kappa,r,a2,t2,icase)
929 !
930  else
931  inv=-1
932  pt1=v(2,node2)
933  pt2=v(2,node1)
934  xflow=-v(1,nodem)*iaxial
935  tt1=v(0,node2)-physcon(1)
936  tt2=v(0,node1)-physcon(1)
937  icase=0
938  call ts_calc(xflow,tt1,pt1,kappa,r,a1,t1,icase)
939  call ts_calc(xflow,tt2,pt2,kappa,r,a2,t2,icase)
940 !
941  endif
942 !
943 ! calculation of the dynamic viscosity
944 !
945  if(lakon(nelem)(2:3).eq.'RE') then
946  icase=0
947  elseif(lakon(nelem)(2:5).eq.'REEX') then
948  if(lakon(nint(prop(index+4)))(2:6).eq.'GAPFA') then
949  icase=0
950  elseif(lakon(nint(prop(index+4)))(2:6).eq.'GAPFI') then
951  icase=1
952  endif
953  endif
954 !
955  if(dabs(dvi).lt.1e-30) then
956  write(*,*) '*ERROR in restrictor: '
957  write(*,*) ' no dynamic viscosity defined'
958  write(*,*) ' dvi= ',dvi
959  call exit(201)
960  endif
961 !
962 ! Reynolds number calculation
963 !
964  if(lakon(nelem)(2:5).eq.'REBR') then
965  d=dsqrt(4d0*a1/pi)
966  reynolds=dabs(xflow)*d/(dvi*a1)
967  else
968  d=prop(index+3)
969  if(a1.le.a2) then
970  reynolds=dabs(xflow)*d/(dvi*a1)
971  else
972  reynolds=dabs(xflow)*d/(dvi*a2)
973  endif
974  endif
975 !
976  if(xflow_oil.lt.1.d-10) then
977  xflow_oil=0.d0
978  endif
979 !
980 ! BEND MILLER with oil
981 !
982  if(lakon(nelem)(2:7).eq.'REBEMI') then
983  if(xflow_oil.ne.0.d0) then
984  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
985  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
986  & v,dvi,cp,r,k_oil,phi,zeta,nshcon,nrhcon,
987  & shcon,rhcon,ntmat_,mi,iaxial)
988 
989  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
990  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
991 !
992  zeta_phi=phi*zeta
993  else
994 !
995  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
996  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
997  phi=1.d0
998  zeta_phi=phi*zeta
999 !
1000  endif
1001  elseif(lakon(nelem)(2:7).eq.'RELOID') then
1002 !
1003 ! long orifice in a wall with oil after Idelchik
1004 !
1005  if(xflow_oil.ne.0.d0) then
1006  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
1007  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
1008  & v,dvi,cp,r,k_oil,phi,zeta,nshcon,nrhcon,
1009  & shcon,rhcon,ntmat_,mi,iaxial)
1010 !
1011  zeta_phi=phi*zeta
1012  else
1013 !
1014  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
1015  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
1016  phi=1.d0
1017  zeta_phi=phi*zeta
1018  endif
1019 !
1020  else
1021 !
1022 ! every other zeta elements with/without oil
1023 !
1024  if(xflow_oil.ne.0.d0) then
1025  call two_phase_flow(tt1,pt1,t1,pt2,xflow,
1026  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
1027  & v,dvi,cp,r,k_oil,phi,zeta,nshcon,nrhcon,
1028  & shcon,rhcon,ntmat_,mi,iaxial)
1029 !
1030  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
1031  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
1032 !
1033  zeta_phi=phi*zeta
1034  else
1035  phi=1.d0
1036  call zeta_calc(nelem,prop,ielprop,lakon,reynolds,zeta,
1037  & isothermal,kon,ipkon,r,kappa,v,mi,iaxial)
1038  zeta_phi=phi*zeta
1039  endif
1040  endif
1041 !
1042  if(zeta.le.0.d0) then
1043  pt1=v(2,node1)
1044  pt2=v(2,node2)
1045  xflow=v(1,nodem)*iaxial
1046  tt1=v(0,node1)
1047  tt2=v(0,node2)
1048 !
1049  endif
1050 !
1051  if(.not.isothermal) then
1052  pt1pt2_crit=(0.5d0*kp1)**(zeta*kdkm1)
1053  else
1054  pt1pt2_crit=0.5d0*(3*kappa-1)**(zeta*kdkm1)
1055  endif
1056  pt1pt2=pt1/pt2
1057 !
1058 ! Mach number calculation
1059 !
1060  m1=dsqrt(2d0/km1*(tt1/t1-1d0))
1061  if((1.d0-m1).le.1e-3) then
1062 
1063  if(zeta.gt.0.d0)then
1064  if(xflow_oil.eq.0.d0) then
1065  call limit_case_calc(a2,pt1,tt2,xflow,zeta,r,kappa,
1066  & pt2_lim,m2)
1067  else
1068  call limit_case_calc(a2,pt1,tt2,xflow,zeta_phi,r,kappa
1069  & ,pt2_lim,m2)
1070  endif
1071  endif
1072  else
1073  m2=dsqrt(2d0/km1*(tt2/t2-1d0))
1074  endif
1075 !
1076  if(iflag.eq.3) then
1077  write(1,*) ''
1078  write(1,55) ' from node ',node1,
1079  & ' to node ', node2,' : air massflow rate = '
1080  & ,xflow,' '
1081  & , ' , oil massflow rate = ',xflow_oil,' '
1082 !
1083  if(lakon(nelem)(4:5).ne.'BR') then
1084 !
1085 ! for restrictors
1086 !
1087  if(inv.eq.1) then
1088  write(1,56)' Inlet node ',node1,' : Tt1= ',tt1,
1089  & ' , Ts1 = ',t1,' , Pt1 = ',pt1,
1090  & ' , M1 = ',m1
1091  write(1,*)' Element ',nelem,lakon(nelem)
1092  write(1,57)' dyn.visc. = ',dvi,' , Re = '
1093  & ,reynolds
1094  write(1,58)' PHI = ',phi,' , ZETA = ',zeta,
1095  & ' , ZETA_PHI = ',phi*zeta
1096  write(1,56)' Outlet node ',node2,' : Tt2 = ',
1097  & tt2,
1098  & ' , Ts2 = ',t2,' , Pt2 = ',pt2,
1099  & ' , M2= ',m2
1100 !
1101  elseif(inv.eq.-1) then
1102  write(1,56)' Inlet node ',node2,' : Tt1= ',tt1,
1103  & ' , Ts1= ',t1,' , Pt1= ',pt1,
1104  & ' , M1= ',m1
1105  write(1,*)' Element ',nelem,lakon(nelem)
1106  write(1,57)' dyn.visc. =',dvi,' , Re = '
1107  & ,reynolds
1108  write(1,58)' PHI = ',phi,' , ZETA = ',zeta,
1109  & ' , ZETA_PHI = ',phi*zeta
1110  write(1,56)' Outlet node ',node1,' : Tt2 = ',
1111  & tt2,
1112  & ' , Ts2 = ',t2,' , Pt2 = ',pt2,
1113  & ' , M2 = ',m2
1114  endif
1115  else
1116 !
1117 ! for branches
1118 !
1119  if(inv.eq.1) then
1120  write(1,56)' Inlet node ',node1,' : Tt1 = ',
1121  & tt1,
1122  & ' , Ts1 = ',t1,' , Pt1 = ',pt1,
1123  & ' , M1 = ',m1
1124  write(1,*)' Element ',nelem,lakon(nelem)
1125  write(1,57)' dyn.visc. = ',dvi,' , Re = '
1126  & ,reynolds
1127  write(1,58)' PHI = ',phi,' , ZETA = ',zeta,
1128  & ' , ZETA_PHI = ',phi*zeta
1129  write(1,56)' Outlet node ',node2,' : Tt2 = ',
1130  & tt2,
1131  & ' , Ts2 = ',t2,' , Pt2 = ',pt2,
1132  & ' , M2 = ',m2
1133 !
1134  elseif(inv.eq.-1) then
1135  write(1,56)' Inlet node ',node2,' : Tt1 = ',
1136  & tt1,
1137  & ', Ts1 = ',t1,' , Pt1 = ',pt1,
1138  & ' , M1 = ',m1
1139  write(1,*)' Element ',nelem,lakon(nelem)
1140  write(1,57)' dyn.visc. =',dvi,' , Re = '
1141  & ,reynolds
1142  write(1,58)' PHI = ',phi,' , ZETA = ',zeta,
1143  &' , ZETA_PHI = ',phi*zeta
1144  write(1,56)' Outlet node ',node1,' : Tt2 = ',
1145  & tt2,
1146  & ' , Ts2 = ',t2,' , Pt2 = ',pt2,
1147  & ' , M2 = ',m2
1148  endif
1149  endif
1150 !
1151  endif
1152  endif
1153 !
1154  55 format(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a)
1155  56 format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
1156  57 format(1x,a,e11.4,a,e11.4)
1157  58 format(1x,a,e11.4,a,e11.4,a,e11.4)
1158 !
1159  xflow=xflow/iaxial
1160  df(3)=df(3)*iaxial
1161 !
1162  return
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
ts_calc
subroutine ts_calc(xflow, Tt, Pt, kappa, r, A, Ts, icase)
Definition: ts_calc.f:20
pt2_lim_calc
subroutine pt2_lim_calc(pt1, a2, a1, kappa, zeta, pt2_lim)
Definition: pt2_lim_calc.f:25
two_phase_flow
subroutine two_phase_flow(Tt1, pt1, T1, pt2, xflow_air, xflow_oil, nelem, lakon, kon, ipkon, ielprop, prop, v, dvi_air, cp, r, k_oil, phi, lambda, nshcon, nrhcon, shcon, rhcon, ntmat_, mi, iaxial)
Definition: two_phase_flow.f:23
e11
static double * e11
Definition: radflowload.c:42
limit_case_calc
subroutine limit_case_calc(a2, pt1, Tt2, xflow, zeta, r, kappa, pt2_lim, M2)
Definition: limit_case_calc.f:21
Hosted by OpenAircraft.com, (Michigan UAV, LLC)