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

Go to the source code of this file.

Functions/Subroutines

subroutine vortex (node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, cp, R, numf, set, mi, ttime, time, iaxial)
 

Function/Subroutine Documentation

◆ vortex()

subroutine vortex ( integer, intent(out)  node1,
integer, intent(out)  node2,
integer, intent(in)  nodem,
integer, intent(in)  nelem,
character*8, dimension(*), intent(in)  lakon,
integer, dimension(*), intent(in)  kon,
integer, dimension(*), intent(in)  ipkon,
integer, dimension(0:3,*), intent(in)  nactdog,
logical, intent(out)  identity,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(out)  prop,
integer, intent(in)  iflag,
real*8, dimension(0:mi(2),*), intent(in)  v,
real*8, intent(out)  xflow,
real*8, intent(out)  f,
integer, dimension(*), intent(out)  nodef,
integer, dimension(*), intent(out)  idirf,
real*8, dimension(*), intent(out)  df,
real*8, intent(in)  cp,
real*8, intent(in)  R,
integer, intent(out)  numf,
character*81, dimension(*), intent(in)  set,
integer, dimension(*), intent(in)  mi,
real*8, intent(in)  ttime,
real*8, intent(in)  time,
integer, intent(in)  iaxial 
)
22 !
23 ! vortex element
24 !
25 ! author: Yannick Muller
26 !
27  implicit none
28 !
29  logical identity
30  character*8 lakon(*)
31  character*81 set(*)
32 !
33  integer nelem,nactdog(0:3,*),node1,node2,nodem,numf,
34  & ielprop(*),nodef(*),idirf(*),index,iflag,iaxial,
35  & inv,ipkon(*),kon(*),t_chang,nelemswirl,mi(*)
36 !
37  real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,r,cp,
38  & p1,p2,t1,t2,km1,pi,ttime,time,r2d,r1d,eta,u1,
39  & c1u,c2u,cinput,r1,r2,omega,k1,ciu,expon,
40  & ui,kr,cte1,cte2,qred_crit,a,xflow_oil
41 !
42  intent(in) nodem,nelem,lakon,kon,ipkon,
43  & nactdog,ielprop,iflag,v,
44  & cp,r,set,mi,ttime,time,iaxial
45 !
46  intent(out) identity,node1,node2,xflow,numf,nodef,idirf,prop,
47  & f,df
48 !
49  if(iflag.eq.0) then
50  identity=.true.
51 !
52  if(nactdog(2,node1).ne.0)then
53  identity=.false.
54  elseif(nactdog(2,node2).ne.0)then
55  identity=.false.
56  elseif(nactdog(1,nodem).ne.0)then
57  identity=.false.
58  endif
59 !
60  elseif(iflag.eq.1)then
61 !
62  xflow=0.d0
63 !
64  elseif(iflag.eq.2)then
65 !
66  numf=4
67  index=ielprop(nelem)
68  kappa=(cp/(cp-r))
69  km1=kappa-1
70  pi=4.d0*datan(1.d0)
71 !
72 ! radius downstream
73  r2d=prop(index+1)
74 !
75 ! radius upstream
76  r1d=prop(index+2)
77 !
78 ! pressure correction factor
79  eta=prop(index+3)
80 !
81  p1=v(2,node1)
82  p2=v(2,node2)
83 !
84  xflow=v(1,nodem)*iaxial
85 !
86  if(xflow.gt.0.d0) then
87  inv=1.d0
88  p1=v(2,node1)
89  p2=v(2,node2)
90  t1=v(0,node1)
91  t2=v(0,node2)
92  r1=r1d
93  r2=r2d
94 !
95  nodef(1)=node1
96  nodef(2)=node1
97  nodef(3)=nodem
98  nodef(4)=node2
99 !
100  elseif(xflow.lt.0.d0) then
101  inv=-1.d0
102  r1=r2d
103  r2=r1d
104  p1=v(2,node2)
105  p2=v(2,node1)
106  t1=v(0,node2)
107  t2=v(0,node1)
108  xflow=-v(1,nodem)*iaxial
109 !
110  nodef(1)=node2
111  nodef(2)=node2
112  nodef(3)=nodem
113  nodef(4)=node1
114 !
115  endif
116 !
117  idirf(1)=2
118  idirf(2)=0
119  idirf(3)=1
120  idirf(4)=2
121 !
122  kappa=(cp/(cp-r))
123 !
124 ! FREE VORTEX
125 !
126  if(lakon(nelem)(4:5).eq.'FR')then
127 !
128 ! rotation induced loss (correction factor)
129  k1= prop(index+4)
130 !
131 ! tangential velocity of the disk at vortex entry
132  u1=prop(index+5)
133 !
134 ! number of the element generating the upstream swirl
135  nelemswirl=nint(prop(index+6))
136 !
137 ! rotation speed (revolution per minutes)
138  omega=prop(index+7)
139 !
140 ! Temperature change
141  t_chang=prop(index+8)
142 !
143  if(omega.gt.0.d0) then
144 !
145 ! rotation speed is given if the swirl comes from a rotating part
146 ! typically the blade of a coverplate
147 !
148 
149 ! C_u is given by radius r1d (see definition of the flow direction)
150 ! C_u related to radius r2d is a function of r1d
151 !
152  if(inv.gt.0) then
153  c1u=omega*r1
154 !
155 ! flow rotation at outlet
156  c2u=c1u*r1/r2
157 !
158  elseif(inv.lt.0) then
159  c2u=omega*r2
160 !
161  c1u=c2u*r2/r1
162  endif
163 !
164  elseif(nelemswirl.gt.0) then
165 ! preswirl nozzle
166  if(lakon(nelemswirl)(2:5).eq.'ORPN') then
167  cinput=prop(ielprop(nelemswirl)+5)
168 ! rotating orifices
169  else if((lakon(nelemswirl)(2:5).eq.'ORMM').or.
170  & (lakon(nelemswirl)(2:5).eq.'ORMA').or.
171  & (lakon(nelemswirl)(2:5).eq.'ORPM').or.
172  & (lakon(nelemswirl)(2:5).eq.'ORPA')) then
173  cinput=prop(ielprop(nelemswirl)+7)
174 ! forced vortex
175  elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
176  cinput=prop(ielprop(nelemswirl)+7)
177 ! free vortex
178  elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
179  cinput=prop(ielprop(nelemswirl)+9)
180 ! Moehring
181  elseif(lakon(nelemswirl)(2:4).eq.'MRG') then
182  cinput=prop(ielprop(nelemswirl)+10)
183 ! RCAVO
184  elseif((lakon(nelemswirl)(2:4).eq.'ROR').or.
185  & (lakon(nelemswirl)(2:4).eq.'ROA'))then
186  cinput=prop(ielprop(nelemswirl)+6)
187 ! RCAVI
188  elseif(lakon(nelemswirl)(2:4).eq.'RCV') then
189  cinput=prop(ielprop(nelemswirl)+5)
190  else
191  write(*,*) '*ERROR in vortex:'
192  write(*,*) ' element',nelemswirl
193  write(*,*) ' referred by element',nelem
194  write(*,*) ' is not a swirl generating element'
195  cinput=0.d0
196  endif
197 !
198  cinput=u1+k1*(cinput-u1)
199 !
200  if(inv.gt.0) then
201  c1u=cinput
202  c2u=c1u*r1/r2
203  elseif(inv.lt.0) then
204  c2u=cinput
205  c1u=c2u*r2/r1
206  endif
207  endif
208 !
209 ! storing the tengential velocity for later use (wirbel cascade)
210  if(inv.gt.0) then
211  prop(index+9)=c2u
212  elseif(inv.lt.0) then
213  prop(index+9)=c1u
214  endif
215 !
216 ! inner rotation
217 !
218  if(r1.lt.r2) then
219  ciu=c1u
220  elseif(r1.ge.r2) then
221  ciu=c2u
222  endif
223 !
224  expon=kappa/km1
225 !
226  if(r2.ge.r1) then
227 !
228  cte1=c1u**2/(2*cp*t1)
229  cte2=1-(r1/r2)**2
230 
231  f=p2/p1-1d0-eta*((1+cte1*cte2)**expon-1d0)
232 !
233  df(1)=-p2/p1**2
234 !
235  df(2)=eta*expon*cte1/t1*cte2*
236  & (1+cte1*cte2)**(expon-1)
237 !
238  df(3)=0
239 !
240  df(4)=1/p1
241 !
242  elseif(r2.lt.r1) then
243 !
244  cte1=c2u**2/(2*cp*t2)
245  cte2=1-(r2/r1)**2
246 !
247  f=p1/p2-1d0-eta*((1+cte1*cte2)**expon-1d0)
248 !
249  df(1)=1/p2
250 !
251  df(2)=eta*expon*cte1/t1*cte2*
252  & (1+cte1*cte2)**(expon-1)
253 !
254  df(3)=0
255 !
256  df(4)=-p1/p2**2
257 !
258  endif
259 !
260 ! FORCED VORTEX
261 !
262  elseif(lakon(nelem)(4:5).eq.'FO') then
263 !
264 ! core swirl ratio
265  kr=prop(index+4)
266 !
267 ! rotation speed (revolution per minutes) of the rotating part
268 ! responsible for the swirl
269  omega=prop(index+5)
270 !
271 ! Temperature change
272  t_chang=prop(index+6)
273 !
274  if(r2.ge.r1) then
275  ui=omega*r1
276  c1u=ui*kr
277  c2u=c1u*r2/r1
278  elseif(r2.lt.r1) then
279  ui=omega*r2
280  c2u=ui*kr
281  c1u=c2u*r1/r2
282  endif
283 !
284 ! storing the tengential velocity for later use (wirbel cascade)
285  if(inv.gt.0) then
286  prop(index+7)=c2u
287  elseif(inv.lt.0) then
288  prop(index+7)=c1u
289  endif
290 !
291  expon=kappa/km1
292 !
293  if(((r2.ge.r1).and.(xflow.gt.0d0))
294  & .or.((r2.lt.r1).and.(xflow.lt.0d0)))then
295 !
296  cte1=(c1u)**2/(2*cp*t1)
297  cte2=(r2/r1)**2-1
298 !
299  f=p2/p1-1-eta*((1+cte1*cte2)**expon-1)
300 !
301 ! pressure node1
302  df(1)=-p2/p1**2
303 !
304 ! temperature node1
305  df(2)=eta*expon*cte1/t1*cte2*(1+cte1*cte2)**(expon-1)
306 !
307 ! massflow nodem
308  df(3)=0
309 !
310 ! pressure node2
311  df(4)=1/p1
312 !
313  elseif(((r2.lt.r1).and.(xflow.gt.0d0))
314  & .or.((r2.gt.r1).and.(xflow.lt.0d0)))then
315  cte1=(c2u)**2/(2*cp*t2)
316  cte2=(r1/r2)**2-1
317 !
318  f=p1/p2-1-eta*((1+cte1*cte2)**expon-1)
319 !
320 ! pressure node1
321  df(1)=1/p2
322 !
323 ! temperature node1
324  df(2)=eta*expon*cte1/t2*cte2*(1+cte1*cte2)**(expon-1)
325 !
326 ! massflow nodem
327  df(3)=0
328 !
329 ! pressure node2
330  df(4)=-p1/p2**2
331 !
332  endif
333  endif
334 !
335 ! outpout
336 !
337  elseif(iflag.eq.3) then
338 !
339  index=ielprop(nelem)
340  kappa=(cp/(cp-r))
341  km1=kappa-1
342  pi=4.d0*datan(1.d0)
343 !
344 ! radius downstream
345  r2d=prop(index+1)
346 !
347 ! radius upstream
348  r1d=prop(index+2)
349 !
350 ! pressure correction factor
351  eta=prop(index+3)
352 !
353  p1=v(2,node1)
354  p2=v(2,node2)
355 !
356  xflow=v(1,nodem)*iaxial
357 !
358  if(xflow.gt.0.d0) then
359  inv=1.d0
360  p1=v(2,node1)
361  p2=v(2,node2)
362  t1=v(0,node1)
363  t2=v(0,node2)
364  r1=r1d
365  r2=r2d
366 !
367  nodef(1)=node1
368  nodef(2)=node1
369  nodef(3)=nodem
370  nodef(4)=node2
371 !
372  elseif(xflow.lt.0.d0) then
373  inv=-1.d0
374  r1=r2d
375  r2=r1d
376  p1=v(2,node2)
377  p2=v(2,node1)
378  t1=v(0,node2)
379  t2=v(0,node1)
380  xflow=v(1,nodem)*iaxial
381 !
382  nodef(1)=node2
383  nodef(2)=node2
384  nodef(3)=nodem
385  nodef(4)=node1
386 !
387  endif
388 !
389  idirf(1)=2
390  idirf(2)=0
391  idirf(3)=1
392  idirf(4)=2
393 !
394  kappa=(cp/(cp-r))
395 !
396 ! FREE VORTEX
397 !
398  if(lakon(nelem)(4:5).eq.'FR')then
399 !
400 ! rotation induced loss (correction factor)
401  k1= prop(index+4)
402 !
403 ! tengential velocity of the disk at vortex entry
404  u1=prop(index+5)
405 !
406 ! number of the element generating the upstream swirl
407  nelemswirl=nint(prop(index+6))
408 !
409 ! rotation speed (revolution per minutes)
410  omega=prop(index+7)
411 !
412 ! Temperature change
413  t_chang=prop(index+8)
414 !
415  if(omega.gt.0.d0) then
416 !
417 ! rotation speed is given if the swirl comes from a rotating part
418 ! typically the blade of a coverplate
419 !
420 ! C_u is given by radius r1d (see definition of the flow direction)
421 ! C_u related to radius r2d is a function of r1d
422 !
423  if(inv.gt.0) then
424  c1u=omega*r1
425 !
426 ! flow rotation at outlet
427  c2u=c1u*r1/r2
428 !
429  elseif(inv.lt.0) then
430  c2u=omega*r2
431 !
432  c1u=c2u*r2/r1
433  endif
434 !
435  elseif(nelemswirl.gt.0) then
436 ! swirl generating element
437 !
438 ! preswirl nozzle
439  if(lakon(nelemswirl)(2:5).eq.'ORPN') then
440  cinput=prop(ielprop(nelemswirl)+5)
441 ! rotating orifices
442  elseif((lakon(nelemswirl)(2:5).eq.'ORMM').or.
443  & (lakon(nelemswirl)(2:5).eq.'ORMA').or.
444  & (lakon(nelemswirl)(2:5).eq.'ORPM').or.
445  & (lakon(nelemswirl)(2:5).eq.'ORPA')) then
446  cinput=prop(ielprop(nelemswirl)+7)
447 ! forced vortex
448  elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
449  cinput=prop(ielprop(nelemswirl)+7)
450 ! free vortex
451  elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
452  cinput=prop(ielprop(nelemswirl)+9)
453 ! Moehring
454  elseif(lakon(nelemswirl)(2:4).eq.'MRG') then
455  cinput=prop(ielprop(nelemswirl)+10)
456 ! RCAVO
457  elseif((lakon(nelemswirl)(2:4).eq.'ROR').or.
458  & (lakon(nelemswirl)(2:4).eq.'ROA'))then
459  cinput=prop(ielprop(nelemswirl)+6)
460 ! RCAVI
461  elseif(lakon(nelemswirl)(2:4).eq.'RCV') then
462  cinput=prop(ielprop(nelemswirl)+5)
463  else
464  write(*,*) '*ERROR in vortex:'
465  write(*,*) ' element',nelemswirl
466  write(*,*) ' referred by element',nelem
467  write(*,*) ' is not a swirl generating element'
468  cinput=0.d0
469  endif
470 !
471  cinput=u1+k1*(cinput-u1)
472 !
473  if(inv.gt.0) then
474  c1u=cinput
475  c2u=c1u*r1/r2
476  elseif(inv.lt.0) then
477  c2u=cinput
478  c1u=c2u*r2/r1
479  endif
480  endif
481 !
482 ! storing the tengential velocity for later use (wirbel cascade)
483  if(inv.gt.0) then
484  prop(index+9)=c2u
485  elseif(inv.lt.0) then
486  prop(index+9)=c1u
487  endif
488 !
489 ! inner rotation
490 !
491  if(r1.lt.r2) then
492  ciu=c1u
493  elseif(r1.ge.r2) then
494  ciu=c2u
495  endif
496 !
497  expon=kappa/km1
498 !
499  if(r2.ge.r1) then
500 !
501  cte1=c1u**2/(2*cp*t1)
502  cte2=1-(r1/r2)**2
503 
504  f=p2/p1-1d0-eta*((1+cte1*cte2)**expon-1d0)
505 !
506  df(1)=-p2/p1**2
507 !
508  df(2)=eta*expon*cte1/t1*cte2*
509  & (1+cte1*cte2)**(expon-1)
510 !
511  df(3)=0
512 !
513  df(4)=1/p1
514 !
515  elseif(r2.lt.r1) then
516 !
517  cte1=c2u**2/(2*cp*t2)
518  cte2=1-(r2/r1)**2
519 !
520  f=p1/p2-1d0-eta*((1+cte1*cte2)**expon-1d0)
521 !
522  df(1)=1/p2
523 !
524  df(2)=eta*expon*cte1/t1*cte2*
525  & (1+cte1*cte2)**(expon-1)
526 !
527  df(3)=0
528 !
529  df(4)=-p1/p2**2
530 !
531  endif
532 !
533 ! FORCED VORTEX
534 !
535  elseif(lakon(nelem)(4:5).eq.'FO') then
536 !
537 ! core swirl ratio
538  kr=prop(index+4)
539 !
540 ! rotation speed (revolution per minutes) of the rotating part
541 ! responsible for the swirl
542  omega=prop(index+5)
543 !
544 ! Temperature change
545  t_chang=prop(index+6)
546 !
547 ! no element generating the upstream swirl
548  nelemswirl=0
549 !
550  if(r2.ge.r1) then
551  ui=omega*r1
552  c1u=ui*kr
553  c2u=c1u*r2/r1
554  elseif(r2.lt.r1) then
555  ui=omega*r2
556  c2u=ui*kr
557  c1u=c2u*r1/r2
558  endif
559 !
560 ! storing the tengential velocity for later use (wirbel cascade)
561  if(inv.gt.0) then
562  prop(index+7)=c2u
563  elseif(inv.lt.0) then
564  prop(index+7)=c1u
565  endif
566 !
567  expon=kappa/km1
568  endif
569 !
570  xflow_oil=0.d0
571 !
572  write(1,*) ''
573  write(1,55) ' from node ',node1,
574  &' to node ', node2,' : air massflow rate = ',xflow,' ',
575  &' , oil massflow rate = ',xflow_oil,' '
576 !
577  if(inv.eq.1) then
578  write(1,56)' Inlet node ',node1,' : Tt1 = ',t1,
579  & ' , Ts1 = ',t1,' , Pt1 = ',p1
580  write(1,*)' Element ',nelem,lakon(nelem)
581  write(1,57)' C1u = ',c1u,
582  &' , C2u = ',c2u
583  write(1,56)' Outlet node ',node2,' : Tt2 = ',t2,
584  & ' , Ts2 = ',t2,' , Pt2 = ',p2
585 !
586  else if(inv.eq.-1) then
587  write(1,56)' Inlet node ',node2,': Tt1 = ',t1,
588  & ' , Ts1 = ',t1,' , Pt1 = ',p1
589  write(1,*)' Element ',nelem,lakon(nelem)
590  write(1,57)' C1u = ',c1u,
591  &' , C2u = ',c2u
592  write(1,56)' Outlet node ',node1,' Tt2 = ',
593  & t2,' , Ts2 = ',t2,' , Pt2 = ',p2
594  endif
595  endif
596 !
597  55 format(1x,a,i6,a,i6,a,e11.4,a,a,e11.4,a)
598  56 format(1x,a,i6,a,e11.4,a,e11.4,a,e11.4,a,e11.4)
599  57 format(1x,a,e11.4,a,e11.4,a)
600 !
601  xflow=xflow/iaxial
602  df(3)=df(3)*iaxial
603 !
604  return
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
r1
static double * r1
Definition: filtermain.c:48
e11
static double * e11
Definition: radflowload.c:42
Hosted by OpenAircraft.com, (Michigan UAV, LLC)