esys-escript

Apply a Neumann boundary condition on a 3D cylinder domain

Asked by whr

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
mydomain = ReadGmsh("cylinder.msh",numDim=3,integrationOrder=2); # read mesh info
x = mydomain.getX()
bx = FunctionOnBoundary(mydomain).getX()
outsideBF = whereZero(sqrt(bx[0]**2+bx[1]**2)-ro,rtol=2.e-2)
# Dirichlet BC, rough loading ends
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.

[1] https://doi.org/10.1016/j.compgeo.2016.01.020

Question information

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

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

Revision history for this message
whr (huanran-wu) said : 		#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 : 		#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 : 		#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 : 		#5

.

Can you help with this problem?

Provide an answer of your own, or ask whr for more information if necessary.

To post a message you must log in.