Stokes Equations

Interface

class pystopt.stokes.StokesInterface(ambient_dim: int)

Generic interface for solutions and representations of Stokes flow.

ambient_dim: int

Dimension of the ambient space.

__init__(ambient_dim: int) → None

The various functionality associated with the StokesInterface is implemented through a singledispatch() setup. Basically all interested parties can implement a subset of the following functions by registering it, e.g. using

@velocity.register(MyType)
def _velocity_interface(self: MyType, *args, **kwargs):
    # compute some symbolic or explicit value related to the velocity
    return u

All the classes in this section implement a subset of this interface with different input arguments.

pystopt.stokes.streamfunction(self: StokesInterface, *args, **kwargs)

pystopt.stokes.velocity(self: StokesInterface, *args, **kwargs)

pystopt.stokes.pressure(self: StokesInterface, *args, **kwargs)

pystopt.stokes.traction(self: StokesInterface, *args, **kwargs)

The traction is the normal component of the stress tensor, i.e. \(\mathbf{n} \cdot \sigma \triangleq n_i \sigma_{ij}\).

pystopt.stokes.tangent_traction(self: StokesInterface, *args, **kwargs)

This is the tangential component of the stress tensor, i.e. \(\mathbf{t}^\alpha \cdot \sigma \triangleq t_i^\alpha \sigma_{ij}\), where \(\alpha \in \{1, \dots, d - 1\}\) is a basis for the tangent space.

pystopt.stokes.normal_velocity_gradient(self: StokesInterface, *args, **kwargs)

This is the normal component of the velocity gradient, i.e. \(\mathbf{n} \cdot \nabla \mathbf{u} \triangleq n_i \nabla_i u_j\).

pystopt.stokes.tangent_velocity_gradient(self: StokesInterface, *args, **kwargs)

This is the tangential component of the velocity gradient, i.e. \(\mathbf{t}^\alpha \cdot \nabla \mathbf{u} \triangleq t_i^\alpha \nabla_i u_j\).

pystopt.stokes.normal_velocity_laplacian(self: StokesInterface, *args, **kwargs)

This is the normal component of the velocity Laplacian, i.e. \(\mathbf{n} \cdot \Delta \mathbf{u} \triangleq n_j \nabla_i \nabla_i u_j\).

pystopt.stokes.normal_velocity_hessian(self: StokesInterface, *args, **kwargs)

This is the normal component of the velocity Hessian, i.e. \(\mathbf{n} \cdot \nabla \nabla \mathbf{u} \cdot \mathbf{n} \triangleq n_i n_k \nabla_i \nabla_k u_j\).

pystopt.stokes.deltaf(self: StokesInterface, *args, **kwargs)

This is the jump in traction at the interface. In non-dimensional variables, for a viscosity ratio \(\lambda\), it is

\[[[\mathbf{n} \cdot \sigma]] = \mathbf{n} \cdot \sigma_+ - \lambda \mathbf{n} \cdot \sigma_+.\]

pystopt.stokes.deltau(self: StokesInterface, *args, **kwargs)

This is the jump in velocity at the interface. In all current applications, this is zero.

Layer Potential Kernels

All of the classes in this section implement pystopt.stokes.StokesInterface with the arguments

@velocity.register(Wrapper)
def _wrapper_velocity(
        self: Wrapper, density, *, qbx_forced_limit=None):
    pass

and similar for the other variables. The only departure from this interface are pystopt.stokes.tangent_traction() and pystopt.stokes.tangent_velocity_gradient() which have

@tangent_traction.register(Wrapper)
def _wrapper_traction(
        self: Wrapper, density, *, tangent_idx, qbx_forced_limit=None):
    pass

to specific which tangent vector in the basis to use. The tangent basis is obtained by pytential.symbolic.primitives.tangential_onb().

class pystopt.stokes.Stresslet(ambient_dim: int)

Wrapper for sumpy.kernel.StressletKernel and pytential.symbolic.stokes.StressletWrapper.

