Neumann Boundary Condition for Elasticity Problem

General
1

Hi All,

I am trying to define a external force T=10 (along x-direction) on a surface of a rectangular geometry (60 * 10 * 5) on a surface with id=15 with actual area (5 * 10 = 50). I have specifically created these id’s in Gmsh and are labelled correctly. The surface id=16 is given a DirichletBC and works perfectly.

  1. But the ds(15) gives me 100. Can anyone help me with this?
  2. There is no error while running the program but the solution is also not converging for this.

I am attaching my files and code for reference.

https://drive.google.com/drive/folders/1UwOmKv_oxyOpJIStliJrya-b2NB5tgGu?usp=sharing

from fenics import *
import numpy as np

mesh = Mesh("Rectangle-v1.xml")

V = VectorFunctionSpace(mesh, 'P', 1)
u = Function(V)
v = TestFunction(V)

fd = MeshFunction('size_t', mesh, "Rectangle-v1_facet_region.xml")

ds = Measure("ds", domain=mesh, subdomain_data=fd)
ds_15 = ds(15)

E0 = 30e3
nu = 0.30
mu = E0/(2*(1 + nu))
lmbda = E0*nu/((1 + nu)*(1 - 2*nu))

def epsilon(u):
    return 0.5*(grad(u) + grad(u).T)

def sigma(u):
    return 2.0*mu*epsilon(u) + lmbda*tr(epsilon(u))*Identity(len(u))

bc_left = DirichletBC(V, Constant((0, 0, 0)), fd, 16)
bcs = [bc_left]

f = Constant((0, 0, 0))
T = Constant((10, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds_15
F = a - L

J=derivative(F, u)

# Solve the problem
u = Function(V)
solve(F == 0, u, J=J)


u_x, u_y, u_z = u.split(deepcopy=True)

file_u_x = File("Test1/u_x_plt.pvd")
file_u_x << u_x

file_u_y = File("Test1/u_y_plt.pvd")
file_u_y << u_y

file_u_z = File("Test1/u_z_plt.pvd")
file_u_z << u_z
2

Without the mesh file, and properly formatted code:

```python
#add code here
```

There is very little people can do to help you

3

Thanks for getting back to my post. I have edited as requested, if anything else needed please do let me know.

4

It should be 100, if you look at the area with id 15 in Paraview:

100_surface
100_surface1677×765 77 KB

As you can see from this figure, x=60, y goes from 0 to 20, z from 0 to 5.

There are no strong boundary conditions in your code:

and as the problem is linear, I would start by solving it with a linear solver:

bc_left = DirichletBC(V, Constant((0, 0, 0)), fd, 16)
bcs = [bc_left]

u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((10, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds_15

# Solve the problem
uh = Function(V)
solve(a==L,uh, bcs=bcs)