Coupled simulation of chemical lasers based on intracavity partially coherent light model and 3 D CFD model

Coupled simulation based on intracavity partially coherent light model and 3D CFD model is firstly achieved in this paper. The dynamic equation of partially coherent intracavity field is derived based on partially coherent light theory. A numerical scheme for the coupled simulation as well as a method for computing the intracavity partially coherent field is given. The presented model explains the formation of the sugar scooping phenomenon, and enables studies on the dependence of the spatial mode spectrum on physical parameters of laser cavity and gain medium. Computational results show that as the flow rate of iodine increases, higher order mode components dominate in the partially coherent field. Results obtained by the proposed model are in good agreement with experimental results. ©2011 Optical Society of America OCIS codes: (140.0140) Lasers and laser optics; (140.1550) Chemical lasers. References and links 1. E. A. Duff and K. A. Truesdell, “Chemical oxygen iodine laser (COIL) technology and development,” Proc. SPIE 5414, 52–68 (2004). 2. D. A. Copeland, C. Warner, and A. H. Bauer, “Simple model for optical extraction from a flowing oxygeniodine medium using a Fabry-Perot resonator,” Proc. SPIE 1224, 474–499 (1990). 3. G. D. Hager, C. A. Helms, K. A. Truesdell, D. Plummer, J. Erkkila, and P. Crowell, “A simplified analytic model for gain saturation and power extraction in the flowing chemical oxygen-iodine laser,” IEEE J. Quantum Electron. 32(9), 1525–1536 (1996). 4. T. J. Madden, “Aspects of 3D chemical oxygen-iodine laser simulation,” Proc. SPIE 5120, 363–375 (2003). 5. R. C. Buggeln, S. Shamroth, A. Lampson, and P. G. Crowell, “Three-dimensional (3-D) Navier-Stokes analysis of the mixing and power extraction in a supersonic chemical oxygen iodine laser (COIL) with transverse I2 injection,” presented at 25th AIAA Plasmadynamics and Lasers Conference, Colorado Springs, CO, 20–23 June, 1994. 6. J. Paschkewitz, J. Shang, J. Miller, and T. Madden, “An assessment of COIL physical property and chemical kinetic modeling methodologies,” presented at 31st AIAA Plasmadynamics and Lasers Conference, Denver, CO, 19–22 June, 2000. 7. M. Endo, T. Masuda, and T. Uchiyama, “Development of hybrid simulation for supersonic chemical oxygeniodine laser,” AIAA J. 45(1), 90–97 (2007). 8. A. I. Lampson, D. N. Plummer, J. Erkkila, and P. G. Crowell, “Chemical oxygen iodine laser (COIL) beam quality predictions using 3-d Navier-Stokes (MINT) and wave optics (OCELOT) codes,” presented at 29th AIAA Plasmadynamics and Lasers Conference, Albuquerque, NM, 15–18 June, 1998. 9. B. D. Barmashenko, “Analysis of lasing in chemical oxygen-iodine lasers with unstable resonators using a geometric-optics model,” Appl. Opt. 48(13), 2542–2550 (2009). 10. D. Yu, F. Sang, Y. Jin, and Y. Sun, “Output beam analysis of an unstable resonator with a large Fresnel number for a chemical oxygen iodine laser,” Opt. Eng. 41(10), 2668–2674 (2002). 11. M. Suzuki, H. Matsueda, and W. Masuda, “Numerical simulation of Q-switched supersonic flow chemical oxygen-iodine laser by solving time-dependent paraxial wave equation,” JSME Int. J. Ser. B 49, 1212–1219 (2006). 12. D. A. Copeland and A. H. Bauer, “Optical saturation and extraction from the chemical oxygen-iodine laser medium,” IEEE J. Quantum Electron. 29(9), 2525–2539 (1993). #155372 $15.00 USD Received 26 Sep 2011; revised 11 Nov 2011; accepted 11 Nov 2011; published 9 Dec 2011 (C) 2011 OSA 19 December 2011 / Vol. 19, No. 27 / OPTICS EXPRESS 26295 13. B. Barmashenko, D. Furman, and S. Rosenwaks, “Analysis of lasing in gas-flow lasers with stable resonators,” Appl. Opt. 37(24), 5697–5705 (1998). 14. T. T. Yang, “Modeling of cw HF chemical lasers with rotational nonequilibrium,” J. Phys. C9, 51–57 (1980). 15. K. Waichman, B. D. Barmashenko, and S. Rosenwaks, “Comparing modeling and measurements of the output power in chemical oxygen-iodine lasers: a stringent test of I2 dissociation mechanisms,” J. Chem. Phys. 133(8), 084301 (2010). 16. A. G. Fox and T. Li, “Resonator modes in a maser interferometer,” Bell Syst. Tech. J. 40, 453–488 (1961). 17. A. E. Siegman and E. A. Sziklas, “Mode calculations in unstable resonators with flowing saturable gain. 1:hermite-gaussian expansion,” Appl. Opt. 13(12), 2775–2791 (1974). 18. E. A. Sziklas and A. E. Siegman, “Mode calculations in unstable resonators with flowing saturable gain. 2: Fast Fourier transform method,” Appl. Opt. 14(8), 1874–1889 (1975). 19. M. Endo, M. Kawakami, K. Nanri, S. Takeda, and T. Fujioka, “Two-dimensional simulation of an unstable resonator with a stable core,” Appl. Opt. 38(15), 3298–3307 (1999). 20. A. Bhowmik, “Closed-cavity solutions with partially coherent fields in the space-frequency domain,” Appl. Opt. 22(21), 3338–3346 (1983). 21. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995). 22. E. Wolf, “New theory of partial coherence in the space-frequency domain. Part II: Steady-state fields and higherorder correlations,” J. Opt. Soc. Am. A 3(1), 76–85 (1986). 23. M. Hishida, N. Azami, K. Iwamoto, W. Masuda, H. Fujii, T. Atsutu, and M. Muro, “Flow and optical fields in a supersonic flow chemical oxygen-iodine laser,” presented at 28th Plasmadynamics and Lasers Conference, Atlanta, GA, 23–25 June, 1997. 24. Y. Huai, S. Jia, and Y. Jin, “Analysis and optimization of mixing process with large eddy simulation: An application to SCOIL,” presented at 40th AIAA Plasmadynamics and Lasers Conference, San Antonio, TX, 22– 25 June, 2009. 25. M. Guizar-Sicairos and J. C. Gutiérrez-Vera, “Coupled mode competition in unstable resonators using the exact cavity equations of motion with dynamic gain,” J. Opt. B Quantum Semiclassical Opt. 7(9), 253–263 (2005). 26. G. F. Calvo, A. Picon, and R. Zambrini, “Measuring the complete transverse spatial mode spectrum of a wave field,” Phys. Rev. Lett. 100(17), 173902 (2008). 27. F. Gori, M. Santarsiero, R. Borghi, and G. Guattari, “Intensity-based modal analysis of partially coherent beams with Hermite-Gaussian modes,” Opt. Lett. 23(13), 989–991 (1998). 28. A. E. Siegman and S. W. Townsend, “Output beam propagation and beam quality from a multimode stableCavity laser,” IEEE J. Quantum Electron. 29(4), 1212–1217 (1993). 29. K. Shimizu and S. Yoshida, “High power chemical oxygen-iodine laser of good beam quality,” in Lasers '89; Proceedings of the International Conference (STS Press, 1990), pp. 218—222. 30. J. Bachar and S. Rosenwaks, “An efficient, small scale chemical oxygen-iodine laser,” Appl. Phys. Lett. 41(1), 16–17 (1982). 31. F. Sang, C. Gu, J. Pang, M. Li, F. Li, Y. Sun, Y. Jin, and Q. Zhuang, “Experimental study of a cw chemical oxygen-iodine laser,” High Power Laser Part. Beams 5, 389–393 (1993) (In Chinese).


