On Jun 29, 2007, at 7:24 PM, ASReeve wrote:
Daniel,
Thanks for the reply. After playing with my code for the last
couple of days, everything seems to be working fine. I'm not sure
what I changed to correct the problem. Despite this, you raised a
few points I'd like to have clarified.
I'm attempting to simulate flow in a vertical profile, with
recharge across the top of the model. I'm treating recharge as a
simple source term. I was under the impression that the 'FixedFlux'
boundary condition only worked for no-flow boundaries. Am I mistaken?
It works if you don't have convection terms. In your case you only
have a diffusion term so you can use it.
Writing the 2-D equation out:
S dh/dt = d/dx { K_x (Z_sat / Z_cell) dh/dx} + d/dz(K_z dh/dz) +-
source
z is vertical, x is horizontal
S is a storage constant (units (1/length)
K is a conductivity (Length/time)
h is hydraulic head (Length)
Z_cell is the height of the cell
Z_sat is the water level relative to the cell (base of cell has 0
relative height), and is related to hydraulic head. I don't think
there is an easy way to linearize the equation.
I think I am correct in saying that this equation can be written in
the sharp interface limit independent of mesh quantities such that,
S dh/dt = d/dx { K_x H(z_water - z) dh/dx} + d/dz(K_z dh/dz) +-
source
where H is a Heaviside step function (H(z) = 1 when z > 0 and H(z) =
0 when z < 0) and z_water is the physical height of the water and
could be a function of x. I think this is a better way to think about
the problem.
The next step is how to represent H on a grid. I believe your
formulation above is the equivalent of a volume of fluid representation.
An alternative approach, often used in the level set method, is to
smooth H across the interface using a sine or cosine function.
Z_sat/Z_cell is a proportion of saturation, and as the water level
in the cell drops, the area that water can flow through drops and
the flux across the cell face drops. The method I'm using never
allows this to go to zero, instead, this ratio is made very small
when the cell dewaters.
In my situation, the source is recharge and only occurs at the top of
the model (where rain enters the system). Why did you suggest using a
'getDivergence' term in my model?
This is another way to create a boundary condition. I think in your
case, FixedFlux will work fine as you
don't have any convection terms.
Is this similar in some way to
assigning a source term,
Yes.
or was my previous e-mail so confusing that
you misunderstood what I was trying to do?
No. I understood what you were trying to do. I was suggesting an
alternative method for including exterior fluxes.
I'm unsure how to specify, in a generic way, where the top row of
cells are without directly indexing the cells. I understand I can
use the 'where' functionality to define parameters in a cell
variable, but don't know how to specify the coordinates for the
bottom and top of a particular row of cells.
You shouldn't be assigning a source only to the top row of cells.
Either use the FixedFlux boundary condition or assign a flux
over faces that is only non-zero on the exterior faces that you are
interested in and add flux.getDivergence() to your equation.
Something like this,
flux = VectorFaceVariable(mesh=mesh)
flux.setValue((0,fluxValue), where=m.getFaceCenters()[...,1] ==
heightOfDomain)
I'm also using indexing to retrieve the cell elvations using:
topBot=mesh.getFaceCenters()[...,1]
That's fine. There you are indexing for an x,y,z coordinate. Those
will always be the same regardless of the mesh you use.
and then use topBot[:nx*nz] to define the 'bottom elevation of each
cell and topBot[nx:nx*(nz+1)] to define the top elevation of each
cell. If this indexing isn't the best way to do this, what should I
be doing?
Shouldn't be indexing whatever you do. Use the Heaviside function.
z = mesh.getCellCenters()[...,1]
H = (z_water - z) > 0
For a structured mesh, I could just assign recharge using something
like
rech.setValue(1e-6/delZ, where=(mesh.getCellCenters[...,1]>=lenZ-
delZ/2.))
but this won't work in unstructured meshes.
Exactly. The point is that you should be able to change your mesh,
know nothing about the cell ordering,
and still solve your problem.
--
Daniel Wheeler
