Problem imposing tangent stress condition


I am solving a compressible Stokes problem with a prescribed pressure on a hexagon with the following BCs on each face i:

  1. =0
  2. ti.( = gi

But in my scheme the 2nd condition is not imposed to a satisfactory degree.

As explained here, I am using penalization for the first condition and projection on the tangent to impose the 2nd condition.

problem stokes([ux,uy],[vx,vy], solver=sparsesolver) = 
        - int2d(Th)( grad(pr)'*[vx, vy] )                         // pressure
        - int2d(Th)( eta*Strain(ux, uy)'*Grad(vx, vy) )           // viscosity
        - int2d(Th)( xi*[ux, uy]'*[vx, vy] )                      // friction 
        - int1d(Th,1)( g1 * [T1[0], T1[1]]'*[vx, vy] )            // ti.( = gi
        - int1d(Th,2)( g2 * [T2[0], T2[1]]'*[vx, vy] )    
        - int1d(Th,3)( g3 * [T3[0], T3[1]]'*[vx, vy] )   
        - int1d(Th,4)( g4 * [T4[0], T4[1]]'*[vx, vy] )    
        - int1d(Th,5)( g5 * [T5[0], T5[1]]'*[vx, vy] )     
        - int1d(Th,6)( g6 * [T6[0], T6[1]]'*[vx, vy] )             
        - int1d(Th)( penal * (N'*[ux, uy]) * (N'*[vx, vy]) );     // penalize u.n

Note: each gi is defined as P1 on a line with tangent Ti (precomputed). The line meshes coincide with the Hexagon boundaries. The issue persists also if g are defined in 2D, and regardless of interpolation order.

The discrepancy between the actual ti.( (after solving for u) and the input gi is what I am struggling to fix. The problem is mostly concentrated near the edges of each line but not only.

I appreciate the help. Thanks.

ffbctest.edp (5.3 KB)

Your code looks correct. I checked that the error in the tangential boundary condition decreases when refining the mesh, linearly in the mesh size. I think that the corners of the domain induce some singularities in the solution, so that the error is not as small as you expect. Nevertheless I think that what you get is correct.

Thank you François,
Could it be that imposing such BCs through sharp corners is somehow ill-posed?
Note that a clearer way to see the problem is setting all g’s to 1.

I think that the problem is well-posed because the conditions are set a.e. on the boundary (at the continuous level).
But the solution u is probably in the Sobolev space H^2 and no more than 2 derivatives.
Looking at the numerical error on the BCs, if you take a test function v which is arbitrary on the boundary satifying v.N=0, and vanishing for the degreees of freedom inside the domain, you see that the boundary error E on the tangent stress verifies \int_\Gamma E.v=\int_\Omega div(\sigma).v
Assume that div(\sigma) is in L^\infty. Then taking v=E on the boundary you get
that v is in L^4 in \Omega (because 2d) estimated by the L^2 norm of E on the boundary. Then ||E|| in L^2 is of the order of (the measure of the cells of \Omega that touch the boundary)^{3/4}, which is h^3/4 .
Thus the error is going to zero slowly as h tends to zero.

When you measure the error numerically you have to normalize your estimate by the length of the boundary, which is about 36.
In the attached file I compute the normalized L^1 error on each of the 6 boundaries
ffbctest.edp (5.7 KB)
I get the errors
0.0839251 0.088971 0.0793675 0.00596786 0.00204029 0.0143004
You see that the normalized measured error is indeed not as bad as you thought!

Thanks very much for this analysis. I will examine the convergence rate and verify that exponent.