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

Go to the source code of this file.

Functions/Subroutines

subroutine initialnet (itg, ieg, ntg, ac, bc, lakon, v, ipkon, kon, nflow, ikboun, nboun, prop, ielprop, nactdog, ndirboun, nodeboun, xbounact, ielmat, ntmat_, shcon, nshcon, physcon, ipiv, nteq, rhcon, nrhcon, ipobody, ibody, xbodyact, co, nbody, network, iin_abs, vold, set, istep, iit, mi, ineighe, ilboun, channel, iaxial, nmpc, labmpc, ipompc, nodempc, coefmpc, ttime, time, iponoel, inoel)
 

Function/Subroutine Documentation

◆ initialnet()

subroutine initialnet ( integer, dimension(*)  itg,
integer, dimension(*)  ieg,
integer  ntg,
real*8, dimension(nteq,nteq)  ac,
real*8, dimension(nteq)  bc,
character*8, dimension(*)  lakon,
real*8, dimension(0:mi(2),*)  v,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
integer  nflow,
integer, dimension(*)  ikboun,
integer  nboun,
real*8, dimension(*)  prop,
integer, dimension(*)  ielprop,
integer, dimension(0:3,*)  nactdog,
integer, dimension(*)  ndirboun,
integer, dimension(*)  nodeboun,
real*8, dimension(*)  xbounact,
integer, dimension(mi(3),*)  ielmat,
integer  ntmat_,
real*8, dimension(0:3,ntmat_,*)  shcon,
integer, dimension(*)  nshcon,
real*8, dimension(*)  physcon,
integer, dimension(*)  ipiv,
integer  nteq,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer, dimension(2,*)  ipobody,
integer, dimension(3,*)  ibody,
real*8, dimension(7,*)  xbodyact,
real*8, dimension(3,*)  co,
integer  nbody,
integer  network,
integer  iin_abs,
real*8, dimension(0:mi(2),*)  vold,
character*81, dimension(*)  set,
integer  istep,
integer  iit,
integer, dimension(*)  mi,
integer, dimension(*)  ineighe,
integer, dimension(*)  ilboun,
integer  channel,
integer  iaxial,
integer  nmpc,
character*20, dimension(*)  labmpc,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
real*8  ttime,
real*8  time,
integer, dimension(*)  iponoel,
integer, dimension(2,*)  inoel 
)
34 !
35  implicit none
36 
37  logical identity,gravity,gaspipe
38 !
39  character*8 lakon(*)
40  character*20 labmpc(*)
41  character*81 set(*)
42 !
43  integer mi(*),ieg(*),nflow,i,j,ntg,ielmat(mi(3),*),ntmat_,
44  & node1,node2,ider,iaxial,nmpc,ipompc(*),nodempc(3,*),
45  & nelem,index,nshcon(*),ipkon(*),kon(*),ikboun(*),nboun,idof,
46  & nodem,idirf(8),nactdog(0:3,*),imat,ielprop(*),id1,id2,
47  & nodef(8),ndirboun(*),nodeboun(*),itg(*),node,kflag,ipiv(*),
48  & nrhs,info,idof1,idof2,nteq,nrhcon(*),ipobody(2,*),ibody(3,*),
49  & nbody,numf,network,iin_abs,icase,index2,index1,nelem1,nelem2,
50  & node11,node21,node12,node22,istep,iit,ineighe(*),iponoel(*),
51  & ilboun(*),idir,channel,inoel(2,*),indexe,iel
52 !
53  real*8 ac(nteq,nteq), bc(nteq),prop(*),shcon(0:3,ntmat_,*),
54  & f,df(8),xflow,xbounact(*),v(0:mi(2),*),cp,r,tg1,
55  & tg2,gastemp,physcon(*),pressmin,dvi,rho,g(3),z1,z2,
56  & rhcon(0:1,ntmat_,*),co(3,*),xbodyact(7,*),kappa,
57  & a,tt,pt,ts,pressmax,constant,vold(0:mi(2),*),
58  & coefmpc(*),ttime,time,xflow360
59 !
60  kflag=1
61  ider=0
62  channel=0
63  gaspipe=.false.
64 !
65 ! applying the boundary conditions
66 !
67  do j=1,nboun
68  v(ndirboun(j),nodeboun(j))=xbounact(j)
69  enddo
70 !
71 ! check for channel elements
72 !
73  do i=1,nflow
74  nelem=ieg(i)
75  if(lakon(nelem)(2:5).eq.'LICH') then
76  channel=1
77  return
78  endif
79  enddo
80 !
81 ! initializing ac and bc
82 !
83  do i=1,nteq
84  do j=1,nteq
85  ac(i,j)=0.d0
86  enddo
87  bc(i)=0.d0
88  enddo
89 !
90 ! for all but purely thermal networks:
91 ! determining the initial pressure
92 ! and identifying the chamber and gas pipe nodes
93 !
94  if(network.gt.2) then
95 !
96  call preinitialnet(ieg,lakon,v,ipkon,kon,nflow,prop,ielprop,
97  & ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,mi,iponoel,inoel,
98  & itg,ntg)
99 !
100 ! determining whether pressure initial conditions
101 ! are provided for all nodes
102 !
103  pressmin=-1.d0
104  pressmax=0.d0
105  constant=1.55d0
106 
107  do i=1,ntg
108  node=itg(i)
109  if(v(2,node).lt.1.d-10) then
110  v(2,node)=0.d0
111  else
112  if(pressmin.lt.0.d0) then
113  pressmin=v(2,node)
114  elseif(v(2,node).lt.pressmin) then
115  pressmin=v(2,node)
116  endif
117 !
118  if(v(2,node).gt.pressmax)then
119  pressmax=v(2,node)
120  endif
121 !
122  endif
123  enddo
124 !
125  if(pressmin.lt.0.d0) then
126  write(*,*)
127  & '*ERROR in initialnet: minimum initial pressure'
128  write(*,*) ' is smaller than zero'
129  call exit(201)
130  endif
131 !
132 ! in nodes in which no initial pressure is given v(2,*)
133 ! is replaced by -n, where n is the number of elements the
134 ! node belongs to: allows to find boundary nodes of the
135 ! network
136 !
137  do i=1,nflow
138  nelem=ieg(i)
139  index=ipkon(nelem)
140  node1=kon(index+1)
141  node2=kon(index+3)
142  call nident(itg,node1,ntg,id1)
143  call nident(itg,node2,ntg,id2)
144 !
145  if(iin_abs.eq.0) then
146  if ((lakon(nelem)(2:5).eq.'GAPF').or.
147  & ((lakon(nelem)(2:3).eq.'RE')
148  & .and.(lakon(nelem)(2:7).ne.'REWAOR')).or.
149  & (lakon(nelem)(2:3).eq.'UP')) then
150 !
151 ! In the case of a element of type GASPIPE or RESTRICTOR
152 ! (except TYPE= RESTRICTOR WALL ORIFICE)
153 ! the number of pipes connected to node 1 and 2
154 ! are computed and stored in ineighe(id1)
155 ! respectively ineighe(id2)
156 !
157  gaspipe=.true.
158  if(node1.ne.0) then
159  if (ineighe(id1).ge.0) then
160 !
161  if(node2.ne.0)then
162  ineighe(id1)=ineighe(id1)+1
163  endif
164  endif
165  endif
166  if(node2.ne.0) then
167  if (ineighe(id2).ge.0) then
168  if(node1.ne.0) then
169  ineighe(id2)=ineighe(id2)+1
170  endif
171  endif
172  endif
173  else
174 !
175 ! for all other elements (different from GASPIPE or
176 ! RESTRICTOR), including RESTRICTOR WALL ORIFICE
177 ! ineighe(idi)=-1
178 ! which means that they are connected to chambers
179 ! i.e. static and total values are equal
180 !
181  if (node1.ne.0) then
182  ineighe(id1)=-1
183  endif
184  if(node2.ne.0) then
185  ineighe(id2)=-1
186  endif
187  endif
188  endif
189 !
190  if((node1.eq.0).or.(node2.eq.0)) cycle
191  if(v(2,node1).lt.1.d-10) then
192  v(2,node1)=v(2,node1)-1.d0
193  endif
194  if(v(2,node2).lt.1.d-10) then
195  v(2,node2)=v(2,node2)-1.d0
196  endif
197  enddo
198 !
199 ! for each end node i: if ineighe(i)<0: chamber
200 ! else: ineighe(i)=number of pipe connections
201 !
202 ! assigning values to the boundary nodes of the network
203 ! (i.e. nodes belonging to only one element)
204 !
205  do i=1,ntg
206  node=itg(i)
207 ! boundary condition nodes
208  if(nactdog(2,node).eq.0) then
209  cycle
210  endif
211 ! initial condition nodes
212  if((nactdog(2,node).ne.0)
213  & .and.(v(2,node).gt.0d0)) then
214  cycle
215  endif
216 ! nodes neither defined as *BOUNDARY nor as *INITIAL CONDITIONS
217  if(abs(v(2,node)+1.d0).lt.1.d-10) then
218 !
219 ! if a nonzero flux is leaving the node increase
220 ! the maximum pressure, else decrease the minimum
221 ! pressure
222 !
223 ! determining the exit element (= element with one
224 ! zero node)
225 !
226  index=iponoel(node)
227  do
228  iel=inoel(1,index)
229  indexe=ipkon(iel)
230  node1=kon(indexe+1)
231  nodem=kon(indexe+2)
232  node2=kon(indexe+3)
233  if((node1.eq.0).or.(node2.eq.0)) exit
234  index=inoel(2,index)
235  enddo
236 !
237  if(((node1.eq.node).and.(v(1,nodem).ge.0.d0)).or.
238  & ((node2.eq.node).and.(v(1,nodem).le.0.d0))) then
239  v(2,node)=0.95d0*pressmin
240  pressmin=0.95d0*pressmin
241  else
242  v(2,node)=1.05d0*pressmax
243  pressmax=1.05d0*pressmax
244  endif
245  endif
246  enddo
247 !
248 ! transforming the pressures using the tangent
249 ! function
250 !
251  do i=1,ntg
252  node=itg(i)
253 ! boundary condition nodes
254  if(nactdog(2,node).eq.0) then
255  v(2,node)=dtan(v(2,node)
256  & *constant/pressmax)
257 ! initial condition nodes
258  elseif(v(2,node).gt.0d0) then
259  v(2,node)=dtan(v(2,node)
260  & *constant/pressmax)
261  endif
262  enddo
263 !
264 ! mass flow and geometry dofs: 1 on the diagonal
265 !
266  do i=1,ntg
267  node=itg(i)
268  do j=1,1
269  if(nactdog(j,node).ne.0) then
270  idof=nactdog(j,node)
271  ac(idof,idof)=1.d0
272  endif
273  enddo
274  do j=3,3
275  if(nactdog(j,node).ne.0) then
276  idof=nactdog(j,node)
277  ac(idof,idof)=1.d0
278  endif
279  enddo
280  enddo
281 !
282 ! pressure dofs
283 !
284  do i=1,nflow
285  nelem=ieg(i)
286  index=ipkon(nelem)
287  node1=kon(index+1)
288  node2=kon(index+3)
289  if((node1.eq.0).or.(node2.eq.0)) cycle
290  idof1=nactdog(2,node1)
291  idof2=nactdog(2,node2)
292  if(idof1.ne.0) then
293  if(v(2,node1).gt.0.d0) then
294 ! initial pressure given in node1
295  ac(idof1,idof1)=1.d0
296  bc(idof1)=v(2,node1)
297  else
298  ac(idof1,idof1)=ac(idof1,idof1)+1.d0
299  if(v(2,node2).gt.0.d0) then
300 ! initial pressure given in node2
301  bc(idof1)=bc(idof1)+v(2,node2)
302  else
303  ac(idof1,idof2)=ac(idof1,idof2)-1.d0
304  endif
305  endif
306  endif
307  if(idof2.ne.0) then
308  if(v(2,node2).gt.0.d0) then
309 ! initial pressure given in node2
310  ac(idof2,idof2)=1.d0
311  bc(idof2)=v(2,node2)
312  else
313  ac(idof2,idof2)=ac(idof2,idof2)+1.d0
314  if(v(2,node1).gt.0.d0) then
315 ! initial pressure given in node1
316  bc(idof2)=bc(idof2)+v(2,node1)
317  else
318  ac(idof2,idof1)=ac(idof2,idof1)-1.d0
319  endif
320  endif
321  endif
322  enddo
323 !
324 ! checking for nodes not belonging to network elements
325 !
326  do i=1,ntg
327  node=itg(i)
328  if(nactdog(2,node).ne.0) then
329  idof=nactdog(2,node)
330  if(ac(idof,idof).eq.0.d0) then
331  ac(idof,idof)=1.d0
332  bc(idof)=v(2,node)
333  endif
334  endif
335  enddo
336  endif
337 !
338 ! temperature conditions for all but purely hydrodynamic
339 ! networks
340 !
341  if(network.ne.4) then
342 !
343 ! temperature dofs
344 !
345  do i=1,nflow
346  nelem=ieg(i)
347  index=ipkon(nelem)
348  node1=kon(index+1)
349  nodem=kon(index+2)
350  node2=kon(index+3)
351  if((node1.eq.0).or.(node2.eq.0)) cycle
352  idof1=nactdog(0,node1)
353  idof2=nactdog(0,node2)
354  if(idof1.ne.0) then
355  if(v(0,node1).gt.0.d0) then
356 ! initial temperature given in node1
357  ac(idof1,idof1)=1.d0
358  bc(idof1)=v(0,node1)
359 !
360 ! temperature taken into account only for incoming
361 ! flux (from the viewpoint of node1).
362 !
363  elseif(v(1,nodem).lt.0.d0) then
364  ac(idof1,idof1)=ac(idof1,idof1)+1.d0
365  if(v(0,node2).gt.0.d0) then
366 ! initial temperature given in node2
367  bc(idof1)=bc(idof1)+v(0,node2)
368  else
369  ac(idof1,idof2)=ac(idof1,idof2)-1.d0
370  endif
371  endif
372  endif
373  if(idof2.ne.0) then
374  if(v(0,node2).gt.0.d0) then
375 ! initial temperature given in node2
376  ac(idof2,idof2)=1.d0
377  bc(idof2)=v(0,node2)
378 !
379 ! temperature taken into account only for incoming
380 ! flux (from the viewpoint of node2). Default (if
381 ! no flux is predefined) is positive flux.
382 !
383  elseif(v(1,nodem).ge.0.d0) then
384  ac(idof2,idof2)=ac(idof2,idof2)+1.d0
385  if(v(0,node1).gt.0.d0) then
386 ! initial temperature given in node1
387  bc(idof2)=bc(idof2)+v(0,node1)
388  else
389  ac(idof2,idof1)=ac(idof2,idof1)-1.d0
390  endif
391  endif
392  endif
393  enddo
394 !
395 ! checking for nodes not belonging to network elements
396 !
397  do i=1,ntg
398  node=itg(i)
399  if(nactdog(0,node).ne.0) then
400  idof=nactdog(0,node)
401  if(ac(idof,idof).eq.0.d0) then
402  ac(idof,idof)=1.d0
403  bc(idof)=v(0,node)
404  endif
405  endif
406  enddo
407  endif
408 !
409 ! additional multiple point constraints
410 !
411  do j=1,nteq
412  if(dabs(ac(j,j)).lt.1.d-20) ac(j,j)=1.d0
413  enddo
414 !
415 ! solving the system
416 !
417  if(nteq.gt.0) then
418  nrhs=1
419  call dgesv(nteq,nrhs,ac,nteq,ipiv,bc,nteq,info)
420  if(info.ne.0) then
421  write(*,*) '*ERROR in initialnet: singular matrix'
422  call exit(201)
423  endif
424  endif
425 !
426 ! storing the initial pressure in v
427 ! (not for purely thermal networks)
428 !
429  if(network.gt.2) then
430  do i=1,ntg
431  node=itg(i)
432  if(nactdog(2,node).eq.0) then
433  v(2,node)=pressmax
434  & *datan(v(2,node))/constant
435  cycle
436  endif
437  v(2,node)=pressmax
438  & *datan(bc(nactdog(2,node)))/constant
439  enddo
440  endif
441 !
442 ! storing the initial temperature in v
443 ! (not for purely hydrodynamic networks)
444 !
445  if(network.ne.4) then
446  do i=1,ntg
447  node=itg(i)
448  if(nactdog(0,node).ne.0) v(0,node)=bc(nactdog(0,node))
449  enddo
450  endif
451 !
452 ! for ELEMENT TYPE BRANCH
453 ! initial pressures in branch 1 and 2 are equal and set to the
454 ! maximum of the two pressures
455 !
456  do i=1,nflow
457  nelem=ieg(i)
458  if(lakon(i)(2:5).ne.'REBR') cycle
459  index=ielprop(nelem)
460 !
461  nelem1=nint(prop(index+2))
462  index1=ipkon(nelem1)
463  node11=kon(index1+1)
464  node21=kon(index1+3)
465 !
466  nelem2=nint(prop(index+3))
467  index2=ipkon(nelem2)
468  node12=kon(index2+1)
469  node22=kon(index2+3)
470 !
471  if(node11.eq.node12) then
472  node1=node21
473  node2=node22
474  elseif(node11.eq.node22) then
475  node1=node21
476  node2=node12
477  elseif(node21.eq.node12) then
478  node1=node11
479  node2=node22
480  elseif(node21.eq.node22) then
481  node1=node11
482  node2=node12
483  endif
484 !
485  if(lakon(i)(6:6).eq.'S') then
486 !
487 ! maximum
488 !
489  if(v(2,node1).ge.v(2,node2)) then
490  v(2,node2)=v(2,node1)
491  else
492  v(2,node1)=v(2,node2)
493  endif
494  elseif(lakon(i)(6:6).eq.'J') then
495 !
496 ! minimum
497 !
498  if(v(2,node1).ge.v(2,node2)) then
499  v(2,node1)=v(2,node2)
500  else
501  v(2,node2)=v(2,node1)
502  endif
503  endif
504  enddo
505 !
506 ! identifying the chamber nodes for purely thermal
507 ! gas networks (needed to determine the static
508 ! temperature which is used for the material properties)
509 !
510 c if(network.le.2) then
511  if((network.le.2).and.(iin_abs.eq.0)) then
512  do i=1,nflow
513  nelem=ieg(i)
514  if((lakon(nelem)(2:3).eq.'LP').or.
515  & (lakon(nelem)(2:3).eq.'LI')) cycle
516  index=ipkon(nelem)
517  node1=kon(index+1)
518  node2=kon(index+3)
519  if(node1.ne.0) then
520  call nident(itg,node1,ntg,id1)
521  ineighe(id1)=-1
522  endif
523  if(node2.ne.0) then
524  call nident(itg,node2,ntg,id2)
525  ineighe(id2)=-1
526  endif
527  enddo
528  endif
529 !
530 ! initialisation of bc
531 !
532  do i=1,nteq
533  bc(i)=0.d0
534  enddo
535 !
536 ! determining the initial mass flow in those nodes for which no
537 ! flux boundary conditions are defined
538 !
539  if(network.gt.2)then
540  do i=1,nflow
541  nelem=ieg(i)
542 !
543  index=ipkon(nelem)
544  node1=kon(index+1)
545  node2=kon(index+3)
546  if((node1.eq.0).or.(node2.eq.0)) cycle
547  nodem=kon(index+2)
548 !
549 ! cycle if both the geometry and the mass flow are known
550 !
551  if((nactdog(1,nodem).eq.0).and.(nactdog(3,nodem).eq.0))
552  & cycle
553 !
554 ! for liquids: determine the gravity vector
555 !
556  if(lakon(nelem)(2:3).eq.'LI') then
557  gravity=.false.
558  do j=1,3
559  g(j)=0.d0
560  enddo
561  if(nbody.gt.0) then
562  index=nelem
563  do
564  j=ipobody(1,index)
565  if(j.eq.0) exit
566  if(ibody(1,j).eq.2) then
567  g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j)
568  g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j)
569  g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j)
570  z1=-g(1)*co(1,node1)-g(2)*co(2,node1)
571  & -g(3)*co(3,node1)
572  z2=-g(1)*co(1,node2)-g(2)*co(2,node2)
573  & -g(3)*co(3,node2)
574  gravity=.true.
575  endif
576  index=ipobody(2,index)
577  if(index.eq.0) exit
578  enddo
579  endif
580  if(.not.gravity) then
581  write(*,*)
582  & '*ERROR in initialnet: no gravity vector'
583  write(*,*)
584  & ' was defined for liquid element',nelem
585  call exit(201)
586  endif
587  endif
588 !
589  tg1=v(0,node1)
590  tg2=v(0,node2)
591 !
592 ! For liquid pipes elements the upstream temperature
593 ! is used to determine the material properties
594 ! For all other elements the averaged value of the
595 ! temperature between inlet and outlet is used
596 !
597  if((lakon(nelem)(2:3).ne.'LP').and.
598  & (lakon(nelem)(2:3).ne.'LI')) then
599  gastemp=(tg1+tg2)/2.d0
600  else
601  xflow=v(1,nodem)
602  if(xflow.gt.0) then
603  gastemp=tg1
604  else
605  gastemp=tg2
606  endif
607  endif
608 !
609  imat=ielmat(1,nelem)
610  call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,cp,
611  & r,dvi,rhcon,nrhcon,rho)
612 !
613 ! If inlet and outlet pressure are equal this leads to a massflow rate
614 ! equal to 0 and in turn possibly to a singular matrix configuration
615 !
616 ! gas
617  if(((lakon(nelem)(2:3).ne.'LI').and.
618  & (v(2,node1).eq.v(2,node2))).or.
619 ! liquid
620  & ((lakon(nelem)(2:3).eq.'LI').and.
621  & (v(2,node1)+z1.eq.v(2,node2)+z2))) then
622 
623 !
624 ! if neither inlet nor outlet pressure are active D.O.F: error-> stop
625 !
626  if((nactdog(2,node1).eq.0)
627  & .and.(nactdog(2,node2).eq.0)) then
628  WRITE(*,*) '**************************************'
629  write(*,*) '*ERROR: in subroutine initialnet.f'
630  write(*,*) ' no pressure gradient'
631  write(*,*) ' in element', nelem
632  write(*,*) ' Inlet and outlet pressures are '
633  write(*,*) ' boundary conditions '
634  write(*,*) ' node1',node1,' pressure',
635  & v(2,node1)
636  write(*,*) ' node2',node2,' pressure',
637  & v(2,node2)
638  call exit(201)
639 !
640 ! if inlet pressure is an active degree of freedom
641 !
642  else if((nactdog(2,node1).ne.0)
643  & .and.(nactdog(2,node2).eq.0))then
644  WRITE(*,*) '**************************************'
645  write(*,*) '*WARNING: in subroutine initialnet.f'
646  write(*,*) ' no pressure gradient'
647  write(*,*) ' in element', nelem
648  write(*,*)
649  & ' Inlet pressure initial condition '
650  write(*,*) ' is changed '
651  write(*,*) ' node1',node1,
652  & ' given initial pressure',v(2,node1)
653  v(2,node1)=1.1*v(2,node1)
654  write(*,*) ' node1',node1,
655  & ' new initial pressure',v(2,node1)
656  write(*,*) ' node2',node2,' pressure',
657  & v(2,node2)
658 !
659 ! if outlet pressure is an active D.O.F.
660 !
661  else if((nactdog(2,node1).eq.0)
662  & .and.(nactdog(2,node2).ne.0))then
663  WRITE(*,*) '**************************************'
664  write(*,*) '*WARNING: in subroutine initialnet.f'
665  write(*,*) ' no pressure gradient'
666  write(*,*) ' in element', nelem
667  write(*,*)
668  & ' Outlet pressure initial condition '
669  write(*,*) ' is changed '
670  write(*,*) ' node1',node1,' pressure'
671  & ,v(2,node1)
672  write(*,*) ' node2',node2,
673  & 'given intial pressure',
674  & v(2,node2)
675  v(2,node2)=0.9*v(2,node2)
676  write(*,*) ' node2',node2,
677  & ' new initial pressure',v(2,node2)
678 !
679 ! if both inlet and outlet pressures are active D.O.F.
680 !
681  else if((nactdog(2,node1).ne.0)
682  & .and.(nactdog(2,node2).ne.0))then
683  WRITE(*,*) '**************************************'
684  write(*,*) '*WARNING: in subroutine initialnet.f'
685  write(*,*) ' no pressure gradient'
686  write(*,*) ' in element', nelem
687  write(*,*) ' Inlet and outlet pressure '
688  write(*,*) ' initial condition are changed '
689  write(*,*) ' node1',node1,
690  & ' given initial pressure',v(2,node1)
691  v(2,node1)=1.05*v(2,node2)
692  write(*,*) ' node1',node1,
693  & ' new intial pressure',v(2,node1)
694  write(*,*) ' node2',node2,
695  & ' given initial pressure',v(2,node2)
696  v(2,node2)=0.95*v(2,node2)
697  write(*,*) ' node2',node2,
698  & ' new intial pressure',v(2,node2)
699  endif
700  endif
701 !
702 ! calculating flux if the flux is an unknown AND there was
703 ! no initial flux defined by the user
704 !
705  if((nactdog(1,nodem).ne.0).and.(v(1,nodem).eq.0.d0)) then
706  call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
707  & nactdog,identity,ielprop,prop,kflag,v,xflow,f,
708  & nodef,idirf,df,cp,r,rho,physcon,g,co,dvi,numf,
709  & vold,set,shcon,nshcon,rhcon,nrhcon,ntmat_,mi,ider,
710  & ttime,time,iaxial)
711  v(1,nodem)=xflow
712  endif
713 !
714  if((lakon(nelem)(2:4).ne.'LIP').and.
715  & (lakon(nelem)(2:3).ne.'VO').and.
716  & (lakon(nelem)(2:4).ne.'ATR').and.
717  & (lakon(nelem)(2:4).ne.'RTA')) then
718  if(v(1,nodem).eq.0d0) then
719  WRITE(*,*) '**************************************'
720  write(*,*) '*ERROR: in subroutine initialnet.f'
721  write(*,*) ' in element', nelem,
722  & lakon(nelem)(1:6)
723  write(*,*) ' mass flow rate value = 0 !'
724  write(*,*) ' node1',node1,' pressure',
725  & v(2,node1)
726  write(*,*) ' node2',node2,' pressure',
727  & v(2,node2)
728  call exit(201)
729  endif
730  if (v(1,nodem)*(v(2,node1)-v(2,node2)).lt.0) then
731  WRITE(*,*) '**************************************'
732  write(*,*) '*WARNING: in subroutine initialnet.f'
733  write(*,*) ' in element', nelem
734  write(*,*)
735  & ' initial mass flow rate value does not'
736  write(*,*) ' correspond to pressure gradient'
737  write(*,*) ' node1',node1,'pressure',
738  & v(2,node1)
739  write(*,*) ' node2',node2,'pressure',
740  & v(2,node2)
741  write(*,*) ' nodem',nodem,'mass flow',
742  & v(1,nodem)
743  write(*,*)
744  & ' the sign of the initial mass flow'
745  write(*,*) ' rate will be changed'
746  endif
747  endif
748  enddo
749  endif
750 !
751  call postinitialnet(ieg,lakon,v,ipkon,kon,nflow,prop,ielprop,
752  & ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,mi,iponoel,inoel,
753  & itg,ntg)
754 !
755 ! calculating the static temperature for nodes belonging to gas pipes
756 ! and restrictors (except RESTRICTOR WALL ORIFICE)
757 !
758  if (gaspipe.and.(iin_abs.eq.0)) then
759 !
760 ! ineighe(i) is set to -1 (= chamber node) if more than 2 pipes
761 ! are connected to the node
762 !
763  do i=1,ntg
764  if(ineighe(i).gt.2) then
765  ineighe(i)=-1
766  write(*,*) '*WARNING :in subroutine initialnet.f'
767  write(*,*) ' more than 2 elements GASPIPE'
768  write(*,*) ' or RESTRICTOR are connected '
769  write(*,*) ' to node',itg(i),'. The common'
770  write(*,*)
771  & ' node is converted into a chamber.'
772  write(*,*) ' Total and static parameters are'
773  write(*,*) ' equal'
774  endif
775  enddo
776 !
777  do i=1,nflow
778  nelem=ieg(i)
779  index=ipkon(nelem)
780  node1=kon(index+1)
781  node2=kon(index+3)
782  call nident(itg,node1,ntg,id1)
783  call nident(itg,node2,ntg,id2)
784 !
785 ! for each end node i:
786 ! if ineighe(i)=-1: chamber
787 ! = 0: middle node
788 ! > 0: number of pipe connections (max. 2 allowed)
789 !
790 ! the number of pipe connections for ineighe(i)>0 is
791 ! replaced by the number of an element the node belongs to:
792 !
793  if(node1.gt.0) then
794  if((ineighe(id1).eq.1).or.
795  & (ineighe(id1).eq.2)) then
796  if(node2.ne.0)then
797  ineighe(id1)=nelem
798  endif
799  endif
800  endif
801 !
802  if(node2.gt.0) then
803  if((ineighe(id2).eq.1).or.
804  & (ineighe(id2).eq.2)) then
805  if (node1.ne.0) then
806  ineighe(id2)=nelem
807  endif
808  endif
809  endif
810  enddo
811 !
812 ! The static temperature is calculated and stored in v(3,node)
813 ! total temperatures are supposed equal (adiabatic pipe)
814 !
815  do i=1,ntg
816  node=itg(i)
817  if(ineighe(i).gt.0) then
818 !
819  nelem=ineighe(i)
820  index=ielprop(nelem)
821  nodem=kon(ipkon(nelem)+2)
822 !
823  imat=ielmat(1,nelem)
824  call materialdata_tg(imat,ntmat_,v(0,node),
825  & shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho)
826  kappa=cp/(cp-r)
827  xflow=v(1,nodem)
828  tt=v(0,node)
829  pt=v(2,node)
830 !
831  if(lakon(nelem)(2:5).eq.'GAPF') then
832  a=prop(index+1)
833  if(lakon(nelem)(2:6).eq.'GAPFA') then
834  icase=0
835  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
836  icase=1
837  endif
838 !
839  elseif(lakon(nelem)(2:3).eq.'RE') then
840  index2=ipkon(nelem)
841  node1=kon(index2+1)
842  node2=kon(index2+3)
843  if(lakon(nelem)(4:5).eq.'EX') then
844  if(lakon(nint(prop(index+4)))(2:6).eq.'GAPFA')
845  & then
846  icase=0
847  elseif
848  & (lakon(nint(prop(index+4)))(2:6).eq.'GAPFI')
849  & then
850  icase=1
851  endif
852  else
853  icase=0
854  endif
855 !
856  if(lakon(nelem)(4:5).eq.'BE') then
857  a=prop(index+1)
858 !
859  elseif(lakon(nelem)(4:5).eq.'BR') then
860  if(lakon(nelem)(4:6).eq.'BRJ') then
861  if(nelem.eq.nint(prop(index+2)))then
862  a=prop(index+5)
863  elseif(nelem.eq.nint(prop(index+3))) then
864  a=prop(index+6)
865  endif
866  elseif(lakon(nelem)(4:6).eq.'BRS') then
867  if(nelem.eq.nint(prop(index+2)))then
868  a=prop(index+5)
869  elseif(nelem.eq.nint(prop(index+3))) then
870  a=prop(index+6)
871  endif
872  endif
873 !
874  else
875  if(node.eq.node1) then
876  a=prop(index+1)
877  elseif(node.eq.node2) then
878  a=prop(index+2)
879  endif
880  endif
881  elseif(lakon(nelem)(2:3).eq.'UP') then
882 !
883 ! user elements whose names start with UP are assumed
884 ! to be Pipe-like (static and total temperatures differ)
885 ! The cross area is assumed to be the first property
886 !
887  a=prop(index+1)
888  icase=0
889  endif
890 !
891  if(v(3,node).eq.0.d0) then
892  xflow360=xflow*iaxial
893  call ts_calc(xflow360,tt,pt,kappa,r,a,ts,icase)
894  v(3,node)=ts
895  endif
896 !
897 ! if the element is not of gaspipe or branch type,
898 ! total and static temperatures are equal for all endnodes
899 !
900  endif
901  enddo
902  endif
903 !
904 ! for chambers the static temperature equals the total
905 ! temperature
906 !
907  do i=1,ntg
908  if(ineighe(i).eq.-1) v(3,itg(i))=v(0,itg(i))
909 c write(*,*) 'initialnet ',i,v(0,itg(i)),
910 c & v(1,itg(i)),v(2,itg(i)),ineighe(i)
911  enddo
912 !
913 c write(*,*) 'initialnet '
914 c do i=1,ntg
915 c write(*,'(i10,3(1x,e11.4))') itg(i),(v(j,itg(i)),j=0,2)
916 c enddo
917 !
918  return
z1
static double * z1
Definition: filtermain.c:48
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
flux
subroutine flux(node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, kflag, v, xflow, f, nodef, idirf, df, cp, R, rho, physcon, g, co, dvi, numf, vold, set, shcon, nshcon, rhcon, nrhcon, ntmat_, mi, ider, ttime, time, iaxial)
Definition: flux.f:24
preinitialnet
subroutine preinitialnet(ieg, lakon, v, ipkon, kon, nflow, prop, ielprop, ielmat, ntmat_, shcon, nshcon, rhcon, nrhcon, mi, iponoel, inoel, itg, ntg)
Definition: preinitialnet.f:22
ts_calc
subroutine ts_calc(xflow, Tt, Pt, kappa, r, A, Ts, icase)
Definition: ts_calc.f:20
dgesv
subroutine dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:58
postinitialnet
subroutine postinitialnet(ieg, lakon, v, ipkon, kon, nflow, prop, ielprop, ielmat, ntmat_, shcon, nshcon, rhcon, nrhcon, mi, iponoel, inoel, itg, ntg)
Definition: postinitialnet.f:22
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
materialdata_tg
subroutine materialdata_tg(imat, ntmat_, t1l, shcon, nshcon, sph, r, dvi, rhcon, nrhcon, rho)
Definition: materialdata_tg.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)