Rotating fields and expressions in fenics

1

Suppose I have an expression like the one of the classic Taylor Green 2D solution of the Navier-Stokes equations:

u0='-sin(pi*x[1])*cos(pi*x[0])*exp(-2.*pi*pi*nu*t)',
u1='sin(pi*x[0])*cos(pi*x[1])*exp(-2.*pi*pi*nu*t)',
p='-(cos(2*pi*x[0])+cos(2*pi*x[1]))*exp(-4.*pi*pi*nu*t)/4.0'

I would like to apply such a field as boundary condition on the boundary of a 2D circle, for example, but rotating it in time. In other words, I want to compose the field above with a rotation (with axis z of course). What is the most efficient and precise way to do this? I tried to define a custom UserExpression (suppose RM_ is the array of a certain rotation matrix 3by3), in order to rotate the coordinates and only then apply the starting Expression to the new coordinates:

RM_ = np.array([[0.70711,  0.70711,  0.0],
                [-0.70711, 0.70711,  0.0],
                [0.0,      0.0,      1.0]])

class RotExpression(UserExpression):
       def eval(self, value, x):
           xa_ = np.array([x[0], x[1], 0.])
           xr_ = RM_ @ xa_     #ROTATE THE COORDINATES, MATRIX VECTOR PRODUCT
           value[0] = xr_[0]
           value[1] = xr_[1]
       def value_shape(self):
           return (2,)

cr_ = RotExpression(degree = 3)
u_ex_rot = Expression(('-(cos(pi*(cr[0]))*sin(pi*(cr[1])))*exp(-2.0*nu*pi*pi*t)',
               ' (cos(pi*(cr[1]))*sin(pi*(cr[0])))*exp(-2.0*nu*pi*pi*t)'),
               cr=cr_, nu=nu, t=t, degree=3)

The domain here is a 2D disc with centre at the origin, so I rotate every point on the boundary by multiplying with a rotation matrix.

And then use this last u_ex_rot expression which I obtain in an usual DirichletBC on the boundary domain. The results I obtain however seem wrong and I would like to rule out that there’s something wrong that I’m doing with this implementation of the boundary conditions.

2

Before you blame the rotational representation, I would solve the problem on a rectangular domain,where you can have no-slip (Zero-conditions) on all boundaries.
Then you can check if you obtain the expected convergence rates.

To obtain a conservative rotation field, you should solve an additional PDE, see this demo where a conservative rotational field is computed with a Crank-Nicolson discretization.

1 Like
3

Thanks @dokken, I am reading the tutorial and it seems quite similar to what I want to do in the first part. I think it makes sense to start first by rotating the domain instead of rotating the field of the boundary conditions. However I don’t understand something:

  • Why do you suggest solving with a rectangular domain? I solved the Navier Stokes equations on both a (fixed) rectangle and a (fixed) circle in 2D, applying the Taylor Green solution as boundary conditions on all of the exterior boundary, and I get a solution which makes sense, and error norms that decade in time (H1 and L2 velocity, L2 pressure). Would there be any difference in taking a rotating rectangle instead of a rotating circle?
  • The first step I would like to do is to solve the NS equations on a domain with rotating points and recover the Taylor Green solution. That is, I want a mesh with an arbitrary movement on a background field and I still recover the solution, which should be the idea of the Arbitrary Lagrangian Eulerian framework. The movement of the mesh should not change the solution compared to the fixed case if the discretisation of the ALE terms is done correctly.
  • As a last step, I would like to achieve another thing: since the domain is rotating rigidly, instead of actually moving the mesh, I could rotate the axes of coordinates. To do this, I think the only thing missing is to change the boundary conditions accordingly, i.e. to rotate the Taylor Green known field, which I have as expression, by composing it with a rotation, which is what I was trying to do above. Do you think it’s feasible?

I hope my question is clear. I will try first to move the mesh following the link you posted.

4

The reason i suggested solving on a fixed domain was to make sure that you solved the equations correctly. As you have done that, you can move on to using a conservative velocity formulation.

5

Dear @dokken, I have added to my previous code the mesh movement procedure from the tutorial. Before solving the NS equations now I move the mesh at each time step using the rotational field from the solution of the PDE discretized with Crank-Nicolson.

The NS solver is the same which worked for the Taylor Green in the fixed case, only with the addition to of the ALE domain velocity in the convective term.
However I see some strange behaviour with the solution quickly degrading, with the velocity losing regularity from the starting Taylor Green solution:

velocity
velocity1467×839 104 KB

And the pressure is even worse, understandably:

pressure
pressure1467×839 104 KB

I see that in this case, the integral of the absolute value of the divergence of velocity over the domain steadily increases at each time step, while the error of the velocity in H1 norm quickly explodes.
With the same procedure I also see that the divergence of the displacement field I calculate, with which I then move the mesh, is not exactly zero. Could this be related to the fields being non divergence free?

6

As this is not a problem I have encountered previously, I am not 100 % certain of what’s wrong in this case.

However, some thoughts:

  • what happens if you reduce the magnitude of the velocity ?
  • how are you solving the Navier Stokes equation? Are you using a splitting scheme or a coupled solver? If you are using a splitting scheme, how do you treat the convective term?
7

@dokken It’s a strange result and I’m not sure what causes it, however:

  • decreasing the value of \omega, even of several orders of magnitude, still doesn’t solve the problem: the solution looks right, at least for the first time steps, and the errors for velocity, H1 and L2, increase less in magnitude each time step, but they are still monotonically increasing with time. The L2 error on the pressure doesn’t seem to be affected that much by the reduction of the rotation rate.
  • I am using Oasis to solve Navier Stokes, therefore a splitting pressure correction scheme. In particular the scheme is the classic incremental pressure correction scheme (IPCS), with an Adam Bashforth discretization of the convective term: \left( \left( \frac{3}{2}\mathbf{u}^0 - \frac{1}{2}\mathbf{u}^{00} - \mathbf{\Omega}\times \mathbf{r}\right) \cdot \nabla \right) \mathbf{u},
    where I added by hand the last term, that is the convective ALE contribution of the mesh velocity, that in this case is a rotation. Calculating this as a finite difference of the displacement divided by \Delta t doesn’t seem to make any difference.
8

At least, you should discretize the mesh velocity in the same way as you have discretize u with A-B semi-implicitly.