1!
2!     CalculiX - A 3-dimensional finite element program
3!     Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine massflow_percent(node1,node2,nodem,nelem,lakon,kon,
20     &        ipkon,nactdog,identity,ielprop,prop,iflag,v,xflow,f,
21     &        nodef,idirf,df,cp,r,physcon,dvi,numf,set,shcon,
22     &        nshcon,rhcon,nrhcon,ntmat_,co,vold,mi,ttime,time,
23     &        iaxial,iplausi)
24!
25!     partial massflow  element
26!
27!     author: Yannick Muller
28!
29      implicit none
30!
31      logical identity
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,
37     &     inv,ipkon(*),kon(*),number,kgas,iaxial,
38     &     nodea,nodeb,mi(*),i,itype,nodemup,
39     &     nrhcon(*),ntmat_,nshcon(*),iplausi
40!
41      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),kappa,R,a,d,xl,
42     &     p1,p2,T1,physcon(*),pi,xflow_oil,T2,co(3,*),vold(0:mi(2),*),
43     &     xflow_sum,percent_xflow,cp,dvi,pt1,pt2,Tt1,Tt2,ttime,time,
44     &     shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*)
45!
46!
47!
48      pi=4.d0*datan(1.d0)
49      index=ielprop(nelem)
50!
51      if(iflag.eq.0) then
52         identity=.true.
53!
54         if(nactdog(2,node1).ne.0)then
55            identity=.false.
56         elseif(nactdog(2,node2).ne.0)then
57            identity=.false.
58         elseif(nactdog(1,nodem).ne.0)then
59            identity=.false.
60         endif
61!
62      elseif(iflag.eq.1)then
63         if(v(1,nodem).ne.0.d0) then
64            xflow=v(1,nodem)
65            return
66         endif
67!
68         percent_xflow=prop(index+1)
69         xflow_sum=0
70!
71         do i=2,10
72            if(nint(prop(index+i)).ne.0) then
73               nodemup=kon(ipkon(nint(prop(index+i)))+2)
74               if(v(1,nodemup).gt.0)then
75                  xflow_sum=xflow_sum+v(1,nodemup)*iaxial
76               endif
77            endif
78         enddo
79!
80         if(xflow_sum.eq.0d0) then
81            xflow_sum=0.001d0
82         endif
83!
84         xflow=xflow_sum*percent_xflow
85!
86      elseif((iflag.eq.2).or.(iflag.eq.3))then
87!
88         percent_xflow=prop(index+1)
89         xflow_sum=0
90         do i=2,10
91            if(nint(prop(index+i)).ne.0) then
92               nodemup=kon(ipkon(nint(prop(index+i)))+2)
93               if(v(1,nodemup).gt.0)then
94                  xflow_sum=xflow_sum+v(1,nodemup)*iaxial
95               endif
96            endif
97         enddo
98
99         if(xflow_sum.eq.0.d0) then
100            xflow_sum=1d-5
101         endif
102!
103         inv=1
104!
105         pt1=v(2,node1)
106         pt2=v(2,node2)
107         xflow=v(1,nodem)*iaxial
108         Tt1=v(0,node1)-physcon(1)
109         Tt2=v(0,node2)-physcon(1)
110!
111         nodef(1)=node1
112         nodef(2)=node1
113         nodef(3)=nodem
114         nodef(4)=node2
115!
116         idirf(1)=2
117         idirf(2)=0
118         idirf(3)=1
119         idirf(4)=2
120!
121         if(iflag.eq.2) then
122            numf=4
123!
124            f=xflow/xflow_sum-percent_xflow
125!
126            df(1)=0
127            df(2)=0
128            df(3)=1/xflow_sum
129            df(4)=0
130!
131!     output
132!
133         elseif(iflag.eq.3)then
134!
135            xflow_oil=0
136!
137            write(1,*) ''
138            write(1,55) ' from node ',node1,
139     &           ' to node ', node2,' :   air massflow rate = '
140     &           ,inv*xflow,
141     &           ', oil massflow rate = ',xflow_oil
142 55         format(1X,A,I6,A,I6,A,e11.4,A,A,e11.4,A)
143!
144            write(1,56)'       Inlet node ',node1,' :   Tt1 = ',Tt1,
145     &           '  , Ts1 = ',Tt1,'  , Pt1 = ',pt1
146!
147            write(1,*)'             Element ',nelem,lakon(nelem)
148            write(1,57)'        Massflow upstream = ',xflow_sum,
149     &        ' [kg/s]'
150            write(1,58)'        Massflow fraction = ', percent_xflow
151            write(1,56)'      Outlet node ',node2,':    Tt2=',Tt2,
152     &           ', Ts2=',Tt2,', Pt2=',pt2
153!
154         endif
155      endif
156!
157 56   format(1X,A,I6,A,e11.4,A,e11.4,A,e11.4,A)
158 57   format(1X,A,e11.4,A)
159 58   format(1X,A,e11.4)
160!
161      xflow=xflow/iaxial
162      df(3)=df(3)*iaxial
163!
164      return
165      end
166