Fourth-order partial differential equation is solved but not satisfied in the bulk

No, I am also imposing boundary conditions on the normal derivative via Nitsche’s method, see

F_N = alpha / r_mesh * ( \
            + (((n_circle( omega ))[i] * omega[i] - omega_r) * ((n_circle( omega ))[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth_circle( omega, c_r ) * (1.0 / r) * ds_r \
            + (((n_circle( omega ))[i] * omega[i] - omega_R) * ((n_circle( omega ))[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth_circle( omega, c_R ) * (1.0 / R) * ds_R \
    )
[...]
F = (F_z + F_omega + F_mu + F_nu) + F_N

where n_circle() is the boundary normal and \omega_i = \nabla_i z. The boundary conditions are correct, I suspect that there is something else…
I observed a similar behavior in a simpler example for the bi-harmonic equation to which you replied, I posted it here, maybe it can help understand the problem in this post too ! :slight_smile:

Thank you for your help