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