Method of lines solution to the transient SBS equations for nanosecond Stokes pulses

The spectral and temporal evolution of distributed sensing based on stimulated Brillouin scattering (SBS) in optical ﬁbers for several-nanosecond Stokes pulses is demonstrated by using the method of lines (MOL) solution of the transient SBS equations. A superbee ﬂux limiter is utilized to avoid numerical damping and dispersion that would otherwise be brought on by the approximation of spatial derivatives associated with the partial differential equations (PDEs). In order to increase computational efﬁciency, an approach is adopted whereby the sparse PDE Jacobian matrix integrator option of the ODE solver(s) is employed. Simulation examples of SBS-based sensing for ﬁbers containing sections with different Brillouin frequencies are presented. To the best of our knowledge, this MOL solution is proposed for the ﬁrst time for modeling of the transient SBS equations for nanosecond Stokes pulses with different waveforms in a SBS based ﬁber optic sensor. [DOI: http://dx.doi.org/10.2971/jeos.2013.13049]


INTRODUCTION
Stimulated Brillouin scattering (SBS) in optical fibers permits us to measure temperature and/or strain on a truly distributed basis, over kilometric ranges with high resolution [1]- [3].It consists in the coupling between two counterpropagating optical waves, the pump and the Stokes wave, and an acoustic wave.In this three-wave mixing process, power is transferred from the pump light wave to the Stokes light wave (that is, to the light wave having a lower frequency) and also to the acoustic wave.The interaction strength depends on the frequency offset between the two light waves and attains its maximum at the so-called Brillouin frequency ν B .As the Brillouin frequency shift changes linearly with temperature and strain, a distributed temperature-strain sensor can be realized using SBS.
The theoretical model used to simulate the sensing in singlemode fiber is based on the wave equations for the pump and Stokes waves with the slowly-varying field amplitudes A p (z, t) and A s (z, t), respectively, which interact nonlinearly by the excitation of the acoustic wave ρ(z, t) [4]: where A p (z, t) and A s (z, t) [W] stand for the amplitudes of pump E p (z, t) and stokes E s (z, t) waves respectively, and n is the refractive index of the medium, c is the speed of light [m/sec], γ e is the is the electrostrictive constant, ρ 0 is the unperturbed density [g/m 3 ] and Γ B is the Brillouin linewidth [MHz] [4].A * represents the complex conjugate of the respective field.In the case of a transient SBS regime for nanosecond optical pulses participating in SBS, the slowly-varying amplitude approximation normally used may no longer be valid for the acoustic field.Therefore, the evolution of acoustic wave equation for the field ρ(z, t) with second-order time derivative follows [5]: Simulation has proven to be a valuable tool in understanding SBS based distributed sensor systems.Accurate simulation models allow for the rapid design and optimization of both sensor infrastructure and signal processing techniques.Visualization of the spatiotemporal behavior of all three fields involved in SBS, at all points through the fiber, and is especially desirable in investigating phenomena whose origins may not be clear from the 'observable' data in an experiment [6]- [8].
If long pulses are assumed, the time dependence in Eq. ( 1) can be neglected, resulting in a simple and accurate solution to the problem.However, in the short pulse regime, the temporal derivatives cannot be neglected and the solution becomes more complex [9].For this solution, high speed and accuracy with the low complexity is highly desirable.In general, for the solution of Eq. ( 1), conventional first order finite difference time domain (FDTD) solution schemes are preferred since they offer simplicity with respect to the form of the solution, although they suffer from several drawbacks, e.g. they require more spatial nodes to achieve the same accuracy as higher order schemes, and this often results in longer computer run-times.
To solve Eq. ( 1), a numerical method based on the Simpson's rule to approximate temporal integrals was introduced by [10].However, they used an implicit method and employed linearization by assuming that E m+1 sn ≈ |E m sn | and this re- placement weakens the coupling of the two laser fields and hence the results deviate from the exact solutions.In Ref. [6] the authors proposed a modified solution which is valid only when the time step size is equal to spatial step size.Moreover, their solution differs considerably from the real solution at locations where abrupt change takes place.Recently, we introduced the MOL solution of the SBS equations with appropriate flux limiter and sparse matrix integrator [11].In this paper, we extent the results in Ref. [11] while computing a system of three second-order hyperbolic PDEs (Eqs.( 1)).
In this paper we introduce, for the first time to the best of our knowledge, the solution of transient SBS equations for nanosecond Stokes pulses using Matlab ODE solvers [12].Bu using this solution, the spectral and temporal evolution of distributed sensing based on stimulated Brillouin scattering (SBS) in optical fibers for several-nanosecond Stokes pulses is demonstrated.We introduce the solution technique and have adapted the software to provide a fast and stable solution by employing sparse matrix techniques and by using an appropriate flux limiter.Briefly, our solution technique combined with the Matlab software abilities provides a powerful method of simulating this SBS problem.
In arriving at a solution, the problem domain is first discretized over a finite grid and organized in a way that provides solutions for the optical and acoustic fields.The resolution of the solution can thus be varied to suit computational requirements.We have also employed an appropriate flux limiter to combat numerical damping/dispersion which avoids spurious oscillations in the solution.We have used software which identifies the sparsity pattern of the solution which is then supplied as an input to the solver.Additionally, we have determined the appropriate Matlab ode solver for this particular application.By use of an appropriate sparse matrix ode solver along with an effective flux limiter, for nanosecond pulse widths, it is shown that the computation can be successfully performed in a reasonable computation time without the damping factor exceeding 3%.Examples are presented, showing the utility of this efficient simulation technique.

The Dynamics of Stimulated Brillouin Scattering (SBS)
In this paper we consider a system whereby continuous wave (CW) light, at frequency, ν p , is injected into the fiber at position z = L. Pulses of light, known as Stokes pulses, are injected into the fiber at z = 0, at frequency, ν s .Due to the SBS, there will be coupling of a counter propagating CW pump wave and a Stokes pulse wave via an induced acoustic wave.An en- hanced interaction between the two beams occurs when the frequency difference of the lasers matches the frequency of the longitudinal acoustic phonons of the optical fiber.There is then a transfer of energy from the high frequency beam (pump) to the low frequency beam (Stokes).The amount of loss of the pump is recorded at z = 0 as a function of the frequency difference in the form of a Brillouin spectrum and thus the power of the CW field intensity over time at z = 0 corresponds to the available time domain information in a real sensing experiment.The maximum loss occurs when the frequency difference of two beams matches the Brillouin frequency of the fiber.This maximum loss depends on the resonance detuning parameter, given by where δν Res (z) is the Brillouin frequency along the fiber and depends on the local strain and temperature.The resonance detuning parameter, δ(z), is the value of the pump-Stokes frequency shift for a resonant interaction at a given point.The resonance detuning parameter, δ(z) is dependent on the temperature and strain along the fiber.This dependence can form the basis of distributed fiber optic strain and temperature sensors.
A distributed Brillouin sensor, based on a Brillouin Optical Time Domain Analysis (BOTDA), is a device where Stokes and pump beams counter-propagate and where frequency and time dependent variations of the Stokes (or pump) intensity is detected.These sensors have been shown to provide the best performances in terms of sensing length, spatial resolution, temperature, and, strain accuracies [2,13].The interaction is a maximum when the optical frequency difference between pump (ν p ) and Stokes (ν s ) corresponds to the Brillouin frequency (ν B ).At the Brillouin frequency, the Stokes beam is amplified at the expense of the pump.ν B is the longitudinal acoustic phonon frequency and is a local material signature.An increase (decrease) in temperature or strain induces a proportional increase (decrease) of ν B .The sensor detects the Brillouin frequency variation as a function of position.The distributed nature of the sensor is achieved by field modulation of the probe beam.The topology of this sensor is shown in Figure 1 [14].

Theoretical Model
By making suitable transformation, Eqs.(1) can be written as in Eqs.
(3) which is a system of three second-order hyperbolic The pump pulse T determines the spatial resolution [15]. PDEs, where E p , E s , and E a are the pump, Stokes, and acoustic fields, respectively, and j represents imaginary number √ −1.Γ = 1/τ ph is the relaxation rate ω = ω p − ω s , ω p , ω s are the acoustic field, pump, and Stokes frequencies, respectively, Ω B is the Brillouin frequency, g = νΓΩg B and g B is the SBS gain factor.Here, τ ph is the phonon lifetime and it is 10ns for silica fibers.

Boundary and Initial Conditions
The CW light from a laser source representing the pump field, Ep(z = L) is injected into the fiber at z = L. Therefore, the boundary condition of the pump wave is: where E p (z = L) is the power of the CW laser.At the opposite end, the Stokes field, E s , is introduced via a pulsed source.We assume a three-section pump pulse of the form (Fig. 1a): where α,β,γ are amplitudes of the corresponding time intervals.Here u(•) represent the Heaviside unit step function.The first section of the pulse of height αA 0 s and width t 0 enters the fiber (z = 0) at t = 0, followed by a βA 0 s -high, T-wide middle section, ending with an infinitely long section of height γA 0 s .In Matlab u(•) can be implemented via the function heaviside function.For the Bright pulse with a CW component waveform the base level of the Stokes pulse is characterized by the extinction ratio ER = 10log(β/α).
After propagating through the fiber, it is assumed that the pump and Stokes pulses exit the fiber without reflection.The acoustic field is assumed to be zero everywhere at t = 0, and the CW light and Stokes leakage intensities turn on everywhere at t = 0. Therefore the initial conditions are [6], E s (z, t = 0) is β for intensity pulse and α for the bright pulse with a CW component, dark pulse waveform and π-phase pulse, respectively.

AC Detection of the pump output in the case of a finite CW component (base) of the Stokes pulse
Distributed sensing based on stimulated Brillouin scattering (SBS) is obtained by temporally resolved Brillouin spectra, i.e. the depleted output pump as a function of the detuning from the Brillouin frequency.The Brillouin frequency gives the strain and/or temperature along the fiber.The timedependent output pump power or the time-dependent deviation from the stationary output power, induced by the Stokes pulse, is measured in the experiment.There are two possibilities of pump loss definition as DC or AC detection, respectively.However, AC detection is most likely to be used in the experiment to measure the small relative variations of the output pump signal in the case of several-nanosecond Stokes pulses.The pump loss defined as AC detection, α AC = E p (z = 0, t = 0) − E p (z = 0, t), each as a function of the time and frequency detuning ∆Ω = Ω − Ω B , represent the Brillouin spectrum.Here E p (z = 0, t) is the time-dependent output pump power at z = 0, E p (z = 0, t = 0) is the output pump power at the initial time moment, defined by the stationary pump-base SBS before the Stokes pulse arrives in the fiber [7].

Finite Difference Time Domain (FDTD) Solution
To obtain a solution to Eqs. (3a)-(3c), we follow the authors in Ref. [9] and substitute the following: into Eqs.(3a)-(3c) and with making transformation z = z 1/ν • γ a to remove ν yields the equations, Where, Φ = Φ a + Φ s − Φ p .The steady state is obtained by setting the time derivatives to zero.Combining Eqs.(8c) and (8f) yields tan Taking the positive part of the amplitudes, from Eq. (8c) it can be concluded that, The basic idea of the method of lines (MOL) is to replace the spatial derivatives with algebraic approximations [16].Since, spatial variables are discretized and time is used as the continuous variable, this approach is called a semidiscretization.This effectively removes the spatial derivatives from the PDE and, since only the initial value independent variable remains, e.g., t, the PDE has been converted to a system of approximating ordinary differential equations (ODEs) that can integrated by standard, well-established numerical algorithms for initial value ODEs.
Eqs. (3a)-(3c) are formulated by making some transformations which are t real • γ a = t trans f ormed , z real • γ a /v = z trans f ormed , where v is the speed of light in the fiber and γ a is 110 • 10 6 .Thus, the problem is discretized on a domain of, The solution is obtained on uniform grids in z and t, z trans f ormed i , t trans f ormed j := z i , t j : For the MOL solution of Eq.( 1b), E p (z, t) can be approximated by using the upwind finite difference (FD) approximations for the first derivatives in z so that E p ≈ (E p i − E p i −1 )/∆z, where i is an index designating a position along a grid in z and ∆z is the spacing in z along the grid (i − 1 designates the upwind direction).For the MOL solution of Eq.(1b), E s (z, t) can also be approximated by upwind FD approximations, i.e., E s ≈ (E s i+1 − E s i )/∆z (i + 1 designates the upwind direction).
The solver integrates the PDEs via the Matlab statement, [t,u]=ode23t(@pde,tout,u0,options); The initial conditions u0, transformed time interval tout, and options, which include tolerances and Jacobian sparsity information, are supplied to the solver.

Application of a Sparse Matrix Integrator
For a large system of ODEs, it is typical that only a few components of y appear in each equation.If component j of y does not appear in component i of f (t, y), then the partial derivative ∂ f i /∂y j is zero.If most of the entries of a matrix are zero, the matrix is said to be sparse.By storing only the nonzero entries of a sparse Jacobian, storage is reduced from the square of the number of equations d to a modest multiple of d.If the Jacobian is sparse then so is the iteration matrix.As with storage, the cost of solving linear systems by elimination can be reduced dramatically by paying attention to zero entries in the matrix.By taking into account the known value of zero for most of the entries in the Jacobian, it is typically possible to approximate all the nonzero entries of several columns at a time.An important special case of a sparse Jacobian is one that has all its nonzero entries located in a band of diagonals [17,18].
If the system of PDEs could be represented purely by banded Jacobians, then a banded integrator would be even more efficient than the sparse integrator since the banded integrator would know in advance where the nonzero elements occur (they all are in the band).It would therefore not have to search for the nonzero elements, and follow the resulting logic to use these nonzero elements, all of which add complexity to a sparse integrator.However, these nonzero elements are located/displayed with the Matlab command spy( ).For the PDE systems of Eqs.(10a)-(10f), some of the nonzero elements are located along the main diagonal as well as the others located outside the band.Therefore, for this case, a sparse integrator is particularly effective since it also operates on out-ofband elements, the so-called outliers, in performing the numerical integration.Thus, this combination of elements in a band along the main diagonal, plus outliers, leads to the efficiency of sparse matrix integration [11].
Since this particular application has a main diagonal, plus outliers, we have used the routine jpattern˙num() to produce a Jacobian map for the sparse matrix option of the ode solver.In this routine, the elements of the Jacobian matrix are computed by finite differences (FDs).This requires a base point, or base points, around which the numerical derivatives are computed.A base point can be any value within the range of the variation of the associated dependent variable.A precise value is not required (but if the ode solver fails, some experimentation with this value may be required).In our application we set the base points ybase(i), equal to the initial condition for each dependent variable which is sufficiently accurate to provide a good solution.
In order to calculate the partial derivatives in the Jacobian matrix by FD approximations, the derivatives at the base point, dyi/dt, are also required.Note that the routine pde() Once the Jacobian matrix has been defined by this call to numjac(), a map of the Jacobian is plotted by the code that follows the call to numjac().With a call to the routine spones(), each non-zero element is replaced by a "1" so as to create a "0-1" map of the Jacobian matrix [11,16]. S=spones(sparse(Jac)); In Figure 2, the sparsity pattern of Eqs.(10a)-( 10f) is plotted by a subsequent call to the routine spy(S).Code for jpattern˙num() is available as a library routine which can be downloaded from Ref. [19].
Note that only 0.185 % of the 6x300 × 6x300 elements of the Jacobian map of Figure 2 are nonzero.This is not uncommon; that is, most of the elements are zero.This condition is the reason for the use of sparse matrix integrators since they will in general process only the nonzero elements and will not expend computer time processing the many zero elements.The processing time with and without jpattern num() is illustrated in Table 1  ing in computer time can be very substantial.It is clear that the additional logic of the sparse matrix integrator (to detect the nonzero elements of the Jacobian matrix and then perform the numerical integration using only these nonzero elements) is well worth the additional complexity of the coding for the sparse matrix.Moreover, for ode23s, using Jpattern num() is the difference between success and failure.

Use of a Flux Limiter
Numerical dissipation is generally observed when solving strongly convective or strongly hyperbolic PDEs.The main cause of the dissipation is low order approximation of the discretization.It should be noted that numerical dissipation produced by first order approximations is to be expected due to truncation of the Taylor series that is the basis for the approximations.In recent years, various high resolution schemes have been developed to obviate this effect with a high degree of accuracy, albeit at the expense of algorithmic and computational complexity.Examples of particularly effective schemes are based upon flux/slope limiters [20] and WENO methods [21].
Flux limiters are used in numerical schemes to solve problems in science and engineering that are described by partial differential equations (PDEs).Their main purpose is to avoid the spurious oscillations (wiggles) that would otherwise occur with high-order spatial discretization due to shocks, discontinuities, or steep gradients in the solution domain.They can be used directly on finite difference schemes for simple applications.Use of flux limiters, together with an appropriate high-resolution scheme, makes the solutions total variation diminishing.Flux limiters are also referred to as slope limiters because they both have the same mathematical form and both have the effect of limiting the solution gradient near shocks or discontinuities.In general, the term flux limiter is used when the limiter acts on system fluxes, and slope limiter is used when the limiter acts on system states.The main idea behind the construction of flux limiter schemes is to limit the spatial derivatives to realistic values; this usually means physically realizable values.They are used in high-resolution schemes for solving problems described by PDE's and only come into operation when sharp wave fronts are present [22].
Considering the semidiscrete scheme below, where, for a finite difference scheme, F(u i+1/2 ) and F(u i−1/2 ) represent flux values on the grid at point x = x i+1/2 and point x = x i−1/2 .If these fluxes can be represented by low-and high resolution schemes, then a flux limiter can switch between these schemes depending upon the gradients close to the particular cell as follows: 2 where, f low low-resolution flux, f high high-resolution flux,φ(r) = flux limiter function, and r represents the ratio of successive gradients on the solution mesh, i.e., The limiter function is constrained to be greater than or equal to zero, i.e., r ≥ 0. Therefore, when the limiter is equal to zero (sharp gradient, opposite slopes, or zero gradient), the flux is represented by a low-resolution scheme.Similarly, when the limiter is equal to 1 (smooth solution), it is represented by a high-resolution scheme.The various limiters listed below have differing switching characteristics and are selected to suit the particular problem and numerical solution scheme.
No particular limiter has been found to work well for all problems, and a particular choice is usually made on a trial-anderror basis [11,22].
van Leer: smart: φ sm (r) = max [0, min(2r, (0.25 + 0.75r), 4)] ; lim superbee: φ sb (r) = max [0, min(2r, 1, min(r, 2))] ; The parameters of the flux limiter, which available as a library routine which can be downloaded from It should be noted that when applying flux limiter solutions it is necessary to have sufficient grid points (say 10 or more) between the occurrence of adjacent steep gradients.If insufficient grid points are used, then smearing of the solution may occur or, in extreme cases, the solution could blow up.
In order to establish the best flux limiter function to suit our particular application we performed numerical experiments to compute the degree of smearing.In Figure 3, The effect of van Leer, Smart and Suberbee limiters in terms of smearing can be seen.
In Figure 3(a), it is clearly seen that the Stokes magnitude is damped due to the numerical dissipation.In Figure 3(b)-(e), the effect of different flux limiters on the Stokes pulse is demonstrated.Damping which can be considered as the amplitude ratio at z = 0 and at z = L is calculated as, (b) 62.2 % for van Leer, (c) 55.44 % for smart (d) 47.74% for superbee flux limiter, respectively (n=100).In Figure 3(e), the damping effect is decreased to only 2.8 % with superbee flux limiter for 300 discretization points.In Figure 3(b)-(e), it is clearly seen that, for this particular application, the superbee limiter is more capable of avoiding numerical damping compared to the performance of van Leer and Smart flux limiters.Since it has superior performance for this application, the superbee flux limiter is preferred.

EXAMPLE
The simulation models an optical fiber, 15 m in length, and the parameter values of the simulation are inferred from the pa- The expected result of this simulation is that the Stokes pulse will propagate down the section of fiber with no loss or dispersion.Since the fiber's natural δν Res (z) is much different from ν p − ν s , there should be little interaction in the 'unstrained' portions of the fiber.Upon entering the 'strained' section, the pulse will interact with the CW field, causing a rise in the acoustic field intensity.The CW field will be depleted, and the pulse will be enhanced.The range of the depletion depends on the extinction ratio and the interaction is stronger for a low extinction ratio.After the Stokes pulse passes through the strained section, the acoustic field will decay away and the CW field will return to its steady state value.The effect of the interaction will take the form of a depleted region of CW light propagating to the end of the fiber [6].
Figure 4 shows the Stokes, CW, and acoustic fields, simulated with and without the superbee flux limiter.As can be seen in Figure 4(a), dissipation and damping is clearly observed.Over time, the Stokes pulse broadens, as well as decreases in intensity.However, the expectation from this simulation is that Stokes pulse will propagate down the section of fiber with no loss or dispersion.It must be noted that, in the 'unstrained' portions of the fiber ∆ = 0 i.e., ν p − ν s = δν Res (z), in the "strained" section of the fiber ∆ = 0, i.e., ν p − ν s = δν Res (z).Plots (c) and (d) show the fields simulated with the superbee flux limiter.Especially in plot (c), with flux limiter, it can be clearly seen that damping is suppressed and the ratio of the amplitudes at z = 0 and at z = L becomes 0.972 which means that Stokes wave lost only 2.8% of its magnitude.Another expectation of the simulation is that, when the pulse enters the strained section there should be a rise in the acoustic field intensity.However, the superbee flux limiter prevents this rise during stabilization.
For the case of a 4-ns Stokes pulse with ER = 15 dB and 50 dB, in Figure 5, AC detection is plotted as a function of the time and frequency detuning ∆Ω = Ω − Ω B .Here t 0 and fiber length is assumed to be 5 ns and 15m, respectively.As can be seen in Figure 5, AC detection shows different characteristics at the time t < t0, before the probe pulse arrives in the fiber, and at t > 2L/ν + t 0 + τ s (159 ns), after it has left the fiber and the related pump depletion has propagated from z = 0 to the output z = L.The AC-detected Brillouin spectrum shows a Lorentzian line shape with the Brillouin natural width ∆ν B = Γ/π ≈ 32MHz resulting from the pumpbase SBS and this corresponds to no signal before pulse arrival (t < t0) and transient relaxation again to zero signal through complicated Brillouin loss signal oscillations with even negative values when out-of-resonance (t > 2L/ν + t 0 + τ s ).
Now we address the SBS-based sensing for longitudinally inhomogeneous distribution of the Brillouin frequency along the fiber.First we consider AC-detected Brillouin spectra of the fiber consisting of two 10-m long sections with Brillouin frequencies 12.800 and 12.840 GHz for the case of a 4-ns Stokes pulse (Figure 6).We verify how the sections and especially the boundary between them are identified AC detection for two different base levels of 15 and 50 dB.The thick dashed curve is drawn at the time t = 157 ns, when the pulse peak center is passing through the middle of the fiber, and this point should be identified by sensing as a boundary between the sections.Additionally, at the moment in time corresponding to the boundary, the spectrum transition from one line to another is clearly visible.However, due to broad pulse spectrum this transition is complicated by the prolonged relaxation of the acoustic field at the frequency related to the previous section, as the pulse has already progressed to the next section.The other parameters are taken from Figure 5. Thick dashed curves are drawn for the output pump at times when the pulse peak center passes through the section's center, and at the frequency shift corresponding to the Brillouin frequency of the section.They denote the time and the frequency, which should be revealed from the detection as a section's spatial and spectral "location".
However, it was observed that, for pulses substantially shorter than the acoustic lifetime in silica, the measured gain profiles abandon their expected spectral broadening and the linewidth was shown to gradually return to the natural value, as determined by the acoustic lifetime [15].This was explained by the fact that in a real experiment the Stokes pulses have a small cw component (base) due to unavoidable leakage of any optical amplitude modulator used to generate nanosecond pulses.
In this case the linewidth could still be close to the natural linewidth, provided the relative contribution of the transient regime on the output pump depletion is weak [7].In our numerical experiment the same trend is observed in Figure 8.

CONCLUSIONS
In this paper we introduce, for the first time to the best of our knowledge, the solution of transient SBS equations for nanosecond Stokes pulses with different waveforms using Matlab ODE solvers (MATLAB, 2012).Bu using this solution the spectral and temporal evolution of the AC detection is simulated.The MOL scheme provided the basis for a solution which was adapted to a suitable mesh size allowing rapid calculation of the fields involved in the SBS interaction.We introduce the solution technique and have adapted the software to provide a fast and stable solution by employing sparse matrix techniques and by using an appropriate flux limiter.The simulation method presented here can be used to foster a rapid understanding of the internal dynamics of the distributed sensing based on SBS, lending insight for both experimental design and interpretation of results.
FIG. 8 FWHM versus Stokes pulse duration for AC detection at t =157 ns, different extinction ratio and the same pump and Stokes pulse parameters as in Figure 5.
Dashed line shows the natural Brillouin linewidth.

FIG. 1
FIG. 1 Diagram of a typical BOTDA distributed sensor system.The energy is transferred from the CW pump to the low frequency Stokes pulse.The amount of loss of the pump is recorded at z = 0 as a function of the frequency difference in the form of the Brillouin spectrum, Thus the power of the CW field intensity over time at z = 0 corresponds to the available time domain information in a real sensing experiment.

FIG. 2
FIG. 2 Sparsity Pattern of the Eqs.(10a)-(10f).The solution is obtained with n = 300 [23], are as follows, y1x=-flux˙function(xl,xu,n,y1,1); y2x= flux˙function(xl,xu,n,y2,-1); The fifth parameter is direction of the flow for the linear advection equation u t + cu x = 0.This argument is required, since the limiter requires the direction of flow or wave propagation.It takes the values +1 or −1 depending on the direction of the upwinding of the FD.The fifth parameter value of −1 for y2x indicates that the pulse moves right to left in z at velocity c = −1 ( i.e., counter propagating in z).The spatial domain in z is defined as xl ≤ z ≤ xu with n points.The interval xu − xl is determined by the fiber length.

FIG. 3
FIG.3The simulation of the pulse, (a) without using flux limiter, with using (b) van Leer, (c) smart (d) superbee flux limiter with n=100 discretization points.In (e) superbee flux limiter is used with n=300 discretization points.

FIG. 4
FIG. 4 Using the simulation of Eqs.(10a)-(10f) with n = 190 discretization points and t 0 = 0 ns, evolution of the Stokes (y1), and Acoustic (y3) fields for a 15 m section of fiber is plotted in Figure 4(a)-(c), (b)-(d), respectively.The frequency difference of Stokes pulse wave was set to be on resonance with a 2 m section in the center of the fiber.Since, in plots (a), (b), flux limiter is not used, the dispersion and attenuation is evident.Plots (c), (d), and show the results when the flux limiter is used.Compared with the low order discretization, with the introduction of the flux limiter, dispersive attenuation (smoothing of sharp edges) along the fiber obviously is eliminated.

Figure 7
Figure7illustrates the SBS-based sensing of a 20-cm section in the middle of a 10-m fiber, where the Brillouin frequency of the section is shifted by 80 MHz relative to that of the rest of the fiber.In Figure7, we have used a 2-ns Stokes pulse with ER = 15 dB (Figure7(a)) and 50 dB (Figure7(b)), respectively.The other parameters are taken from Figure5.Thick dashed curves are drawn for the output pump at times when the pulse peak center passes through the section's center, and at the frequency shift corresponding to the Brillouin frequency of the section.They denote the time and the frequency, which should be revealed from the detection as a section's spatial and spectral "location".

FIG. 5
FIG. 5 Brillouin spectra resolved in the time domain by AC detection for the pump power P p = 5 mW and the Stokes pulse with duration τ s = 4 ns, t 0 = 5 ns, peak power P s = 10 mW, extinction ratio (a) ER = 15 dB (b) ER = 50 dB, (c) ER = 50 dB with back view of Brillouin spectra, for the fiber length L = 15 m and n=190 discretization points.