Computation of the "Marsigli Flow"

What you see in this example is a numerical simulation of
Marsigli flow, which has been known since the 17th century (see Gill's book
"Atmosphere-Ocean Dynamics"):
It seems that when Marsigli went to Constantinople in 1679 he was told
about a well-known undercurrent in the Bosphorous: "... for
the fisherman of the towns on the Bosphorous say that the whole stream
does not flow in the direction of Byzantium, but while the upper current
which we can see plainly does flow in this direction, the deep water of
the abyss, as it is called, moves in a direction exactly opposite to
that of the upper current and so flows continuously against the current
which is seen". That is, the undercurrent water flows toward the Black Sea
from the Mediterranean. Marsigli reasoned that the effect was due to density
differences: water from the Black Sea is lighter than water from the
Mediterranean. The lower density of the Black Sea can be attributed to
lower salinity resulting from river runoff. He then performed a laboratory
experiment: A container is initially divided in two by a partition.
The left side contained water taken from the undercurrent in the Bosphorous,
while the right side contained dyed water having the density of surface
water in the Black Sea. The experiment was to put two holes in the partition
to observe the resulting flow. The flow through the lower hole
was in the direction of the undercurrent in the Bosphorous, while the flow
through the upper hole was in the direction of the surface flow.
Such process can be simulated by
Boussinesq flow with two initially piecewise constant temperatures in
an insulated box [0,8] X [0,1]. The above figure shows the visualization
of temperature profiles on the resolution 2048 X 256 at a sequence of times:
t=2, 4, 6, 8. The partition was
located at x=4. The temperature is chosen to be 1.5 at the left
half, which indicates the lower density, 1 at the right half, which
indicates the higher density. (By Boussinesq assumption, the density
difference can be converted into temperature difference with the
reverse ratio). The whole flow was at rest at t=0. A no-slip
boundary condition was imposed for the velocity and adiabatic boundary
condition was imposed for the temperature.
The physical parameters are chosen as: the Reynolds number
Re=5000, the Prandtl number Pr=1, and the Richardson
number (which corresponds to the gravity effect) Ri=4.
Like Riemann shock-tube problem, once the
partition was removed, the flow was driven by the gravity force. The
results indicated clearly the appearance of an upper current flow,
which moved from the left side to the right side, and an undercurrent
flow, which moved in the opposite direction. It coincided with the
phenomenon observed by Marsigli. Consequently, a sharp interface was
formed between the two currents. In other words, two currents with
different moving directions were separated by an interface. Strong
shear flow and vortex sheet came into being along the interface. This
vortex sheet exhibited the Kelvin-Helmholtz instability. As a result,
at t=2, two symmetric vortices and the rolling up structures were
formed. As the time goes on, more and more rolling-up structures were
generated and swelled.
The computational method is based on the fourth order scheme coupled
with 4-th order Runge-Kutta time stepping. Briley's
formula is used as the boundary condition for the vorticity.
The computations are repeated by using two resolutions:
2048 X 256 and 4096 X 512.
To see the details more clearly, we plot the temperature
and vorticity in a zooming region of [ 2.5, 3.5] X [ 0,1 ] at
t=6, on the resolution of 4096 X 512, in the following two pictures.
The computation results of temperature and vorticity
at time t=6 ona y=0.5 cut between two resolutions: 2048 X 256,
4096 X 512 are compared, which match perfectly well.
The flow behavior depends on many physics parameters, such as Reynolds number,
the ratio of densities, etc.

