Equivolumetric Evolution of Planar Curves

and Applied Analysis 3 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (a) 0.2 0.25 0.3 0.35 0.4 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4


Introduction
Many physical phenomena such as flame propagation [1], crystal growth [2,3], and wear of rolling stones [4] can be described by curvature dependent front evolutions.For a given planar curve  0 (), the curvature dependent normal evolution (, ) is defined by the equation   (, ) =  ( (, )) n (, ) , where (, ) is the signed curvature of the front (, ) and n(, ) is the unit normal vector at (, ).The speed  of this evolution is a function in the curvature .
For a flame propagation model, it is natural to choose  as a decreasing function, since the region near a concave point, which is surrounded by the flame front, burns faster than the region near a convex point.A simple choice of a decreasing speed function is the linear model [1], in which  is a linear function of  with a negative slope.A drawback of the linear model is that the speed function  becomes negative for large curvature values.In this case, the flame front may move backward, which means that the region once burnt turns back into unburnt state.This is against physical laws.Moreover, the evidence of nonlinear relation between the curvature and the flame propagation speed has been reported [5].
The main purpose of this paper is to introduce a flame propagation model, which reflects the nonlinear contribution of the curvature to the normal flame propagation speed from a geometric viewpoint.We may assume that the heat energy produced by a segment of a flame front is proportional to the length of the segment, and this heat energy propagates along the normal direction of the flame front.So we propose a model of the evolution speed function, in which the infinitesimal area of the burning region is proportional to the length of the flame front.
This type of the curvature dependent evolution naturally fits into the problem domain of the level set method.After the seminal work by Osher and Sethian [6], the usefulness of the level set method is proved in many applications.One may consult [7,8] for comprehensive discussion on various numerical techniques using the level set method and their applications.
The rest of the paper is organized as follows.In Section 2, we present the backgrounds and motivation of this paper.Section 3 provides a rigorous definition and basic properties of the equivolumetric evolution.In Section 4, we formulate a level set method approach to the equivolumetric evolution.We also present some examples of numerical simulation of the equivolumetric evolution.Finally, Section 5 concludes the paper.

Backgrounds
2.1.Flame Propagation Models.Among many flame propagation models [9], the linear model, proposed by Markstein [1], is the most basic and widely used one.In this model, the flame speed  depends linearly on the signed curvature  of the flame front in such a way that where  0 is the speed of flat flame front and D is the mass diffusivity coefficient.So the flame near a concave point moves faster than the flat flame, and it moves slower near a convex point.
Markstein's linear model is a good description of the curvature dependency of the flame propagation speed for moderate ranges of the curvature of the flame front.However, if the absolute value of the curvature is large, the model generates undesirable behavior of the flame front.Suppose the curvature  is larger than 1/D; then the flame propagation speed  becomes negative, which means the flame moves backward changing the burnt region into the unburnt region.Hence, the model itself violates the entropy condition.This phenomenon can be easily detected near the sharp corner of the initial flame front.
Figure 1(a) illustrates the motion of the flame front under the linear model when the initial front is a rectangle.The rectangular initial flame front soon becomes a smooth curve and then moves outward approaching to a circular shape asymptotically.Figure 1(b) shows the motion of the front at a small region near a corner of the initial flame front with finer time steps.It is clearly noticeable that the flame front moves backward in a short period of time before it moves outward.
Echekki and Chen [5] provided the evidence of the nonlinear relation between the curvature of the flame front and the normal flame propagation speed using direct numerical simulations.The governing transport equation for the mass fraction  of the deficient reactant can be written as where u is the convective velocity,  is the density, and ω is the reaction rate.From this equation, they showed that the normal speed of the flame front could be expressed as where n is the unit normal vector of the flame front, which can be obtained by So the normal speed of flame front could be expressed as the sum of three terms: reaction, normal diffusion, and curvature.Although the explicit contribution of the curvature to the normal speed is linear, the other terms also depend on the curvature in a nonlinear manner.Their simulation results indicate that the normal speed remains positive for the region of the positive curvature, and the contribution of the curvature is much larger than the mass diffusivity coefficient D for the region of the negative curvature.