The wrapper handles the application of the Stresslet kernel \(T_{ijk}\), the associated pressure kernel \(\Pi_{ij}\) and associated stress kernel \(\Theta_{ijkl}\) to a given density. The notation and definitions can be found in [Pozrikidis1992].

ambient_dim

class pystopt.stokes.Stokeslet(ambient_dim)

Wrapper for sumpy.kernel.StokesletKernel and pytential.symbolic.stokes.StokesletWrapper.

The wrapper handles the application of the Stokeslet kernel \(G_{ij}\), the associated pressure kernel \(P_{i}\) and associated stress kernel \(T_{ijk}\) to a given density. The notation and definitions can be found in [Pozrikidis1992].

class pystopt.stokes.StokesletJumpRelations(ambient_dim: int)

Jump relations for all the layer potential operators constructed in Stokeslet.

class pystopt.stokes.StokesKernel(ambient_dim: int, base_kernel: Any)

Bases: StokesInterface

class pystopt.stokes.StokesJumpRelations(ambient_dim: int)

Bases: StokesInterface

Generic Stokes jump relations.

Layer Potentials

The TwoPhaseStokesRepresentation implements pystopt.stokes.StokesInterface with the following arguments

@velocity.register(TwoPhaseStokesRepresentation)
def _velocity_repr(self, density, *, side=None, qbx_forced_limit=None):
    pass

and similar for the other variables. It also implements the functionality for pystopt.operators.RepresentationInterface.

class pystopt.stokes.TwoPhaseStokesRepresentation(knl: StokesKernel, jmp: StokesJumpRelations, bc: TwoPhaseStokesBoundaryCondition, kernel_arguments: Dict[str, Any] | None = None)

Bases: RepresentationInterface, StokesInterface

__init__(knl: StokesKernel, jmp: StokesJumpRelations, bc: TwoPhaseStokesBoundaryCondition, kernel_arguments: Dict[str, Any] | None = None) → None

Parameters:

• knl – a Stokes kernel-like Stokeslet.

• jmp – jump relations for the kernels in knl.

• bc – farfield boundary conditions.

class pystopt.stokes.TwoPhaseSingleLayerStokesRepresentation(bc: TwoPhaseStokesBoundaryCondition, kernel_arguments: Dict[str, Any] | None = None)

Bases: TwoPhaseStokesRepresentation

Constructs a complete representation (layer potentials and jumps) for a single-layer (Stokeslet) representation of the velocity

\[u_i(\mathbf{x}) = \int_\Sigma G_{ij}(\mathbf{x}, \mathbf{y}) q_j(\mathbf{y}) \,\mathrm{d}S.\]

and all the variables described by Stokeslet. The construction largely corresponds to Section 5.3 from [Pozrikidis1992].

Solving

The StokesSolution class implements pystopt.stokes.StokesInterface functions with the arguments

@velocity.register(StokesSolution)
def _solution_velocity(self: StokesSolution, *,
        side=None, qbx_forced_limit=None, dofdesc=None):
    pass

and similar for the other variables.

class pystopt.stokes.StokesSolution(actx: ArrayContext, source: pytential.symbolic.primitives.DOFDescriptor, places: pytential.GeometryCollection, density: ndarray, op: TwoPhaseStokesRepresentation, context: Dict[str, Any], density_name: str = 'q')

Bases: object

Represents the result of solving a Stokes layer potential, e.g. by single_layer_solve().

density

Solved density for the operator op.

op

class pystopt.stokes.TwoPhaseSingleLayerStokesSolver(sym_op: TwoPhaseStokesRepresentation, *, viscosity_ratio, context=None, gmres_arguments=None)

Simple wrapper around single_layer_solve().

__init__(sym_op: TwoPhaseStokesRepresentation, *, viscosity_ratio, context=None, gmres_arguments=None)

__call__(actx, places, *, source=None)

Solve the Stokes integral equation on the given source.

Returns:

a StokesSolution.

