File pressflow.c

RCS Header: /cvsroot/rheoplast/pressflow.c,v 1.51 2006/03/08 11:50:00 hazelsct Exp

This provides rheoplast with all of the functions for modeling fluid flow using the velocity-pressure formulation.

The incompressible Navier-Stokes equations in velocity-pressure form are annoying because of spurrious modes in the pressure. The standard finite difference way around this is to use a "staggered mesh", in which the x-velocity mesh is offset from the pressure mesh by one-half grid spacing in the x-direction, y-velocity is offset in the y-direction, etc. This is a bit of a pain in a general parallel code like RheoPlast, particularly with regard to symmetry boundary conditions (hence the presence of the vortflow module), which must be specially designed for the staggered mesh. On the other hand, the shearstrain module also needs staggered meshes, so the infrastructure is necessary anyway; and velocity-pressure is more efficient in 3-D.

Included Files

Preprocessor definitions

#define __FUNCT__ "pressflow_first_setup"

#define __FUNCT__ "pressflow_labels_initcond"

#define u( point )

#define v( point )

#define w( point )

#define p( point )

#define omega( point )

#define ufunc( point )

#define vfunc( point )

#define wfunc( point )

#define pfunc( 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 DIRAC( width, jump, ycoord )

#define __FUNCT__ "pressflow_temp_parameters_line"

#define __FUNCT__ "pressflow_temp_parameters_boundary_line"

#define __FUNCT__ "pressflow_interior_line_function"

#define interp_function( x )

#define __FUNCT__ "pressflow_boundary_line_function"

#define __FUNCT__ "pressflow_interior_line_jacobian"

Local Variables

sftype
static SFType sftype

Global Function pressflow_boundary_line_function()

This calculates the time derivatives and constraint functions for the velocity-pressure form of the Navier-Stokes equations for an interior line. Continuity and x-motion equations are swapped to avoid zeroes on the Jacobian diagonal.

void pressflow_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 pressparm struct typedef, from which this gets needed parameters.
int side
Side on which to calculate the function values.
void pressflow_interior_line_function It returns nothing.
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.
The stagnation boundary condition is identical to the initial condition, with some special handling for pressure boundary conditions.

Global Function pressflow_first_setup()

The basic setup, setting the number of solved and temporary field variables and the stencil width.

void pressflow_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 pressparm structure this inserts parameters from the command line.

Global Function pressflow_interior_line_function()

This calculates the time derivatives and constraint functions for the velocity-pressure form of the Navier-Stokes equations for an interior line. Continuity and x-motion equations are swapped to avoid zeroes on the Jacobian diagonal.

void pressflow_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 pressparm struct typedef, from which this gets needed parameters.
The equations solved are given in dimensional form with non-unit density, starting with the equations of motion:
  • p: du/dt = -u.grad u - 1/rho dp/dx - div taux + Fx.

  • v: dv/dt = -u.grad v - 1/rho dp/dy - div tauy + Fy,
  • w: dw/dt = -u.grad w - 1/rho dp/dy - div tauy + Fy.
    Next there are two options for the pressure equation. If PRESSURE_MOMENTUM_DIVERGENCE is not defined, then we use the equation of continuity:
  • Div (u) = 0.
    The viscous stress divergence for a uniform-viscosity Newtonian fluid is the Laplacian of the velocity.
    The elastic shear strain term is documented in shearstrain.c, and its divergence is used here. The symmetry of the shear strain tensor, and the incompressibility condition (diagonals of shear strain sum to zero) reduce this to a function of just the xx- and xy-components of shear strain.
    This is also applied for Cahn-Hilliard with C taking the place of phi.
    With shear strain and phase field, we use an interpolation function p(phi) = 3phi2 - 2phi3 of the solid strength to weight the elastic and viscous stresses according to
    tau = p(phi) tauel + (1-p(phi)) tauvisc.

    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, if statsolid is set, then add the force function developed by Tonhardt and Amberg: (insert function here).
    Otherwise, add the body force due to interface curvature as described by Jacqmin given by 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 (density-normalized) x-driving force sinusoidal in y and starting at time t=t0 such that it can be 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.

    Global Function pressflow_interior_line_jacobian()

    This calculates the Jacobian of the equations corresponding to the velocity-pressure Navier-Stokes variables.

    void pressflow_interior_line_jacobian ( PetscScalar* x, PetscScalar* temp, Mat J, int points, int gxm, int gym, int firstrow, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data )

    PetscScalar* x
    The field variables from which to evaluate the Jacobian.
    PetscScalar* temp
    Array of temporary field variables.
    Mat J
    Where to put the evaluated Jacobian.
    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.
    int firstrow
    The matrix row number corresponding to the first point in the line.
    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 chparm struct typedef, from which this gets needed parameters.

    Global Function pressflow_labels_initcond()

    This sets up the field variable labels, uses command-line options to set the parameters in the pressparm struct typedef including the maximum stable explicit timestep size, and builds the initial condition for the velocity-pressure variables.

    void pressflow_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 pressparm 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.
    Membrane formation simulations use the Schmidt number (Sc, default 1) and a dimensionless force parameter (Fp, default 109) to scale the viscous and interfacial curvature fore terms respectively.
    The initial condition is zero velocity with (optional) stagnation flow inward in the x-direction and outward in the y-direction (and z direction in 3-D). With symmetry boundary conditions, the stagnation point is at the origin, otherwise it is at the center of the grid. Because of the staggered mesh, the extra delta x/2 shift needed for the vortflow module is unnecessary here.

    Global Function pressflow_temp_parameters_boundary_line()

    This calculates vorticity along the boundary.

    void pressflow_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 )

    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 pressparm struct typedef, from which this gets needed parameters.
    int side
    Side on which to calculate the function values.

    Global Function pressflow_temp_parameters_line()

    Though not needed for this module, the vorticity is helpful for others, so it is calculated here.

    void pressflow_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 ` +html+ <i>x</i>-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 pressparm struct typedef, from which this gets needed parameters.
    Vorticity is calculated at the corner of the cell, where it is best defined.