Hello
I do a small simulation with FEniCS and plot the results with Paraview, I notice the figure doesn’t make sense. I have posted this problem on fenics support google group, but no responds yet;
Sorry, I have to re-post here for the answer.
Thanks a lot
Suess
Please repost the full code and images from that post, as many users are not members of the Google group. Discourse is where such questions belong.
Maybe @bleyerj can shed some light on this, as he is the author of the original problem. I am unfortunately not very familiar with Orthotropic linear elasticity.
Hi,
yes sorry. It is a silly typo in the code. The compliance matrix S is missing some minus sign on the extra diagonal terms, as correctly written in the formula in the numerical tour. Without that a positive Poisson ratio, yields transverse expansion rather than contraction under uniaxial tension.
I will correct it on the numerical tours website.
Thanks for pointing this out.
Hi Dokken,
I did small simulation with FEniCS, just stretching a 2D plane to the right. I modified an example program called orthotropic_elasticity.py, it worked and I sent the result (displacement) to Paraview. I opened the graph, I got the displacement on the X direction was good, positive on the right and zero on the left, but when I noticed the displacement on Y direction, it was wrong (I think). It’s positive on top and negative on bottom, seems the 2D plane expands vertically when I stretch it horizontally. Why? Is it a bug in FEniCS / Paraview? I’m confused, please help me.
I’m attaching the graph and the script. Thanks in advance.
Suess!
PS. I’d like to attach the script, but it’s not allowed, so put it here:
"""
This script is modification of the orthotropic_elasticity.py, obtained from:
https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/orthotropic_elasticity.py.html
A script to test 2D plane under tension.
"""
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
L = 1.
N = 20 # mesh density
domain = Rectangle(Point(0.,0.), Point(L, L))
mesh = generate_mesh(domain, N)
plot(mesh)
plt.show()
# Stress Strain
Ex, Ey, nuxy, Gxy = 100., 10., 0.3, 5.
S = as_matrix([[1./Ex,nuxy/Ex,0.],[nuxy/Ex,1./Ey,0.],[0.,0.,1./Gxy]])
C = inv(S)
def eps(v):
return sym(grad(v))
def strain2voigt(e):
"""e is a 2nd-order tensor, returns its Voigt vectorial representation"""
return as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s):
"""
s is a stress-like vector (no 2 factor on last component)
returns its tensorial representation
"""
return as_tensor([[s[0], s[2]],
[s[2], s[1]]])
def sigma(v):
return voigt2stress(dot(C, strain2voigt(eps(v))))
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0) and on_boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L) and on_boundary
# exterior facets MeshFunction
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Left().mark(facets, 1)
Right().mark(facets, 2)
ds = Measure('ds', subdomain_data=facets)
# Define function space
V = VectorFunctionSpace(mesh, 'Lagrange', 1)
# Define variational problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = Function(V, name='Displacement')
a = inner(sigma(du), eps(u_))*dx
# uniform traction on the right
T = Constant((10,0))
l = dot(T, u_)*ds(2)
# boundary condition
bc = DirichletBC(V, Constant((0,0)), facets, 1)
solve(a == l, u, bc)
# Save solution in VTK format
file = File("displ_xy.pvd");
file << u
Hi Bleyerj
So, is this the reason that the program gives strange result on Y displacement?
Yes, it’s as if you had used a real Poisson ratio of -0.3. instead of 0.3, so expansion rather than contraction.