pystopt.stokes.single_layer_solve(actx, places, op, *, lambdas, sources=None, context=None, gmres_arguments=None)

Solves a two-phase Stokes problem.

Parameters:

• places – a GeometryCollection or any input that can be passed to its constructor.

• op – a TwoPhaseStokesRepresentation instance that describes the integral equation.

• sources – a list of identifiers in places on which to solve.

• lambdas – a list of viscosity ratios for each source domain.

Returns:

a StokesSolution.

Boundary Conditions

All subclasses of TwoPhaseStokesBoundaryCondition implement pystopt.stokes.StokesInterface with the arguments

@svelocity.register(TwoPhaseStokesBoundaryCondition)
def _bc_velocity(self: TwoPhaseStokesBoundaryCondition):
    pass

and similar for the remamining variables.

class pystopt.stokes.TractionJumpCondition(ambient_dim: int)

ambient_dim

__call__()

Symbolic expression of the jump condition.

class pystopt.stokes.YoungLaplaceJumpCondition(ambient_dim: int, capillary_number_name: str = 'ca')

Bases: TractionJumpCondition

Constant surface tension jump condition. It is given by

\[[\![ \mathbf{f} ]\!] = \frac{1}{\mathrm{Ca}} \kappa \mathbf{n}.\]

capillary_number_name

Name of the symbolic variable representing the Capillary number.

capillary_number

A Variable.

__init__(ambient_dim: int, capillary_number_name: str = 'ca') → None

class pystopt.stokes.VariableYoungLaplaceJumpCondition(ambient_dim: int, capillary_number_name: str = 'ca', surface_tension_name: str = 'gamma')

Bases: YoungLaplaceJumpCondition

Variable surface tension jump condition. It is given by

\[[\![ \mathbf{f} ]\!] = \frac{1}{\mathrm{Ca}} (\gamma \kappa \mathbf{n} + \nabla_\Sigma \gamma),\]

where \(\nabla_\Sigma\) is the tangential gradient on the surface \(\Sigma\).

surface_tension_name

Name of the symbolic variable representing the surface tension coefficient.

surface_tension

A Variable.

__init__(ambient_dim: int, capillary_number_name: str = 'ca', surface_tension_name: str = 'gamma') → None

class pystopt.stokes.GravityJumpCondition(ambient_dim: int, direction: int, bond_number_name: str = 'bo')

Bases: TractionJumpCondition

Gravity-based jump condition. It is given by

\[[\![ \mathbf{f} ]\!] = \frac{1}{\mathrm{Bo}} x_d \mathbf{n}.\]

direction

Integer indicating the axis direction in which gravity actx.

bond_number_name

Name of the symbolic variable representing the Bond number.

bond_number

A Variable.

__init__(ambient_dim: int, direction: int, bond_number_name: str = 'bo') → None

class pystopt.stokes.TwoPhaseStokesBoundaryCondition(ambient_dim: int, deltafs: Tuple[TractionJumpCondition, ...])

Bases: object

Boundary conditions for a generic two-phase Stokes flow.

ambient_dim

deltafs

A tuple of TractionJumpCondition that are summed together to obtain the full jump condition.

__init__(ambient_dim: int, deltafs: Tuple[TractionJumpCondition, ...]) → None

class pystopt.stokes.CapillaryBoundaryCondition(ambient_dim: int, deltafs: Tuple[TractionJumpCondition, ...])

Bases: TwoPhaseStokesBoundaryCondition

Boundary conditions for a two-phase Stokes flow in the presence of constant surface tension at the interface.

capillary_number

Symbolic variable representing the Capillary number.

class pystopt.stokes.FarfieldCapillaryBoundaryCondition(ambient_dim: int, deltafs: Tuple[TractionJumpCondition, ...], uinf: ndarray, pinf: pytential.symbolic.primitives.var, finf: ndarray, tinf: ndarray, dudninf: ndarray, dudtinf: ndarray, ddudninf: pytential.symbolic.primitives.var, dduinf: pytential.symbolic.primitives.var)

