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!     this subroutine solves iteratively the following equation
20!     to determine the pressure for which section A2 is critical
21!
22!     author: Yannick Muller
23!
24      subroutine pt2_lim_calc (pt1,a2,a1,kappa,zeta,pt2_lim)
25!
26      implicit none
27!
28      integer i
29!
30      real*8 pt1,a2,a1,kappa,pt2_lim,x,zeta,f,df,expon1,
31     &     expon2,expon3,cte,a2a1,kp1,km1,delta_x,fact1,fact2,term
32!
33      x=0.999d0
34!
35!     x belongs to interval [0;1]
36!
37!     modified 25.11.2007
38!     since Pt1/Pt2=(1+0.5(kappa)-M)**(zeta*kappa)/(kappa-1)
39!     and for zeta1 elements type M_crit=M1=1
40!     and for zeta2 elements type M_crit=M2 =1
41!     it is not necessary to iteratively solve the flow equation.
42!     Instead the previous equation is solved to find pt2_crit
43!
44      if(zeta.ge.0d0) then
45         kp1=kappa+1.d0
46         km1=kappa-1.d0
47         a2a1=a2/a1
48         expon1=-0.5d0*kp1/(zeta*kappa)
49         expon2=-0.5d0*kp1/km1
50         cte=a2a1*(0.5d0*kp1)**expon2
51         expon3=-km1/(zeta*kappa)
52         i=0
53!
54!
55         do
56            i=i+1
57!
58            f=x**(-1.d0)-cte*x**(expon1)
59     &           *(2.d0/km1*(x**expon3-1.d0))**(-0.5d0)
60!
61            df=-1.d0/X**2-cte*(x**expon1
62     &           *(2.d0/km1*(x**expon3-1.d0))**(-0.5d0))
63     &           *(expon1/X-1.d0/km1*expon3*x**(expon3-1.d0)
64     &           *(2.d0/km1*(x**expon3-1.d0))**(-1.d0))
65
66            delta_x=-f/df
67!
68            if(( dabs(delta_x/x).le.1.d-8)
69     &           .or.(dabs(delta_x/1d0).le.1.d-10)) then
70!
71               pt2_lim=pt1*X
72!
73               exit
74            endif
75            if(i.gt.25)then
76                pt2_lim=pt1/(1.d0+0.5d0*km1)**(zeta*kappa/km1)
77                exit
78             endif
79!
80            x=delta_x+x
81!
82         enddo
83!
84      else
85!
86         do
87            kp1=kappa+1.d0
88            km1=kappa-1.d0
89            a2a1=a2/a1
90            expon1=kp1/(zeta*kappa)
91            expon2=km1/(zeta*kappa)
92            expon3=kp1/km1
93            cte=a2a1**2*(0.5d0*kp1)**(-expon3)*(2.d0/km1)**(-1.d0)
94            fact1=x**(-expon1)
95            fact2=x**(-expon2)
96            term=fact2-1.d0
97!
98            f=x**(-2.d0)-cte*fact1*term**(-1.d0)
99!
100            df=-2.d0*x**(-3.d0)-cte*(x**(-expon1-1.d0)*term**(-1.d0))
101     &           *(-expon1+expon2*(X**(-expon2))*fact2*term**(-1.d0))
102!
103            delta_x=-f/df
104!
105            if(( dabs(delta_x/x).le.1.d-8)
106     &           .or.(dabs(delta_x/1d0).le.1.d-10)) then
107               pt2_lim=pt1*X
108               exit
109            endif
110!
111            x=delta_x+x
112!
113         enddo
114!
115      endif
116
117      return
118      end
119