Introduction
The operational mechanism of gas chemical lasers, e.g. chemical oxygen iodine lasers (COIL) involves the entangled processes of optical oscillation, fluid mixing and chemical reactions [1].To get a better understanding of these processes and improve the laser design, reliable models have to be built.Early studies included only a 1D flowing process, a few chemical reactions that directly related to pumping and lasing, and simplified Fabry-Perot cavity [2,3].Numerous improved models were developed in the following works.For example, 3D computational fluid dynamic (CFD) studies based on the Navier-Stokes equation were carried out for analyzing the mixing process of the molecular iodine and the excited state molecular oxygen [4,5].Chemical reaction packages involving 21 or 35 reactions were also used to give a more comprehensive vision of the energy transfer process taking place inside the reaction chamber (the nozzle and cavity) [6,7].While the near/far field intensity profile and beam quality are of interest, optical computation is needed to be associated with fluid dynamic and chemical kinetic computations [8].
Both geometrical optics and wave optics models have been applied for the optical computation of COIL [9][10][11].For the study of unstable cavities, numerous models faithfully preserve the cavity geometry [8][9][10].In cases of stable resonators, however, the Fabry-Perot cavity and the roof-top cavity are often used as the approximation model [3,[12][13][14].The Fabry-Perot model was successful in estimating the output power and the chemical efficiency with respect to the outcoupling rate and the flow condition.However, certain difficulty arises (C) 2011 OSA when using the Fabry-Perot model to predict the intensity pattern.The upstream-downstream optical coupling effect in COIL with stable cavity, known as the sugar scooping phenomenon, cannot be explained by the Fabry-Perot model [3,[13][14].In order to simulate this phenomenon, Yang adopted the roof-top cavity model which forces the optical field to flip as it is reflected by one of the mirrors [14].Hence, both the Fabry-Perot model and the roof-top model change the geometry of cavity mirrors and obscure the optical mechanism of the device.Another model, namely the constant intracavity intensity model, assumes the intracavity intensity to be uniformly distributed over the resonator cross section [13].Recently this model was coupled with 3D CFD model of the gas flow under various I 2 dissociation mechanism assumptions [15].The predicted output power results were compared with experimental results to examine the I 2 dissociation mechanism assumptions.As shown in [13] the constant intracavity intensity model predicts accurate value of the output power from a stable resonator.However, it neglects the optical propagating mechanism and therefore cannot give reliable information about the output intensity distribution and divergence angle.
It should be noted that many iterative methods based on Fourier optics, including the Fox-Li method [16], the Hermite-Gaussian expansion method [17] and the Fast Fourier transform method [18], work well for unstable resonators which oscillate only at the lowest transverse mode.However, if these methods are directly used on stable resonators with very large Fresnel numbers, the propagating wave expands and converges repetitively making the computational convergence difficult [19].A stable cavity of COIL is often with large Fresnel number (>500) and many transverse modes oscillating simultaneously.Bhowmik noticed the partially coherent nature of the fields, and developed an improved Fox-Li type iteration by starting with a random noise distribution [20].The random noise distribution was supposed to represent the partially coherence field which could fill the cavity width during whole propagation inside the cavity.This method leads to a quasi-steady solution and the final optical intensity is statistically averaged among many rounds.Later Endo developed Bhowmik's method to two-dimensional cases [19].These studies demonstrate the partial coherent nature of the oscillating field and provide a better understanding of the oscillating process.Therefore, a more formalized representation of intracavity field of stable cavities with large Fresnel numbers based on the partially coherent theory, and the coupled computational study of the partially coherent optical field and the chemical reacting flow field are expected.
In this paper, the dynamic equation of the partially coherent intracavity field is established based on the partially coherent light theory [21,22], a numerical method for solving the equation is given, and a numerical scheme for the optical-fluidic-chemical coupled computation is discussed.To our knowledge, a coupled simulation based on the intracavity partially coherent light model and a 3D CFD model is achieved for the first time.The coupled simulation provides critical information about COIL operation and reveals the influence of operational parameters on the formation of the partially coherent beam.

Motion equation of partially coherent light inside loaded stable cavity
According to the partially coherent light theory [22], a stationary optical field of any state of coherence can be represented as an ensemble of monochromatic wave fields.Each member of the ensemble is in the form of ( ) ( ), where ( ) n φ r are eigenfunctions of the Helmholtz equation, r is the Cartesian coordinates of N dimension, n stands for N ordered positive integers ( ) In this paper where loaded stable cavities are studied, three dimensional Cartesian coordinates are established with z axis along the optical axis.And the intracavity partially coherent field consists of the where ( ) are chosen to be the Hermite-Gaussian modes of bare cavity so that they satisfy the orthonormalized relationship [17] ( ) ( ) * , , d .
Intensity of the partially coherent field is expressed as the average intensity of the ensemble members, and by using Eq. ( 3) and Eq. ( 2), equals the uncorrelated superposition of the Hermite-Gaussian modes: , , , .
Here ( ) ( ) ( ) where 2 ∇ ρ is the transverse Laplacian and k is the wave number.Substituting Eq. (3) into Eq.( 6) and using the orthogonality relationship [Eq.( 4)], and noticing ( ) the equations for the expansion coefficients ( ) n a z are obtained as ( ) With some further algebra on Eq. ( 8) we obtain where ( )

Re
denotes the real part of a complex.By having ensemble average on both sides of Eq. ( 10) and noticing that nn nn G G ′ << for all n n′ ≠ , the motion equation of the forward propagating partially coherent light inside the loaded stable cavity is given by We for all n n′ ≠ because in practical chemical lasers, the loaded gain is a slowly varying function with its value near the threshold gain [2,3].And from the orthonormalized relationship [Eq.( 4)] it can be seen that as the gain distribution ( ) , g z ρ tends to a uniform distribution the nn G ′ approaches to zero.Similarly, the motion equation of the backward propagating partially coherent field is given by Equations ( 11), ( 12) and ( 5) constitute the equations for describing the intracavity partially coherent fields of loaded stable cavities and are coupled with a 3D CFD based dynamic gain model in this study.

Numerical scheme of the coupled simulation
The motion equations given in Sec. 2 are working together with a 3D CFD model to perform the coupled simulation of the operational mechanisms of COIL.The schematic of the computation process is shown in Fig. 1.The coupled numerical simulation of intracavity partially coherent light and supersonic flowing chemical reacting media is achieved with an iterative process with the optical computation and the CFD computation executed in turns and providing intermittent data to each other.The CFD computation provides the gain distribution for the optical computation, and the optical computation returns the optical intensity to the CFD codes.This process is repeated until both the gain and the intensity converge to a steady distribution.
Governing equations of the 3D CFD computation include the conservation equations of mass, momentum, energy and species with appended source items representing effects of the chemical reactions and stimulated emission.These equations can be found in [23] and are solved by the finite volume method.Second order central difference schemes are used for spatial discretization of the viscid terms.Roe's flux difference splitting method, which is second order accurate in space, is applied to evaluate the convective flux.In this paper we specifically lay stress on the optical computation methods since the CFD schemes adopted are similar to those in previous works [23].

Computation method of intracavity partially coherent fields
Assuming the cavity length is L .The outcoupling mirror M 1 and the reflecting mirror M 2 are located at plane 1 z z = and plane 2 z z = , respectively.The radii of curvature of mirror M 1 and M 2 are 1 R and 2 R , respectively.For solving Eqs. ( 11), (12), and ( 5), the Hermite- Gaussian modes are firstly computed based on the following formulas [17] ( ) , , , , , n y z φ has equations similar to Eqs. ( 14)-( 16) but just replacing 1 n with 2 n and x with y .The complex radius of curvature ( ) q z and the ( ) ( ) ( ) where ( ) R z and ( ) z ω are the wave front radius and spot radius at plane z , respectively.
( ) Generally, when the mode ( ) given by Eqs. ( 13)-( 19) is reflected by a spherical mirror, e.g. by M 2 , ( ) q z and ( ) Q z of the reflected mode will change into ( ) and Thus the reflected modes need to be computed again.However, by setting the reference plane 0 z z = and the ( ) 0 q z according to Eqs. ( 23)-( 25), we have ( ) ( ) . Under this condition, the modes satisfy the self-consistency criteria of the cavity.It means that the forward and backward propagating modes have exactly the same profile and can both be expressed by Eqs. ( 13)- (19).It provides convenience for numerical realization since these modes keep unchanged after reflection and therefore need to be computed only once: ( ) ( ) given back to the CFD codes to update the gain data ( ) and therefore forms a close-loop coupled simulation.

Calculation and discussion
The code was applied to the kilowatts COIL experimental platform in Key Laboratory of Chemical Lasers, Chinese Academy of Sciences.The device is similar to the one reported in our previous work [24].However, some optical and aerodynamical parameters are updated.In the current study, the device is equipped with a stable resonator that consists of a plane outcoupling mirror M 1 and a spherical total reflector M 2 .The cavity length is L = 0.6 m.The transmittance of M 1 is 3% α = . The radius of curvature of M 2 is 10.0 m.The aperture is a rectangle with side lengths of y a = 30 mm and x a = 60 mm.The COIL uses nitrogen as the carrier gas, uses subsonic mixing of the primary and secondary flows.The total primary and secondary flowrates are 130 mmol/s and 35 mmol/s, respectively.The molar flowrate ratio of chlorine to nitrogen in the primary flow is 1:1, and of iodine to nitrogen in the secondary flow is 1:21.The flowing direction is along the x axis.Flow conditions for CFD calculation are listed in Table 1.And the chemical reactions used in simulation can be found in [7].From Eqs. ( 15)-( 17) we know that the spot size of the fundamental Hermite-Gaussian mode on the reference plane is 0 ω = 0.0997 mm.And the reference plane is located at the output plane 1 0 z z = = . A practical question is how many modes are needed in computation.
We determine the minimum amount of modes by ensuring the expansion basis cover the whole aperture region.For the x n th order 1D Hermite-Gaussian mode, its plot exhibits a quasi-periodic behavior with x n zero crossings, and the function amplitude falls to zero very rapidly outside the point x x n ω = [17].Therefore, the number of modes that must be involved in computation is estimated by max 0 .
And for the same reason we have max 0 .
For the computation we adopt 1000 x y n n = = .A good check on the accuracy of the computed modes is how well the orthonormalized relationship Eq. ( 4) is satisfied numerically [25].The computed modes satisfy Eq. ( 4) within an error of 10 −4 for all computed modes.The coupled simulation experiences 35 iterations before convergence.The output intensity is shown in Fig. 2. The intensity distribution has peaks in both the upstream and downstream areas.This is in agreement with our experiments and some previous publications [13].This effect, known as the sugar scooping phenomenon [13], can be explained by the presented model in the following way.According to the model described in Sec.consists of a large amount of Hermite-Gaussian modes, and the intensity is determined by the uncorrelated summation of each mode.Noticing that the Hermite-Gaussian modes are symmetric, so the peak at downstream is formed by the superposition of Hermite-Gaussian modes which also create the peak at upstream.Therefore, the partially coherent nature of the intracavity field is responsible for the formation of the sugar scooping phenomenon.The computed gain is shown in Fig. 3.And the small signal gain is shown in Fig. 4 for comparison.Simulation result indicates that as the gaseous media enters the cavity aperture, the saturated gain drops dramatically due to the power extraction.However, within the aperture the saturated gain decreases very slowly along x direction (flow direction) and keeps almost unchanged along y direction.This proves the validity of the assumption The slowly varying character of the saturated gain conforms to previous studies on COIL [2,3].
The expansion coefficient ( ) λ represents the intensity of the ( ) the partially coherent field and is referred to as the spatial mode spectrum [26].Previous studies have focused on the reverse problem of determining the spatial mode spectrum from the intensity distribution of certain partially coherent fields [27].The method introduced in this paper, however, enables a dynamic study for predicting the spatial mode spectrum from physical parameters of the laser cavity and the gain medium.spatial mode spectrum of the output field.Furthermore, the dependence of the spatial mode spectrum on the flow condition can be studied.In some previous works [19,20], the partially coherent field is computed by averaging on the round-trip coherent fields starting from a noise distribution.This method leads to a quasisteady solution and demands many rounds of computation for the averaging.The proposed model, however, includes the averaging inherently in the motion equation.The expansion coefficients (spatial mode spectrum)

Conclusion
A coupled simulation of the intracavity partially coherent light model and the 3D CFD model of the chemical reacting flowing lasing media is realized for the first time.The dynamic equation of partially coherent intracavity field is deduced based on the partially coherent light theory.The numerical scheme for the coupled simulation as well as the method for computing the intracavity partially coherent field is given.The coupled computation gives results on the output intensity, the saturated gain, the filled factor, and the spatial mode spectrum.Computational results indicate that the formation of the sugar scooping phenomenon is closely related to the partially coherent nature of the intracavity field.The presented model enables studies on the dependence of the spatial mode spectrum on physical parameters of the laser cavity and the gain medium.Computational results also show that as the flow rate of iodine increases, higher order mode components dominate in the partially coherent field.Experimental results proved the validity of the proposed model.The output power, divergence angle and filled factor predicted by the proposed model are in good agreement with the experimental results.

#
155372 -$15.00USD Received 26 Sep 2011; revised 11 Nov 2011; accepted 11 Nov 2011; published 9 Dec 2011 deduct the motion equation to which the expansion coefficients stable cavities.In the following part of this paper, discrimination among them is necessary.Let us consider the forward propagating field ( )

Fig. 5 .
Fig. 5. Normalized spatial mode spectrum of the output field (mass fraction of molecular Iodine is 30%).

Figure 6
gives the normalized spatial mode spectrum under various iodine flow rates.The results indicate that higher order mode dominates as the flow rate of iodine increases.Results of the proposed model are compared with results of the Fabry-Perot approximation model, the roof-top model, and the experiment.The calculated output power of the proposed model, the Fabry-Perot model and the roof-top model are 1516 W, 1490 W and 1551 W, respectively.The experimentally measured output power is 1530 W. The calculated output power of the proposed model is in good agreement with the experimental result.And all three models predict close values of the output power.However, the divergence angle predicted by the proposed model is quite different from results of the Fabry-Perot model and the roof-top model.Divergence angles (in the y direction) predicted by the Fabry-Perot model and the roof-top model are 9.0 × 10 −5 rad and 8.7 × 10 −5 rad, respectively.Divergence angle (in the y direction) predicted by the proposed model is 9 mrad.Such a huge difference among the models is caused by the following reason.In the Fabry-Perot model and the roof-top model, curvatures of the cavity mirrors are neglected.Divergence of the beam is mainly caused by the diffraction of the aperture and the ununiformity of the active media.The diffraction limited divergence angle is characterized by #155372 -$15.00USD Received 26 Sep 2011; revised 11 Nov 2011; accepted 11 Nov 2011; published 9 Dec 2011 (C) 2011 OSA 19 December 2011 / Vol. 19, No. 27 / OPTICS EXPRESS 26304

Fig. 7 .
Fig. 7. Calculated intensity distribution by (a) the Fabry-Perot model and (b) the roof-top model.

Fig. 8 .
Fig. 8.Comparison of intensity distribution along x axis.(a) The experimental result.(b) The proposed model.(c) The Fabry-Perot model.(d) The roof-top model.
quantities by definition.So in the numerical computation no averaging among many rounds is needed.And the proposed model leads to steadier and smoother results compared with those previous works.#155372 -$15.00USD Received 26 Sep 2011; revised 11 Nov 2011; accepted 11 Nov 2011; published 9 Dec 2011 (C) 2011 OSA 19 December 2011 / Vol. 19, No. 27 / OPTICS EXPRESS 26306