Bases: CapillaryBoundaryCondition

Farfield boundary conditions for constant surface tension Stokes flow.

pystopt.stokes.get_farfield_boundary_from_name(ambient_dim: int, name: str, **kwargs: Any) → FarfieldCapillaryBoundaryCondition

pystopt.stokes.get_extensional_farfield(ambient_dim: int, *, alpha: float = 0.5, pinf: float = 0.0, axis: int = 2, capillary_number_name: str = 'ca') → FarfieldCapillaryBoundaryCondition

Extensional farfield flow with constant surface tension

\[\begin{split}\mathbf{u}_\infty = \begin{cases} (\alpha x, -\alpha y), & \quad d = 2, \\ (\alpha x, \alpha y, -2 \alpha z), & \quad d = 3, \end{cases}\end{split}\]

where the location of the \(-2 \alpha\) can be controlled by the axis argument. Setting axis to \(-1\) extends the 2D case to 3D by padding with zero.

Parameters:

• alpha – strain rate.

• pinf – constant farfield pressure.

• axis – can define the plane in which the extensional flow acts for a three-dimensional flow, but is ignored in two dimensions.

Returns:

a FarfieldCapillaryBoundaryCondition.

pystopt.stokes.get_uniform_farfield(ambient_dim: int, *, uinf: float | ndarray = 1.0, pinf: float = 0.0, capillary_number_name: str = 'ca') → FarfieldCapillaryBoundaryCondition

Uniform farfield flow with constant surface tension.

\[\mathbf{u}^\infty = (U, V, W)\]

Parameters:

• uinf – an ndarray of size ambient_dim or a single number, in which case the velocity field is set to be \((U_\infty, 0, 0)\).

• pinf – constant farfield pressure.

Returns:

a FarfieldCapillaryBoundaryCondition.

pystopt.stokes.get_bezier_uniform_farfield(ambient_dim: int, points_or_order: int | Tuple[ndarray, ...], *, tmax: float | pytential.symbolic.primitives.var, pinf: float = 0.0, tau: str | pytential.symbolic.primitives.var = 'tau', capillary_number_name: str = 'ca') → FarfieldCapillaryBoundaryCondition

Time-dependent uniform farfield flow with constant surface tension.

The control points in ps define a curve

\[\mathbf{p}(\tau) = \sum \mathbf{p}_i \phi_i(\tau),\]

where \(\phi_i(\tau)\), for the normalized time variable \(\tau \in [0, 1]\), are the Bézier curve basis functions for the respective order. The velocity field is then given by

\[\mathbf{u}^\infty(\tau) = \frac{\mathrm{d} \mathbf{p}}{\mathrm{d} \tau}.\]

Parameters:

• points_or_order – a set of control points used to form the Bézier curve. The number of points also defines the order of the curve. If this is a single integer, the points are defined symbolically and named "p{i}".

• pinf – constant farfield pressure.

• tau – variable name for the parametrization of the Bézier curve. This is expected to be normalized to \([0, 1]\) by the user when the expression is eventually evaluated.

Returns:

a FarfieldCapillaryBoundaryCondition.

pystopt.stokes.get_shear_farfield(ambient_dim: int, *, alpha: float = 1.0, pinf: float = 0.0, capillary_number_name: str = 'ca') → FarfieldCapillaryBoundaryCondition

Shear farfield flow with constant surface tension.

\[\mathbf{u}^\infty = (\alpha y, 0, 0)\]

Parameters:

alpha – shear strength.

pystopt.stokes.get_solid_body_rotation_farfield(ambient_dim: int, *, omega: float | ndarray = 1.0, pinf: float = 0.0, capillary_number_name: str = 'ca') → FarfieldCapillaryBoundaryCondition

Solid body rotation flow with constant surface tension.

\[\mathbf{u}^\infty = \mathbf{\omega} \times \mathbf{x},\]

