Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ZN_ssWSB.f90
Go to the documentation of this file.
1 ! *********************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4 ! * *
5 ! * Illinois Rocstar LLC *
6 ! * Champaign, IL *
7 ! * www.illinoisrocstar.com *
8 ! * sales@illinoisrocstar.com *
9 ! * *
10 ! * License: See LICENSE file in top level of distribution package or *
11 ! * http://opensource.org/licenses/NCSA *
12 ! *********************************************************************
13 ! *********************************************************************
14 ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15 ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16 ! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17 ! * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18 ! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19 ! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20 ! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21 ! * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22 ! *********************************************************************
23 ! ---------------------------------------------------------------------------
24 !
25 ! SUBROUTINE : ZN_ssWSB
26 !
27 ! purpose : calculate steady state solution
28 !
29 ! Authors : K.C. Tang and M. Q. Brewster
30 !
31 ! Creation Date : Sep. 15, 2000
32 !
33 ! Modifications :
34 !
35 ! No. Date Authors Description
36 ! 1 09/03/02 K.C. Tang Global variables G_ZN
37 !
38 !
39 ! ---------------------------------------------------------------------------
40 !
41 ! arguments :
42 !
43 !
44 ! G_ZN : Global variables for Rocburn_1D_ZN
45 ! P : pressure (atm)
46 ! qr : radiant flux (cal/cm^2/s)
47 ! To : initial temperature (K)
48 ! rb : burning rate (cm/s)
49 ! Ts : surface temperature (K)
50 ! Tf : flame temperature
51 ! fr : fraction of radiation absorbed below surface reaction
52 ! zone
53 !
54 ! ---------------------------------------------------------------------------
55 !
56  SUBROUTINE zn_sswsb(G_ZN, P, qr, To, rhoc, rb, Ts, Tf, fr, Tn)
57 
58  USE m_rocburn_zn_global_data
59 
60  IMPLICIT NONE
61  include 'mpif.h'
62 
63 
64  TYPE(g_burn_1d), POINTER :: g_zn
65  REAL(DBL), INTENT(IN) :: p, qr, to, rhoc
66  REAL(DBL), INTENT(OUT) :: rb, ts, tf, fr
67  REAL(DBL), INTENT(OUT) :: tn(:)
68 
69 
70  REAL(DBL), parameter :: tol=1.0d-5
71 
72 
73 !
74 ! ============================================================================
75 !
76 ! delcare local variables
77 !
78  REAL(DBL) :: j(2,2),xx(2),e(2)
79  REAL(DBL) :: xreac
80  REAL(DBL) :: term0, term1, term2, term3, term4, term5
81  REAL(DBL) :: xcd, da, xg, xcond, frj
82 
83  INTEGER :: iter, i
84 !
85 ! ============================================================================
86 !
87 
88 !
89 ! Generate initial conditions
90 !
91 !
92 !
93 
94  IF(g_zn%Model_combustion == 1) THEN
95 
96 ! Model_combustion = 1
97 !
98 ! steady state WSB homogeneouse propellant combustion model
99 !
100 ! guess initial values for doublebase propellant
101 ! rb and Ts
102 !
103 
104  rb=g_zn%a_p*(p**g_zn%n_p) ! cm/s
105  ts=550.0 ! K
106  xx(1)=rb*rhoc ! g/(cm^2*s)
107  xx(2)=ts ! K
108 
109  DO iter=0,g_zn%itermax
110 !
111 ! set up residual matrix and Jacobian derivative matrix for
112 ! Ibiricu and Williams condensed phase expression
113 !
114 
115  xreac=2.0*rhoc*g_zn%alfac*g_zn%R*xx(2)/xx(1)/g_zn%Ec ! cm
116  term0=1.0-g_zn%Ka*xreac ! -
117 
118  term1=g_zn%Ac*g_zn%R*(xx(2)*rhoc)**2*g_zn%alfac*g_zn%C* &
119  dexp(-g_zn%Ec/g_zn%R/xx(2))/g_zn%Ec ! cal*g/(cm^4*s^2)
120 
121  fr=dexp(-xreac*g_zn%Ka)
122  term2=g_zn%C*(xx(2)-to)-g_zn%Qc/2.0-fr*qr/xx(1) ! cal/g
123  term3=term1/term2 ! g^2/(cm^4*s^2)
124 
125  e(1)=xx(1)**2-term3 ! g^2/(cm^4*s^2)
126  j(1,1)=2.0*xx(1)+term3/term2*fr*qr/xx(1)**2 &
127  *term0 ! g/(cm^2*s)
128  j(1,2)=-xx(1)**2*((2.0d0+g_zn%Ec/g_zn%R/xx(2))/xx(2) - &
129  (g_zn%C+g_zn%Ka*xreac*qr*fr/xx(1)/xx(2))/term2) ! g^2/(cm^4*s^2*K)
130 
131 !
132 ! set up residual matrix and Jacobian derivative matrix for
133 ! Ward-Son-Brewster gas phase expression
134 
135  xcd=g_zn%lamg/g_zn%C/xx(1) ! cm
136  da=4.0*g_zn%Bg*g_zn%C/g_zn%lamg*(p*g_zn%MW*xcd/g_zn%R)**2 ! -
137  term4=(1.0+da)**0.5 ! -
138  xg=2.0*xcd/(term4-1.0) ! cm
139  term5=g_zn%Qc-g_zn%C*(xx(2)-to) ! cal/g
140  tf=(qr/xx(1)+g_zn%Qc+g_zn%Qg)/g_zn%C+to ! K
141 
142  e(2)=(tf-xx(2))/xg+(qr+xx(1)*term5)/g_zn%lamg ! K/cm
143  j(2,1)=-qr/xg/xx(1)**2/g_zn%C- &
144  (tf-xx(2))/xg/xx(1)/term4+term5/g_zn%lamg ! (K*cm*s)/g
145  j(2,2)=-1.0/xg-1.0/xcd ! 1/cm
146 
147 !
148 ! update guesses for m and Ts using Newton-Rhapson
149 ! & guard against a negative (erroneous) rb solution
150 !
151  CALL newton(xx,j,e)
152  IF(xx(1).LT.0.0) THEN
153 
154  print *,' ROCBURN_ZN: rank=',g_zn%rank
155  print *,' rb negative xx(1)=',xx(1)
156  xx(1)=0.0d0
157 
158  END IF
159 
160 
161  IF(abs(e(1)).LT.tol .AND. abs(e(2)).LT.tol) THEN
162 
163 !
164 ! reset final converged values
165 !
166  rb=xx(1)/rhoc ! cm/s
167  ts=xx(2) ! K
168  xreac=2.0*g_zn%alfac*g_zn%R*ts/rb/g_zn%Ec ! 1/cm
169  fr=dexp(-xreac*g_zn%Ka) ! -
170  xcd=g_zn%lamg/g_zn%C/xx(1) ! cm
171  da=4.0*g_zn%Bg*g_zn%C/g_zn%lamg*(p*g_zn%MW*xcd/g_zn%R)**2 ! -
172  term4=(1.0+da)**0.5 ! -
173  xg=2.0*xcd/(term4-1.0) ! cm
174  tf=(qr/xx(1)+g_zn%Qc+g_zn%Qg)/g_zn%C+to ! K
175 
176 !
177 ! Sets ICs using I&W S.S. T profile
178 !
179 
180  xcond=g_zn%alfac/rb
181  xreac=xcond*g_zn%R*ts*2.0/g_zn%Ec
182  fr=dexp(-xreac*g_zn%Ka)
183  frj=fr*qr/rhoc/rb/g_zn%C/(ts-to)/(1.-g_zn%Ka*xcond)
184  DO i=1,g_zn%nx
185  tn(i)=to+(ts-to)*((1.0-frj)* &
186  dexp(g_zn%x(i)/xcond)+ frj*exp(g_zn%x(i)*g_zn%Ka))
187  END DO ! i
188 
189  IF(abs(tn(g_zn%nx)-to)/to .GT. 0.001) THEN
190 
191  WRITE(*,*) 'ROCBURN_ZN: rank=',g_zn%rank
192  WRITE(*,*) ' Error: xmax not large enough'
193  WRITE(*,*) ' To boundary condition not satisfied'
194  CALL mpi_abort( mpi_comm_world, -1)
195  stop
196 
197  ELSE
198 
199  tn(g_zn%nx) =to
200 
201  END IF
202 
203  RETURN
204  END IF
205  END DO ! iter
206 
207  print *,' ROCBURN_ZN: rank=',g_zn%rank
208  print *,' ROCBURN_ZN: Error-steady state WSB' &
209  ,' in ignition_combustion'
210  print *,' ROCBURN_ZN: Maximum # of ',g_zn%itermax &
211  , ' iterations exceeded'
212  CALL mpi_abort( mpi_comm_world, -1)
213  stop
214 !
215 ! end of steady state WSB homogeneouse propellant combustion model
216 !
217  END IF
218 
219 
220  IF(g_zn%Model_combustion == 2) THEN
221 
222 ! Model_combustion = 2
223 !
224 ! ZN approach for composite propellant with empirical formulatin
225 ! for the gas phase reaction law
226 !
227 ! initial guess of rb and Ts
228 !
229 
230  rb=g_zn%a_p*p**g_zn%n_p ! cm/s
231  ts=727.002 ! K
232  xx(1)=rb*rhoc ! g/(cm^2*s)
233  xx(2)=ts ! K
234 
235  DO iter=0,g_zn%itermax
236 
237 !
238 ! set up residual matrix and Jacobian derivative matrix for
239 ! Ibiricu and Williams condensed phase expression
240 !
241  xreac=2.0*rhoc*g_zn%alfac*g_zn%R*xx(2)/xx(1)/g_zn%Ec ! cm
242  term0=1.0-g_zn%Ka*xreac ! -
243  term1=g_zn%Ac*g_zn%R*(xx(2)*rhoc)**2*g_zn%alfac*g_zn%C* &
244  dexp(-g_zn%Ec/g_zn%R/xx(2))/g_zn%Ec ! cal*g/(cm^4*s^2)
245 
246  fr=dexp(-xreac*g_zn%Ka)
247  term2=g_zn%C*(xx(2)-to)-g_zn%Qc/2.0-fr*qr/xx(1) ! cal/g
248  term3=term1/term2 ! g^2/(cm^4*s^2)
249 
250  e(1)=xx(1)**2-term3 ! g^2/(cm^4*s^2)
251  j(1,1)=2.0*xx(1)+term3/term2*fr*qr/xx(1)**2 &
252  *term0 ! g/(cm^2*s)
253  j(1,2)=-xx(1)**2*((2.0d0+g_zn%Ec/g_zn%R/xx(2))/xx(2) - &
254  (g_zn%C+g_zn%Ka*xreac*qr*fr/xx(1)/xx(2))/term2) ! g^2/(cm^4*s^2*K)
255  e(2)=g_zn%a_T*((xx(2) - to - qr/(xx(1)*g_zn%C) )/ &
256  (xx(2)-300.0))**g_zn%n_T - rhoc*g_zn%a_p*p**g_zn%n_p/xx(1)
257  j(2,1)=rhoc*g_zn%a_p*p**g_zn%n_p/(xx(1)*xx(1))
258  j(2,2)=g_zn%n_T*g_zn%a_T*((xx(2)-to)/(xx(2)-300.0))**(g_zn%n_T-1.0) &
259  *(to-300.0)/(xx(2)-300.0)**2
260 !
261 ! update guesses for m and Ts using Newton-Rhapson
262 ! & guard against a negative (erroneous) rb solution
263 !
264  CALL newton(xx,j,e)
265  IF(xx(1).lt.0.0) THEN
266 
267  print *,' ROCBURN_ZN: rank=',g_zn%rank
268  print *,' rb negative xx(1)=',xx(1)
269  xx(1)=0.0d0
270 
271  END IF
272 
273  IF(abs(e(1)).LT.tol .AND. abs(e(2)).LT.tol) THEN
274 
275 !
276 ! reset final converged values
277 !
278  rb=xx(1)/rhoc ! cm/sec
279  ts=xx(2) ! K
280  xreac=2.0*g_zn%alfac*g_zn%R*ts/rb/g_zn%Ec ! 1/cm
281  fr=dexp(-xreac*g_zn%Ka) ! -
282  tf=(qr/xx(1)+g_zn%Qc+g_zn%Qg)/g_zn%C+to ! K
283 
284 !
285 ! Sets ICs using I&W S.S. T profile
286 !
287 
288  xcond=g_zn%alfac/rb
289  xreac=xcond*g_zn%R*ts*2.0/g_zn%Ec
290  fr=dexp(-xreac*g_zn%Ka)
291  frj=fr*qr/rhoc/rb/g_zn%C/(ts-to)/(1.-g_zn%Ka*xcond)
292  DO i=1,g_zn%nx
293  tn(i)=to+(ts-to)*((1.0-frj)* &
294  dexp(g_zn%x(i)/xcond)+frj*exp(g_zn%x(i)*g_zn%Ka))
295  END DO ! i
296 
297  IF(abs(tn(g_zn%nx)-to)/to .gt. 0.001) THEN
298 
299  WRITE(*,*) ' ROCBURN: rank=',g_zn%rank
300  WRITE(*,*) ' ROCBURN: Error: xmax' &
301  ,' not large enough'
302  WRITE(*,*) ' ROCBURN: To boundary' &
303  ,' condition not satisfied'
304  CALL mpi_abort( mpi_comm_world, -1)
305  stop
306 
307  ELSE
308 
309  tn(g_zn%nx) =to
310 
311  END IF
312 
313  RETURN
314 
315  END IF
316 
317  END DO ! iter
318 
319  print *,' ROCBURN_ZN: rank=',g_zn%rank
320  print *,' ROCBURN_ZN: Error-steady state ZN' &
321  ,' in ignition_combustion'
322  print *,' ROCBURN_ZN: Maximum # of ',g_zn%itermax &
323  , ' iterations exceeded'
324  CALL mpi_abort( mpi_comm_world, -1)
325  stop
326 
327 !
328 ! end of steady state ZN phenomenological combustion model
329 !
330  END IF
331 
332  IF(g_zn%Model_combustion == 3) THEN
333 
334 !
335 ! Model_combustion = 3
336 !
337 ! rb=a*P**n
338 !
339 ! initial guess of rb and Ts
340 !
341 
342  rb=g_zn%a_p*p**g_zn%n_p ! cm/s
343  tn =1000.0d0
344  ts =600.00d0
345 
346  RETURN
347 
348  END IF
349 
350  WRITE(*,*) ' ROCBURN_ZN: rank=',g_zn%rank
351  WRITE(*,*) ' ROCBURN_ZN: Error-ssWSB:'
352  WRITE(*,*) ' No appropriate combustion model'
353  WRITE(*,*) ' ROCBURN_ZN: Model_combustion=', g_zn%Model_combustion
354  CALL mpi_abort( mpi_comm_world, -1)
355  stop
356 
357  RETURN
358 
359  CONTAINS
360 
361 ! ------------------------------------------------------------------------
362 ! INTERNAL PROCEDURES
363 ! ------------------------------------------------------------------------
364 
365 
366 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
367 ! c
368 ! subroutine : Newton c
369 ! c
370 ! c
371 ! Author: c
372 ! c
373 ! Paul Loner c
374 ! c
375 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
376 ! c
377  SUBROUTINE newton(x,J,e)
378 
379  IMPLICIT NONE
380 
381  REAL(DBL) :: x(2),j(2,2),e(2),detj,delx(2),dum
382  integer i
383 
384  detj=j(1,1)*j(2,2)-j(2,1)*j(1,2)
385 
386  if(detj.eq.0) then
387  write(*,*) 'ROCBUNR_ZN Error: singular matrix encountered-could not invert'
388  stop
389  endif
390 !
391 ! invert the 2x2 Jacobian matrix
392 !
393  dum=j(1,1)
394  j(1,1)=j(2,2)
395  j(2,2)=dum
396  j(1,1)=j(1,1)/detj
397  j(1,2)=-j(1,2)/detj
398  j(2,1)=-j(2,1)/detj
399  j(2,2)=j(2,2)/detj
400 !
401 ! multiply J^-1()*e() and update guesses
402 !
403  do 10 i=1,2
404  delx(i)=j(i,1)*e(1)+j(i,2)*e(2)
405  x(i)=x(i)-delx(i)
406  10 continue
407 
408  RETURN
409  END SUBROUTINE newton
410 
411 ! -------------------------------------------------------------------
412 ! END OF INTERNAL PROCEDURES
413 ! -------------------------------------------------------------------
414 
415  END SUBROUTINE zn_sswsb
416 
417 
418 
419 
420 
421 
d
const NT & d
Definition: rational_rotation.h:73
zn_sswsb
subroutine zn_sswsb(G_ZN, P, qr, To, rhoc, rb, Ts, Tf, fr, Tn)
Definition: ZN_ssWSB.f90:56
m_rocburn_zn_global_data
Definition: 02.burn_ZN_global_data.f90:23
i
blockLoc i
Definition: read.cpp:79
x
void int int REAL * x
Definition: read.cpp:74
m_rocburn_zn_global_data::g_burn_1d
Definition: 02.burn_ZN_global_data.f90:48
newton
subroutine newton(x, J, e)
Definition: ZN_calc_burning_rate.f90:432
Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277
j
j indices j
Definition: Indexing.h:6