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