Equivolumetric Offsets.
The flame propagation model in this paper is motivated by the equivolumetric offsets, which is introduced in [10].The conventional offset curve, which is also called the equidistant offset curve, of a planar curve is the equidistant parallel curve of the original curve.For a given parameterized planar curve  with the parameter , its equidistant offset curve   of offset distance  is given by where n() is the unit normal vector field of .The equidistance offset curve is of great importance in the computer aided design (CAD) and the computer aided manufacturing (CAM).In particular, in the pocket machining application of CAM, the equidistant offset is the basic building block of the tool path generation for the contour parallel milling [11].
When we use the milling machine to cut a pocket with the boundary  out of a metal block, the center of the tool should trace the equidistant offset curve   where the offset distance  is the radius of the tool.
In the pocket machining practice, one of the most important engineering constraints is the limitation of the cutting force, which is proportional to the material removal rate.The machining scheme with the constant feedrate, which is the speed of the tool tip, ends up with the variable material removal rate due to the curvature effect.Thus, a variable feedrate scheme [12] is used to achieve the constant material removal rate.In the equidistance offset based machining, it is impossible to control both the feedrate and the material removal rate simultaneously.We can control both the feedrate and the material removal rate as constants by employing the equivolumetric offsets.
The equivolumetric offset γ of a given planar curve  is a variable distance offset whose offset distance is a function of the curvature ; that is, where the function (()) is elaborately chosen to compensate the volume variation due to the curvature effect.Then, for any interval [, ] in the parameter domain, the area between ([, ]) and γ([, ]) is proportional to the length of ([, ]).It is shown in [10] that the function () is given by a solution of a quadratic equation whose coefficients involve the curvature .The concept of equivolumetric offsets has been generalized to surfaces and tubular solids [13,14].

