# Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation with non-convex flow function.

код для вставкиСкачатьAppl. Num. Anal. Comp. Math. 1, No. 1, 113 – 127 (2004) / DOI 10.1002/anac.200310010 Characteristic directions approach to solving scalar onedimensional nonlinear advection equation with non-convex ow function D. N. Bokov∗ Russia, 456770, Snezhinsk, Chelyabinsk reg., Russian Federal Nuclear Center - VNIITF, P. O. Box 245 Received 30 June 2003, revised 30 October 2003, accepted 2 December 2003 Published online 15 March 2004 Key words Characteristics, advection equation, non-convex flow function, discontinuity decay Subject classi cation 76M20, 76R99 A concept of the characteristic technique used to obtain a generalized solution of the scalar one-dimensional nonlinear advection equation with the non-convex flow function is presented. Two grids: characteristic and Eulerian are used to obtain numerical solution. A characteristic grid is adaptive both to the properties of the initial distribution function and to the properties of the boundary condition function. This allows: development of the algorithm for obtaining a numerical solution on characteristic grid using the properties of the solution of nonlinear advection equation in smooth region; to reproduce spatial location and solution value at the discontinuity points and extreme points at the accuracy determined by interpolation and approximation of initial values and boundary condition functions. For the non-convex flow function, algorithms are proposed for the definition of the sequence of Riemann problems (strong discontinuity) and for their solving. Refined expressions are derived for the velocity of a strong non-stationary discontinuity. Construction of the solution with satisfying of integral preservation law for the non-convex flow function is presented. c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 1 Cauchy problem for the nonlinear advection equation Consider a nonlinear, homogeneous first-order partial differential equation: ∂u ∂F (u) + = 0, ∂t ∂x t > 0, x ∈ (−∞, ∞). (1) Solution to the Cauchy problem for this equation is function u(t, x), satisfying equation (1) at t > 0 and a specified function u(0, x) = u0 (x) at t = 0. Correctness of the Cauchy problem statement, the theorems of existence and uniqueness of the solution to continuous function u(t, x) are considered in [1]. The statement of the Cauchy problem for a class of discontinuous functions u(t, x) and two-times continuously - differentiable F (u(t, x)) is considered in [2]. Equation (1) can be written down as: ∂u ∂u + c(u) = 0, ∂t ∂x (2) where: c(u) = ∗ ∂F (u) . ∂u Corresponding author: e-mail: d.n.bokov@vniitf.ru, Phone: +7 35172 547 30, Fax: +7 35172 551 18 c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 114 D. N. Bokov: Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation . . . The solution of (1) in the class of continuous functions is reduced to the solution of the system of two ordinary first-order differential equations: dx = c(u), dt (3) du = 0, dt with the initial data: u(0, x) = u0 (x), x(0) = x0 . Line specified by the equation dx = c(u(t, x))dt, is a characteristic of the (2). Some authors understand the characteristic as a trajectory and value of u(t, x) on this trajectory [2]. In the class of discontinuous functions u(t, x) there are trajectories x = x(t), on which function u(t, x) has a discontinuity: u(t, x − 0) − u(t, x + 0) ≤ M < ∞. If point (t1 , x1 ) is on the discontinuity trajectory, there exist at least two characteristics, whose projections on a plane (t, x) pass through a point (t1 , x1 ) and all points of these characteristics at t < t1 belong to the solution. Besides this, the curvilinear integral F (u(x, t))dt − u(t, x)dx along the closed contour, formed by a straight line t = 0 and the projections of these two characteristics, is equal to zero [2]. Properties of the generalized solution of the Cauchy problem [2]: 1. Trajectory x = x(t) of the discontinuity points of the solution u(t, x) is single-valued and continuous function of t. 2. At each discontinuity point of the solution u(t, x) there are final limits u(t, x − 0) and u(t, x + 0). 3. At each point (t, x(t)) along a discontinuity trajectory the following relation is valid: lim t→t1 >t F (u(t, x + 0)) − F (u(t, x − 0)) x(t1 ) − x(t) = . t1 − t u(t, x + 0) − u(t, x − 0) (4) 2 Generalized solution of the nonlinear advection equation. Numerical algorithm The algorithm of a numerical solution is based on ”the use of scheme the exact solutions of elementary problems” [3]. The basic idea of characteristic directions approach consists in the use of the solution properties of the nonlinear advection equation: 1. Along each characteristic the solution has a constant value. This value depends on: • initial distribution u0 (x) in the area of smoothness; • solution to the Riemann problem; • boundary condition, when the characteristic is directed inwards the range of definition of the solution. 2. Along each characteristic (for the homogeneous equation) the velocity c(u(t, x)) remains constant, therefore the characteristic is a straight line, i.e. the system (3) is integrated precisely at any time step. 3. Discontinuities of the solution result from: • intersection two adjoining characteristics in the area of the solution smoothness; • interaction of discontinuities . 4. Discontinuities of the solution can be set in: • initial distribution function; • boundary condition function. 5. The spatial coordinate of the discontinuity is determined by integrating the discontinuity velocity, determined from a Hugoniot condition (4) [2]. 6. To determine the solution value at the discontinuity, property 2 section 1 is used. The values of the solution are limited by the solution values at the characteristics arriving at the discontinuity point. 7. The solution to the problem of arbitrary discontinuity breakdown in the case of alternating Fuu (u) uses the property of solution discontinuity stability for function Fuu (u) of constant sign [4] [5]. 2.1 Description of grids. Numerical approximation At the time moment t the numerical solution is constructed on two independent grid sets: • Eulerian grid for the result presentation - set of motionless points with coordinates xj , j = 0 . . . J. c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Appl. Num. Anal. Comp. Math. 1, No. 1 (2004) / www.interscience.wiley.com 115 • Characteristic grid - set of mobile points with coordinates xi , i = 0 . . . I. The complete solution is a union of the solutions on the above sets. The Eulerian grid is constructed, proceeding from the requirements to the result presentation, and does not depend on functions of initial distribution of the solution u0 (x) and initial distribution of the movement velocity c0 (x). The characteristic grid is constructed in the following way. Function of initial distribution of the solution u0 (x) in the area of smoothness is approximated at the given accuracy εu by a piecewise-linear function. Function of the velocity distribution c( x) in the smooth area is also approximated at the given accuracy εc by a piecewise-linear function. In both cases the approximating functions contain the extremum points of each function. Discontinuities, being the borders of the smooth areas, are included in the approximating functions as two points with identical spatial coordinates and known values of unilateral limits. This approximation can be made in different ways. The points of the approximating functions, if united in the ascending order of argument x, form a characteristic grid, adaptive to the properties of initial distribution function of the solution and function of initial distribution of the velocity of the characteristics. This allows the development of the algorithm to determine the numerical solution in the characteristic grid, proceeding from the properties of the solution of the nonlinear advection equation with no use of the points of the Eulerian grid. Linear interpolation of the piecewise-linear solution function in the characteristic grid determines the solution values at the points of the Eulerian grid. If a boundary Cauchy problem is considered, an additional relationship appears (X is a spatial coordinate of the boundary): u(t, X) = f (t), t > 0. Thus, the problem range of definition changes: c(u(t, x)) < 0, x ∈ (−∞, X], c(u(t, x)) ≥ 0, x ∈ [X, +∞). In order to correctly include the boundary condition effect, the following is done. The boundary condition function f (t) is approximated at the given accuracy εf by a piecewise-linear function, which contains extremum points and discontinuity points. The discontinuity is described by two nodes of the approximating function with equal time t∗ and values equal to f (t∗ − 0) and f (t∗ + 0). The nodes tk of the approximating function form a time grid. The complete solution shows at the approximation accuracy of initial distributions of the solution, and the velocity of the characteristics and boundary condition function: • spatial coordinates of discontinuities; • solution values at the extremum and discontinuity points. 2.2 Calculation of strong discontinuity Classification of the solutions on strong discontinuity [u− , u+ ] for the function of constant sign Fuu (u) is considered in [5]. The discontinuity is steady, if c− ≥ c+ . In a characteristic grid discontinuity is described by two adjacent points (fronts) with identical spatial coordinates and values equal to u− and u+ , respectively. If c− < c+ , the strong discontinuity is unstable and breaks up in centered rarefaction wave. In this case new points are added to the characteristic grid between the discontinuity fronts xnl and xnl+1 (the border characteristics) in the following way. Two sets of points (item 2.1) are considered which approximate functions F (u) and c(u). We unite these sets and choose only S points, which are within [u− , u+ ]. Assign the same spatial coordinate to these points: xnl = xnl+1 = · · · = xnl+S = xnl+S+1 . The added characteristic points form a fan of the centered rarefaction wave. In the case of alternating function Fuu (u) we deal with the Riemann problem of arbitrary discontinuity breakdown [6]. Find all roots of the equation Fuu (u) = 0 in the interval (u− , u+ ). Considering only those roots that are the extreme points for function c(u), order them in order of decreasing (increasing) argument u: u− > u1 > · · · > um > u+ (u− < u1 < · · · < um < u+ ). These roots divide the interval (u− , u+ ) into areas where function Fuu (u) has constant sign. Within each interval valid are the statements for single-signed function Fuu (u). Intervals where Fuu (u) is of constant sign are considered one after another. Thus, arbitrary discontinuity breaks into a series of alternating steady strong discontinuities and centered rarefaction waves, see Fig.1. The segments [u− , u1 ], [u2 , u3 ] and [u4 , u+ ] correspond to the steady strong discontinuities; [u1 , u2 ] and [u3 , u4 ] correspond to the centered rarefaction waves. Following [4], [6] and [7], take advantage of function F (u) diagram. Each steady strong c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 116 D. N. Bokov: Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation . . . F(u) 3,0 2,5 A 2,0 1,5 C B 1,0 0,5 0 u + u4 2 u3 4 u2 6 10 u- u1 8 u Fig. 1 Qualitative picture of the breakdown of an arbitrary discontinuity. discontinuity [ui , ui+1 ] ∈ [u− , u+ ] , can be geometrically presented as a chord, connecting points F (ui ) and F (ui+1 ) ( Fig.1 - chords A, B and C). To obtain a unique solution, the chord is iteratively moved towards a tangent to the function F (u) plot. If the steady strong discontinuity is the first (last) area where function Fuu (u) is single-signed, the tangent is built from the left (right) point of the interval, see Fig.2 (tangents A , B and C ). Designate arguments of the points of contact as follows: a subscript corresponds to the chord, and a superscript corresponds to the left (right) value of the interval. If two adjacent tangents intersect (see tangents B and C in Fig.2), the areas of these steady strong discontinuities, together with area between them, form a single area of the steady strong discontinuity, and a new tangent is constructed in it (tangent D in Fig.2). If two adjacent tangents have no intersections (tangents A and D in Fig.2), area is a centered rarefaction wave. The resultant configuration consists of the centered rarefaction waves and steady strong discontinuities, which are joined through the condition of contact. On the characteristic corresponding to the condition of contact, the velocities of the characteristic and the discontinuity are equal. 2.3 Construction of the numerical solution n n+1 - next moment of time; τ = tn+1 − tn Introduce the following designations: t - current moment of time; t ∗ time step; LIN T (y (φp , yp ), (φp+1 , yp+1 ))- operator of linear interpolation, defines value φ∗ at the point with coordinate y ∗ within the interval [yp , yp+1 ]. If tn = tk (tk - point of the time grid), a new point is added to the characteristic grid with the coordinates and parameters of the solution at the boundary at time moment tn . If tk is a discontinuity of the boundary condition function, the algorithm described in subsection 2.2 is used with u− = f (tn − 0), u+ = f (tn + 0). As tn+1 cannot be equal to the argument value in the time grid, the solution at the boundary point in this case is determined as un+1 = LIN T (tn+1 (f (tk ), tk ), (f (tk+1 ), tk+1 )). The boundary condition (including the area of its influence) is included by adding points to the characteristic grid at tk moments of time. In the area of smoothness the solution in the point of the characteristic grid remains constant along the characteristic coming from it (see section 1). Therefore, to obtain the solution at the next time c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Appl. Num. Anal. Comp. Math. 1, No. 1 (2004) / www.interscience.wiley.com F(u) 117 + - uD 3,0 uD 2,5 A’ D 2,0 B’ C’ A C 1,5 B 1,0 0,5 u 0 2 + u = uA + - uA uB 4 + 6 + uC uB - 8 10 - - uC = u Fig. 2 Qualitative picture of unique solution. moment it is enough to find spatial coordinate of the characteristic: xn+1 = xni + τ λi , i = 0 . . . I; i un+1 = uni , i = 0 . . . I. i If at time moment tn+1 two characteristics xl and xl+1 (not being points of discontinuity) intersect, the following xn+1 ) ≤ 0 , hence: cl ≥ cl+1 . Find the moment of intersection condition is true for them: (xnl+1 − xnl )(xn+1 l+1 − l ∗ n n with respect to the step τ : τ = (xl+1 − xl ) (cl − cl+1 ). After that run the calculation at τ = τ ∗. Values ul and ul+1 form a steady (cl ≥ cl+1 ) strong discontinuity. The trajectory of steady strong discontinuity is specified by the equation: x∗l = x∗l+1 = Dν τ + xnl+1 . How to determine the value of discontinuity velocity Dν and find the solution at discontinuity (u∗l and u∗l+1 ) is described in subsection 2.4. After the solution at discontinuity is obtained, a correction is made of the characteristic grid points corresponding to the steady strong discontinuity: = x∗l ; xn+1 l un+1 = u∗l ; l ∗ xn+1 l+1 = xl+1 ; ∗ un+1 l+1 = ul+1 . ∗ The interaction of the characteristics with the steady strong discontinuity consists in the following. If xn+1 l−1 = xl , n+1 ∗ point xn+1 l−1 is removed from characteristic grid after calculation of steady strong discontinuity. If xl−1 > xl , the time step τ is corrected so that the characteristic xl−1 reaches the steady strong discontinuity at time moment tn+1 : τ = (x∗l − xnl−1 ) (cl−1 − D). The same is done with characteristic xl+2 . If when the characteristic interacts with steady strong discontinuity, the condition of contact is fulfilled at the other front of discontinuity, then at the next time step a new characteristic is run from this front and with this front parameters. Such addition of the points to the characteristic grid allows a correct description of how the discontinuity interacts with the characteristics from the area of the solution smoothness. The interaction of steady strong discontinuities among themselves reduces to the Riemann problem within the interval corresponding to the parameters on the edge fronts of the interacting discontinuities. In case of the linear advection equation F (u) = cu, c ≡ const the velocities of the characteristics and discontinuities coincide and are equal c. All discontinuities are steady. The values at the nodes of piecewise-linear function are transferred along the x axis unchanged. Thus, the exact solution ui is c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 118 D. N. Bokov: Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation . . . obtained on the characteristic grid. The values un+1 at points of Eulerian grid xj are obtained with the accuracy j of linear interpolation over the points of the characteristic grid with the boundary points included: n+1 = LIN T (xj (un+1 , xn+1 ), (un+1 xn+1 ≤ xj ≤ xn+1 j = 0 . . . J. un+1 j i i i+1 , xi+1 )), i i+1 , 2.4 Steady strong discontinuity: velocity and solution To determine discontinuity velocity Dν , calculate a curvilinear integral F (u(t, x))dt − u(t, x)dx = 0 ABCDA over the closed contour formed by a straight line t = tn and projections of the characteristics, coming to point (tn+1 , x∗l+1 ) [2], see Fig.3. Values at the discontinuity are determined as follows: xl∗ = x l∗+1 t n +1 t n ∗ ul∗ A ul +1 B x L , uL xln , uln xln+1 , uln+1 D xR , u R C Fig. 3 Closed contour to determine the discontinuity velocity. n+1 n+1 n ∗ u∗l = LIN T (x∗l (unl−1 , xn+1 )), xn+1 ; l−1 ), (ul , xl l−1 ≤ xl ≤ xl n+1 n+1 n ∗ xn+1 u∗l+1 = LIN T (x∗l+1 (unl+1 , xn+1 l+1 ), (ul+2 , xl+2 )), l+1 ≤ xl+1 ≤ xl+2 . (5) (6) n Points (xL , uL ) and (xR , uR ) belong to the solution at t and characteristics running from them come to point (tn+1 , x∗l+1 ). Having calculated the curvilinear integral, obtain the relationship: Dv = 2(F (u∗l+1 ) − F (u∗l )) + c(u∗l+1 )(unl+1 − u∗l+1 ) − c(u∗l )(unl − u∗l ) . (u∗l+1 + unl+1 ) − (u∗l + unl ) (7) This relationship coincides with the Hugoniot condition both in case of the stationary discontinuity (un+1 = unl ; l n+1 ul+1 = unl+1 ) and in the limiting case, when τ → 0. To determine Dν , an iterative process is built using relationships (5), (6) and initial approximation D0 = (F (unl+1 ) − F (unl )) (unl+1 − unl ). c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Appl. Num. Anal. Comp. Math. 1, No. 1 (2004) / www.interscience.wiley.com 119 As mentioned in subsection 2.2 after the solution to the problem of arbitrary discontinuity breakdown is obtained, a configuration is formed with discontinuities, each having the satisfied condition of contact at one of its fronts (or on both fronts). The condition of contact is equivalent to the fact that the velocity of the characteristic coming from some point at this front is equal to the velocity of the discontinuity itself, i.e. for this front there are no characteristics coming onto it, see Fig.4 and Fig.5. Hence, for such discontinuities relationship (5) or t n+1 xl∗ = x l∗+1 uln +1 B uln +1 t n+1 ul∗+1 D C n l n l x ,u n l +1 x xR , u R n l +1 ,u ul∗+1 A B A tn xl∗ = x l∗+1 Fig. 4 Contact condition at the left front. tn D C x L , uL n l n l x ,u xln+1 , uln+1 Fig. 5 Contact condition at the right front. relationship (6) is replaced by the condition of contact and relationship (7) becomes like this: Dv = 2(F (u∗l+1 ) − F (unl )) + c(u∗l+1 )(unl+1 − u∗l+1 ) − c(unl )(u∗l − unl ) ; (u∗l+1 + unl+1 ) − (u∗l + unl ) (8) Dv = 2(F (unl+1 ) − F (u∗l )) + c(unl+1 )(u∗l+1 − unl+1 ) − c(u∗l )(unl − u∗l ) . (u∗l+1 + unl+1 ) − (u∗l + unl ) (9) or The initial approximation remains the same. As a result, obtain: D - velocity of the steady strong discontinuity, u∗l and u∗l+1 - values of the solution at the fronts of the steady strong discontinuity. If the condition of contact is fulfilled, then D = c(u∗l ) or D = c(u∗l+1 ). 2.5 Selection of the time step To choose the time step, time moments of the following events are considered with respect to tn : 1. t1 - intersection of two adjacent characteristics in the area of solution smooth; 2. t2 - intersection of the characteristic with the trajectory of strong discontinuity; 3. t3 - intersection of trajectories of two strong discontinuities; 4. t4 - time moment corresponding to the time grid point of the boundary condition function; 5. t5 - time moment of the result presentation. Time step is selected as τ = min (ti ). i=1,...,5 3 Conservatism A generalized solution of the advection equation needs the zero value of the contour integral taken along the closed contour [6]. Consider the following features of the proposed approach. Characteristic. The contour integral taken along the characteristic is equal to zero at any time step τ [6]. Hence, the integral will be equal to zero when integration is over time [t0 , T ] (range of definition of the time problem) along the characteristic trajectory. Steady strong discontinuity. The contour integral taken along the discontinuity trajectory is equal to zero [6]. c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 120 D. N. Bokov: Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation . . . This property was used to determine the discontinuity velocity. Hence, the integral will be equal to zero when integration is over time [t0 , T ] along the discontinuity trajectory. Consider the contour integral in the area of solution smoothness, formed by two adjacent characteristics at two consecutive time moments, see Fig.6. Consider how the law of conservation is fulfilled for contour ABCDA. t xin+1+1 , u i +1 xin +1 , ui n +1 B A C tn D xin+1 , u i +1 xin , ui Fig. 6 Closed contour in the area of solution smoothness. The contour integral can be presented as: B = ABCDA C + A D + B A + C . (10) D Consider each integral separately. n C t F (u(t, x))dt − u(t, x)dx = B Fi dt − ui dx = −τ Fi + ui ci τ. xB tn+1 n+1 t A F (u(t, x))dt − u(t, x)dx = D xC xA Fi+1 dt − tn ui+1 dx = τ Fi+1 − ui+1 ci+1 τ. xD Since sections [A, B] and [C, D] are on the appropriate identical time layers, the integrals over time will be equal to zero, i.e.: B xB F (u(t, x))dt − u(t, x)dx = − A u(t, x)dx, xA D xD F (u(t, x))dt − u(t, x)dx = − C u(t, x)dx. xC c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Appl. Num. Anal. Comp. Math. 1, No. 1 (2004) / www.interscience.wiley.com 121 Suppose, that distribution u(tn , x), x ∈ [xi , xi+1 ] is linear then: xD u(t, x)dx = 0.5(ui + ui+1 )(xni+1 − xni ) xC is an exact relationship. Assume, that distribution u(tn+1 , x), x ∈ [xi , xi+1 ] is also linear. Hence: xD u(t, x)dx = 0.5(ui + ui+1 )(xni+1 − xni ), xC or, writing this relationship through xn , obtain: xB u(t, x)dx = 0.5τ (ui + ui+1 )(ci+1 − ci ) + 0.5(ui + ui+1 )(xni+1 − xni ). xA Then after substitution of above formulas in formula (10) it will be written in the following form: = 0.5τ {(Fi+1 + ci+1 (ui − ui+1 ) − Fi ) − (Fi + ci (ui+1 − ui ) − Fi+1 )}. (11) ABCDA Take a Taylor series expansion of functions F (ui ) and F )ui+1 ) in the vicinity of points ui+1 and ui , respectively, the reminder term being in Lagrangian form: F (ui+1 ) = Fi + ci (ui+1 − ui ) + 0.5Fuu (ũi )(ui+1 − ui )2 , ũi ∈ [ui , ui+1 ]. 2 F (ui ) = Fi+1 + ci+1 (ui − ui+1 ) + 0.5Fuu (ũi+1 )(ui − ui+1 ) , ũi+1 ∈ [ui , ui+1 ]. Then the equality (11) will be as follows: = 0.25τ (Fuu (ũi ) − Fuu (ũi+1 ))(ui+1 − ui )2 . ABCDA Consequence 1. If ui = ui+1 , then ABCDA constant at any flow function F (u). Consequence 2. If Fuu (u) = const, then (12) = 0, i.e. constant distribution u(tn , x), x ∈ [xi , xi+1 ] remains ABCDA = 0, i.e. linear initial distribution u(0, x), x ∈ [xi , xi+1 ] remains linear at any moment of time. Consequence 3. Fuu (u) = const. From relation (12) follows, that linear distribution u(tn , x), x ∈ [xi , xi+1 ] becomes nonlinear at tn+1 moment of time. Approximate nonlinear distribution u(tn+1 , x), x ∈ [xi , xi+1 ] by two segments of the straight lines running from points xi and xi+1 and intersecting at point n+1 (x̄n+1 i+0.5 , ūi+0.5 ), n+1 x̄n+1 , xn+1 i+0.5 ∈ [xi i+1 ] so that to compensate deviation from linear distribution, see Fig.7. Suppose, that n+1 (x̄n+1 + xn+1 i+0.5 = 0.5(xi i+1 ). Then, proceeding from equality of contour integral to zero, obtain a formula to determine ūn+1 i+0.5 : ūn+1 i+0.5 = 0.5(ui + ui+1 ) + xn+1 i+1 δτ , − xn+1 i c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 122 D. N. Bokov: Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation . . . u x n +1 i + 0.5 ,u xin++11 , uin+1 n +1 i + 0.5 xin +1 , uin xin+1 , u in+1 xin , uin x n i + 0.5 ,u n i + 0.5 x Fig. 7 A qualitative picture of nonlinear distribution approximation. where δ = 0.5{(Fi+1 + ci+1 (ui − ui+1 ) − Fi ) − (Fi + ci (ui+1 − ui ) − Fi+1 )}. From a condition of conservation of the solution monotony, require that the following condition is met: xn+1 i+1 δτ < 0.5(ui+1 − ui ). − xn+1 i Consider: xn+1 i+1 δτ δ δτ < ; n+1 = (xn n ci+1 − ci − xi i+1 − xi ) + τ (ci+1 − ci ) 0.25(Fuu (u∗i+1 ) − Fuu (u∗i ))(ui+1 − ui )2 Fuuu (u∗∗ ) δ (ui+1 − ui )2 . = ∼ 0.25 ci+1 − ci ci+1 − ci Fuu (u∗ ) Then, Fuu (u) = 0 in the interval (ui , ui+1 ), as the value of the solution u, at which Fuu (u) = 0, is a point of the characteristic grid. Since Fuu (u) is a continuously-differentiable function within the interval (ui , ui+1 ), then Fuuu (u) is limited within it. Hence, if the following is valid F (u∗∗ ) uuu (13) 0.5 ui+1 − ui < 1, ∗ Fuu (u ) the solution remains monotonous. Whether relation (13) is valid is checked when the characteristic grid is formed. n+1 If relation (13) is violated, point (x̄n+1 i+0.5 , ūi+0.5 ) is added to the characteristic grid. From the above it follows that the contour integral taken along the characteristic grid, complemented by the central points (x̄, ū), will be equal to zero, that is a necessary condition of the generalized solution [8]. Consider a case, when one of the contour characteristics is front of a steady strong discontinuity, see Fig.8 and Fig.9. For the case shown in Fig.8 the initial contour ABCDA is divided into two contours: contour ABCD A, which is calculated as described above, and contour AD DA, which is a part of the closed contour equal to zero (see subsection 2.4). For the c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Appl. Num. Anal. Comp. Math. 1, No. 1 (2004) / www.interscience.wiley.com A B t n+1 tn t n+1 123 B A’ A tn C D’ D Fig. 8 CB - trajectory of the characteristic. DA - trajectory of the steady strong discontinuity. C D Fig. 9 CB - trajectory of the characteristic. DA - trajectory of the steady strong discontinuity. case presented in Fig.9 the initial contour ABCDA is divided into two contours: contour A BCDA, which is calculated as described above, and contour AA DA, which is a part of the closed contour equal to zero (see subsection 2.4). The other strong discontinuity front is treated similarly. If the time moment comes, when the value of local disbalance (disbalance within one interval of the characteristic grid) exceeds the permitted value εbal , the central point of this interval is added into the characteristic grid. A generalized numerical solution at Fuu (u) = const is obtained by uniting the points of the characteristic grid and points (x̄, ū) at each interval. It guarantees conservatism within the whole range of definition of the problem. Analysis of the numerical results on test problems showed, that the value of disbalance over the contour of the characteristic grid is only 10−5 %. 4 Conclusions The paper proposes the numerical algorithm of obtaining the generalized solution of one-dimensional scalar nonlinear advection equation with non-convex flow function. The algorithms are described of obtaining the solution at the characteristic grid points in the area of the solution smoothness, on the steady strong discontinuities and any point of space within the range of definition of the problem. It is shown how to determine the strong discontinuity velocity at the finite time step. For non-convex flow function the algorithms are proposed to determine the configuration (in case of strong discontinuity breakdown) and the solution at the configuration elements. The algorithms of constructing a characteristic grid are presented which allow description of initial distribution function at the given accuracy. It is shown how to construct the grid where the integral law of conservation for the non-convex flow function is valid. The method presented gives approaches to the solution of the following problems: 1. obtain the solution in case of ”thickening of peculiar features”, i.e. when there are several discontinuities of different types at one point in space; 2. conservation of extreme points of the solution distribution function; 3. conservation of the solution properties in asymptotics, i.e. during evolution of the solution on the large number of scales with respect to the initial distribution [9]; 4. grid adaptation in compliance with the solution properties. Unlike the finite-difference methods, the numerical solution obtained with the proposed method has the following properties: 1. decision does not depend on the difference of the scales of the spatial grid; 2. time step is determined by the interaction of the solution elements and does not depend on the spatial grid; 3. there is no loss of resolution, i.e. it is an impossible situation that the spatial distribution of the solution elements is smaller than the spatial step of the grid; 4. in the area of smoothness, addition of ”new” points to the characteristic grid does not change the solution value at the ”old” points, but improves the approximation accuracy of the solution distribution profile. c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 124 D. N. Bokov: Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation . . . For a linear flow function the technique considered gives an exact solution on the characteristic grid. For a nonlinear flow function the solution accuracy is determined by the accuracy of approximation of the initial solution distribution and the integral law of conservation is fulfilled on the extended grid. 5 Test problems Test 1. Problem statement is taken from [10]. A non-viscous Burgers equation is considered: ∂u ∂F (u) + = 0, ∂t ∂x where F (u) = 0.5u2 . Initial distribution: u(0, x) = sin(πx + π), x ∈ [−1, 1]; 0, x∈ / [−1, 1]. To approximate the initial distribution, 32 points were selected uniformly located along the spatial coordinate x were chosen. Fig.10 shows the profiles of the solution distribution at four moments of time. The steady strong discontinuity is formed at t = 0.3183. The discontinuity reaches a maximal amplitude at t = 0.5. u 1,0 0,5 0,0 t=0 t = 0.3183 t = 0.5 t=1 -0,5 -1,0 x -1,0 -0,8 -0,6 -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 Fig. 10 Profiles of the solution distribution for different time moments. Test 2. Problem statement is taken from [11]. The Buckley-Leverett equation is considered: ∂u ∂F (u) + = 0, ∂t ∂x where F (u) = 4u2 . 4u2 + (1 − u)2 Initial distribution: u(0, x) = 1, x ∈ [−0.5, 0]; 0, x ∈ / [−0.5, 0]. c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Appl. Num. Anal. Comp. Math. 1, No. 1 (2004) / www.interscience.wiley.com 125 u 1,0 t=0 t = 0.1 t = 0.3 t = 0.45 0,8 0,6 0,4 0,2 0,0 x -0,5 0,0 0,5 1,0 Fig. 11 Profiles of the solution distribution for different time moments. u 2,0 1,5 1,0 0,5 0,0 t=0 t = 0.5 t = 1.0 t = 1.65 t = 2.3 -0,5 -1,0 -1,5 -2,0 x -1 0 1 2 3 4 5 6 7 8 9 Fig. 12 Profiles of distribution u(t, x) at 5 time moments prior to the first interaction of the strong discontinuities. Function F (u) has one inflection point within the interval u ∈ [0, 1]. To approximate the initial distribution, 6 points are used corresponding to the borders of areas where u(0, x) is specified. Calculation of the discontinuity breakdown at point x = −0.5 gives a centered rarefaction wave with an upward convexity [11], containing 12 points of the characteristic grid and steady strong discontinuity. Calculation of the discontinuity breakdown at c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 126 D. N. Bokov: Characteristic directions approach to solving scalar one-dimensional nonlinear advection equation . . . u 2,0 A 1,5 B C 1,0 C B B t = 2.315 t = 5.0 t = 8.079 0,5 0,0 -0,5 x -4,5 -3,0 -1,5 0,0 1,5 3,0 4,5 Fig. 13 Profiles of distribution u(t, x) at 3 time moments. At time moments t = 2.315 and t = 8.079 the steady strong discontinuities interacted with each other. point x = 0 gives a centered rarefaction wave with a downward convexity [11], containing 12 points of the characteristic grid and steady strong discontinuity. Fig.11 shows the profiles of distribution u(t, x) at 4 moments of time. Test 3. Problem statement is taken from [12]. A nonlinear advection equation is considered: ∂u ∂F (u) + = 0, ∂t ∂x where F (u) = 0.25(u2 − 1)(u2 − 4). Initial distribution: 2, u(0, x) = x < 0; − 2, x ∈ [0, 2]; 2, x > 2. The function F (u) has two inflection points within the interval u ∈ [−2, 2]. To approximate the initial distribution, 6 points are used corresponding to the borders of areas where u(0, x) is specified. Calculation of the discontinuity breakdown at point x = 0 gives two steady strong discontinuities separated by a centered rarefaction wave. Fan of the centered rarefaction wave contains 12 points of the characteristic grid. Calculation of the discontinuity breakdown at point x = 2 gives two centered rarefaction waves separated with a contact discontinuity. Fans of centered rarefaction waves contain 12 points of the characteristic grid each. Fig.12 shows the profiles of distribution u(t, x) at 5 time moments prior to the first interaction of the strong discontinuities. Fig.13 presents the profiles of distribution u(t, x) at 3 time moments. At time moments t = 2.315 and t = 8.079 the steady strong discontinuities interacted with each other. In Fig.13 an area between the characteristics marked as A and B is a centered rarefaction wave formed due to discontinuity breakdown. An area between the characteristics marked as B and C is a rarefaction wave c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Appl. Num. Anal. Comp. Math. 1, No. 1 (2004) / www.interscience.wiley.com 127 u 2,0 1,9 1,8 t = 14.25 t = 74.25 t = 134.25 t = 194.25 t = 254.25 t = 314.25 1,7 1,6 x 0 200 400 600 800 1000 Fig. 14 Behavior of the solution u(t, x) at t → ∞. resulting from the interaction of steady strong discontinuity with the characteristics from the area of the solution smoothness before the discontinuity front. Fig.14 shows behavior of the solution u(t, x) at t → ∞. References [1] I.G.Petrovsky, Lectures about the equations with partial derivatives, Moscow, GITTL, 1950. (in Russian). [2] O.A.Oleinik, About a Cauchy problem for the nonlinear equations in a class of discontinuity functions, DAN, XCV(3), 451(1954). (in Russian). [3] S.K.Godunov, Memoirs about differential schemes. Novosibirsk, Scientific book, 1997. (in Russian). [4] I.M.Gelfand, Some problems of the theory of quasi-linear equations, UMN, XIV(2(86)), 87(1959). (in Russian). [5] B.L.Rozhdestvensky and N.N.Yanenko, Systems of quasi-linear equations and their appendix to gas dynamics, Moscow, Nauka, 1978. (in Russian). [6] O.A.Oleinik, About uniqueness and stability of the generalized solution of a Cauchy problem for quasi-linear equation, UMN, XIV(2(86)), 165(1959). (in Russian). [7] W.B.Lindquist, The scalar Riemann problem in two spatial dimensions: piecewise smoothness of solutions and its breakdown, SIAM Journal of Mathematical Analysis, 17(5), 1178(1986). [8] B.L.Rozhdestvensky, Discontinuity solutions of systems of quasi-linear equations of a hyperbolic type, UMN, XV(6(96)), 59(1960). (in Russian). [9] A.A.Samarsky, Differential approximation of convective advection with spatial splitting of time derivative, MM, 10(1), 86(1998). (in Russian). [10] A.Kurganov and E.Tadmor, New high-resolution central schemes for nonlinear conservation laws and convectiondiffusion equations, Journal of Computational Physics, 160, 241(2000). [11] C.-W.Shu and S.Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics, 77(2), 439(1988). [12] A.Harten, B.Engquist, S.Osher and S.R.Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes, III, Journal of Computational Physics, 71(2), 231(1987). c 2004 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

1/--страниц