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

Go to the source code of this file.

Functions/Subroutines

subroutine initialchannel (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, ttime, time, iaxial)
 

Function/Subroutine Documentation

◆ initialchannel()

subroutine initialchannel ( 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,
real*8  ttime,
real*8  time,
integer  iaxial 
)
33 !
34  implicit none
35 !
36  logical identity,calcinitialpressure,gravity,gaspipe
37 !
38  character*8 lakon(*)
39  character*81 set(*)
40 !
41  integer mi(*),ieg(*),nflow,i,j,ntg,ielmat(mi(3),*),ntmat_,id,
42  & node1,node2,
43  & nelem,index,nshcon(*),ipkon(*),kon(*),ikboun(*),nboun,idof,
44  & nodem,idirf(8),nactdog(0:3,*),imat,ielprop(*),id1,id2,
45  & nodef(8),ndirboun(*),nodeboun(*),itg(*),node,kflag,ipiv(*),
46  & nrhs,info,idof1,idof2,nteq,nrhcon(*),ipobody(2,*),ibody(3,*),
47  & nbody,numf,network,iin_abs,icase,index2,index1,nelem1,nelem2,
48  & node11,node21,node12,node22,istep,iit,ineighe(*),
49  & ilboun(*),nelemup,k,node2up,ider,iaxial
50 !
51  real*8 ac(nteq,nteq), bc(nteq),prop(*),shcon(0:3,ntmat_,*),
52  & f,df(8),xflow,xbounact(*),v(0:mi(2),*),cp,r,tg1,
53  & tg2,gastemp,physcon(*),pressmin,dvi,rho,g(3),z1,z2,
54  & rhcon(0:1,ntmat_,*),co(3,*),xbodyact(7,*),kappa,
55  & a,tt,pt,ts,pressmax,constant,vold(0:mi(2),*),href,
56  & ttime,time
57 !
58  kflag=1
59 !
60 ! applying the boundary conditions
61 !
62  do j=1,nboun
63  v(ndirboun(j),nodeboun(j))=xbounact(j)
64  enddo
65 !
66 ! determining the initial pressure (not for purely thermal networks)
67 ! and identifying the chamber and gas pipe nodes (only for gas
68 ! networks)
69 !
70  if(network.gt.2) then
71 !
72 ! determining whether pressure initial conditions
73 ! are provided for all nodes
74 !
75  pressmin=-1.d0
76  pressmax=0.d0
77  constant=1.55d0
78 
79  do i=1,ntg
80  node=itg(i)
81  if(v(2,node).lt.1.d-10) then
82  v(2,node)=0.d0
83  else
84  if(pressmin.lt.0.d0) then
85  pressmin=v(2,node)
86  elseif(v(2,node).lt.pressmin) then
87  pressmin=v(2,node)
88  endif
89 !
90  if(v(2,node).gt.pressmax)then
91  pressmax=v(2,node)
92  endif
93 !
94  endif
95  enddo
96 !
97  if(pressmin.lt.0.d0) then
98  write(*,*)
99  & '*ERROR in initialchannel: minimum initial pressure'
100  write(*,*) ' is smaller than zero'
101  call exit(201)
102  endif
103 !
104 ! in nodes in which no initial pressure is given v(2,*)
105 ! is replaced by -n, where n is the number of elements the
106 ! node belongs to: allows to find boundary nodes of the
107 ! network
108 !
109  gaspipe=.false.
110  calcinitialpressure=.false.
111 !
112  do i=1,nflow
113  nelem=ieg(i)
114  index=ipkon(nelem)
115  node1=kon(index+1)
116  node2=kon(index+3)
117  call nident(itg,node1,ntg,id1)
118  call nident(itg,node2,ntg,id2)
119 !
120  if (((lakon(nelem)(1:5).eq.'DGAPF').and.(iin_abs.eq.0))
121  & .or.((lakon(nelem)(1:3).eq.'DRE')
122  & .and.(lakon(nelem)(1:7).ne.'DREWAOR')
123  & .and.(iin_abs.eq.0))) then
124 !
125 ! In the case of a element of type GASPIPE or RESTRICTOR
126 ! (except TYPE= RESTRICTOR WALL ORIFICE)
127 ! the number of pipes connected to node 1 and 2
128 ! are computed and stored in ineighe(id1)
129 ! respectively ineighe(id2)
130 !
131  gaspipe=.true.
132  if(node1.ne.0) then
133  if (ineighe(id1).ge.0) then
134 !
135  if(node2.ne.0)then
136  ineighe(id1)=ineighe(id1)+1
137  endif
138  endif
139  endif
140  if(node2.ne.0) then
141  if (ineighe(id2).ge.0) then
142  if(node1.ne.0) then
143  ineighe(id2)=ineighe(id2)+1
144  endif
145  endif
146  endif
147  else
148  if(iin_abs.eq.0) then
149 !
150 ! for all other elements (different from GASPIPE or
151 ! RESTRICTOR), including RESTRICTOR WALL ORIFICE
152 ! ineighe(idi)=-1
153 ! which means that they are connected to chambers
154 ! i.e. static and total values are equal
155 !
156  if (node1.ne.0) then
157  ineighe(id1)=-1
158  endif
159  if(node2.ne.0) then
160  ineighe(id2)=-1
161  endif
162  endif
163  endif
164 !
165  if((node1.eq.0).or.(node2.eq.0)) cycle
166  if(v(2,node1).lt.1.d-10) then
167  v(2,node1)=v(2,node1)-1.d0
168  calcinitialpressure=.true.
169  endif
170  if(v(2,node2).lt.1.d-10) then
171  v(2,node2)=v(2,node2)-1.d0
172  calcinitialpressure=.true.
173  endif
174  enddo
175 !
176 ! for each end node i: if ineighe(i)<0: chamber
177 ! else: ineighe(i)=number of pipe connections
178 !
179  else
180 !
181 ! identifying the chamber nodes for purely thermal
182 ! gas networks (needed to determine the static
183 ! temperature which is used for the material properties)
184 !
185  do i=1,nflow
186  nelem=ieg(i)
187  if((lakon(nelem)(2:3).eq.'LP').or.
188  & (lakon(nelem)(2:3).eq.'LI')) cycle
189  index=ipkon(nelem)
190  node1=kon(index+1)
191  node2=kon(index+3)
192  if(node1.ne.0) then
193  call nident(itg,node1,ntg,id1)
194  ineighe(id1)=-1
195  endif
196  if(node2.ne.0) then
197  call nident(itg,node2,ntg,id2)
198  ineighe(id2)=-1
199  endif
200  enddo
201  endif
202 !
203 ! temperature initial conditions
204 !
205  do i=1,ntg
206  node=itg(i)
207  if (nactdog(0,node).eq.0) cycle
208  if (v(0,node)-physcon(1).lt.1.d-10) then
209  write(*,*)
210  & '*WARNING in initialchannel : the initial temperature f
211  &or node',node
212  write(*,*)
213  & 'is O Kelvin or less; the default is taken (293 K)'
214  write(*,*)
215  v(0,node)=293.d0-physcon(1)
216  endif
217  enddo
218 !
219 ! initialisation of bc
220 !
221  do i=1,nteq
222  bc(i)=0.d0
223  enddo
224 !
225 ! determining the initial mass flow in those nodes for which no
226 ! flux boundary conditions are defined
227 ! liquid channels are treated separately
228 !
229  if(network.gt.2) then
230 !
231 ! calculate the initial mass flow
232 !
233 ! check whether the mass flow is given as a boundary condition
234 !
235  do j=1,nflow
236  nelem=ieg(j)
237  index=ipkon(nelem)
238  nodem=kon(index+2)
239  if(nactdog(1,nodem).eq.0) then
240  idof=8*(nodem-1)+1
241  call nident(ikboun,idof,nboun,id)
242  if(id.gt.0) then
243  if(ikboun(id).eq.idof) then
244  xflow=xbounact(ilboun(id))
245  if(dabs(xflow).gt.1.d-30) exit
246  endif
247  endif
248  endif
249  enddo
250 !
251  if(dabs(xflow).gt.1.d-30) then
252 !
253 ! if nonzero: set all mass flow to this value
254 !
255  do j=1,nflow
256  nelem=ieg(j)
257  index=ipkon(nelem)
258  nodem=kon(index+2)
259  if(nactdog(1,nodem).ne.0) v(1,nodem)=xflow
260  enddo
261  else
262 !
263 ! calculate the mass flow: look for a sluice gate or weir
264 !
265  do j=1,nflow
266  nelem=ieg(j)
267  if((lakon(nelem)(6:7).ne.'SG').and.
268  & (lakon(nelem)(6:7).ne.'WE')) cycle
269  index=ipkon(nelem)
270  node1=kon(index+1)
271  node2=kon(index+3)
272  if((node1.eq.0).or.(node2.eq.0)) cycle
273  nodem=kon(index+2)
274 !
275 ! determine the gravity vector
276 !
277  gravity=.false.
278  do k=1,3
279  g(k)=0.d0
280  enddo
281  if(nbody.gt.0) then
282  index=nelem
283  do
284  k=ipobody(1,index)
285  if(k.eq.0) exit
286  if(ibody(1,k).eq.2) then
287  g(1)=g(1)+xbodyact(1,k)*xbodyact(2,k)
288  g(2)=g(2)+xbodyact(1,k)*xbodyact(3,k)
289  g(3)=g(3)+xbodyact(1,k)*xbodyact(4,k)
290  gravity=.true.
291  endif
292  index=ipobody(2,index)
293  if(index.eq.0) exit
294  enddo
295  endif
296  if(.not.gravity) then
297  write(*,*)
298  & '*ERROR in initialchannel: no gravity vector'
299  write(*,*) ' was defined for liquid element',
300  & nelem
301  call exit(201)
302  endif
303 !
304  tg1=v(0,node1)
305  tg2=v(0,node2)
306  gastemp=(tg1+tg2)/2.d0
307  imat=ielmat(1,nelem)
308  call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,
309  & cp,r,dvi,rhcon,nrhcon,rho)
310 !
311  call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
312  & nactdog,identity,ielprop,prop,kflag,v,xflow,f,
313  & nodef,idirf,df,cp,r,rho,physcon,g,co,dvi,numf,
314  & vold,set,shcon,nshcon,rhcon,nrhcon,ntmat_,mi,ider,
315  & ttime,time,iaxial)
316 !
317  if(dabs(xflow).gt.1.d-30) exit
318  enddo
319 !
320  if(dabs(xflow).gt.1.d-30) then
321 !
322 ! if nonzero: set all mass flow to this value
323 !
324  do j=1,nflow
325  nelem=ieg(j)
326  index=ipkon(nelem)
327  nodem=kon(index+2)
328  if(nactdog(1,nodem).ne.0) v(1,nodem)=xflow
329  enddo
330  else
331  write(*,*) '*ERROR in initialchannel: initial mass flow'
332  write(*,*) ' cannot be determined'
333  call exit(201)
334  endif
335  endif
336 !
337 ! calculate the depth
338 !
339  if(calcinitialpressure) then
340 !
341 ! determine the streamdown depth for sluice gates,
342 ! weirs and discontinuous slopes
343 !
344  do j=1,nflow
345  nelem=ieg(j)
346  if((lakon(nelem)(6:7).ne.'SG').and.
347  & (lakon(nelem)(6:7).ne.'WE').and.
348  & (lakon(nelem)(6:7).ne.'DS')) cycle
349  index=ipkon(nelem)
350  node1=kon(index+1)
351  node2=kon(index+3)
352  if((node1.eq.0).or.(node2.eq.0)) cycle
353  nodem=kon(index+2)
354 !
355 ! determine the gravity vector
356 !
357  gravity=.false.
358  do k=1,3
359  g(k)=0.d0
360  enddo
361  if(nbody.gt.0) then
362  index=nelem
363  do
364  k=ipobody(1,index)
365  if(k.eq.0) exit
366  if(ibody(1,k).eq.2) then
367  g(1)=g(1)+xbodyact(1,k)*xbodyact(2,k)
368  g(2)=g(2)+xbodyact(1,k)*xbodyact(3,k)
369  g(3)=g(3)+xbodyact(1,k)*xbodyact(4,k)
370  gravity=.true.
371  endif
372  index=ipobody(2,index)
373  if(index.eq.0) exit
374  enddo
375  endif
376  if(.not.gravity) then
377  write(*,*)
378  & '*ERROR in initialchannel: no gravity vector'
379  write(*,*) ' was defined for liquid element',
380  & nelem
381  call exit(201)
382  endif
383 !
384  tg1=v(0,node1)
385  tg2=v(0,node2)
386  gastemp=(tg1+tg2)/2.d0
387  imat=ielmat(1,nelem)
388  call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,
389  & cp,r,dvi,rhcon,nrhcon,rho)
390 !
391  call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
392  & nactdog,identity,ielprop,prop,kflag,v,xflow,f,
393  & nodef,idirf,df,cp,r,rho,physcon,g,co,dvi,numf,
394  & vold,set,shcon,nshcon,rhcon,nrhcon,ntmat_,mi,ider,
395  & ttime,time,iaxial)
396 !
397  enddo
398 !
399 ! for all other elements the depth is taken to be
400 ! 0.9 of the depth in the downstream node of the
401 ! streamup reference element
402 !
403  do j=1,nflow
404  nelem=ieg(j)
405  if((lakon(nelem)(6:7).eq.'SG').or.
406  & (lakon(nelem)(6:7).eq.'WE').or.
407  & (lakon(nelem)(6:7).eq.'DS')) cycle
408 !
409  index=ipkon(nelem)
410  node1=kon(index+1)
411  node2=kon(index+3)
412  if((node1.eq.0).or.(node2.eq.0)) cycle
413 !
414  index=ielprop(nelem)
415  nelemup=nint(prop(index+6))
416  node2up=kon(ipkon(nelemup)+3)
417  href=0.9d0*v(2,node2up)
418  if(nactdog(2,node1).ne.0)
419  & v(2,node1)=href
420  if(nactdog(2,node2).ne.0)
421  & v(2,node2)=href
422  enddo
423 !
424 ! reapplying the boundary conditions (the depth of the
425 ! sluice gate may have changed if it exceeded the critical
426 ! value
427 !
428  do j=1,nboun
429  v(ndirboun(j),nodeboun(j))=xbounact(j)
430  enddo
431  endif
432 !
433 ! calculating the static temperature for nodes belonging to gas pipes
434 ! and restrictors (except RESTRICTOR WALL ORIFICE)
435 !
436  endif
437 !
438 ! for chambers the static temperature equals the total
439 ! temperature
440 !
441  do i=1,ntg
442  if(ineighe(i).eq.-1) v(3,itg(i))=v(0,itg(i))
443  enddo
444 !
445  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
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)