Help with a weak form

1

Hi,

I am trying to solve the following problem and wanted to check and see if my formulation is correct.

-\mu \nabla^2 {\bf v}+\nabla p- (\mu/3) \nabla div({\bf v}) = 0
div({\bf v}) = f
Q\hat{n} = - \kappa \hat{n}

Where \hat{n} is the unit outward normal vector and Q is the stress tensor given by Q = \mu (\nabla {\bf v} + \nabla {\bf v}^T) -(p + (2 \mu/3) div({\bf v})) (so that -\nabla Q = 0 gives you the first equation above).
Now I take a vector test function \xi and a scalar test function q and write the following weak form:
\int_{\Omega} \mu \nabla{\bf v} \cdot \nabla \xi \ dx - \int_{\Omega} p \ div(\xi) \ dx + (\mu/3) \int_{\Omega} div({\bf v}) \ div(\xi) \ dx + \int_{\Omega} div({\bf v}) \ q \ dx - \int_{\Omega} f \ q \ dx +\\ \int_{\partial \Omega} \kappa \hat{n} \cdot \xi \ dS=0

Is this weak form correct or am I missing something? Can this be implemented exactly like this in FEniCS with a mixed pressure-velocity space?

Thank you in advance.

2

If you want -\boldsymbol{\kappa} to have the interpretation of flux \mathbf{Q}, you want to directly integrate the flux divergence form of the problem by parts:

\int_\Omega\mathbf{Q}:\nabla\boldsymbol{\xi}\,d\mathbf{x} + \int_{\partial\Omega}\boldsymbol{\kappa}\hat{\mathbf{n}}\cdot\boldsymbol{\xi}\,d\mathbf{s} + \int_\Omega(\nabla\cdot\mathbf{v} - f)q\,d\mathbf{x}= 0\quad\forall\boldsymbol{\xi},q\text{ ,}

i.e., skip the simplification

\nabla\cdot(2\nabla^s\mathbf{v} - (2/3)(\nabla\cdot\mathbf{v})\mathbf{I}) = \nabla^2\mathbf{v} + (1/3)\nabla(\nabla\cdot\mathbf{v})

in the strong form before deriving the weak form. With the weak form given in the post, you would enforce

-\boldsymbol{\kappa}\hat{\mathbf{n}} = \left(\mu\nabla\mathbf{v} - p\mathbf{I} + (\mu/3)(\nabla\cdot\mathbf{v})\mathbf{I}\right)\hat{\mathbf{n}}~~~\neq~~~\mathbf{Q}\hat{\mathbf{n}}\text{ .}
1 Like
3

I see. Thank you. Can this be directly implemented in Fenics using the following weak form?
mu*inner(grad(v)+grad(v).T, grad(xi))*dx +p*div(xi)*dx+(2*mu/3)*div(v)*div(xi)*dx +dot(kappa*NORMAL,xi)*ds(1)+ div(v)*q*dx-f*q*dx
I know that I need a mixed function space for pressure and velocity and their respective test functions. But is there more to it?

4

Comparing with the formula for \mathbf{Q}, there are sign errors in the 2nd and 3rd terms, but, otherwise, yes, that looks like how you would translate it into UFL. You could also use an intermediate definition of Q like

I = Identity(len(v))
Q = 2*mu*sym(grad(v)) - (p + (2*mu/3)*div(v))*I
residual = inner(Q,grad(xi))*dx + dot(kappa*NORMAL,xi)*ds(1) + (div(v)-f)*q*dx
a = lhs(residual)
L = rhs(residual)

to avoid typos. The definitions of the bilinear form a and linear form L are assuming that v,p = TrialFunctions(V) for a mixed space V, so you can use the solve syntax for a linear problem, e.g.,

w = Function(V)
solve(a==L,w,bcs)
1 Like
5

Thank you so much. It was really helpful.