File vortflow.c

RCS Header: /cvsroot/rheoplast/vortflow.c,v 1.85 2006/03/06 18:15:03 wanida 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 u( point )

#define v( point )

#define omega( point )

#define vis( 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 C2( point )

#define C3( point )

#define Mu2( point )

#define Mu3( point )

#define __FUNCT__ "vortflow_temp_parameters_line"

#define __FUNCT__ "vortflow_temp_parameters_boundary_line"

#define __FUNCT__ "vortflow_interior_line_function"

#define interp_function( x )

This calculates the time derivatives and constraint functions for the velocity-vorticity form of the Navier-Stokes equations for a boundary line.
void vortflow_boundary_line_function It returns nothing.
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.
int side Side on which to calculate the function values.

#define __FUNCT__ "vortflow_boundary_line_function"

Global Function vortflow_boundary_line_function()

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_boundary_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, int side )

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.
int side
 
void vortflow_interior_line_function It returns nothing.

Global Function vortflow_first_setup()

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.

Global Function vortflow_interior_line_function()

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.
Dynamics of the elastic shear strain term are documented in shearstrain.c. The curl of the shear stress divergence, which is necessary here, is given by an equation you can find in the LaTeX-based documentation. The symmetry of the shear strain tensor, and the incompressibility condition (diagonals of shear strain sum to zero) reduce this easily to a function of just the xx- and xy-components of shear strain. Because the shear strain is on a staggered mesh, this becomes a somewhat complex thing whose second derivatives are calculated on a stencil four nodes across.
With shear strain and phase field, we use an interpolation function of the order parameter to weight the elastic and viscous stresses according to
tau = p(phi) tauel + (1-p(phi)) tauvisc.
This interpolation function is most simply expressed as the integral of (phi-0.1)(1-phi) scaled so that it is 0 for phi<0.1 and 1 for phi>1 (details in the LaTeX-based documentation).
This is also applied for Cahn-Hilliard with C taking the place of phi.
With shear strain and no phase field, we can use this for incompressible viscoelastic mechanics by applying both the elastic and viscous stresses.
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().
For the chternary, since it's basically Cahn-Hilliard, we add a similar curvature body force given by the sum of negative Ci times the gradient of mui.

Global 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.

Global Function vortflow_temp_parameters_boundary_line()

void vortflow_temp_parameters_boundary_line ( PetscScalar* x, PetscScalar* temp, int points, int gxm, int gym, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data, int side )

Global Function vortflow_temp_parameters_line()

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.