close

Вход

Забыли?

вход по аккаунту

?

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
Документ
Категория
Без категории
Просмотров
4
Размер файла
219 Кб
Теги
flow, characteristics, dimensions, nonlinear, advection, convex, approach, equations, non, scala, one, directional, function, solving
1/--страниц
Пожаловаться на содержимое документа