Skip to content

Momentum equation

momentum_equation

Derived terms and associated equations for the Stokes system.

All terms are considered as if they were on the right-hand side of the equation, leading to the following UFL expression returned by the residual method:

This sign convention ensures compatibility with Thetis's time integrators. In general, however, we like to think about the terms as they are on the left-hand side. Therefore, in the residual methods below, we first sum the terms in the variable F as if they were on the left-hand side, i.e.

and then return -F.

viscosity_term(eq, trial)

Viscosity term in the momentum equation.

Using the symmetric interior penalty method (Epshteyn & Rivière, 2007), the weak form becomes

where σ is a penalty parameter.

Epshteyn, Y., & Rivière, B. (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2), 843-872.

Source code in g-adopt/gadopt/momentum_equation.py 
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def viscosity_term(
    eq: Equation, trial: Argument | ufl.indexed.Indexed | Function
) -> Form:
    r"""Viscosity term $-nabla * (mu nabla u)$ in the momentum equation.

    Using the symmetric interior penalty method (Epshteyn & Rivière, 2007), the weak
    form becomes

    $$
    {:( -int_Omega nabla * (mu grad u) phi dx , = , int_Omega mu (grad phi) * (grad u) dx ),
      ( , - , int_(cc"I" uu cc"I"_v) "jump"(phi bb n) * "avg"(mu grad u) dS
          -   int_(cc"I" uu cc"I"_v) "jump"(u bb n) * "avg"(mu grad phi) dS ),
      ( , + , int_(cc"I" uu cc"I"_v) sigma "avg"(mu) "jump"(u bb n) * "jump"(phi bb n) dS )
    :}
    $$

    where σ is a penalty parameter.

    Epshteyn, Y., & Rivière, B. (2007).
    Estimation of penalty parameters for symmetric interior penalty Galerkin methods.
    Journal of Computational and Applied Mathematics, 206(2), 843-872.
    """
    dim = eq.mesh.geometric_dimension()
    identity = Identity(dim)
    compressible_stress = eq.approximation.compressible

    mu = eq.approximation.mu
    stress = eq.stress
    F = inner(nabla_grad(eq.test), stress) * eq.dx

    sigma = interior_penalty_factor(eq)
    sigma *= FacetArea(eq.mesh) / avg(CellVolume(eq.mesh))
    if not is_continuous(eq.trial_space):
        trial_tensor_jump = tensor_jump(eq.n, trial) + tensor_jump(trial, eq.n)
        if compressible_stress:
            trial_tensor_jump -= 2 / 3 * identity * jump(trial, eq.n)

        F += (
            sigma
            * inner(tensor_jump(eq.n, eq.test), avg(mu) * trial_tensor_jump)
            * eq.dS
        )
        F -= inner(avg(mu * nabla_grad(eq.test)), trial_tensor_jump) * eq.dS
        F -= inner(tensor_jump(eq.n, eq.test), avg(stress)) * eq.dS

    # NOTE: Unspecified boundaries result in free stress (i.e. free in all directions).
    # NOTE: "un" can be combined with "stress" provided the stress component is
    # tangential (e.g. no normal flow with wind)
    for bc_id, bc in eq.bcs.items():
        if "u" in bc and any(bc_type in bc for bc_type in ["stress", "un"]):
            raise ValueError(
                '"stress" or "un" cannot be specified if "u" is already given.'
            )
        if "normal_stress" in bc and any(bc_type in bc for bc_type in ["u", "un"]):
            raise ValueError(
                '"u" or "un" cannot be specified if "normal_stress" is already given.'
            )

        if "u" in bc:
            trial_tensor_jump = outer(eq.n, trial - bc["u"])
            if compressible_stress:
                trial_tensor_jump -= (
                    2 / 3 * identity * (dot(eq.n, trial) - dot(eq.n, bc["u"]))
                )
            trial_tensor_jump += transpose(trial_tensor_jump)
            # Terms below are similar to the above terms for the DG dS integrals.
            F += (
                2
                * sigma
                * inner(outer(eq.n, eq.test), mu * trial_tensor_jump)
                * eq.ds(bc_id)
            )
            F -= inner(mu * nabla_grad(eq.test), trial_tensor_jump) * eq.ds(bc_id)
            F -= inner(outer(eq.n, eq.test), stress) * eq.ds(bc_id)

        if "un" in bc:
            trial_tensor_jump = outer(eq.n, eq.n) * (dot(eq.n, trial) - bc["un"])
            if compressible_stress:
                trial_tensor_jump -= 2 / 3 * identity * (dot(eq.n, trial) - bc["un"])
            trial_tensor_jump += transpose(trial_tensor_jump)
            # Terms below are similar to the above terms for the DG dS integrals.
            F += (
                2
                * sigma
                * inner(outer(eq.n, eq.test), mu * trial_tensor_jump)
                * eq.ds(bc_id)
            )
            F -= inner(mu * nabla_grad(eq.test), trial_tensor_jump) * eq.ds(bc_id)
            # We only keep the normal part of stress; the tangential part is assumed to
            # be zero stress (i.e. free slip) or prescribed via "stress".
            F -= dot(eq.n, eq.test) * dot(eq.n, dot(stress, eq.n)) * eq.ds(bc_id)

        if "stress" in bc:  # a momentum flux, a.k.a. "force"
            # Here we need only the third term because we assume jump_u = 0
            # (u_ext = trial) and stress = n . (mu . stress_tensor).
            F -= dot(eq.test, bc["stress"]) * eq.ds(bc_id)

        if "normal_stress" in bc:
            F += dot(eq.test, bc["normal_stress"] * eq.n) * eq.ds(bc_id)

    return -F