I consider the following situation.
I have a domain comprising of a free flow part Z^+=(0,1)x(0,4) and a porous-medium part Z^-= (0,1)x(0,4)\P, where solid inclusions are present and P denotes the set of solid inclusions.
Both subdomains are divided by the interface S=(0,1)x{0}.
The overall flow domain is denoted by Z=Z^+ \cup S \cup Z^-.
Now I consider two problems defined on the domain Z.
First a Stokes problem with a jump condition across the interface.
The solution to this problem is {u,p}.
What I know: u stabilizes exponentially fast to a constant C for x_2 → 4.
The second problem is then defined as follows:
problem eta([ue], [uhe], init=reuseMatrix, solver=sparsesolver) =
int2d(Th)(
dx(ue)*dx(uhe)
+ dy(ue)dy(uhe)
)
+ int2d(Th)((u - C(y>0))uhe)
+ int1d(Th,interface)( -int2d(Th)((u)(y<0)) * uhe )
/// we set the following constraint to define eta uniquely
/// but the exact constant for us constant does not play a role
+ int1d(Th, bottom)(ue * uhe)
;
Since for later calculations I need that the third derivative of ue exists I need to use P3 FEM at least.
However, when I use P3 FEM I observe that when I compute dx(ue) on the left and right boundary there appear some sort of numerical singularities (see picture attached).
But indeed, for the second problem itself it is sufficient to use P2 FEM.
So I did this and using P2 the singularities for dx(ue) do not appear.
Has somebody an idea where these problems come from?
Can you post the full code as an attachment? I guess if you just want to
get around the 3rd derivative issue, you can introduce another
variable as in mixed formulation, y2=dx(y1) and this may have some
benefits or not.
I would guess that has something to do with the periodic boundary conditions.
I have not used those enough to comment but I guess there are some
issues there. I can’t run your code with so many nodes in the adaptmesh
but may play with it later. In any case the FF experts may have more interest if
you mention the periodic thing. I take it your porous medium them is
periodic and not random. Do you expect physics due to this ? Resonances under
some flow conditions? I guess analytically it would be suitable for FFT analysis
probably. I guess it is possible under some conditions you get singularities.
Also 'sum rules" may be important as I found out playing with a coloumb example
with a net charge per cell In any case, if you change the forcing functions curious if
the peaks move or go away.
I’m not sure of the issues with adaptmesh and periodic boundary condtions
but
I guess given the physics I would try to include 2 periods in one cell,
I think this is called a supercell, and see if you get a singularity inside the
compsite cell. If you suddenly reduce the flow rate that fluid has to go
somewhere and the singularities look like that could be where What
happens with the flow in the x direction which I guess is v?
You have u and v constrained to zero at the bottom.
I had to reduce the vertex count in adaptmesh and could not
reproduce the peaks but it sounds like they make some sense.
I guess you could play with a position dependent height/depth or
density and viscosity. If the equations don’t have a meaningful
steady state result the output may vary with details of the model.
I’m not sure if there are any "“long range forces”’ that allow distant
cells to interact as you see in electromagnetics but maybe you don’t
need them here
Thanks for your answer.
I run the code now without using adaptmesh but I had the same problems with the singularities on the left and right boundary (they just appeared in different positions).
I did not get the idea with the supercell you meant. Could you please explain again?
For the first problem (called stokes in my code), the problem is well-posed and the constraint setting u=v=0 on the bottom is also physically consistent. For this problem also the solution looks fine.
Singularities appear only when considering the gradient of ue.
But as stated before, this is only the case when I use P3 elements (that I need for later purposes) to solve the problem eta. When I use P2 elements the gradient of ue looks fine.
If the derivative is positive on one side and negative on the other I guess
you get a cusp at the boundary. I don’t remember the details but depending on
the boundary conditions it looked like fluid could be “reflect” from the obstacles
and the location of singularity looked like a good place for it to pile up.
IIRC, you had two obstacles one above the other with periodic boundary conditions
in the X direction. You could create 4 obstacles in a a square pattern and see if you
get the singularity on the middle of the mesh as well as the periodic boundaries.
This may also let you perturb their locations slightly if you are worried about
stuff like that.
In any case, if you change the forcing functions curious if