RCS Header: /cvsroot/rheoplast/vortflow.c,v 1.65 2004/08/23 16:47:51 hazelsct Exp
This provides rheoplast with all of the functions for modeling fluid flow
using the velocity-vorticity formulation (which seems easier than
velocity-pressure).
The incompressible Navier-Stokes equations in velocity-pressure form are a
pain in the neck to solve, because one must worry about spurrious modes in
the pressure. Sure, there are ways around it, like staggered meshes in
finite difference and Taylor-Hood elements in finite elements. But there are
alternate forms which don't require such tricks, including
velocity-vorticity, which is particularly useful in this phase field code
because the vorticity gives the rotation rate for the order parameter vector
and the elastic strain tensor.
For the velocity-vorticity formulation in cartesian coordinates, our
variables will be
u and v for x- and y-direction velocities,
and omega for
vorticity defined by
omega = delxu =
dv/dx - du/dy
The incompressible Navier-Stokes equations start with continuity:
Equation not yet translated, see LaTeX-based docs.
If we differentiate that with respect to
x, that becomes equivalent to
Equation not yet translated, see LaTeX-based docs.
The two middle terms are
d omega/dy,
so we can rewrite this as
Equation not yet translated, see LaTeX-based docs.
We can also differentiate
the continuity equation above with respect to y, and
through a similar manipulation end up with
Equation not yet translated, see LaTeX-based docs.
Next we turn to the incompressible equations of motion:
Equations not yet translated, see LaTeX-based docs.
Now we just subtract the
y-derivative of the x-momentum equation from the
x-derivative of the y-momentum equation.
The left side gives:
Equations not yet translated, see LaTeX-based docs.
The right side is a bit simpler:
Equations not yet translated, see LaTeX-based docs.
So when we put
the left and right sides
together, we get:
Equation not yet translated, see LaTeX-based docs.
and for a uniform-density uniform-viscosity fluid, this simplifies to
Equation not yet translated, see LaTeX-based docs.
The x and y vorticity equations of continuity, and the
vorticity equation of motion, together
comprise the velocity-vorticity form of the incompressible Navier-Stokes
equations. It is interesting to note that this form consists of two
equations of continuity and one of motion.
The velocity-vorticity formulation is used here in
RheoPlast
for a couple of reasons. It has the numerical advantage of no zeroes on the
diagonal, unlike velocity-pressure whose "pressure equation", which is zero
divergence of velocity, has no pressure in it. Also unlike
velocity-pressure, it has no spurious modes in difference equations with all
of the variables computed at each node, so we don't need a staggered mesh for
stability (though we do need a staggered mesh for shear strain when that is
included). Finally, the velocity and vorticity happen to be just the
parameters needed for the convective and rotational terms in the
vector-valued phase field and shear strain tensor field equations.
It remains to be seen whether these advantages will be sufficient in three
dimensions to justify the additional two fields and equations required for
the three components of the vorticity vector field, vs. the scalar pressure
field. If memory is not a problem, that should be the case.
Included Files
Preprocessor definitions
#define __FUNCT__ "vortflow_first_setup"
#define __FUNCT__ "vortflow_labels_initcond"
#define __FUNCT__ "vortflow_temp_parameters_line"
#define u( point )
#define v( point )
#define omega( point )
#define ufunc( point )
#define vfunc( point )
#define omegafunc( point )
#define phi( point )
#define gxx( point )
#define gxy( point )
#define C( point )
#define mu( point )
#define phi2( point )
#define phi3( point )
#define mu2( point )
#define mu3( point )
#define __FUNCT__ "vortflow_interior_line_function"
The basic setup, setting the number of solved and temporary field
variables and the stencil width.
void vortflow_first_setup ( PetscTruth threedee, int* vars, int* tempvars, int* stencilwid, AppCtx* data )
- PetscTruth threedee
- Request support for 3-D.
- int* vars
- Pointer to the number of solved field variables.
- int* tempvars
- Pointer to the number of temporary field variables.
- int* stencilwid
- Pointer to the stencil width.
- AppCtx* data
- Pointer to the
AppCtx
struct typedef, into whose
vortparm
structure this inserts parameters from the command line.
This calculates the time derivatives and constraint functions for the
velocity-vorticity form of the Navier-Stokes equations for an interior line.
Note the
shearstrain
module interaction programmed here;
void vortflow_interior_line_function ( PetscScalar* x, PetscScalar* func, PetscScalar* temp, PetscTruth** mixed_constraints, int points, int gxm, int gym, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data )
- PetscScalar* x
- The field variables from which to evaluate the function.
- PetscScalar* func
- Where to put the evaluated function.
- PetscScalar* temp
- Array of temporary field variables.
- PetscTruth** mixed_constraints
- Arrays of boolean variables indicating
constraint equations in mixed timestep-constraint fields.
- int points
- Number of points to evaluate at.
- int gxm
- The
x-width
of the ``local'' vector's array, including shadow nodes, for the
y-increment.
- int gym
- The
y-width
of the ``local'' vector's array, including shadow nodes, for the
z-increment.
- PetscScalar xmin
- First node
x-coordinate.
- PetscScalar xmax
- Last node plus one
x-coordinate.
- PetscScalar ycoord
- This line
y-coordinate.
- PetscScalar zcoord
- This line
z-coordinate.
- PetscScalar time
- Current simulation time.
- AppCtx* data
- Pointer to the main simulation parameter structure, which
includes the
vortparm
struct typedef, from which this gets needed parameters.
The
u equation is the constraint in the x-continuity equation above.
above, divided by deltax_m2 to bring the Jacobian in line for small
deltax.
Likewise, the
v equation is the constraint in the y-continuity equation
above, divided by deltay_m2 to bring the Jacobian in line for small
deltax.
The vorticity equation is the time-derivative in
the vorticity equation of motion above.
Here it is implemented in three parts: the convective terms, the stress
terms (viscous and, if shearstrain is present, elastic), and the body
force terms, including interface curvature from the Cahn-Hilliard
module.
For Cahn-Hilliard, we add the body force due to interface
curvature as described by Jacqmin \cite{jacqderive}:
negative C times the gradient of mu.
For the membrane, since it's basically Cahn-Hilliard, we add a
similar curvature body force given by
the sum of negative phii times the
gradient of mui.
I'm adding a
x-driving force sinusoidal in y and starting at
time t=t0
such that it is zero during any initial explicit timesteps. A force
amplitude of
4pi2 (meaning force curl amplitude of
8pi3)
gives a velocity amplitude of one, which should be good for what we're
looking for. The
-sineforcet0 and -sineforcemax
command-line switches control the force starting time and amplitude
respectively, and the latter is normalized by density in the function
vortflow_labels_initcond().
This uses options to set the parameters in the
vortparm
struct typedef, sets up the field variable labels, maximum stable explicit
deltat, and initial condition for the velocity-vorticity variables.
void vortflow_labels_initcond ( PetscScalar* globalarray, int nx, int ny, int nz, int xm, int ym, int zm, int xs, int ys, int zs, int vars, AppCtx* data, PetscScalar* max_explicit_deltat )
- PetscScalar* globalarray
- The global field array.
- int nx
- Overall
x-width
of the global array.
- int ny
- Overall
y-width
of the global array.
- int nz
- Overall
z-width
of the global array.
- int xm
- The
x-width
of the local part of the array.
- int ym
- The
y-width
of the local part of the array.
- int zm
- The
z-width
of the local part of the array.
- int xs
- The (integer)
x-coordinate
of the start of the local part of the array.
- int ys
- The (integer)
y-coordinate
of the start of the local part of the array.
- int zs
- The (integer)
z-coordinate
of the start of the local part of the array.
- int vars
- Total number of field variables to be solved.
- AppCtx* data
- Pointer to the
AppCtx
struct typedef, whose
vortparm
structure this uses for various purposes.
- PetscScalar* max_explicit_deltat
- Pointer to the largest allowable explicit
timestep size for this equation, which this function can set/modify.
There are no temporary field variables for velocity-vorticity flow, so this
does nothing.
void vortflow_temp_parameters_line ( PetscScalar* x, PetscScalar* temp, int points, int gxm, int gym, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data )
- PetscScalar* x
- Array with the "real" field variables.
- PetscScalar* temp
- Array with the temporary field variables.
- int points
- Number of points at which to calculate the temporary variables.
- int gxm
- The
x-width
of the ``local'' vector's array, including shadow nodes, for the
y-increment.
- int gym
- The
y-width
of the ``local'' vector's array, including shadow nodes, for the
z-increment.
- PetscScalar xmin
- First node
x-coordinate.
- PetscScalar xmax
- Last node plus one
x-coordinate.
- PetscScalar ycoord
- This line
y-coordinate.
- PetscScalar zcoord
- This line
z-coordinate.
- PetscScalar time
- Current simulation time.
- AppCtx* data
- Pointer to the main simulation parameter structure, which
includes the
vortparm
struct typedef, from which this gets needed parameters.