Commit 6879e9ac authored by Jeremy BLEYER's avatar Jeremy BLEYER

Plane strain and axisymmetric computations in small strain von Mises

parent 65db9759
...@@ -3,8 +3,6 @@ import mfront_wrapper as mf ...@@ -3,8 +3,6 @@ import mfront_wrapper as mf
import numpy as np import numpy as np
import ufl import ufl
hypothesis = "axisymmetric" # axisymmetric
Re, Ri = 1.3, 1. # external/internal radius Re, Ri = 1.3, 1. # external/internal radius
# elastic parameters # elastic parameters
E = 70e3 E = 70e3
...@@ -20,47 +18,51 @@ mesh = Mesh("../meshes/thick_cylinder.xml") ...@@ -20,47 +18,51 @@ mesh = Mesh("../meshes/thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "../meshes/thick_cylinder_facet_region.xml") facets = MeshFunction("size_t", mesh, "../meshes/thick_cylinder_facet_region.xml")
ds = Measure('ds', subdomain_data=facets) ds = Measure('ds', subdomain_data=facets)
if hypothesis == "plane_strain":
V = VectorFunctionSpace(mesh, "CG", 2)
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
n = FacetNormal(mesh)
for hypothesis in ["plane_strain", "axisymmetric"]:
u = Function(V, name="Displacement")
if hypothesis == "plane_strain":
print("Expansion of a thick cylinder in plane strain")
q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0) q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0)
measure = 1 measure = 1
elif hypothesis == "axisymmetric": elif hypothesis == "axisymmetric":
x = SpatialCoordinate(mesh) x = SpatialCoordinate(mesh)
q_lim = float(2*ln(Re/Ri)*sig0) q_lim = float(2*ln(Re/Ri)*sig0)
measure = 2*pi*abs(x[0]) measure = 2*pi*abs(x[0])
print("---------------------------------------------")
loading = Expression("-q*t", q=q_lim, t=0, degree=2)
V = VectorFunctionSpace(mesh, "CG", 2) mat_prop = {"YoungModulus": E,
u = Function(V, name="Displacement")
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
n = FacetNormal(mesh)
loading = Expression("-q*t", q=q_lim, t=0, degree=2)
mat_prop = {"YoungModulus": E,
"PoissonRatio": nu, "PoissonRatio": nu,
"HardeningSlope": H, "HardeningSlope": H,
"YieldStrength": sig0} "YieldStrength": sig0}
material = mf.MFrontNonlinearMaterial('../materials/src/libBehaviour.so', material = mf.MFrontNonlinearMaterial('../materials/src/libBehaviour.so',
'IsotropicLinearHardeningPlasticity', 'IsotropicLinearHardeningPlasticity',
hypothesis=hypothesis, hypothesis=hypothesis,
material_properties=mat_prop) material_properties=mat_prop)
problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4) problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4)
problem.set_loading(loading*dot(n, u)*measure*ds(4)) problem.set_loading(loading*dot(n, u)*measure*ds(4))
problem.bc = bc problem.bc = bc
p = problem.get_state_variable(name="EquivalentPlasticStrain") p = problem.get_state_variable(name="EquivalentPlasticStrain")
assert (ufl.shape(p) == ()) assert (ufl.shape(p) == ())
epsel = problem.get_state_variable(name="ElasticStrain") epsel = problem.get_state_variable(name="ElasticStrain")
assert (ufl.shape(epsel)==(4, )) assert (ufl.shape(epsel)==(4, ))
file_results = XDMFFile("results/plasticity_results.xdmf") file_results = XDMFFile("results/plasticity_{}_results.xdmf".format(hypothesis))
file_results.parameters["flush_output"] = True file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0) P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic_strain") p_avg = Function(P0, name="Plastic_strain")
Nincr = 40 Nincr = 40
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5 load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2)) results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps): for (i, t) in enumerate(load_steps):
loading.t = t loading.t = t
problem.solve(u.vector()) problem.solve(u.vector())
...@@ -71,8 +73,8 @@ for (i, t) in enumerate(load_steps): ...@@ -71,8 +73,8 @@ for (i, t) in enumerate(load_steps):
results[i+1, :] = (u(Ri, 0)[0], t) results[i+1, :] = (u(Ri, 0)[0], t)
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
plt.plot(results[:, 0], results[:, 1], "-o") plt.plot(results[:, 0], results[:, 1], "-o")
plt.xlabel("Displacement of inner boundary") plt.xlabel("Displacement of inner boundary")
plt.ylabel(r"Applied pressure $q/q_{lim}$") plt.ylabel(r"Applied pressure $q/q_{lim}$")
plt.show() plt.show()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment