small_strain_vonMises_plasticity.py 2.63 KB
 Jeremy BLEYER committed Feb 21, 2020 1 2 3 ``````from dolfin import * import mfront_wrapper as mf import numpy as np `````` Jeremy BLEYER committed Feb 24, 2020 4 ``````import ufl `````` Jeremy BLEYER committed Feb 21, 2020 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) `````` Jeremy BLEYER committed Mar 06, 2020 16 ``````mesh = Mesh("meshes/thick_cylinder.xml") `````` Jeremy BLEYER committed Feb 21, 2020 17 `````` `````` Jeremy BLEYER committed Mar 06, 2020 18 ``````facets = MeshFunction("size_t", mesh, "meshes/thick_cylinder_facet_region.xml") `````` Jeremy BLEYER committed Feb 21, 2020 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) `````` Jeremy BLEYER committed Feb 24, 2020 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} `````` Jeremy BLEYER committed Mar 06, 2020 43 `````` material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so", `````` Jeremy BLEYER committed Feb 24, 2020 44 `````` "IsotropicLinearHardeningPlasticity", `````` Jeremy BLEYER committed Feb 24, 2020 45 46 47 48 `````` hypothesis=hypothesis, material_properties=mat_prop) problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=4) problem.bc = bc `````` Jeremy BLEYER committed Feb 24, 2020 49 `````` problem.set_loading(loading*dot(n, u)*measure*ds(4)) `````` Jeremy BLEYER committed Feb 21, 2020 50 `````` `````` Jeremy BLEYER committed Mar 08, 2020 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, )) `````` Jeremy BLEYER committed Feb 21, 2020 55 `````` `````` Jeremy BLEYER committed Feb 24, 2020 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") `````` Jeremy BLEYER committed Feb 21, 2020 61 `````` `````` Jeremy BLEYER committed Feb 24, 2020 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()) `````` Jeremy BLEYER committed Feb 21, 2020 68 `````` `````` Jeremy BLEYER committed Feb 24, 2020 69 `````` file_results.write(u, t) `````` Jeremy BLEYER committed Feb 21, 2020 70 `````` `````` Jeremy BLEYER committed Mar 08, 2020 71 72 73 `````` p_avg.assign(project(p, P0)) file_results.write(p_avg, t) results[i+1, :] = (u(Ri, 0)[0], t) `````` Jeremy BLEYER committed Feb 21, 2020 74 75 `````` `````` Jeremy BLEYER committed Feb 24, 2020 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()``````