logarithmic_strain_plasticity.py 1.84 KB
Newer Older
1
import matplotlib.pyplot as plt
2 3 4
from dolfin import *
import mfront_wrapper as mf
import numpy as np
5
import ufl
6 7 8 9 10 11

length, width , height = 1., 40e-3, 100e-3
nb_elt_p, nb_elt_l = 10, 30
mesh = BoxMesh(Point(0, -width/2, -height/2.), Point(length, width/2, height/2.), nb_elt_l, 2, nb_elt_p)

V = VectorFunctionSpace(mesh, "CG", 2)
12
u = Function(V, name="Displacement")
13 14 15 16 17 18 19 20 21

def left(x, on_boundary):
    return near(x[0], 0) and on_boundary
def right(x, on_boundary):
    return near(x[0], length) and on_boundary

bc = [DirichletBC(V, Constant((0.,)*3), left),
      DirichletBC(V.sub(0), Constant(0.), right)]

Jeremy BLEYER's avatar
Jeremy BLEYER committed
22
file_results = XDMFFile("results/finie_strain_plasticity.xdmf")
23 24 25 26
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

selfweight = Expression(("0", "0", "-t*qmax"), t=0., qmax = 50e6, degree=0)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
27

28
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
29
                                      "LogarithmicStrainPlasticity")
30
problem = mf.MFrontNonlinearProblem(u, material, bcs=bc)
31 32
problem.set_loading(dot(selfweight, u)*dx)

33 34
p = problem.get_state_variable("EquivalentPlasticStrain")
epsel = problem.get_state_variable("ElasticStrain")
35 36
print(ufl.shape(epsel))

37 38 39
prm = problem.solver.parameters
prm["absolute_tolerance"] = 1e-6
prm["relative_tolerance"] = 1e-6
40
prm["linear_solver"] = "mumps"
41 42 43

P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
Jeremy BLEYER's avatar
Jeremy BLEYER committed
44

45 46 47 48 49 50
Nincr = 30
load_steps = np.linspace(0., 1., Nincr+1)
results = np.zeros((Nincr+1, 3))
for (i, t) in enumerate(load_steps[1:]):
    selfweight.t = t

Jeremy BLEYER's avatar
Jeremy BLEYER committed
51 52 53
    # problem.solve(u.vector())

    results[i+1, 0] = u(length, 0, 0)
54 55 56
    results[i+1, 1] = t

    file_results.write(u, t)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
57 58
    p_avg.assign(project(p, V0))
    file_results.write(p_avg, t)
59 60 61 62 63

plt.figure()
plt.plot(results[:, 0], results[:, 1], "-o")
plt.xlabel(r"Displacement")
plt.ylabel("Load")
Jeremy BLEYER's avatar
Jeremy BLEYER committed
64