small_strain_vonMises_plasticity.py 2.63 KB
Newer Older
1 2 3
from dolfin import *
import mfront_wrapper as mf
import numpy as np
4
import ufl
5 6 7 8 9 10 11 12 13 14 15

Re, Ri = 1.3, 1.   # external/internal radius
# elastic parameters
E = 70e3
nu = 0.3
# yield strength
sig0 = 250.
Et = E/100.
# hardening slope
H = E*Et/(E-Et)

16
mesh = Mesh("meshes/thick_cylinder.xml")
17

18
facets = MeshFunction("size_t", mesh, "meshes/thick_cylinder_facet_region.xml")
19 20 21 22 23 24 25
ds = Measure('ds', subdomain_data=facets)


V = VectorFunctionSpace(mesh, "CG", 2)
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
n = FacetNormal(mesh)

26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
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)
        measure = 1
    elif hypothesis == "axisymmetric":
        x = SpatialCoordinate(mesh)
        q_lim = float(2*ln(Re/Ri)*sig0)
        measure = 2*pi*abs(x[0])
    print("---------------------------------------------")
    loading = Expression("-q*t", q=q_lim, t=0, degree=2)

    mat_prop = {"YoungModulus": E,
               "PoissonRatio": nu,
               "HardeningSlope": H,
               "YieldStrength": sig0}
43
    material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
44
                                          "IsotropicLinearHardeningPlasticity",
45 46 47 48
                                          hypothesis=hypothesis,
                                          material_properties=mat_prop)
    problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4)
    problem.bc = bc
49
    problem.set_loading(loading*dot(n, u)*measure*ds(4))
50

51 52 53 54
    p = problem.get_state_variable("EquivalentPlasticStrain")
    assert (ufl.shape(p) == ())
    epsel = problem.get_state_variable("ElasticStrain")
    assert (ufl.shape(epsel)==(4, ))
55

56 57 58 59 60
    file_results = XDMFFile("results/plasticity_{}_results.xdmf".format(hypothesis))
    file_results.parameters["flush_output"] = True
    file_results.parameters["functions_share_mesh"] = True
    P0 = FunctionSpace(mesh, "DG", 0)
    p_avg = Function(P0, name="Plastic_strain")
61

62 63 64 65 66 67
    Nincr = 40
    load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
    results = np.zeros((Nincr+1, 2))
    for (i, t) in enumerate(load_steps):
        loading.t = t
        problem.solve(u.vector())
68

69
        file_results.write(u, t)
70

71 72 73
        p_avg.assign(project(p, P0))
        file_results.write(p_avg, t)
        results[i+1, :] = (u(Ri, 0)[0], t)
74 75


76 77 78 79 80
    import matplotlib.pyplot as plt
    plt.plot(results[:, 0], results[:, 1], "-o")
    plt.xlabel("Displacement of inner boundary")
    plt.ylabel(r"Applied pressure $q/q_{lim}$")
    plt.show()