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 carbon_seal(node1,node2,nodem,nelem,lakon,
20     &     nactdog,identity,ielprop,prop,iflag,v,xflow,f,
21     &     nodef,idirf,df,R,physcon,dvi,numf,set,mi,ttime,time,
22     &     iaxial,iplausi)
23!
24!     carbon seal element calculated with Richter method
25!      Richter "Rohrhydraulik", Springer ,1971,p. 175
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,mi(*),iaxial,iplausi
38!
39      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),R,d,dl,
40     &     p1,p2,T1,physcon(*),dvi,pi,s,T2,ttime,time
41!
42!
43!
44      if(iflag.eq.0) then
45         identity=.true.
46!
47         if(nactdog(2,node1).ne.0)then
48            identity=.false.
49         elseif(nactdog(2,node2).ne.0)then
50            identity=.false.
51         elseif(nactdog(1,nodem).ne.0)then
52            identity=.false.
53         endif
54!
55      elseif(iflag.eq.1)then
56         if(v(1,nodem).ne.0.d0) then
57            xflow=v(1,nodem)
58            return
59         endif
60!
61         index=ielprop(nelem)
62         d=prop(index+1)
63         s=prop(index+2)
64         dl=prop(index+3)
65         pi=4.d0*datan(1.d0)
66!
67         p1=v(2,node1)
68         p2=v(2,node2)
69         if(p1.ge.p2) then
70            inv=1
71            T1=v(0,node1)-physcon(1)
72         else
73            inv=-1
74            p1=v(2,node2)
75            p2=v(2,node1)
76            T1=v(0,node2)-physcon(1)
77         endif
78!
79         if(lakon(nelem)(2:6).eq.'CARBS') then
80!
81!     gapflow
82!     Richter "Rohrhydraulik", Springer ,1971,p. 175
83!
84            xflow=inv*Pi*d*s**3*(p1**2-p2**2)/(24.d0*R*T1*dvi*dl)
85
86         elseif(lakon(nelem)(2:6).ne.'CARBS') then
87            write(*,*) '*WARNING in Carbon_seal.f'
88            write(*,*) 'unable to perform carbon seal calculation'
89            write(*,*) 'check input file'
90         endif
91!
92      elseif(iflag.eq.2)then
93!
94         numf=4
95         p1=v(2,node1)
96         p2=v(2,node2)
97         if(p1.ge.p2) then
98            inv=1
99            xflow=v(1,nodem)*iaxial
100            T1=v(0,node1)-physcon(1)
101            nodef(1)=node1
102            nodef(2)=node1
103            nodef(3)=nodem
104            nodef(4)=node2
105         else
106            inv=-1
107            p1=v(2,node2)
108            p2=v(2,node1)
109            xflow=-v(1,nodem)*iaxial
110            T1=v(0,node2)-physcon(1)
111            nodef(1)=node2
112            nodef(2)=node2
113            nodef(3)=nodem
114            nodef(4)=node1
115         endif
116!
117         idirf(1)=2
118         idirf(2)=0
119         idirf(3)=1
120         idirf(4)=2
121!
122         index=ielprop(nelem)
123         d=prop(index+1)
124         s=prop(index+2)
125         dl=prop(index+3)
126         pi=4.d0*datan(1.d0)
127
128!
129         if(lakon(nelem)(2:6).eq.'CARBS') then
130!
131            f=xflow*T1-pi*d*s**3*(p1**2-p2**2)/(24.d0*R*dvi*dl)
132!
133            df(1)=-(pi*d*s**3*p1)/(12.d0*R*dvi*dl)
134            df(2)=xflow
135            df(3)=T1
136            df(4)=(pi*d*s**3*p2)/(12.d0*R*dvi*dl)
137!
138         endif
139
140      elseif(iflag.eq.3) then
141         p1=v(2,node1)
142         p2=v(2,node2)
143         if(p1.ge.p2) then
144            inv=1
145            xflow=v(1,nodem)*iaxial
146            T1=v(0,node1)-physcon(1)
147            T2=v(0,node2)-physcon(1)
148            nodef(1)=node1
149            nodef(2)=node1
150            nodef(3)=nodem
151            nodef(4)=node2
152         else
153            inv=-1
154            p1=v(2,node2)
155            p2=v(2,node1)
156            xflow=-v(1,nodem)*iaxial
157            T1=v(0,node2)-physcon(1)
158            T2=v(0,node1)-physcon(1)
159            nodef(1)=node2
160            nodef(2)=node2
161            nodef(3)=nodem
162            nodef(4)=node1
163         endif
164
165         write(1,*) ''
166         write(1,55) ' from node',node1,
167     &' to node', node2,':   air massflow rate=',xflow
168 55      FORMAT(1X,A,I6,A,I6,A,e11.4,A,A,e11.4,A)
169
170         if(inv.eq.1) then
171            write(1,56)'       Inlet node ',node1,':   Tt1=',T1,
172     &           ' , Ts1=',T1,' , Pt1=',p1
173
174            write(1,*)'             Element',nelem,lakon(nelem)
175
176            write(1,56)'      Outlet node ',node2,':   Tt2=',T2,
177     &           ' , Ts2=',T2,' , Pt2=',p2
178!
179         else if(inv.eq.-1) then
180            write(1,56)'       Inlet node ',node2,':    Tt1=',T1,
181     &           ' , Ts1=',T1,' , Pt1=',p1
182     &
183            write(1,*)'             Element',nelem,lakon(nelem)
184
185            write(1,56)'      Outlet node ',node1,':    Tt2=',T2,
186     &           ' , Ts2=',T2,' , Pt2=',p2
187
188         endif
189
190 56      FORMAT(1X,A,I6,A,e11.4,A,e11.4,A,e11.4,A)
191      endif
192!
193      xflow=xflow/iaxial
194      df(3)=df(3)*iaxial
195!
196      return
197      end
198
199
200