where \(\mathbf{omega}\) is the axis of the solid body rotation. If a scalar is given, the axis is set to \((0, 0, \omega)\). For the 2D case, we omega can only be a scalar and the rotation is always in the \(x-y\) plane.

Parameters:

omega – magnitude or axis of rotation.

pystopt.stokes.get_helical_farfield(ambient_dim: int, *, tmax: float | pytential.symbolic.primitives.var, height: float | pytential.symbolic.primitives.var = 1.0, omega: float | pytential.symbolic.primitives.var = 6.283185307179586, pinf: float | pytential.symbolic.primitives.var = 0.0, capillary_number_name: str = 'ca') → FarfieldCapillaryBoundaryCondition

Analytical Solutions

All subclasses of AnalyticStokesSolution implement pystopt.stokes.StokesInterface with the arguments

@svelocity.register(AnalyticStokesSolution)
def _velocity_solution(self: AnalyticStokesSolution, side: int):
    pass

and similar for the remamining variables. In the case of the jumps, the side can be omitted for obvious reasons.

class pystopt.stokes.SympySolution(pressure: Expr, velocity: ndarray, traction: ndarray, normal_velocity_gradient: ndarray, tangent_traction: ndarray, tangent_velocity_gradient: ndarray, normal_velocity_laplacian: Expr, normal_velocity_hessian: Expr, streamfunction: Expr = None)

__init__(pressure: Expr, velocity: ndarray, traction: ndarray, normal_velocity_gradient: ndarray, tangent_traction: ndarray, tangent_velocity_gradient: ndarray, normal_velocity_laplacian: Expr, normal_velocity_hessian: Expr, streamfunction: Expr = None) → None

class pystopt.stokes.AnalyticStokesSolution(ambient_dim: int, *, viscosity_ratio_name: str = 'viscosity_ratio', capillary_number_name: str = 'ca')

Bases: StokesInterface

Analytical solution to a two-phase Stokes flow.

viscosity_ratio_name: str

A name for the viscosity ratio symbolic variable.

capillary_number_name: str

A name for the capillary number symbolic variable.

property farfield: TwoPhaseStokesBoundaryCondition

The farfield boundary conditions for this solution.

property as_sympy: PymbolicToSympyMapper

A mapper used to translate an expression to sympy.

__getitem__(side: int) → SympySolution

Parameters:

side – interior (-1) or exterior (+1) solution.

Returns:

a SympySolution with the solution on the respective side.

class pystopt.stokes.SphericalStokesSolution(*, radius=1.0, viscosity_ratio_name='viscosity_ratio', capillary_number_name='ca')

Bases: AnalyticStokesSolution

Two-phase Stokes solutions on a sphere.

property normal: ndarray

Exterior normal vector on the sphere:

\[\mathbf{n} = \frac{\mathbf{x}}{R}.\]

property summed_curvature: float

Summed curvature of the sphere (sum of principal curvatures)

\[\kappa = \frac{d - 1}{R}.\]

class pystopt.stokes.PolarStokesSolution(*, radius: float = 1.0, uinf: float = 1.0, pinf: float = 0.0, alpha: float | None = None, viscosity_ratio_name: str = 'viscosity_ratio', capillary_number_name: str = 'ca')

Bases: AnalyticStokesSolution

Two-phase Stokes solutions on a circle.

radius: float

Radius of the circle representing the droplet interface.

uinf: float

Farfield velocity field.

pinf: float

Farfield pressure value.

alpha: float

Additional parameter used to determine a unique solution. A special value of this parameter leads to a steady state solution.

property normal: ndarray

Exterior normal vector on the circle.

property tangent: ndarray

Tangent vector on the circle (counterclockwise).

property curvature: float

Curvature of the circle (constant based on the radius).

class pystopt.stokes.HadamardRybczynski(*, radius=1.0, alpha=None, uinf=1.0, pinf=0.0, capillary_number_name='ca', viscosity_ratio_name='viscosity_ratio')

Bases: SphericalStokesSolution

Hadamard-Rybczynski solution.