Equivolumetric Evolution of Planar Curves
For a given initial curve  0 (), which is a simple piecewise smooth closed planar curve with the parameter domain [0, ], thenormal evolution (, ) with the speed function  is characterized by where n(, ) is the outward unit normal vector of the front (, ).So the point (, ) on the front moves along the normal direction of the front with the speed (, ).Let t(, ) denote the unit tangent vector field of the front, which is t(, ) =   (, )/‖  (, )‖.The signed curvature (, ) of the front (, ) is defined by where (, ) = ((, ), (, )).Since we are dealing with theoutward normal vector n(, ), the relations between the tangent vector and the normal vector are The speed function (, ) may depend on the position of the front and the curvature of the front.If the speed function (, ) is a constant, then the corresponding normal evolution is equivalent to the one parameter family of equidistant offsets of the initial front.Instead of dealing with the speed function at each point on the front, we here consider the volumetric speed of a segment of the front.For a fixed time  0 , let us consider a segment ([, ],  0 ) where [, ] ⊂ [0, ].The evolution velocities (,  0 )n(,  0 ) attached to each point of the segment form a vector field defined on ([, ],  0 ).This vector field is naturally embedded in the normal bundle of the front and the normal bundle of the front is again immersed in the plane.So the velocity vector field corresponding to the segment ([, ],  0 ) occupies some region in the plane.This region can be parameterized by for  ≤  ≤  and 0 ≤ V ≤ 1.
Definition 1.For a given normal evolution (, ) with the positive speed function (, ), the volumetric speed of the segment ([, ],  0 ) at  =  0 is defined by the area of the region occupied by the velocity vector field (,  0 )n(,  0 ) for  ≤  ≤  in the plane.
Figure 2(a) illustrates the velocity field of a constant speed normal evolution for the elliptical initial front.Two curve segments on the ellipse of the same arc-length are indicated by thick curves, and the corresponding volumetric speeds are the areas of the shaded regions.Although they have the same arc-lengths, their volumetric speeds are different due to the curvature effect.
The signed curvature (, ) of the evolution curve is negative at a concave point.If the speed function (, ) is bigger than the absolute value of the radius of curvature near a concave point, then the velocity vector field goes beyond the focal point of the normal bundle.In this case, the area element has negative sign and the ambiguity between oriented and nonoriented areas arises.To avoid this ambiguity, we restrict our discussion to the case in which the velocity vector field does not reach the focal point.If (, ) ≥ 0, then there is not a restriction on (, ), since the normal bundle does not have a focal point.If (, ) < 0, then the speed function should satisfy (, ) ≤ −1/(, ).These conditions can be written as Remark 2. The volumetric speed should not be confused with the infinitesimal volume variation of the normal evolution.For a small time interval [ 0 ,  0 + Δ], the normal evolution corresponding to the segment ([, ],  0 ) sweeps the region which is the total flux of the velocity field (,  0 )n(,  0 ) on the segment ([, ],  0 ).So the infinitesimal volume variation does not capture the curvature of the front.This is the main reason for us to use the volumetric speed measured in the normal bundle.
We now define the equivolumetric evolution of a planar curve.
Definition 3. The normal evolution (, ) of the initial curve (, 0) with the positive speed function (, ) is the equivolumetric evolution with the volumetric ratio  0 if the volumetric speed of a front segment ([, ], ) is proportional to the length of the segment with the proportionality constant  0 , for any time  and an arbitrary interval Since the volumetric speed of a segment ([, ], ) is affected by the curvature, the speed function of the equivolumetric evolution should be a function of the curvature to compensate the curvature effect.The following theorem provides an explicit expression of the speed function of the equivolumetric evolution.Theorem 4. If (, ) is the positive speed function of the equivolumetric evolution with the volumetric ratio  0 and (, ) satisfies condition (12), then (, ) is the function of the curvature (, ) given by provided that the curvature is bounded below by −1/2 0 .
So the speed function () is continuous at  = 0 although it is expressed as a special case in (14).Figure 3 illustrates the shape of the equivolumetric speed function () for a specific volumetric ratio.This curve illustrates a great similarity with the simulation result in [5], which shows the convex decreasing property of the flame speed as a function in the curvature.
For the equivolumetric evolution to be a well-posed problem, we need the variation diminishing property [15].The total variation Var() of the front at time  is defined by and the variation should decrease during the evolution.It has been shown in [15] that the variation diminishing property is true if   (0) ≤ 0. The speed function of the equivolumetric evolution clearly satisfies this condition since The speed function () given in (14) does not describe the evolution of all possible front because it is defined only for the curvature  greater than −1/2 0 .Thus, we need to extend the speed function to the curvature less than −1/2 0 in order to perform the evolution for arbitrary initial fronts.If  < −1/2 0 , then we cannot impose the equivolumetric condition, since (17) does not have real roots for .But we still want the extended speed function to be a continuous decreasing function in the curvature.One possible extension is the reflection of () in ( 14) with respect to the point (−1/2 0 , 2 0 ).The extended speed function F() is given by and F(0) =  0 and F(−1/ 0 ) = 3 0 .Although the denominators in the above expression are zeros at  = 0 or at  = −1/ 0 , the extended speed function is a continuous function on the whole range of the curvature.So we express the extended speed function for only two cases as in (21) for the rest of the paper.

