Neumann boundary conditions and different geometries

[JanSuntajs]

Hello,

while solving some linear elasticity problems on a sphere, I've encountered issues probably related to the specification of the Neumann boundary conditions at the sphere's outer surface.

In formulation of the problem, I follow the Gridap tutorial on linear elasticity. I use a spherical mesh generated using the Open Cascade Kernel in GridapGmsh. As I typically use the no-displacement boundary condition at the sphere's center, I've assigned a physical domain named Center to a single point at the sphere's origin. The other physical domain is named Casing and consists of the sphere's outer surface. As per the gotchas section of GridapGmsh I took special care to also include the 0-d and 1-d entities to the Casing physical region. The code seems to work just fine if I specify the zero-displacement Dirichlet boundary condition at the Center and some other Dirichlet-type boundary condition at the Casing. However, it fails to produce meaningful results once I try to apply the Neumann boundary conditions at the Casing. Typically, the magnitude of the displacement vector would be abnormally large or/and would not exhibit the symmetry properties of the applied external pressure. Below I illustrate this with two MWEs for the above mentioned cases.

The MWE for the working case with Dirichlet BCs at both physical regions:

using Gridap
using Gridap.Geometry

model = DiscreteModelFromFile("./sphere_shape_1.0_lc_0.05_.msh")  #unzip and save the appended .msh file
writevtk(model, "model")

# prepare the vector-valued FE space
order = 1
reffe = ReferenceFE(lagrangian, VectorValue{3, Float64}, order)
V0 = TestFESpace(model,reffe;
  conformity=:H1,
  dirichlet_tags=["Center", "Casing"],
  dirichlet_masks=[(true, true, true), (true, true,true)])

# specify the BC - displacement in the radial direction, modulated by some periodic function
g_ = 0.001
absx(x) = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
g1(x) = g_ * cos(x[3] * pi / 2) * cos(x[2] * pi / 2) * cos(x[1] * pi/2)* VectorValue(x[1], x[2], x[3]) / absx(x)
g0(x) = VectorValue(0., 0., 0.)

U = TrialFESpace(V0, [g0, g1])

# constitutive law
const E_ = 70.0e9
const ν = 0.33
const λ = (E_*ν)/((1+ν)*(1-2*ν))
const μ = E_/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

# weak form
degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
l(v) = 0

op = AffineFEOperator(a,l,U,V0)
uh = solve(op)

writevtk(Ω, "sphere_dirichlet_bc",cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ∘ε(uh)])

In thise case, the results seem to be reasonable, as
shown by these ParaView plots:
.

Now, let's turn to the combined case with Dirichlet and Neumann BC. In that case, the code is as follows:

using Gridap
using Gridap.Geometry

# constants and constitutive relations
const E_ = 70.0e9
const ν = 0.33
const λ = (E_*ν)/((1+ν)*(1-2*ν))
const μ = E_/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

# load the model
model = DiscreteModelFromFile("./sphere_shape_1.0_lc_0.05_.msh")  #unzip and save the appended .msh file
writevtk(model, "model")

# prepare the vector-valued FE space
order = 1
reffe = ReferenceFE(lagrangian, VectorValue{3, Float64}, order)
V0 = TestFESpace(model,reffe;
  conformity=:H1,
  dirichlet_tags=["Center"],
  dirichlet_masks=[(true, true, true)])

g0(x) = VectorValue(0., 0., 0.)

U = TrialFESpace(V0, [g0])

# prepare the Neumann surfaces and related quantities
degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

neumann_tags = ["Casing"]
Γ = BoundaryTriangulation(model, tags=neumann_tags)
dΓ = Measure(Γ, degree)
n_Γ = get_normal_vector(Γ)
gvn(x) =  E_ * 1e-004 * cos(x[3] * pi / 2) * cos(x[2] * pi / 2) * cos(x[1] * pi/2)


a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
l(v) = ∫(gvn * (n_Γ ⋅ v)) * dΓ

op = AffineFEOperator(a,l,U,V0)
ls = BackslashSolver()
solver = LinearFESolver(ls)
uh = solve(solver, op)

writevtk(Ω,"sphere_dirichlet_neumann_bc",cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ∘ε(uh)])

This approach yields the following (wrong) results:

Are there any glaring errors in my above approach? If so, what should I do to fix them? I assume this is somehow related to the specification of the boundary conditions; for instance, if I solve the same problem on a cylindrical domain, applying the Dirichlet BC at its bases and Neumann BC along the casing, the problem does not appear.
Thank you for your response.
Jan

[arnold-pdev]

I don't think I can help you, but wanted to say that Neumann boundary conditions are not named for John von Neumann, but German mathematician Carl Neumann.

[JanSuntajs]

I don't think I can help you, but wanted to say that Neumann boundary conditions are not named for John von Neumann, but German mathematician Carl Neumann.

Thanks, edited :)

[JanSuntajs]

I have since managed to pinpoint the culprit and solve my problem; it would say it is not that much of a Gridap problem, more a FEM related one. If only Neumann conditions are specified at the sphere's outer shell, the solutions are still symmetric against rigid rotations.

This (with properly tuned $\alpha$ effectively solved my problem:

degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

neumann_tags = ["Casing"]
Γ = BoundaryTriangulation(model, tags=neumann_tags)
dΓ = Measure(Γ, degree)
n_Γ = get_normal_vector(Γ)

α = 0.1
a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ + ∫(α * (curl(u)) ⋅ (curl(v)))*dΩ # stabilization term
l(v) = ∫(pbound * (n_Γ ⋅ v)) * dΓ

Regards,
Jan