multiphase_model.py 2.87 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Feb  1 08:49:11 2019

@author: bleyerj
"""
from dolfin import *
import mfront_wrapper as mf
import numpy as np
from ufl import diag

width = 1.
height = 1.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
15
mesh = RectangleMesh(Point(0., 0.), Point(width, height), 10, 10)
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

Ve = VectorElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([Ve, Ve]))
u = Function(V, name="Displacements")
(u1, u2) = split(u)

def bottom(x, on_boundary):
    return near(x[1], 0) and on_boundary
def top(x, on_boundary):
    return near(x[1], height) and on_boundary
def left(x, on_boundary):
    return near(x[0], 0) and on_boundary


bc = [DirichletBC(V.sub(0).sub(0), Constant(0), left),
      DirichletBC(V.sub(1).sub(0), Constant(0), left),
      DirichletBC(V.sub(0).sub(1), Constant(0), bottom),
      DirichletBC(V.sub(1).sub(1), Constant(0), bottom),
      DirichletBC(V.sub(0).sub(1), Constant(-1), top),
      DirichletBC(V.sub(1).sub(1), Constant(-1), top)]

facets = MeshFunction("size_t", mesh, 1)
ds = Measure("ds", subdomain_data=facets)


Jeremy BLEYER's avatar
Jeremy BLEYER committed
41 42 43 44 45 46 47
mat_prop = {"MatrixYoungModulus": 70e3,
            "MatrixPoissonRatio": 0.2,
            "FiberYoungModulus": 200e3,
            "FiberPoissonRatio": 0.3,
            "FiberVolumeFraction": 0.1,
            "Size": 0.1}

48
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
49 50 51
                                      "MultiphaseModel",
                                      hypothesis="plane_strain",
                                      material_properties=mat_prop)
52 53 54 55 56 57 58 59 60

problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
problem.bc = bc
problem.register_gradient("MatrixStrain", sym(grad(u1)), symmetric=True)
problem.register_gradient("InclusionStrain", sym(grad(u2)), symmetric=True)
problem.register_gradient("RelativeDisplacement", diag(u2-u1), symmetric=True)

problem.solve(u.vector())

Jeremy BLEYER's avatar
Jeremy BLEYER committed
61
x = np.linspace(width/2, width, 21)
62 63
import matplotlib.pyplot as plt
plt.figure()
Jeremy BLEYER's avatar
Jeremy BLEYER committed
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "ob", label=r"$u_x$ (matrix, size={})".format(mat_prop["Size"]))
plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "or", label=r"$u_x$ (fiber, size={})".format(mat_prop["Size"]))

# mat_prop["Size"] = 0.02

# problem2 = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
# problem2.bc = bc
# problem2.register_gradient("MatrixStrain", sym(grad(u1)), symmetric=True)
# problem2.register_gradient("InclusionStrain", sym(grad(u2)), symmetric=True)
# problem2.register_gradient("RelativeDisplacement", diag(u2-u1), symmetric=True)
# u = Function(V)
# problem2.solve(u.vector())

# plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "sb", markerfacecolor="None", label=r"$u_x$ (matrix, size={})".format(mat_prop["Size"]))
# plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "sr", markerfacecolor="None", label=r"$u_x$ (fiber, size={})".format(mat_prop["Size"]))

80 81
plt.legend()
plt.show()