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

Go to the source code of this file.

Functions/Subroutines

subroutine mafillnet (itg, ieg, ntg, ac, nload, sideload, nelemload, xloadact, lakon, ntmat_, v, shcon, nshcon, ipkon, kon, co, nflow, iinc, istep, dtime, ttime, time, ielmat, nteq, prop, ielprop, nactdog, nacteq, physcon, rhcon, nrhcon, ipobody, ibody, xbodyact, nbody, vold, xloadold, reltime, nmethod, set, mi, nmpc, nodempc, ipompc, coefmpc, labmpc, iaxial, cocon, ncocon, iponoel, inoel)
 

Function/Subroutine Documentation

◆ mafillnet()

subroutine mafillnet ( integer, dimension(*)  itg,
integer, dimension(*)  ieg,
integer  ntg,
real*8, dimension(nteq,*)  ac,
integer  nload,
character*20, dimension(*)  sideload,
integer, dimension(2,*)  nelemload,
real*8, dimension(2,*)  xloadact,
character*8, dimension(*)  lakon,
integer  ntmat_,
real*8, dimension(0:mi(2),*)  v,
real*8, dimension(0:3,ntmat_,*)  shcon,
integer, dimension(*)  nshcon,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
real*8, dimension(3,*)  co,
integer  nflow,
integer  iinc,
integer  istep,
real*8  dtime,
real*8  ttime,
real*8  time,
integer, dimension(mi(3),*)  ielmat,
integer  nteq,
real*8, dimension(*)  prop,
integer, dimension(*)  ielprop,
integer, dimension(0:3,*)  nactdog,
integer, dimension(0:3,*)  nacteq,
real*8, dimension(*)  physcon,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer, dimension(2,*)  ipobody,
integer, dimension(3,*)  ibody,
real*8, dimension(7,*)  xbodyact,
integer  nbody,
real*8, dimension(0:mi(2),*)  vold,
real*8, dimension(2,*)  xloadold,
real*8  reltime,
integer  nmethod,
character*81, dimension(*)  set,
integer, dimension(*)  mi,
integer  nmpc,
integer, dimension(3,*)  nodempc,
integer, dimension(*)  ipompc,
real*8, dimension(*)  coefmpc,
character*20, dimension(*)  labmpc,
integer  iaxial,
real*8, dimension(0:6,ntmat_,*)  cocon,
integer, dimension(2,*)  ncocon,
integer, dimension(*)  iponoel,
integer, dimension(2,*)  inoel 
)
28 !
29  implicit none
30 !
31  logical identity,crit
32 !
33  character*8 lakonl,lakon(*)
34  character*20 sideload(*),labmpc(*)
35  character*81 set(*)
36 !
37  integer mi(*),itg(*),ieg(*),ntg,nteq,nflow,nload,inv,
38  & ielmat(mi(3),*),iaxial,ncocon(2,*),iponoel(*),inoel(2,*),
39  & nelemload(2,*),nope,nopes,mint2d,i,j,k,l,iflag,
40  & node,imat,ntmat_,id,ifaceq(8,6),ifacet(6,4),
41  & ifacew(8,5),node1,node2,nshcon(*),nelem,ig,index,konl(20),
42  & ipkon(*),kon(*),idof,iinc,ibody(3,*),istep,jltyp,nfield,
43  & ipobody(2,*),nodem,ieq,kflag,nrhcon(*),numf,k_oil,
44  & idofp1,idofp2,idofm,idoft1,idoft2,idoft,nactdog(0:3,*),
45  & nacteq(0:3,*),ielprop(*),nodef(8),idirf(8),nbody,
46  & nmethod,icase,nmpc,nodempc(3,*),ipompc(*),idir,ider
47 !
48  real*8 ac(nteq,*),xloadact(2,*),cp,h(2),physcon(*),dvi,
49  & xl2(3,8),coords(3),dxsj2,temp,xi,et,weight,xsj2(3),
50  & gastemp,v(0:mi(2),*),shcon(0:3,ntmat_,*),co(3,*),shp2(7,8),
51  & ftot,field,prop(*),f,df(8),tg1,tg2,r,rho,tl2(8),dl,
52  & dtime,ttime,time,areaj,xflow,tvar(2),g(3),coefmpc(*),
53  & rhcon(0:1,ntmat_,*),xbodyact(7,*),sinktemp,ts1,ts2,xs2(3,7),
54  & pt1,pt2,vold(0:mi(2),*),xloadold(2,*),reltime,pi,d,
55  & cocon(0:6,ntmat_,*),xflow360,xflow_oil,ks,form_fact,
56  & qred_crit,reynolds,pt2zpt1_c,qred1,pt2zpt1,phi,lambda,
57  & m1,m2,bb,cc,km1,kp1,a,kappa,ff1,ff2,tt1,tt2,t1,t2,m1_c,
58  & heatnod,heatfac
59 !
60  include "gauss.f"
61 !
62  data ifaceq /4,3,2,1,11,10,9,12,
63  & 5,6,7,8,13,14,15,16,
64  & 1,2,6,5,9,18,13,17,
65  & 2,3,7,6,10,19,14,18,
66  & 3,4,8,7,11,20,15,19,
67  & 4,1,5,8,12,17,16,20/
68  data ifacet /1,3,2,7,6,5,
69  & 1,2,4,5,9,8,
70  & 2,3,4,6,10,9,
71  & 1,4,3,8,10,7/
72  data ifacew /1,3,2,9,8,7,0,0,
73  & 4,5,6,10,11,12,0,0,
74  & 1,2,5,4,7,14,10,13,
75  & 2,3,6,5,8,15,11,14,
76  & 4,6,3,1,12,15,9,13/
77 !
78  data iflag /2/
79 !
80  kflag=2
81  ider=1
82 !
83  pi=4.d0*datan(1.d0)
84  tvar(1)=time
85  tvar(2)=ttime+time
86 !
87 ! reinitialisation of the Ac matrix
88 !
89  do i=1,nteq
90  do j=1,nteq
91  ac(i,j)=0.d0
92  enddo
93  enddo
94 !
95 ! solving for the gas temperatures in forced convection
96 !
97  ftot=0.d0
98 !
99 ! element contribution.
100 !
101 
102  do i=1,nflow
103  nelem=ieg(i)
104  index=ipkon(nelem)
105  node1=kon(index+1)
106  nodem=kon(index+2)
107  node2=kon(index+3)
108 !
109  xflow=v(1,nodem)
110 !
111 ! next section is for liquids only (no gas)
112 !
113  if((lakon(nelem)(2:3).ne.'LP').and.
114  & (lakon(nelem)(2:3).ne.'LI')) then
115  if(node1.eq.0) then
116  tg1=v(0,node2)
117  tg2=tg1
118  ts1=v(3,node2)
119  ts2=ts1
120  elseif(node2.eq.0) then
121  tg1=v(0,node1)
122  tg2=tg1
123  ts1=v(3,node1)
124  ts2=ts1
125  else
126  tg1=v(0,node1)
127  tg2=v(0,node2)
128  ts1=v(3,node1)
129  ts2=v(3,node2)
130  endif
131 !
132  gastemp=(ts1+ts2)/2.d0
133 !
134 ! for liquid pipe element only the upstream temperature is used to
135 ! determine thematerial properties
136 !
137  else
138 
139  if(xflow.gt.0) then
140  if(node1.eq.0) then
141  gastemp=v(0,node2)
142  else
143  gastemp=v(0,node1)
144  endif
145  else
146  if(node2.eq.0) then
147  gastemp=v(0,node1)
148  else
149  gastemp=v(0,node2)
150  endif
151  endif
152 !
153  if(node1.eq.0) then
154  tg2=v(0,node2)
155  tg1=tg2
156  elseif(node2.eq.0) then
157  tg1=v(0,node1)
158  tg2=tg1
159  else
160  tg1=v(0,node1)
161  tg2=v(0,node2)
162  endif
163  endif
164 !
165  imat=ielmat(1,nelem)
166 !
167  call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,cp,r,dvi,
168  & rhcon,nrhcon,rho)
169 !
170  kappa=(cp/(cp-r))
171 !
172  if(node1.ne.0) then
173  idoft1=nactdog(0,node1)
174  idofp1=nactdog(2,node1)
175  else
176  idoft1=0
177  idofp1=0
178  endif
179  if(node2.ne.0) then
180  idoft2=nactdog(0,node2)
181  idofp2=nactdog(2,node2)
182  else
183  idoft2=0
184  idofp2=0
185  endif
186  idofm=nactdog(1,nodem)
187 !
188 ! Definitions of the constant for isothermal flow elements
189 !
190  if(lakon(nelem)(2:6).eq.'GAPFI') then
191  if((node1.ne.0).and.(node2.ne.0)) then
192  icase=1
193 !
194 ! properties
195 !
196  index=ielprop(nelem)
197  a=prop(index+1)
198  d=prop(index+2)
199  dl=prop(index+3)
200  ks=prop(index+4)
201  form_fact=prop(index+5)
202  xflow_oil=prop(index+6)
203  k_oil=nint(prop(index+7))
204 !
205  km1=kappa-1.d0
206  kp1=kappa+1.d0
207 !
208 c xflow360=xflow*iaxial
209 !
210 ! orientation of the element based on the direction of the
211 ! mass flow; determination of the total pressures and total
212 ! temperatures (K)
213 !
214 ! xflow360 is the absolute value of the mass flow for
215 ! 360 degrees
216 !
217 c if(xflow.gt.0.d0) then
218  if(v(2,node1).gt.v(2,node2)) then
219  inv=1
220  xflow360=xflow*iaxial
221  pt1=v(2,node1)
222  pt2=v(2,node2)
223  tt1=v(0,node1)-physcon(1)
224  tt2=v(0,node2)-physcon(1)
225  else
226  inv=-1
227  xflow360=-xflow*iaxial
228  pt1=v(2,node2)
229  pt2=v(2,node1)
230  tt1=v(0,node2)-physcon(1)
231  tt2=v(0,node1)-physcon(1)
232 !
233  idoft1=nactdog(0,node2)
234  idofp1=nactdog(2,node2)
235  idoft2=nactdog(0,node1)
236  idofp2=nactdog(2,node1)
237  endif
238 !
239 ! calculation of the static temperatures T1 and T2
240 ! (K)
241 !
242  call ts_calc(xflow360,tt1,pt1,kappa,r,a,t1,icase)
243 c call ts_calc(xflow360,Tt2,Pt2,kappa,r,a,T2,icase)
244  endif
245  endif
246 !
247  if(node1.ne.0) then
248 !
249 ! energy equation contribution node1
250 !
251  if(nacteq(0,node1).ne.0) then
252  ieq=nacteq(0,node1)
253  if((xflow.le.0.d0).and.(nacteq(3,node1).eq.0))then
254 !
255 ! incoming flow in node1 and genuine energy conservation equation
256 !
257  if(idoft1.ne.0) then
258  ac(ieq,idoft1)=ac(ieq,idoft1)-cp*xflow
259  endif
260 !
261  if(idoft2.ne.0)then
262  ac(ieq,idoft2)=ac(ieq,idoft2)+cp*xflow
263  elseif(node2.eq.0) then
264  write(*,*) '*WARNING in mafillnet'
265  write(*,*) ' temperature at inlet'
266  write(*,*) ' node',node1,' unknown'
267 c call exit(201)
268  endif
269 !
270  if(idofm.ne.0) then
271  ac(ieq,idofm)=ac(ieq,idofm)-cp*(tg1-tg2)
272  endif
273 !
274  elseif((node2.ne.0).and.(nacteq(3,node1).eq.node2)) then
275 !
276 ! isothermal element
277 !
278  pt2zpt1=pt2/pt1
279 !
280 ! calculation of the dynamic viscosity
281 !
282  if(dabs(dvi).lt.1e-30) then
283  write(*,*) '*ERROR in gaspipe_fanno: '
284  write(*,*) ' no dynamic viscosity defined'
285  write(*,*) ' dvi= ',dvi
286  call exit(201)
287  endif
288 !
289  reynolds=xflow360*d/(dvi*a)
290  if(reynolds.lt.1) then
291  reynolds=1.d0
292  endif
293 !
294 ! calculation of the friction coefficient
295 !
296  if(xflow_oil.ne.0d0) then
297 !
298 ! two-phase-flow
299 !
300  call two_phase_flow(tt1,pt1,t1,pt2,xflow360,
301  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
302  & v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
303  & shcon,rhcon,ntmat_,mi)
304  lambda=lambda*phi
305  else
306 !
307 ! for pure air
308 !
309  phi=1.d0
310  call friction_coefficient(dl,d,ks,reynolds,
311  & form_fact,lambda)
312  endif
313 !
314 ! calculating the critical conditions
315 !
316  call pt2zpt1_crit(pt2,pt1,tt1,lambda,kappa,r,dl,d,
317  & pt2zpt1_c,qred_crit,crit,icase,m1_c)
318 !
319  qred1=xflow360*dsqrt(tt1)/(a*pt1)
320 !
321 ! check whether flow is critical
322 ! assigning the physical correct sign to xflow360
323 !
324  if(crit) then
325 !
326 ! critical flow
327 !
328  xflow360=inv*qred_crit*a*pt1/dsqrt(tt1)
329  m1=dsqrt(2/km1*((tt1/t1)-1.d0))
330 c M1=min(M1,0.999d0/dsqrt(kappa))
331  else
332 !
333 ! subcritical flow
334 !
335  if(qred1.gt.qred_crit) then
336  xflow360=inv*qred_crit*a*pt1/dsqrt(tt1)
337  m1=m1_c
338  else
339  xflow360=inv*xflow360
340  m1=dsqrt(2/km1*((tt1/t1)-1.d0))
341  endif
342  call ts_calc(xflow360,tt2,pt2,kappa,r,a,t2,icase)
343  m2=dsqrt(2/km1*((tt2/t2)-1.d0))
344  endif
345 !
346  bb=km1/2.d0
347  cc=-kp1/(2.d0*km1)
348  ff1=(2.d0*bb*m1**2)/(1.d0+bb*m1**2*(1.d0+2.d0*cc))
349 !
350  if(.not.crit) then
351  ff2=(2.d0*bb*m2**2)/(1.d0+bb*m2**2*(1.d0+2.d0*cc))
352  if(idofp1.ne.0) then
353  ac(ieq,idofp1)=ff1*t1/pt1
354  endif
355  if(idoft1.ne.0) then
356  ac(ieq,idoft1)=(1.d0-ff1/2.d0)/(1.d0+bb*m1**2)
357  endif
358  if(idofm.ne.0) then
359  ac(ieq,idofm)=(ff2*t2-ff1*t1)/xflow360
360  endif
361  if(idofp2.ne.0) then
362  ac(ieq,idofp2)=-ff2*t2/pt2
363  endif
364  if(idoft2.ne.0) then
365  ac(ieq,idoft2)=(ff2/2.d0-1.d0)/(1.d0+bb*m2**2)
366  endif
367  else
368  if(idofp1.ne.0) then
369  ac(ieq,idofp1)=ff1*t1/pt1
370  endif
371  if(idoft1.ne.0) then
372  ac(ieq,idoft1)=(1.d0-ff1/2.d0)/(1.d0+bb*m1**2)
373  endif
374  if(idofm.ne.0) then
375  ac(ieq,idofm)=-ff1*t1/xflow360
376  endif
377  if(idoft2.ne.0) then
378  ac(ieq,idoft2)=-1.d0/(1.d0+bb/kappa)
379  endif
380  endif
381 !
382  endif
383  endif
384 !
385 ! mass equation contribution node1
386 !
387  if (nacteq(1,node1).ne.0) then
388  ieq=nacteq(1,node1)
389  if (idofm.ne.0) then
390  ac(ieq,idofm)=1.d0
391  endif
392  endif
393  endif
394 !
395  if(node2.ne.0) then
396 !
397 ! energy equation contribution node2
398 !
399  if(nacteq(0,node2).ne.0) then
400  ieq=nacteq(0,node2)
401  if((xflow.ge.0d0).and.(nacteq(3,node2).eq.0))then
402 !
403 ! adiabatic element
404 !
405  if(idoft1.ne.0)then
406  ac(ieq,idoft1)=ac(ieq,idoft1)-cp*xflow
407  elseif(node1.eq.0) then
408  write(*,*) '*WARNING in mafillnet'
409  write(*,*) ' temperature at inlet'
410  write(*,*) ' node',node2,' unknown'
411 c call exit(201)
412  endif
413 !
414  if(idoft2.ne.0) then
415  ac(ieq,idoft2)=ac(ieq,idoft2)+cp*xflow
416  endif
417 !
418  if(idofm.ne.0) then
419  ac(ieq,idofm)=ac(ieq,idofm)+cp*(tg2-tg1)
420  endif
421 !
422  elseif((node1.ne.0).and.(nacteq(3,node2).eq.node1))then
423 !
424 ! isothermal element
425 !
426  pt2zpt1=pt2/pt1
427 !
428 ! calculation of the dynamic viscosity
429 !
430  if(dabs(dvi).lt.1e-30) then
431  write(*,*) '*ERROR in gaspipe_fanno: '
432  write(*,*) ' no dynamic viscosity defined'
433  write(*,*) ' dvi= ',dvi
434  call exit(201)
435  endif
436 !
437  reynolds=xflow360*d/(dvi*a)
438  if(reynolds.lt.1) then
439  reynolds=1.d0
440  endif
441 !
442 ! calculation of the friction coefficient
443 !
444  if(xflow_oil.ne.0d0) then
445 !
446 ! two-phase-flow
447 !
448  call two_phase_flow(tt1,pt1,t1,pt2,xflow360,
449  & xflow_oil,nelem,lakon,kon,ipkon,ielprop,prop,
450  & v,dvi,cp,r,k_oil,phi,lambda,nshcon,nrhcon,
451  & shcon,rhcon,ntmat_,mi)
452  lambda=lambda*phi
453  else
454 !
455 ! for pure air
456 !
457  phi=1.d0
458  call friction_coefficient(dl,d,ks,reynolds,
459  & form_fact,lambda)
460  endif
461 !
462 ! calculating the critical conditions
463 !
464  call pt2zpt1_crit(pt2,pt1,tt1,lambda,kappa,r,dl,d,
465  & pt2zpt1_c,qred_crit,crit,icase,m1_c)
466 !
467  qred1=xflow360*dsqrt(tt1)/(a*pt1)
468 !
469 ! check whether flow is critical
470 ! assigning the physical correct sign to xflow360
471 !
472  if(crit) then
473 !
474 ! critical flow
475 !
476  xflow360=inv*qred_crit*a*pt1/dsqrt(tt1)
477  m1=dsqrt(2/km1*((tt1/t1)-1.d0))
478 c M1=min(M1,0.999d0/dsqrt(kappa))
479  else
480 !
481 ! subcritical flow
482 !
483  if(qred1.gt.qred_crit) then
484  xflow360=inv*qred_crit*a*pt1/dsqrt(tt1)
485  m1=m1_c
486  else
487  xflow360=inv*xflow360
488  m1=dsqrt(2/km1*((tt1/t1)-1.d0))
489  endif
490  call ts_calc(xflow360,tt2,pt2,kappa,r,a,t2,icase)
491  m2=dsqrt(2/km1*((tt2/t2)-1.d0))
492  endif
493 !
494  bb=km1/2.d0
495  cc=-kp1/(2.d0*km1)
496  ff1=(2.d0*bb*m1**2)/(1.d0+bb*m1**2*(1.d0+2.d0*cc))
497 !
498  if(.not.crit) then
499  ff2=(2.d0*bb*m2**2)/(1.d0+bb*m2**2*(1.d0+2.d0*cc))
500  if(idofp1.ne.0) then
501  ac(ieq,idofp1)=ff1*t1/pt1
502  endif
503  if(idoft1.ne.0) then
504  ac(ieq,idoft1)=(1.d0-ff1/2.d0)/(1.d0+bb*m1**2)
505  endif
506  if(idofm.ne.0) then
507  ac(ieq,idofm)=(ff2*t2-ff1*t1)/xflow360
508  endif
509  if(idofp2.ne.0) then
510  ac(ieq,idofp2)=-ff2*t2/pt2
511  endif
512  if(idoft2.ne.0) then
513  ac(ieq,idoft2)=(ff2/2.d0-1.d0)/(1.d0+bb*m2**2)
514  endif
515  else
516  if(idofp1.ne.0) then
517  ac(ieq,idofp1)=ff1*t1/pt1
518  endif
519  if(idoft1.ne.0) then
520  ac(ieq,idoft1)=(1.d0-ff1/2.d0)/(1.d0+bb*m1**2)
521  endif
522  if(idofm.ne.0) then
523  ac(ieq,idofm)=-ff1*t1/xflow360
524  endif
525  if(idoft2.ne.0) then
526  ac(ieq,idoft2)=-1.d0/(1.d0+bb/kappa)
527  endif
528  endif
529 !
530  endif
531  endif
532 !
533 ! mass equation contribution node2
534 !
535  if (nacteq(1,node2).ne.0) then
536  ieq=nacteq(1,node2)
537  if(idofm.ne.0)then
538  ac(ieq,idofm)=-1.d0
539  endif
540  endif
541  endif
542 !
543 ! element equation
544 !
545  if (nacteq(2,nodem).ne.0) then
546  ieq=nacteq(2,nodem)
547 !
548 ! for liquids: determine the gravity vector
549 !
550  if(lakon(nelem)(2:3).eq.'LI') then
551  do j=1,3
552  g(j)=0.d0
553  enddo
554  if(nbody.gt.0) then
555  index=nelem
556  do
557  j=ipobody(1,index)
558  if(j.eq.0) exit
559  if(ibody(1,j).eq.2) then
560  g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j)
561  g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j)
562  g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j)
563  endif
564  index=ipobody(2,index)
565  if(index.eq.0) exit
566  enddo
567  endif
568  endif
569 !
570  call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
571  & nactdog,identity,ielprop,prop,kflag,v,xflow,f,
572  & nodef,idirf,df,cp,r,rho,physcon,g,co,dvi,numf,
573  & vold,set,shcon,nshcon,rhcon,nrhcon,ntmat_,mi,ider,
574  & ttime,time,iaxial)
575 !
576  do k=1,numf
577  idof=nactdog(idirf(k),nodef(k))
578  if(idof.ne.0)then
579  ac(ieq,idof)=df(k)
580  endif
581  enddo
582  endif
583  enddo
584 !
585 ! convection with the walls
586 !
587  do i=1,nload
588  if(sideload(i)(3:4).eq.'FC') then
589  nelem=nelemload(1,i)
590  index=ipkon(nelem)
591  if(index.lt.0) cycle
592  lakonl=lakon(nelem)
593  node=nelemload(2,i)
594  ieq=nacteq(0,node)
595  if(ieq.eq.0) then
596  cycle
597  endif
598 !
599  call nident(itg,node,ntg,id)
600 !
601 ! calculate the area
602 !
603  read(sideload(i)(2:2),'(i1)') ig
604 !
605 ! number of nodes and integration points in the face
606 !
607  if(lakonl(4:4).eq.'2') then
608  nope=20
609  nopes=8
610  elseif(lakonl(4:4).eq.'8') then
611  nope=8
612  nopes=4
613  elseif(lakonl(4:5).eq.'10') then
614  nope=10
615  nopes=6
616  elseif(lakonl(4:4).eq.'4') then
617  nope=4
618  nopes=3
619  elseif(lakonl(4:5).eq.'15') then
620  nope=15
621  elseif(lakonl(4:4).eq.'6') then
622  nope=6
623  else
624  cycle
625  endif
626 !
627  if(lakonl(4:5).eq.'8R') then
628  mint2d=1
629  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
630  & then
631  if(lakonl(7:7).eq.'A') then
632  mint2d=2
633  else
634  mint2d=4
635  endif
636  elseif(lakonl(4:4).eq.'2') then
637  mint2d=9
638  elseif(lakonl(4:5).eq.'10') then
639  mint2d=3
640  elseif(lakonl(4:4).eq.'4') then
641  mint2d=1
642  endif
643 !
644  if(lakonl(4:4).eq.'6') then
645  mint2d=1
646  if(ig.le.2) then
647  nopes=3
648  else
649  nopes=4
650  endif
651  endif
652  if(lakonl(4:5).eq.'15') then
653  if(ig.le.2) then
654  mint2d=3
655  nopes=6
656  else
657  mint2d=4
658  nopes=8
659  endif
660  endif
661 !
662 ! connectivity of the element
663 !
664  index=ipkon(nelem)
665  if(index.lt.0) then
666  write(*,*) '*ERROR in mafillnet: element ',nelem
667  write(*,*) ' is not defined'
668  call exit(201)
669  endif
670  do k=1,nope
671  konl(k)=kon(index+k)
672  enddo
673 !
674 ! coordinates of the nodes belonging to the face
675 !
676  if((nope.eq.20).or.(nope.eq.8)) then
677  do k=1,nopes
678  tl2(k)=v(0,konl(ifaceq(k,ig)))
679  do j=1,3
680  xl2(j,k)=co(j,konl(ifaceq(k,ig)))+
681  & v(j,konl(ifaceq(k,ig)))
682  enddo
683  enddo
684  elseif((nope.eq.10).or.(nope.eq.4)) then
685  do k=1,nopes
686  tl2(k)=v(0,konl(ifacet(k,ig)))
687  do j=1,3
688  xl2(j,k)=co(j,konl(ifacet(k,ig)))+
689  & v(j,konl(ifacet(k,ig)))
690  enddo
691  enddo
692  else
693  do k=1,nopes
694  tl2(k)=v(0,konl(ifacew(k,ig)))
695  do j=1,3
696  xl2(j,k)=co(j,konl(ifacew(k,ig)))+
697  & v(j,konl(ifacew(k,ig)))
698  enddo
699  enddo
700  endif
701 !
702 ! integration to obtain the area and the mean
703 ! temperature
704 !
705  do l=1,mint2d
706  if((lakonl(4:5).eq.'8R').or.
707  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
708  xi=gauss2d1(1,l)
709  et=gauss2d1(2,l)
710  weight=weight2d1(l)
711  elseif((lakonl(4:4).eq.'8').or.
712  & (lakonl(4:6).eq.'20R').or.
713  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
714  xi=gauss2d2(1,l)
715  et=gauss2d2(2,l)
716  weight=weight2d2(l)
717  elseif(lakonl(4:4).eq.'2') then
718  xi=gauss2d3(1,l)
719  et=gauss2d3(2,l)
720  weight=weight2d3(l)
721  elseif((lakonl(4:5).eq.'10').or.
722  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
723  xi=gauss2d5(1,l)
724  et=gauss2d5(2,l)
725  weight=weight2d5(l)
726  elseif((lakonl(4:4).eq.'4').or.
727  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
728  xi=gauss2d4(1,l)
729  et=gauss2d4(2,l)
730  weight=weight2d4(l)
731  endif
732 !
733  if(nopes.eq.8) then
734  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
735  elseif(nopes.eq.4) then
736  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
737  elseif(nopes.eq.6) then
738  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
739  else
740  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
741  endif
742 !
743  dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
744  & xsj2(3)*xsj2(3))
745  areaj=dxsj2*weight
746 !
747  temp=0.d0
748  do k=1,3
749  coords(k)=0.d0
750  enddo
751  do j=1,nopes
752  temp=temp+tl2(j)*shp2(4,j)
753  do k=1,3
754  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
755  enddo
756  enddo
757 !
758  if(sideload(i)(5:6).ne.'NU') then
759  h(1)=xloadact(1,i)
760  else
761  read(sideload(i)(2:2),'(i1)') jltyp
762  jltyp=jltyp+10
763  sinktemp=v(0,node)
764  call film(h,sinktemp,temp,istep,
765  & iinc,tvar,nelem,l,coords,jltyp,field,nfield,
766  & sideload(i),node,areaj,v,mi,ipkon,kon,lakon,
767  & iponoel,inoel,ielprop,prop,ielmat,shcon,nshcon,
768  & rhcon,nrhcon,ntmat_,cocon,ncocon,
769  & ipobody,xbodyact,ibody,heatnod,heatfac)
770 c if(nmethod.eq.1) h(1)=xloadold(1,i)+
771 c & (h(1)-xloadold(1,i))*reltime
772  endif
773 !
774  idoft=nactdog(0,node)
775  if(idoft.gt.0) then
776  if(lakonl(5:7).eq.'0RA') then
777  ac(ieq,idoft)=ac(ieq,idoft)+2.d0*h(1)*dxsj2*weight
778  else
779  ac(ieq,idoft)=ac(ieq,idoft)+h(1)*dxsj2*weight
780  endif
781  endif
782  enddo
783  endif
784  enddo
785 !
786 ! additional multiple point constraints
787 !
788  j=nteq+1
789  do i=nmpc,1,-1
790  if(labmpc(i)(1:7).ne.'NETWORK') cycle
791  j=j-1
792  if(labmpc(i)(8:8).eq.' ') then
793  index=ipompc(i)
794 !
795 ! linear equation
796 !
797  do
798  node=nodempc(1,index)
799  idir=nodempc(2,index)
800  if(nactdog(idir,node).ne.0) then
801  ac(j,nactdog(idir,node))=coefmpc(index)
802  endif
803  index=nodempc(3,index)
804  if(index.eq.0) exit
805  enddo
806  else
807 !
808 ! user-defined network equation
809 !
810  call networkmpc_lhs(i,ipompc,nodempc,coefmpc,labmpc,
811  & v,nactdog,ac,j,mi,nteq)
812  endif
813  enddo
814 !
815 c write(*,*) nteq,' mafillnet '
816 c do i=1,nteq
817 c write(*,'(17(1x,e11.4))') (ac(i,j),j=1,nteq)
818 c enddo
819 !
820  return
pt2zpt1_crit
subroutine pt2zpt1_crit(pt2, pt1, Tt1, lambda, kappa, r, l, d, pt2zpt1_c, Qred1_crit, crit, icase, M1)
Definition: pt2zpt1_crit.f:35
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
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
ts_calc
subroutine ts_calc(xflow, Tt, Pt, kappa, r, A, Ts, icase)
Definition: ts_calc.f:20
friction_coefficient
subroutine friction_coefficient(l, d, ks, reynolds, form_fact, lambda)
Definition: friction_coefficient.f:26
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
networkmpc_lhs
subroutine networkmpc_lhs(i, ipompc, nodempc, coefmpc, labmpc, v, nactdog, ac, j, mi, nteq)
Definition: networkmpc_lhs.f:23
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
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
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
film
subroutine film(h, sink, temp, kstep, kinc, time, noel, npt, coords, jltyp, field, nfield, loadtype, node, area, vold, mi, ipkon, kon, lakon, iponoel, inoel, ielprop, prop, ielmat, shcon, nshcon, rhcon, nrhcon, ntmat_, cocon, ncocon, ipobody, xbody, ibody, heatnod, heatfac)
Definition: film.f:24
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)