Level Set Formulation and Numerical Results
The equivolumetric evolution fits naturally into the framework of the level set method [7,8].For the given initial front  0 (), let  0 (, ) be the initial level set function defined on R 2 which satisfies the following: (i)  0 () is the zero level set of  0 (, ); (ii)  0 (, ) is negative inside the initial front; and (iii)  0 (, ) is positive outside the initial front.Then the normal evolution of  0 () with the speed () can be obtained by the zero level set of the solution (, , ) of the partial differential equation with the initial condition  (, , 0) =  0 (, ) .
The curvature  of the front can be computed directly from  as the divergence of the unit gradient vector field of , which is For the equivolumetric evolution with the extended speed function (21), the level set PDE is given by Since the speed function of the equivolumetric evolution has a nonzero value  0 at  = 0, the equivolumetric evolution has both the convection and the diffusion.As mentioned in [6], to apply the general level set scheme, we need to separate the Figure 4: (a) The equivolumetric evolution starting from a rectangular initial flame is illustrated.The initial front is a square centered at the origin whose length of a side is 0.6.The evolution speed is given by ( 21) with the volumetric ratio  0 = 0.04, which is the same as the speed of planar flame in Figure 1.The sequence of fronts is plotted from  = 0 to  = 10 with Δ = 1.(b) This plot shows an enlargement of the motion of the front near the corner marked as a dotted box on (a).This is the motion of the front from  = 0 to  = 1 with Δ = 0.1.convection term and the diffusion term of the speed function, which yields The convection term  0 ‖∇‖ on the left should be approximated by the upwind difference scheme [7], and the diffusion term on the right should be approximated by the central difference scheme.We can employ the simple forward Euler scheme for the time evolution.We now provide the results of the numerical experiments of the equivolumetric evolution for various initial fronts.
Example 5. We choose the rectangular initial front for the comparison of the equivolumetric evolution with the linear combustion model discussed in Section 2. With the same rectangular initial front as in Figure 1, the result obtained by the level set method for the equivolumetric evolution is shown in Figure 4.In this computation, we used the signed distance function as the initial level set function.As expected, the equivolumetric evolution does not show any backward evolution near the corner.So this flame propagation model satisfies the entropy condition.Another noticeable difference between the linear model and the equivolumetric evolution is the duration of the existence of the flat front.Although the flame fronts approach to circular shapes asymptotically in both models, the equivolumetric evolution preserves the flat region of the front for longer period of time.It is because the backward evolution of the linear model causes more rapid smoothing of the front curvature.Example 6.We consider a seven-pointed star which is similar to the example presented in [6].The initial front is given by  () = 0.5 + 0.2 sin (7) , 0 ≤  ≤ 2,

Conclusion
We proposed a nonlinear flame propagation model, whose volumetric speed is proportional to the length of the flame front.We present an explicit expression of the curvature dependent speed function and the numerical results using the level set method.

Figure 1 :
Figure 1: (a) The front evolution starting from a rectangular initial front is illustrated.The initial front is a square centered at the origin whose length of a side is 0.6.The front evolution speed is  0 (1 − D) with  0 = 0.04 and D = 0.05.The sequence of fronts is plotted from  = 0 to  = 10 with Δ = 1.(b) This figure shows an enlargement of the motion of the front near the corner marked as a dotted box on (a).This is the motion of the front from  = 0 to  = 1 with Δ = 0.1.

Figure 2 :
Figure 2: The velocity fields of two types of normal evolution for the elliptical initial front are shown.(a) The velocity field of a constant speed normal evolution is plotted.Two segments corresponding to the shaded regions have the same arc-length.Due to the curvature effect, their areas are different.(b) The velocity field of an equivolumetric evolution is plotted.Two regions occupied by the velocity field corresponding to the segments of the same arc-length are the same.

Figure 3 :
Figure 3: The curvature dependent speed function of the equivolumetric evolution with the volumetric ratio  0 = 1 is shown.
in the polar coordinate.The initial level set function is chosen as (, ) =  0.5 + 0.2 sin (7) − 1,(28)where (, ) is the polar coordinate corresponding to (, ).The whole domain[−1, 1] × [−1, 1] is discretized by 201 mesh points per side.Figure5(a) illustrates the front evolution under the linear model with  0 = 0.04 and D = 0.05 and Figure5(b) illustrates the equivolumetric evolution with the volumetric ratio  0 = 0.04.The difference between these two simulation results is particularly noticeable near the point of high curvature.If we compare the first few evolution curves near the tips of the star-shaped initial front in Figures5(a) and 5(b), then we can notice that the linear model shows slower propagation than the equivolumetric evolution.This is because the linear model suggests very small or even negative evolution speed for the high curvature region.On the other hand, the equivolumetric evolution shows much regular behavior near the tips of the initial front.

Figure 5 :
Figure 5: (a) The front evolution under the linear model starting from a seven-pointed star shaped initial front is illustrated.The evolution speed is  0 (1 − D) with  0 = 0.04 and D = 0.05.The sequence of fronts is plotted from  = 0 to  = 10 with Δ = 1.(b) The equivolumetric evolution with the volumetric ratio  0 = 0.04 for the same initial front is plotted at the same time steps as in (a).