# Apply a Neumann boundary condition on a 3D cylinder domain

Dear all,

I'm working on a 3D cylinder domain (see figure 5 in [1] https://www.dropbox.com/s/ghz0oddt35dx7x9/mesh.png?dl=0). The boundary condition will be:
(1) perfect rough boundary for the top and bottom (no horizontal sliding), the top will move downwards to simulate the loading process;
(2) constant confining pressure for the side of the cylinder (here is the problem).

I'm using escript 4.2 and the following code to apply the boundary condition:

vel = -0.1e-3; confining=-4.e6;
ro=0.025; lz = 0.1; # sample dimension
x = mydomain.getX()
bx = FunctionOnBoundary(mydomain).getX()
outsideBF = whereZero(sqrt(bx[0]**2+bx[1]**2)-ro,rtol=2.e-2)
Dbc = whereZero(x[2])*[1,1,1]+whereZero(x[2]-lz)*[1,1,1]
Vbc = whereZero(x[2])*[0,0,0]+whereZero(x[2]-lz)*[0,0,vel]
# Neumann BC, constant lateral confining
Nbc = outsideBF*confining*mydomain.getNormal()

From my results, the applied confining pressure is smaller than I expected (around -3.2e6 https://www.dropbox.com/s/96agcz0n2r6dxh6/stress_x.png?dl=0).

Is there any suggestions for me? Thanks in advance.

## Question information

Language:
English Edit question
Status:
For:
esys-escript Edit question
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Lutz Gross (l-gross) said on 2018-04-09: #1

How did you visualize or calculate the " applied confining pressure"?

 Revision history for this message whr (huanran-wu) said on 2018-04-10: #2

I use saveVTK to save the stress state of the gauss point.

If the confining pressure is applied properly, I expect all the gauss points will have stress along x[0] and x[1] direction close to -4e6. At least, for the elements close to the boundary and almost perpendicular to x[0] direction, they will have stress along x[0] direction stress[0][0] close to -4e6.

Am I right?

 Revision history for this message Lutz Gross (l-gross) said on 2018-04-10: #3

For VTK the stresses are averaged over the elements which can be seen as approximations at the element centers.
This moves the location of inspection even more away from boundary. Have you tried what happened when you refine the mesh?
Unfortunately, escript does not provide an easy mechanism to recover the stress on the surface.

 Revision history for this message whr (huanran-wu) said on 2018-04-11: #4

Yeah, that's what I expect.
Actually, I'm conducting a FEMxDEM simulation and the initial state of gauss points (DEM packing consists of discrete particles) is of -4.e6 stress along x[0] and x[1]. So I expect at the initial state, with the boundary condition I apply, the stress for the element will be the same -4.e6.

The finer mesh for this cylinder problem is not possible because of the computational cost is too high.

Now, I see the problem might be the shear stress due to the coarse mesh. Actually, I have archived reasonable initial state in a 2D case with much finer mesh (80 along the hoop direction).

Thanks anyway.

 Revision history for this message Joel Fenwick (j-fenwick1) said on 2018-04-11: #5

.