Problems with the two-phase flow problem calculated on a hyperball (Step.21)

[simon doersam]

Hello,

In Step.21 (The two phase flow problem) I replaced the domain "Hyper_cube" with "Hyper_ball". 
The numerical solution is unrealistic and appears to be numerical problems at all points,
where the boundary intersects with one of the following planes:
E1: = {x \ in R ^ 3, x1 = x2},
E2: = {x \ in R ^ 3, x1 = x3},
E3: = {x \ in R ^ 3, x2 = x3}.

The pressure and the viscosity are partially negative at these points.

Can someone help me and explain why this is happening?

I am happy for any hint, literature or explanation.

Best wishes,
Simon

[Wolfgang Bangerth]

On 8/26/20 7:05 AM, 'simon doersam' via deal.II User Group wrote:
> In Step.21 (The two phase flow problem) I replaced the domain "Hyper_cube"
> with "Hyper_ball".
> The numerical solution is unrealistic and appears to be numerical problems at
> all points,
> where the boundary intersects with one of the following planes: E1: = {x \ in
> R ^ 3, x1 = x2}, E2: = {x \ in R ^ 3, x1 = x3}, E3: = {x \ in R ^ 3, x2 = x3}.
>
> The pressure and the viscosity are partially negative at these points.
>
> Can someone help me and explain why this is happening?

You also need to change the definition of boundary values:

template <int dim>
class SaturationBoundaryValues : public Function<dim>
{
public:
SaturationBoundaryValues()
: Function<dim>(1)
{}

virtual double value(const Point<dim> &p,
const unsigned int /*component*/ = 0) const override
{
if (p[0] == 0)
return 1;
else
return 0;
}
};

This simply doesn't make any sense with the geometry you are now using.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[spass...@googlemail.com]

Thank you for your respone. Sorry, I did not think about the meaning of the modification in step.21. It was naive. But in my work, I solve only the Darcy equation (v velocity, p pressure) on a more complicated mesh, which is similar to a Hyper_Ball. I have the same problems. So I tried to find a “simple case” with the same behavior to understand the problem.

For example I solve the Darcy equation on a "HyperCube" with constant Dirichlet condition for the pressure (p=p_D on partial Omega). Then I get the trivial solution v = 0, p = p_D (=2000).
But if I use a “Hyper_Shell” or “Hyper_Ball” instead of a “Hyper_Cube”, then there are problems(see attachments) and I don't get the trivial solution. Why?

From the data, I have only an average pressure. So I choose the pressure on a part of the boundary constant. Certainly this is not entirely true, but it doesn't seem that this is my problem. I suppose the grid related to the Raviart-Thomas elements is the problem. But I will also think again about my boundary conditions.

Best, Simon

[Wolfgang Bangerth]

Simon,

> visit0000.pngThank you for your respone. Sorry, I did not think about the

> meaning of the modification in step.21. It was naive. But in my work, I solve
> only the Darcy equation (v velocity, p pressure) on a more complicated mesh,

> whichis similar to a Hyper_Ball. I have the same problems. So I tried to find
> a “simple case” with the same behavior to understandthe problem.

>
> For example I solve the Darcy equation on a "HyperCube" with constant
> Dirichlet condition for the pressure (p=p_D on partial Omega). Then I get the
> trivial solution v = 0, p = p_D (=2000).
> But if I use a “Hyper_Shell” or “Hyper_Ball” instead of a “Hyper_Cube”, then
> there are problems(see attachments) and I don't get the trivial solution. Why?

The question is: How exactly do you prescribe the boundary values? There are
several steps associated with getting boundary values right:
* You need to generate a geometry
* If you have curved boundaries, you also need to make sure that the different
parts of the boundary have the appropriate manifold objects associated with
them to ensure that the mesh is refined in the correct way
* You need to associate different boundary indicators to different parts of
the boundary
* You need to write a function that for a given points returns the appropriate
boundary value.

I don't know how you implement all of these steps, but you should double check.