Small step-22 modification gives unexpected result

[Patrick Eastham]

Hello all,

I'm very new to deal.II and I've made some small modifications to the step-22 tutorial in an effort to validate (for myself) how to simulate a simple channel-flow problem. I've run into an issue that might be a bug, but more likely I believe I'm missing something (hopefully obvious) about how boundary constraints are handled. Below I detail what I have changed from step-22, as well as include the code and a screenshot of the output which shows the issue. For a one-sentence summary: I'm trying to do channel flow (left-to-right) while prescribing velocity boundary conditions on the inflow (left) and walls (top,bottom), and homogeneous neumann on the outflow (right); the result is generally what is expected except for the solution near the outflow boundary.

From step-22, I have made two modifications: (1) boundary id's changed so that top, left, and bottom have boundaryid=1 and right boundary has boundaryid=0, and (2) implement exact solution function for simple 2D channel flow, to be implement on boundaryid=1. This should give flow from left to right. Nothing else about the program is changed.

As can be seen from the screenshot, there is error near the outflow where the flow has non-zero vertical component (it should be zero) and the pressure is not constant on the outflow boundary. The velocity values on the boundaries do satisfy the boundary conditions. I've tried (a) using a direct solver and (b) implementing boundary conditions after assembly and both of those changes (not included in attached code) gave the same solution.

I've tried looking through tutorials, as well as other similar questions relate to boundary conditions on this forum, but nothing has worked so far. Any help would be greatly appreciated!

Best,

Patrick

[Praveen C]

Some things I noticed in your code.

Your velocity seems to be from right to left since you put a minus sign

    if (component == 0)

      return -umax*p[1]*(1.0 + p[1]);

but you say you put neumann bc on right boundary.

You are using equality test to identify boundary faces

        if (cell->face(f)->center()[dim - 1] == 0)    // top boundary

          cell->face(f)->set_all_boundary_ids(1);

This may or may not work correctly, it is better to put some tolerance.

Best
praveen

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/f459c6dd-2b01-4163-83f1-fff30e3db3f0%40googlegroups.com.
<step-22_channel.cc><solution-00.vtk><solution-00_screenshot.png>

[Praveen C]

Sorry, I did not see your y is from [-1,0] so your velocity bc is correct.

Since the flow is driven by pressure gradient, it may be better to impose the pressure gradient as a body force term.

Ideally you would want to use periodic bc along x.

Best

praveen

[Patrick Eastham]

Praveen,

Thanks for the suggestion on adding tolerance to assigning boundary conditions -- that could come back to bite me in the future.

I also appreciate your other suggestions, although I'd like to get this working without those conditions. I know that the FEM problem of assigning dirichlet velocity BC on 3 inflow and walls, and neumann conditions on the outflow is well-posed and leads to the exact solution (as the exact velocity is quadratic and exact pressure is linear, just like my element family) as I've done this same problem in a different FEM code.

It seems like the cell at the intersection of the dirichlet and neumann conditions is not being handled correctly based on how I have assigned constraints/boundary conditions. Has anyone else had any experience with this?

Best,

Patrick

To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.

[Praveen C]

The step-22 code uses the stress form of the NS equations. If you use the form

-mu*Laplace(u) + grad(p) = 0

and use the same dirichlet+neumann bc, then it seems to give the correct solution. Please check it.

These are the lines to be changed

    std::vector<Tensor<2, dim>> symgrad_phi_u(dofs_per_cell);

                symgrad_phi_u[k] =

                  fe_values[velocities].gradient(k, q);

local_matrix(i, j) +=

                      (scalar_product(symgrad_phi_u[i], symgrad_phi_u[j]) // (1)

                       - div_phi_u[i] * phi_p[j]                 // (2)

                       - phi_p[i] * div_phi_u[j])                // (3)

                      * fe_values.JxW(q);                        // * dx

Best

praveen

[Patrick Eastham]

You're absolutely right! That solves the problem completely. I greatly appreciate your help Praveen!