The Hadamard-Rybczynski solution is an axisymmetric solution in a uniform farfield flow. It allows stream function solutions of the form

\[\psi_\pm = \left( \frac{A_\pm}{r} + B_\pm r + C_\pm r^2 + D_\pm r^4 \right) \sin^2 \theta,\]

where the constants are determined (as a one-parameter family) by the boundary conditions

• mass conservation at the interface \([\![ \mathbf{u} ]\!] = 0\).

• jump in normal traction at the interface \([\![ \mathbf{f} \cdot \mathbf{n} ]\!] = \kappa / \mathrm{Ca}\).

• jump in tangential traction at the interface \([\![ \mathbf{f} \cdot \mathbf{t} ]\!] = 0\).

• finite solutions at the origin and at infinity.

For details, see Equation 3-7 in [Clift1978] and the surrounding discussion.

__init__(*, radius=1.0, alpha=None, uinf=1.0, pinf=0.0, capillary_number_name='ca', viscosity_ratio_name='viscosity_ratio')

Parameters:

• alpha – parameter for the solution. If None, the classic solution from [Clift1978] is computed.

• uinf – magnitude of the \(x\) component of farfield velocity field.

• pinf – farfield pressure field.

class pystopt.stokes.Taylor(*, radius=1.0, alpha=1.0, pinf=0.0, capillary_number_name='ca', viscosity_ratio_name='viscosity_ratio')

Bases: SphericalStokesSolution

Taylor solution.

This is a fully three-dimensional solution in an extensional flow. For details see [Taylor1932].

__init__(*, radius=1.0, alpha=1.0, pinf=0.0, capillary_number_name='ca', viscosity_ratio_name='viscosity_ratio')

Parameters:

• alpha – shear strength.

• pinf – farfield pressure field.

Adjoint Boundary Conditions

pystopt.stokes.make_static_adjoint_boundary_conditions(ambient_dim: int, cost: ShapeFunctional) → AdjointBoundaryCondition

Construct a AdjointBoundaryCondition based on the static jump conditions AdjointTractionJumpCondition.

pystopt.stokes.make_unsteady_adjoint_boundary_conditions(ambient_dim: int, cost: ShapeFunctional, *, xstar_name: str = 'xstar') → AdjointBoundaryCondition

Construct a AdjointBoundaryCondition based on the quasi-static jump conditions UnsteadyAdjointTractionJumpCondition.

class pystopt.stokes.AdjointTractionJumpCondition(ambient_dim: int, cost: ShapeFunctional)

Bases: TractionJumpCondition

Defines adjoint boundary conditions for the static Stokes system.

In the static case, the jump conditions are entirely defined by the cost function. In particular, we have that

\[[\![ \mathbf{f}^* ]\!] = \nabla_{\mathbf{u}} J.\]

cost

class pystopt.stokes.UnsteadyAdjointTractionJumpCondition(ambient_dim: int, cost: ShapeFunctional, xstar_name: str = 'xstar')

Bases: AdjointTractionJumpCondition

Defines adjoint boundary conditions for the quasi-static Stokes system.

In the quasi-static case, the jump conditions are defined by the cost function and the adjoint variable constraining the forward kinematic condition. We have that

\[[\![ \mathbf{f}^* ]\!] = x^* \mathbf{n} + \nabla_{\mathbf{u}} J,\]

under the assumptions that the forward kinematic condition is given by

\[\dot{\mathbf{x}} = (\mathbf{u} \cdot \mathbf{n}) \mathbf{n}.\]

xstar_name

xstar

class pystopt.stokes.AdjointBoundaryCondition(ambient_dim: int, deltafs: Tuple[TractionJumpCondition, ...])

Bases: TwoPhaseStokesBoundaryCondition

Adjoint boundary conditions for a static Stokes system.

Implements the interface of StokesInterface with the signature:

sti.velocity(adjoint_bc)

ambient_dim

cost

__init__(ambient_dim: int, deltafs: Tuple[TractionJumpCondition, ...]) → None