### (Deleted) Elastic Membrane Problems

48

views

0

To solve a elastic problem i generated a Cylindrical mesh using Gmsh and imported as XML file in fenics.

1st question) Does the xml file, converted with dolfin keep exact the same mesh ?

I use Jupyter inside a Docker container and i converted the file from inside . I say this because i wrote a code for

a simple elastic problem as in the tutorial but it seems it has a few problems .

from fenics import *

from dolfin import *

from mshr import *

import numpy as np

#Geometries and Meshing with GMSH

mesh = Mesh("me.xml");

plot(mesh)

#Elasticity

V = VectorFunctionSpace(mesh, "P", 1)

#Constants

d = 799

r = d/2

mu = 1

rho = 1

beta = 1.25

lambda_ = beta

g = 9.81

#Def boundary condition, with tolerance

tol = 1E-14

def clamped_boundary(x, on_boundary):

return on_boundary and np.abs(x[0]**2 + x[1]**2 > r**2)

bc = DirichletBC(V, (0,0,0), clamped_boundary)

# Define Strain

def epsilon(u):

return 0.5*(grad(u)+grad(u).T)

def sigma(u):

return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem

u = TrialFunction(V)

v = TestFunction(V)

f = Constant((0, 0, -rho*g))

d = u.geometric_dimension() #SPACE DIMENSION

# print(d)

T = Constant((0,0,0))

#GOVERNING EQUATION

a = inner(sigma(u), epsilon(v))*dx

L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution

u = Function(V)

solve(a == L, u, bc)

vtkfile = File("displacement.pvd")

vtkfile << (u)

#plot stresses

s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) # deviatoric stress

von_Mises = sqrt(3./2*inner(s, s))

V = FunctionSpace(mesh, "Lagrange", 2)

von_Mises = project(von_Mises, V)

plot(von_Mises, title="Stresses")

vtkfile = File("stresses.pvd")

vtkfile << (von_Mises)

This is the code and these are the results:

Geometry and mesh imported from gmsh and resaved as vtk file, it looks completely fine but...

1st question) Does the xml file, converted with dolfin keep exact the same mesh ?

I use Jupyter inside a Docker container and i converted the file from inside . I say this because i wrote a code for

a simple elastic problem as in the tutorial but it seems it has a few problems .

from fenics import *

from dolfin import *

from mshr import *

import numpy as np

#Geometries and Meshing with GMSH

mesh = Mesh("me.xml");

plot(mesh)

#Elasticity

V = VectorFunctionSpace(mesh, "P", 1)

#Constants

d = 799

r = d/2

mu = 1

rho = 1

beta = 1.25

lambda_ = beta

g = 9.81

#Def boundary condition, with tolerance

tol = 1E-14

def clamped_boundary(x, on_boundary):

return on_boundary and np.abs(x[0]**2 + x[1]**2 > r**2)

bc = DirichletBC(V, (0,0,0), clamped_boundary)

# Define Strain

def epsilon(u):

return 0.5*(grad(u)+grad(u).T)

def sigma(u):

return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem

u = TrialFunction(V)

v = TestFunction(V)

f = Constant((0, 0, -rho*g))

d = u.geometric_dimension() #SPACE DIMENSION

# print(d)

T = Constant((0,0,0))

#GOVERNING EQUATION

a = inner(sigma(u), epsilon(v))*dx

L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution

u = Function(V)

solve(a == L, u, bc)

vtkfile = File("displacement.pvd")

vtkfile << (u)

#plot stresses

s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) # deviatoric stress

von_Mises = sqrt(3./2*inner(s, s))

V = FunctionSpace(mesh, "Lagrange", 2)

von_Mises = project(von_Mises, V)

plot(von_Mises, title="Stresses")

vtkfile = File("stresses.pvd")

vtkfile << (von_Mises)

This is the code and these are the results:

Geometry and mesh imported from gmsh and resaved as vtk file, it looks completely fine but...

Displecement

Stresses and here it starts to show problems, it could be still fine but ...

And this is the displecemnt once i show WarpByVector, and it's not that i m watching only one axis , if i rotate it still keep being a spike like if the elements fall in themself. When if i use WarpByVector in the Linear Elastic Beam problem in the tutorial everything works fine.

I'm thinking there are some problem with the mesh but i'm not sure about it.

Thanks to everyone

Community: FEniCS Project

### 1 Answer

1

Your displacement is of the order 1e+7, that's huge.

The diameter of your disc is 799 meters while the Lamé constants are 1 Pa (should be in the range of 1e7-1e9)

You can set the scale factor really small in warp by vector. But I'd start using reasonable values for the constants.

Also:

If you use piecewise linear function space for the displacements, why would you project the mises stress onto a piecewise quadratic space? Naturally, the Mises stress should be of the same order as your strains, in this case come from a piecewise constant space:

Maybe this would lift the problem of the awkward oscillations in your Mises plot.

The diameter of your disc is 799 meters while the Lamé constants are 1 Pa (should be in the range of 1e7-1e9)

You can set the scale factor really small in warp by vector. But I'd start using reasonable values for the constants.

Also:

If you use piecewise linear function space for the displacements, why would you project the mises stress onto a piecewise quadratic space? Naturally, the Mises stress should be of the same order as your strains, in this case come from a piecewise constant space:

`DGSpace = FunctionSpace(mesh, FiniteElement('DG', mesh.ufl_cell, 0))`

Maybe this would lift the problem of the awkward oscillations in your Mises plot.

Please login to add an answer/comment or follow this question.

The thread is closed. No new answer/comment may be added.