DTIC AD1005707: Sources of Bias in Particle-Mesh Methods for PDF Models for Turbulent Flows

24 Pages·0.25 MB·English

by  Defense Technical Information Center

SOURCES OF BIAS IN PARTICLE(cid:2)MESH METHODS FOR
PDF MODELS FOR TURBULENT FLOWS
by
Jun Xu and Stephen B(cid:3) Pope
FDA (cid:4)(cid:5)(cid:2)(cid:6)(cid:7) January(cid:8) (cid:7)(cid:4)(cid:4)(cid:5)
Abstract
Numerical errors(cid:8) in particular bias(cid:8) in PDF(cid:2)based particle(cid:2)mesh methods for turbulence
modeling have been explored(cid:3) It is shown that bias decreases linearly with the increase of
the number of particles(cid:8) but increases with grid re(cid:9)nement(cid:3) The (cid:10)uctuations in mean (cid:9)elds
which are fed back into the coe(cid:11)cients of stochastic di(cid:12)erential equations are attributed to
be the sources of bias(cid:3) The Frozen Coe(cid:2)cient approach has been proposed and adopted to
pinpoint the sources of bias in detail(cid:3) These results provide guidelines for improving the
numerical accuracy of PDF models for turbulent (cid:10)ows(cid:3)
(cid:7)
INTRODUCTION
For complex turbulent reactive (cid:10)ows(cid:8) probability density function (cid:13)PDF(cid:14) methods have
been well developed and o(cid:12)er great potential (cid:15)(cid:16)(cid:17) (cid:15)(cid:18)(cid:17)(cid:3) In the application of these methods(cid:8) a
particle(cid:2)mesh method is used to solve a modeled equation for the evolution of a PDF(cid:8) e(cid:3)g(cid:3)(cid:8)
the joint velocity(cid:2)frequency(cid:2)composition PDF(cid:3) This is a Monte(cid:2)Carlo method in which the
PDF is represented by a large set of particles distributed in the physical space(cid:3) This space is
divided into a number of cells in order to estimate the mean (cid:9)elds such as the mean velocity
as a function of position(cid:3) Such a method has been implemented in the PDF(cid:19)DV code (cid:15)(cid:20)(cid:17)
which is a Fortran code to calculate the properties of statistically two(cid:2)dimensional (cid:13)plane or
axi(cid:2)symmetric(cid:14) turbulent reactive (cid:10)ows using the joint velocity(cid:2)frequency(cid:2)composition PDF
model(cid:3) Inthepresentwork(cid:8) numerical errors(cid:8) inparticularbias(cid:8) inPDF(cid:19)DVareinvestigated(cid:3)
ThenumericalaccuracyandconvergenceofPDFmethodshavebeendiscussed byPope(cid:15)(cid:21)(cid:17)
andWeltonet(cid:3)al(cid:3) (cid:15)(cid:4)(cid:17)(cid:3) Usually(cid:8)weakconvergence(cid:8) i(cid:3)e(cid:3)(cid:8)theconvergence ofexpectationsinstead
of the PDF(cid:8) is sought for PDF methods(cid:3) Four di(cid:12)erent types of numerical errors have been
identi(cid:9)ed by considering estimating a mean quantity(cid:22) statistical error(cid:8) spatial discretization
error(cid:8) temporal discretization error and bias(cid:3) The convergence of numerical solutions to the
modeledequations(cid:8) whicharestochasticdi(cid:12)erential equations(cid:13)SDE(cid:14)(cid:8)requiresthatnumerical
errors vanish as the particle number per cell N tends to in(cid:9)nity(cid:8) and as the time step (cid:23)t and
(cid:0)(cid:0)(cid:2)
grid size h approach to zero(cid:3) It has been shown that statistical error (cid:13)SE(cid:14) scales as (cid:7)(cid:2)N
(cid:15)(cid:21)(cid:17)(cid:3) It is reasonable to postulate that temporal error and spatial error(cid:8) behave as that in
(cid:9)nite di(cid:12)erence method(cid:8) i(cid:3)e(cid:3)(cid:8) both vanish as (cid:23)t and h tend to zero for (cid:9)xed N(cid:3)
However(cid:8) in early experiences with the PDF(cid:19)DV code it has been observed that there
appeared to be relatively large bias(cid:3) The bias is the deterministic error resulting from using
a (cid:9)nite number of particles(cid:3) Its features in the PDF methods are not very clear yet(cid:3) This
study is devoted to understanding the behavior and the sources of bias in PDF(cid:19)DV(cid:8) which
will provide one of the guidelines for improving the accuracy of the code(cid:3)
The description of bias in PDF(cid:19)DV is presented in the next section(cid:8) and then the strate(cid:2)
gies for pinpointing the sources of bias is described in the following section(cid:3) Two test cases
are used to search for the sources of bias(cid:22) stationary homogeneous turbulence and Couette
(cid:10)ow(cid:3) Results for these cases are discussed(cid:3) Conclusions are drawn in the (cid:9)nal section(cid:3)
(cid:19)
DESCRIPTION OF BIAS IN PDF(cid:2)DV
Before presenting the behavior of bias in PDF(cid:19)DV(cid:8) the model equations and some numerical
techniques are introduced(cid:3) A Lagrangian approach is taken both in the modeling and in the
numerical methodwhich isaMonte(cid:2)Carloparticle(cid:2)meshmethod(cid:3) The(cid:10)uid(cid:2)particlepossesses
(cid:3) (cid:3)
properties(cid:22) velocity u (cid:13)t(cid:14) and turbulence relaxation rate (cid:13)turbulence frequency(cid:14) (cid:3) (cid:13)t(cid:14) (cid:3)
(cid:0)
These properties are then modelled by the corresponding stochastic processes u (cid:13)t(cid:14) and
(cid:0)
(cid:3) (cid:13)t(cid:14) based on the idea of the stochastic Lagrangian modeling approach (cid:15)(cid:5)(cid:17)(cid:3) (cid:13)If combustion
(cid:0)
is considered(cid:8) the stochastic models for scalar (cid:9)elds are needed(cid:3)(cid:14) In PDF(cid:19)DV(cid:8) u (cid:13)t(cid:14) evolves
according to the simpli(cid:9)ed Langevin equation (cid:13)SLM(cid:8) (cid:15)(cid:5)(cid:17)(cid:14)
(cid:0) (cid:7) (cid:26) (cid:0) (cid:0)(cid:0)(cid:2)
du (cid:13)t(cid:14) (cid:24) (cid:0)rhpidt(cid:0) (cid:25) C(cid:4) (cid:27)(cid:13)u (cid:13)t(cid:14)(cid:0)hUi(cid:14)dt(cid:25)(cid:13)C(cid:4)k(cid:27)(cid:14) dW(cid:4) (cid:13)(cid:7)(cid:14)
(cid:19) (cid:16)
(cid:0) (cid:2)
(cid:0)
where W is Wiener process(cid:8) (cid:3) (cid:13)t(cid:14) is solved by the Ito stochastic di(cid:12)erential equation (cid:15)(cid:19)(cid:17)
(cid:0) (cid:0) (cid:0) (cid:2) (cid:0) (cid:0)(cid:0)(cid:2)
d(cid:3) (cid:24) (cid:0)C(cid:5)(cid:13)(cid:3) (cid:0)h(cid:3)i(cid:14)(cid:27)dt(cid:0)h(cid:3)i(cid:3) S(cid:2)dt(cid:25) (cid:19)C(cid:5)(cid:5) (cid:27)(cid:3) h(cid:3)i dW(cid:6) (cid:13)(cid:19)(cid:14)
(cid:3) (cid:4)
The non(cid:2)dimensional source term S(cid:2) is de(cid:9)ned as(cid:22)
(cid:2)
S(cid:2) (cid:24) C(cid:2) (cid:0)C(cid:0)SijSij(cid:2)h(cid:3)i (cid:4) (cid:13)(cid:26)(cid:14)
where Sij is the mean rate of strain(cid:3) For all other terms(cid:8) model constants and the de(cid:9)nition
for conditional mean of turbulence frequency (cid:27) in (cid:13)(cid:7)(cid:14) and (cid:13)(cid:19)(cid:14)(cid:8) refer to (cid:15)(cid:19)(cid:17) (cid:15)(cid:5)(cid:17)(cid:3)
In PDF(cid:19)DV(cid:8) the expectation of a random variable(cid:8) e(cid:3)g(cid:3) hUi(cid:8) is approximated by an
ensemble mean estimated by the cloud(cid:2)in(cid:2)cell method(cid:3) Because the number of samples(cid:8) i(cid:3)e(cid:3)(cid:8)
the number of particles per cell N is (cid:9)nite(cid:8) the ensemble mean itself is a random variable(cid:3) Of
course(cid:8) it also depends on time step (cid:23)t and grid size h(cid:3) Several numerical techniques(cid:8) e(cid:3)g(cid:3)(cid:8)
variance reduction (cid:13)VR(cid:14) and time averaging (cid:13)TAV(cid:14)(cid:8) are adopted to reduce the statistical
(cid:10)uctuations in the mean (cid:9)elds (cid:15)(cid:20)(cid:17)(cid:3)
The bias BQ of a statistics hQi(cid:8) e(cid:3)g(cid:3) hUi(cid:8) is the deterministic error caused by N being
(cid:9)nite(cid:3) Using hQiN(cid:3)(cid:6)t(cid:3)h to represent the ensemble(cid:2)average of Q calculated by (cid:9)nite N(cid:8) (cid:23)t and
h(cid:8) BQ can be written as(cid:8)
BQ (cid:24) hhQiN(cid:3)(cid:6)t(cid:3)hi(cid:0)hQi(cid:2)(cid:3)(cid:6)t(cid:3)h(cid:6) (cid:13)(cid:16)(cid:14)
(cid:26)
1.00 1.30
(a) (b)
0.95 1.25
T T
〉N 〉N
〉U ω〉
〈〈 〈〈
0.90 1.20
0.85 1.15
0.000 0.002 0.004 0.006 0.008 0.010 0.000 0.002 0.004 0.006 0.008 0.010
1/N 1/N
Figure (cid:7)(cid:22) Time(cid:2)averaged mean values vs(cid:3) N in Couette (cid:10)ow(cid:22) (cid:13)a(cid:14) Mean Velocity(cid:28) (cid:13)b(cid:14) Mean
h(cid:3)i(cid:3) Grids in y direction Ng (cid:24) (cid:19)(cid:7)(cid:8) both VR and TAV are o(cid:12)(cid:8) y(cid:2)H (cid:24) (cid:6)(cid:6)(cid:6)(cid:26)(cid:26)(cid:3)
(cid:3)(cid:0)
A simple analysis by Pope (cid:15)(cid:21)(cid:17) suggests that bias BQ scales as N (cid:8) so BQ can be expressed
as(cid:22)
b
BQ (cid:24) (cid:4) (cid:13)(cid:18)(cid:14)
N
where b may be a function of other factors a(cid:12)ecting bias(cid:8) e(cid:3)g(cid:3)(cid:8) grid size and other numerical
techniques in the code(cid:3)
As an example(cid:8) Couette (cid:10)ow (cid:13)for (cid:10)ow descriptions refer to appendix A(cid:14) is calculated
using PDF(cid:19)DV(cid:3) With the same grid size(cid:8) time step and numerical techniques(cid:8) the time
averages hhUiNiT and hh(cid:3)iNiT of one point(cid:8) which are equivalent with the ensemble mean
of hUiN and h(cid:3)iN respectively according to the ergodic assumption(cid:8) are shown in Figure (cid:7)
as a function of the number of particles per cell(cid:3) Since temporal error and spatial error are
independent of particle number(cid:8) the linear relationship in the (cid:9)gure implies that bias scales
as (cid:13)(cid:18)(cid:14)(cid:3) This behavior of bias is common in PDF(cid:19)DV(cid:3) The scale of bias as (cid:13)(cid:18)(cid:14) ensures that
bias approaches zero as N goes to in(cid:9)nity and thus leads to convergence of the scheme with
respect to particle number per cell(cid:3)
On the other hand(cid:8) the slopes of the lines in Figure (cid:7) provide the value of b in (cid:13)(cid:18)(cid:14) which
determines the magnitude of bias(cid:3) For speci(cid:9)c (cid:23)t and h (cid:13)and for all other aspects of the
(cid:16)
102 102
(a) (b)
〉U 101 ω〉 101
b〈 b〈
100 100
10-2 10-1 10-2 10-1
1/Ng 1/Ng
Figure (cid:19)(cid:22) Bias vs(cid:3) grid size in Couette (cid:10)ow(cid:22) (cid:13)a(cid:14) Mean Velocity(cid:8) the slope is (cid:0)(cid:6)(cid:6)(cid:18)(cid:18)(cid:28) (cid:13)b(cid:14) Mean
h(cid:3)i(cid:8) the slope is (cid:0)(cid:6)(cid:6)(cid:20)(cid:6)(cid:3) Both VR and TAV are o(cid:12)(cid:8) y(cid:2)H (cid:24) (cid:6)(cid:6)(cid:6)(cid:26)(cid:26)(cid:3)
numerical method (cid:9)xed(cid:14)(cid:8) we have
(cid:7) (cid:7)
hQiN(cid:0)(cid:3)(cid:6)t(cid:3)h(cid:0)hQiN(cid:2)(cid:3)(cid:6)t(cid:3)h (cid:24) BQN(cid:0) (cid:0)BQN(cid:2) (cid:24) b (cid:0) (cid:4) (cid:13)(cid:21)(cid:14)
N(cid:0) N(cid:2)
(cid:0) (cid:2)
A formula for b is thus obtained
N(cid:2)N(cid:0)
b(cid:13)(cid:23)t(cid:4)h(cid:14) (cid:24) (cid:13)hQiN(cid:0)(cid:3)(cid:6)t(cid:3)h (cid:0)hQiN(cid:2)(cid:3)(cid:6)t(cid:3)h(cid:14)(cid:6) (cid:13)(cid:5)(cid:14)
N(cid:2) (cid:0)N(cid:0)
Consequently(cid:8) using two di(cid:12)erent particle numbers per cell(cid:8) b and BQ can be calculated for
a speci(cid:9)c grid size(cid:3) By (cid:9)xing the time step and the numerical techniques(cid:8) b for di(cid:12)erent
grid sizes in Couette (cid:10)ow is obtained by (cid:13)(cid:5)(cid:14) and is plotted against (cid:7)(cid:2)Ng in Figure (cid:19)(cid:8) where
Ng is the number of cells or grids in the domain(cid:3) Since nonuniform grids are used in this
calculation(cid:8) (cid:7)(cid:2)Ng is used here to represent the averaged grid size(cid:3) Figure (cid:19) shows that b
increases with grid size(cid:3) This behavior is observed in homogeneous stationary turbulence
(cid:13)described in appendix A(cid:14) as well(cid:3) In this case(cid:8) bias is the deviation of turbulence energy
and turbulence frequency from the stationary solutions(cid:3) As shown in Figure (cid:26)(cid:8) increasing
the number of cells results in larger bias(cid:3)
The fact that the bias increases with grid size is of great concern to PDF(cid:19)DV(cid:3) If b
explodes faster with grid size than the convergence of bias due to increasing N(cid:8) bias will not
(cid:18)
1.50 0.02
(a) (b)
1.40
k 1.30 ω〉
as of 〈s of 0.01
Bi 1.20 Bia
1.10
1.00 0.00
0 2 4 6 0 2 4 6
Δy Δy
Figure (cid:26)(cid:22) Bias vs(cid:3) grid size in Homogeneous Stationary Turbulence(cid:22) (cid:13)a(cid:14) Turbulence Energy
k(cid:28) (cid:13)b(cid:14) Mean h(cid:3)i(cid:3) Time step (cid:23)t (cid:24) (cid:6)(cid:6)(cid:7)s(cid:8) for (cid:16)(cid:6)(cid:6)(cid:6) steps(cid:3)
vanish as h approaches zero and N goes to in(cid:9)nity(cid:3) In the view of numerical computations(cid:8)
to minimize the total numerical error(cid:8) one might want to increase Ng to reduce the spatial
discretization error(cid:8) which unfortunately causes larger bias(cid:3) Therefore(cid:8) N must increase as
Ng increases in order to prevent the bias from exploding and to get the total error converged(cid:8)
which thus leads to unacceptably high computational cost(cid:3)
STRATEGIES FOR EXPLORING BIAS
Theoretically(cid:8) the statistical (cid:10)uctuations in the coe(cid:11)cients of thestochastic di(cid:12)erential equa(cid:2)
tions(cid:8) e(cid:3)g(cid:3)(cid:8) the Langevin equation(cid:8) are the major sources of bias(cid:3) The numerical solutions to
(cid:13)(cid:7)(cid:14) and (cid:13)(cid:19)(cid:14) are essentially di(cid:12)erent from the following standard problem(cid:22) given coe(cid:11)cients
a(cid:13)x(cid:4)t(cid:14) and b(cid:13)x(cid:4)t(cid:14)(cid:8) an initial condition X(cid:13)(cid:6)(cid:14) (cid:24) x(cid:4)(cid:8) and a stopping time T (cid:7) (cid:6)(cid:8) integrate the
stochastic di(cid:12)erential equation
dX(cid:13)t(cid:14) (cid:24) a(cid:13)X(cid:13)t(cid:14)(cid:4)t(cid:14)dt(cid:25)b(cid:13)X(cid:13)t(cid:14)(cid:4)t(cid:14)dW(cid:13)t(cid:14)(cid:4) (cid:13)(cid:20)(cid:14)
which has been well studied (cid:15)(cid:26)(cid:17)(cid:3) This is because in (cid:13)(cid:7)(cid:14) and (cid:13)(cid:19)(cid:14)(cid:8) the coe(cid:11)cients depend on
the mean of a function of the process(cid:8) e(cid:3)g(cid:3)(cid:8) hUi(cid:8) which needs to be estimated in numerical
(cid:21)
Sources of Bias
Fixed 〈U〉 〈ω〉 k
Randomness of Mean Fields
SDE’s: Coefficients Numerical Techniques
Fixed Fixed Fixed Fixed
〈U〉 〈ω〉 〈U〉 〈U〉
〈ω〉 k
Vel. Cor.
Drift Term Sω Diffusion Terms Estimate VR, TAV Position Cor.
〈U〉 〈ω〉 〈ω〉 〈uu〉 etc, k 2nd Moments Pres. Algo.
...
∂〈U〉
i
Fixed
∂x
j
Grid Size Dependence of Bias
Figure (cid:16)(cid:22) Strategies for exploring bias in PDF(cid:19)DV
computation and inevitably carries numerical (cid:10)uctuations(cid:3) The mean (cid:9)elds in PDF(cid:19)DV are
calculated by cloud(cid:2)in(cid:2)cell method and fed back into the coe(cid:11)cients(cid:3) This feedback causes
bias(cid:3) The two questions to be addressed are(cid:22)
(cid:7)(cid:3) Which coe(cid:11)cients associated with mean (cid:9)elds yield bias (cid:29)
(cid:19)(cid:3) Why does bias increase with (cid:9)ner grids (cid:29)
For the PDF(cid:19)DV code(cid:8) all mean (cid:9)elds(cid:8) terms in the SDE(cid:30)s and numerical techniques
which may be related to the above two issues are sketched in Figure (cid:16)(cid:3)
If(cid:8) instead(cid:8) the coe(cid:11)cients were non(cid:2)random(cid:8) independent of particle properties(cid:8) the
computed particle properties at any time would be independent(cid:8) identically distributed(cid:8)
(cid:5)
and independent of N(cid:3) It follows then there would be no bias(cid:3) Therefore(cid:8) if the estimate
of mean (cid:9)eld is replaced by a non(cid:2)random input (cid:13)say (cid:9)xed or frozen hUi(cid:14)(cid:8) there is no
correspondingstatistical(cid:10)uctuationinthatmean(cid:9)eldandrelatedtermsinthe SDE(cid:30)ssothat
the contribution of statistical errors in the mean (cid:9)eld to bias will be prevented(cid:3) According
to this Frozen Coe(cid:2)cient approach(cid:8) a mean (cid:9)eld or a term in the SDE(cid:30)s is justi(cid:9)ed as a
source of bias if freezing it results in no bias(cid:3) The things to be tested(cid:8) which are possibly
the sources of bias(cid:8) and associated methods(cid:8) are also shown in Figure (cid:16)(cid:3)
NUMERICAL TESTS AND RESULTS
The approach of Frozen Coe(cid:2)cient is accomplished step by step through numerical tests(cid:3)
Because the (cid:10)uctuations in the mean (cid:9)elds are suspected to be the major sources of bias(cid:8)
(cid:9)rst of all the calculations are made by (cid:9)xing mean (cid:9)elds in all coe(cid:11)cients of the SDE(cid:30)s to
check the behavior of bias(cid:3) If it disappears(cid:8) then we can conclude that bias arises entirely
from the (cid:10)uctuations in the coe(cid:11)cients and no more calculations are needed to search for
other sources of bias due to the numerical techniques(cid:3) Several cases are calculated and
compared to the base case(cid:28) namely(cid:8) the result from the original code without VR and TAV(cid:3)
The conditions for di(cid:12)erent cases are presented in Table I(cid:3)
(cid:0) Sources of Bias(cid:2) General Views
Inthecalculationofcase (cid:7)(cid:8)itis nottheestimatesfromcloud(cid:2)in(cid:2)cellmethodthatareusedbut
constant values for mean (cid:9)elds(cid:3) No bias is expected if bias is only related to the statistical
(cid:10)uctuations in the mean (cid:9)elds(cid:3) The result compared to the base case (cid:13)case (cid:6)(cid:14) is shown
in Table I(cid:3) Obviously(cid:8) bias of velocity and h(cid:3)i are much smaller than case (cid:6) and almost
zero after (cid:9)xing mean (cid:9)elds globally(cid:3) Therefore bias of these two variables is because of
the (cid:10)uctuations in the estimates of mean (cid:9)elds(cid:3) However(cid:8) because the ensemble mean of
velocity is used to estimate second moments(cid:8) bias of turbulent energy is increased instead(cid:8)
which implies that the (cid:10)uctuations in the (cid:9)rst moments are the source of bias for the second
moments(cid:3) This will be discussed further in other cases(cid:3)
(cid:20)
(cid:3) Source of Bias for Velocity
Fixing Mean Fields for Coe(cid:2)cients in Eq(cid:3)(cid:4)(cid:5)(cid:6)
To clarify further that bias of velocity comes from Eq(cid:3)(cid:13)(cid:7)(cid:14)(cid:8) case (cid:19) is calculated by (cid:9)xing
all mean (cid:9)elds in this equation(cid:3) As in Table I(cid:8) this case shows that the bias of the mean
velocity is in the same order as case (cid:7)(cid:8) which indicates that the bias of velocity is due to
the (cid:10)uctuations in the mean (cid:9)elds fed back into the coe(cid:11)cients of Eq(cid:3)(cid:13)(cid:7)(cid:14)(cid:3) This is consistent
with the analysis in the previous section(cid:3) To distinguish the in(cid:10)uence on bias of drift term
from di(cid:12)usion term in Eq(cid:3)(cid:13)(cid:7)(cid:14)(cid:8) the following cases are designed(cid:3)
Fixing Mean Velocity for Drift Term
The mean velocity is a critical variable in Eq(cid:3)(cid:13)(cid:7)(cid:14) because the particle velocity is forced to
relax to it(cid:8) and the particle velocity(cid:8) in turn(cid:8) is used to estimate the mean velocity(cid:3) This
means that there must be a very strong interaction between them(cid:3) As shown in Table I
(cid:13)case (cid:26)(cid:14)(cid:8) the bias of velocity resulting from (cid:9)xing the mean velocity for the drift term in
Eq(cid:3)(cid:13)(cid:7)(cid:14) is very small(cid:3) Hence(cid:8) the (cid:10)uctuation in mean velocity is a major source for bias of
velocity(cid:3) One more interesting observation is that the bias of h(cid:3)i becomes very small in this
case as well(cid:3) It seems that the behavior of the bias of h(cid:3)i is dominated by the mean velocity
although Eq(cid:3)(cid:13)(cid:19)(cid:14) is the same in the form as Eq(cid:3)(cid:13)(cid:7)(cid:14)(cid:3) This will be further discussed later(cid:3)
It may be noticed that in case (cid:26) there is a huge bias of k(cid:3) However(cid:8) this phenomenon
should not be paid too much attention(cid:3) When the drift term is frozen while the di(cid:12)usion
(cid:5) (cid:0)
term is not(cid:8) in the equation for k derived from Eq(cid:3)(cid:13)(cid:7)(cid:14)(cid:8) (cid:0)(cid:7)C(cid:4)(cid:27)(cid:13)u (cid:0)hUi(cid:14)dt in the drift term
will not still balance with the di(cid:12)usion term(cid:8) which may cause a large bias of k(cid:3)
Fixing (cid:27) and k
In order to determine the e(cid:12)ect of (cid:10)uctuations in the di(cid:12)usion coe(cid:11)cients(cid:8) a calculation was
attempted in which (cid:27) and k were (cid:9)xed only in the di(cid:12)usion term in Eq(cid:3)(cid:13)(cid:7)(cid:14)(cid:3) It turns out that
the solution is not stable(cid:3) The explanation could be that di(cid:12)erent values of (cid:27) in the drift
and di(cid:12)usion terms give rise to an instability of numerical solutions to the equation(cid:3) This
case is then modi(cid:9)ed to let (cid:27) take (cid:9)xed values both for the drift term and the di(cid:12)usion term