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