50
The sea is assumed horizontally homogeneous, so that the governing equations of the
problem read
au + f e x u = ~ (K au) ,
at
z
az U az
(22)
ap = ~(K ap) .
at
az c az
(23)
At the initial instant, the velocity u is zero and the density is a linear function of z, p(t=O,z) =
(_pr/" o 2/g)z + p(t=O,z=O), where No is the constant initial Brunt-Viiisala frequency. The fluid
is set into motion by a constant stress 't applied at the sea surface: at z = 0, Ku (au/az) = ~/p.
No mass or heat flux is imposed at the sea surface, implying that Kc (ap/az) = 0 at z = O. The
surface boundary conditions for the turbulent variables are (Blumberg and Mellor, 1987) q =
16.6 1i3 u and Z = 0, where u denotes the friction velocity, defined by u 2 = I~I/p. We take
u/ = 1O~ m 2 s-2 and N o
2 =* 10-4 s-2.
•
Discrete model. Equations (18), (19), (22) and (23) are solved numerically by a standard
finite volume technique. The only difficulty pertains to the turbulent variables: obtaining, from
the numerical algorithm, negative values of lor lz would be devastating. Thus, the equations
of the turbulence closure model (18)-(19) must be discretized in such a way that the positivity
of q2 and lz be guaranteed.
To gain some insight into the problem of ensuring the positivity of the numerical solution of
an evolution equation, let us first examine a simple model, dljf/dt = -D(ljf), where D, the sink
term, is always positive. With an explicit time stepping, ljI'+ 1 = ljI' - tJ.t D( ljI') - where ljI' =
ljf(ntJ.t), n = 0, 1,2 ... -, it is required that ljI' be larger than tJ.tD(ljI') for ljI'+1 to be positive.
This condition may not always be verified. Patankar (1980) suggested an alternative
discretization, consisting in evaluating the destruction term in a pseudo-implicit way, viz
D(vr")vr"+I/vr". Accordingly, the modified scheme is vr"+1 = vr"/(1 + MD(vr")/vr"),
guaranteeing that ljI'+1 is always positive - provided ljI' is positive. This scheme is simple,
consistent and robust - in the sense that it is able to provide positive results, whatever the
value of ljI', D( ljI') and tJ.t. Consequently, Patankar (1980) pointed out, in a footnote, that "this
seemingly minor topic may tum out to contain the most valuable information" in his book.
When F is neglected, the equations of the closure scheme are of the form
aljf = P _ D + ~ (K aljf) ,
at
az
1{1 az
(24)
where ljfstands for lor lz; P (~O) and D(~O) denote appropriate source or sink terms.
It is desirable that the numerical scheme be able to cope with values of the dimensionless
ratio K ';t/!J.z 2 - where !J.z is the mesh size - that may be larger than 1/2. Consequently, to
avoid numerical instability, an implicit discretization of the diffusion term must be implemented
(e.g. Hirsch, 1988).
Using Patankar's (1980) trick, opting for a conservative and implicit discretization of the
diffusion operator, setting 8~l/2 = (Kl!ctl/2 tJ.t/ !J.z2, the discretized counterpart of (24) reads
- 8:-1/2 ~~i + (8:-1/2 + 1 + tJ.t D;/ ~ + 8;+1/2) ~+1 - 8k: 1/2 ~:~
The sea is assumed horizontally homogeneous, so that the governing equations of the
problem read
au + f e x u = ~ (K au) ,
at
z
az U az
(22)
ap = ~(K ap) .
at
az c az
(23)
At the initial instant, the velocity u is zero and the density is a linear function of z, p(t=O,z) =
(_pr/" o 2/g)z + p(t=O,z=O), where No is the constant initial Brunt-Viiisala frequency. The fluid
is set into motion by a constant stress 't applied at the sea surface: at z = 0, Ku (au/az) = ~/p.
No mass or heat flux is imposed at the sea surface, implying that Kc (ap/az) = 0 at z = O. The
surface boundary conditions for the turbulent variables are (Blumberg and Mellor, 1987) q =
16.6 1i3 u and Z = 0, where u denotes the friction velocity, defined by u 2 = I~I/p. We take
u/ = 1O~ m 2 s-2 and N o
2 =* 10-4 s-2.
•
Discrete model. Equations (18), (19), (22) and (23) are solved numerically by a standard
finite volume technique. The only difficulty pertains to the turbulent variables: obtaining, from
the numerical algorithm, negative values of lor lz would be devastating. Thus, the equations
of the turbulence closure model (18)-(19) must be discretized in such a way that the positivity
of q2 and lz be guaranteed.
To gain some insight into the problem of ensuring the positivity of the numerical solution of
an evolution equation, let us first examine a simple model, dljf/dt = -D(ljf), where D, the sink
term, is always positive. With an explicit time stepping, ljI'+ 1 = ljI' - tJ.t D( ljI') - where ljI' =
ljf(ntJ.t), n = 0, 1,2 ... -, it is required that ljI' be larger than tJ.tD(ljI') for ljI'+1 to be positive.
This condition may not always be verified. Patankar (1980) suggested an alternative
discretization, consisting in evaluating the destruction term in a pseudo-implicit way, viz
D(vr")vr"+I/vr". Accordingly, the modified scheme is vr"+1 = vr"/(1 + MD(vr")/vr"),
guaranteeing that ljI'+1 is always positive - provided ljI' is positive. This scheme is simple,
consistent and robust - in the sense that it is able to provide positive results, whatever the
value of ljI', D( ljI') and tJ.t. Consequently, Patankar (1980) pointed out, in a footnote, that "this
seemingly minor topic may tum out to contain the most valuable information" in his book.
When F is neglected, the equations of the closure scheme are of the form
aljf = P _ D + ~ (K aljf) ,
at
az
1{1 az
(24)
where ljfstands for lor lz; P (~O) and D(~O) denote appropriate source or sink terms.
It is desirable that the numerical scheme be able to cope with values of the dimensionless
ratio K ';t/!J.z 2 - where !J.z is the mesh size - that may be larger than 1/2. Consequently, to
avoid numerical instability, an implicit discretization of the diffusion term must be implemented
(e.g. Hirsch, 1988).
Using Patankar's (1980) trick, opting for a conservative and implicit discretization of the
diffusion operator, setting 8~l/2 = (Kl!ctl/2 tJ.t/ !J.z2, the discretized counterpart of (24) reads
- 8:-1/2 ~~i + (8:-1/2 + 1 + tJ.t D;/ ~ + 8;+1/2) ~+1 - 8k: 1/2 ~:~