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 characteristic(node1,node2,nodem,nelem,lakon, 20 & kon,ipkon, 21 & nactdog,identity,ielprop,prop,iflag,v,xflow,f, 22 & nodef,idirf,df,cp,r,physcon,dvi,numf,set, 23 & mi,ttime,time,iaxial,iplausi) 24! 25! This subroutine is used to enables the processing of empiric 26! given under the form 27! massflow*dsqrt(T1)/Pt1=f((Pt1-Pt2)/Pt1) and T1=T2 28! characteristics the subroutine proceeds using 29! linear interpolation to estimate the values for the whole characteristic 30! note that the characteristic is implicitely containing the point (0,0) 31! 32! author: Yannick Muller 33! 34 implicit none 35! 36 logical identity 37! 38 character*8 lakon(*) 39 character*81 set(*) 40! 41 integer nelem,nactdog(0:3,*),node1,node2,nodem,kon(*),ipkon(*), 42 & ielprop(*),nodef(*),idirf(*),index,iflag, 43 & inv,id,numf,npu,i,mi(*),iaxial,iplausi 44! 45 real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),cp,r,dvi, 46 & p1,p2,physcon(*),ttime,time,xmach,kappa, 47 & xpu(100),ypu(100),Qred,p1mp2zp1,T1,scal,T2 48! 49! 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).or.(iflag.eq.2)) then 63 if(iflag.eq.1) then 64 if(v(1,nodem).ne.0.d0) then 65 xflow=v(1,nodem) 66 return 67 endif 68 endif 69! 70 index=ielprop(nelem) 71! 72 npu=nint(prop(index+2)) 73 scal=prop(index+1) 74 75 do i=1, 100 76 xpu(i)=0 77 ypu(i)=0 78 enddo 79! 80 do i=1,npu 81 xpu(i)=prop(index+2*i+1) 82 ypu(i)=prop(index+2*i+2) 83 enddo 84! 85 p1=v(2,node1) 86 p2=v(2,node2) 87! 88 if(p1.ge.p2) then 89 inv=1 90 T1=v(0,node1)-physcon(1) 91 else 92 inv=-1 93 p1=v(2,node2) 94 p2=v(2,node1) 95 T1=v(0,node2)-physcon(1) 96 endif 97! 98 p1mp2zp1=(p1-p2)/p1 99! 100 if(iflag.eq.1) then 101 102 call ident(xpu,p1mp2zp1,npu,id) 103 if(id.eq.0) then 104 Qred=scal*ypu(1) 105 xflow=inv*Qred*p1/dsqrt(T1) 106 elseif(id.ge.npu) then 107 Qred=scal*ypu(npu) 108 xflow=inv*Qred*p1/dsqrt(T1) 109 else 110 Qred=scal*(ypu(id)+(ypu(id+1)-ypu(id)) 111 & *(p1mp2zp1-xpu(id))/(xpu(id+1)-xpu(id))) 112 xflow=inv*Qred*p1/dsqrt(T1) 113 endif 114! 115 elseif (iflag.eq.2) then 116 numf=4 117! 118 p1=v(2,node1) 119 p2=v(2,node2) 120 xflow=v(1,nodem)*iaxial 121! 122 if (p1.ge.p2) then 123! 124 inv=1 125 xflow=v(1,nodem)*iaxial 126 T1=v(0,node1)-physcon(1) 127 nodef(1)=node1 128 nodef(2)=node1 129 nodef(3)=nodem 130 nodef(4)=node2 131! 132 else 133! 134 inv=-1 135 p1=v(2,node2) 136 p2=v(2,node1) 137 T1=v(0,node2)-physcon(1) 138 xflow=-v(1,nodem)*iaxial 139 nodef(1)=node2 140 nodef(2)=node2 141 nodef(3)=nodem 142 nodef(4)=node1 143 endif 144! 145 idirf(1)=2 146 idirf(2)=0 147 idirf(3)=1 148 idirf(4)=2 149! 150 df(2)=xflow/(2.d0*p1*dsqrt(T1)) 151 df(3)=inv*dsqrt(T1)/p1 152! 153 call ident(xpu,p1mp2zp1,npu,id) 154! 155 if(id.eq.0) then 156 f=dabs(xflow)/p1*dsqrt(T1)-scal*ypu(1) 157 df(4)=0.01d0 158 df(1)=-xflow*dsqrt(T1)/p1**2 159! 160 elseif(id.ge.npu) then 161 f=dabs(xflow)/p1*dsqrt(T1)-scal*ypu(npu) 162 df(4)=0.01d0 163 df(1)=-xflow*dsqrt(T1)/p1**2 164! 165 else 166 f=dabs(xflow)/p1*dsqrt(T1)-(scal*ypu(id) 167 & +scal*(ypu(id+1)-ypu(id)) 168 & *(p1mp2zp1-xpu(id))/(xpu(id+1)-xpu(id))) 169! 170 df(4)=scal*(ypu(id+1)-ypu(id))/(xpu(id+1)-xpu(id))*1/p1 171! 172 df(1)=-xflow*dsqrt(T1)/p1**2-(p2/p1**2) 173 & *(scal*(ypu(id+1)-ypu(id))/(xpu(id+1)-xpu(id))) 174 endif 175 endif 176 177 elseif(iflag.eq.3) then 178 p1=v(2,node1) 179 p2=v(2,node2) 180 xflow=v(1,nodem)*iaxial 181 kappa=(cp/(cp-r)) 182 xmach=dsqrt(((p1/p2)**((kappa-1.d0)/kappa)-1.d0)*2.d0/ 183 & (kappa-1.d0)) 184! 185 if (p1.ge.p2) then 186! 187 inv=1 188 xflow=v(1,nodem)*iaxial 189 T1=v(0,node1)-physcon(1) 190 T2=v(0,node2)-physcon(1) 191 nodef(1)=node1 192 nodef(2)=node1 193 nodef(3)=nodem 194 nodef(4)=node2 195! 196 else 197! 198 inv=-1 199 p1=v(2,node2) 200 p2=v(2,node1) 201 T1=v(0,node2)-physcon(1) 202 T2=v(0,node1)-physcon(1) 203 xflow=-v(1,nodem)*iaxial 204 nodef(1)=node2 205 nodef(2)=node2 206 nodef(3)=nodem 207 nodef(4)=node1 208 endif 209! 210 write(1,*) '' 211 write(1,55) ' from node',node1, 212 & ' to node', node2,': air massflow rate=',xflow 213! 214 55 FORMAT(1X,A,I6,A,I6,A,e11.4,A,A,e11.4,A) 215! 216 if(inv.eq.1) then 217 write(1,56)' Inlet node ',node1,': Tt1=',T1, 218 & ', Ts1=',T1,', Pt1=',p1 219 220 write(1,*)' Element ',nelem,lakon(nelem) 221 write(1,57) 'M = ',xmach 222! 223 write(1,56)' Outlet node ',node2,': Tt2=',T2, 224 & ', Ts2=',T2,', Pt2=',p2 225! 226 else if(inv.eq.-1) then 227 write(1,56)' Inlet node ',node2,': Tt1=',T1, 228 & ', Ts1=',T1,', Pt1=',p1 229 & 230 write(1,*)' Element ',nelem,lakon(nelem) 231 write(1,57) 'M = ',xmach 232! 233 write(1,56)' Outlet node ',node1,': Tt2=',T2, 234 & ', Ts2=',T2,', Pt2=',p2 235! 236 endif 237! 238 56 FORMAT(1X,A,I6,A,e11.4,A,e11.4,A,e11.4,A) 239 57 format(40x,a,e11.4) 240! 241 endif 242! 243 xflow=xflow/iaxial 244 df(3)=df(3)*iaxial 245! 246 return 247 end 248