Updated phase field demo with converging tension-compression split model

demos/phase_field.py

-
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
@@ -8,34 +7,98 @@ Created on Fri Feb  1 08:49:11 2019
 """
 from dolfin import *
 import mfront_wrapper as mf
+from mfront_wrapper.utils import local_project
+from ufl import Max
+import matplotlib.pyplot as plt
+from mshr import *
+import mgis.behaviour as mgis_bv
 
-N = 100
-mesh = UnitSquareMesh(N, 2*N, "crossed")
+N = 75
+domain = Rectangle(Point(0, 0), Point(1, 1)) - Circle(Point(0.5, 0.5), 0.2, 40)
+mesh = generate_mesh(domain, N)
+plot(mesh, linewidth=0.2)
+plt.show()
 
+Vu = VectorFunctionSpace(mesh, "CG", 1)
+def top(x, on_boundary):
+    return near(x[1], 1) and on_boundary
+def internal(x, on_boundary):
+    return near((x[0]-0.5)**2+(x[1]-0.5)**2, 0.2**2, 0.05) and on_boundary
+Uimp = Expression(("0", "t"), t=0, degree=0)
+bcu = [DirichletBC(Vu, Constant((0, 0)), internal),
+       DirichletBC(Vu, Uimp, top)]
 
-V = FunctionSpace(mesh, "CG", 1)
-def crack(x, on_boundary):
-    return near(x[0], 0.5, 0.2) and near(x[1], 0.5)
-def border(x, on_boundary):
-    return on_boundary
+Vd = FunctionSpace(mesh, "CG", 1)
+bcd = DirichletBC(Vd, Constant(0.), internal)
 
-bc = DirichletBC(V, Constant(1.), crack)
+u  = Function(Vu, name="Displacement")
+d = Function(Vd, name="Damage")
+dold = Function(Vd, name="Previous damage")
 
-d = Function(V, name="Damage")
 
-material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
-                                      "PhaseField",
+material_u = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
+                                      "PhaseFieldDisplacement",
                                       hypothesis="plane_strain",
-                                      material_properties={"RegularizationLength": 0.02})
+                                      material_properties={"YoungModulus": 200.,
+                                                            "PoissonRatio": 0.2})
+problem_u = mf.MFrontNonlinearProblem(u, material_u, quadrature_degree=0)
+problem_u.bc = bcu
+problem_u.integration_type = mgis_bv.IntegrationType.IntegrationWithTangentOperator
+d0 = problem_u.get_state_variable("Damage")
+psi =problem_u.get_state_variable("DamageDrivingForce")
+prm_u = problem_u.solver.parameters
+prm_u["report"] = False
+
+def solve_u(t):
+    Uimp.t = t
+    d0.interpolate(d)
+    problem_u.affect_state_variables()
+    problem_u.solve(u.vector())
+
+Gc = 1.
+l0 = 0.02
+material_d = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
+                                      "PhaseFieldDamage",
+                                      hypothesis="plane_strain",
+                                      material_properties={"RegularizationLength": l0,
+                                                            "FractureEnergy": Gc})
+
+problem_d = mf.MFrontNonlinearProblem(d, material_d, quadrature_degree=0)
+problem_d.register_gradient("Damage", d)
+problem_d.register_gradient("DamageGradient", grad(d))
+H = problem_d.get_state_variable("HistoryFunction")
+prm_d = problem_d.solver.parameters
+prm_d["report"] = False
+
+# positive part
+ppos = lambda x: (x+abs(x))/2
+model = "without_threshold"
+if model == "without_threshold":
+    newH = psi
+elif model == "with_threshold":
+    newH = ppos(psi-Gc/l0)
+def solve_d(bc=[]):
+    H.assign(local_project(Max(H, newH), H.function_space(), problem_d.dx))
+    problem_d.affect_state_variables()
+    problem_d.bc = bcd
+    dold.assign(d)
+    problem_d.solve(d.vector())
+
+tol = 1e-2
+import numpy as np
+for (i, t) in enumerate(np.linspace(0, 2e-1,21)[1:]):
+    res = 1.
+    j = 1
+    print("Time step:", i)
+    while res > tol:
+        solve_d(bcd)
+        solve_u(t)
+        res = np.max(d.vector().get_local()-dold.vector().get_local())
+        print("   Iteration {}: {}".format(j, res))
+        j += 1
+    p=plot(d, vmin=0, vmax=1)
+    plt.colorbar(p)
+    plt.savefig("phase_field_{:04d}.png".format(i), dpi=400)
+    plt.show()
 
-problem = mf.MFrontNonlinearProblem(d, material, quadrature_degree=0)
-problem.bc = bc
-problem.register_gradient("Damage", d)
-problem.register_gradient("DamageGradient", grad(d))
-problem.solve(d.vector())
-#TODO: manage state variables communication with external computation
 
-import matplotlib.pyplot as plt
-p=plot(d)
-plt.colorbar(p)
-plt.show()

demos/phase_field_coupled.py

-#!/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
-from mfront_wrapper.utils import local_project
-from ufl import Max
-import matplotlib.pyplot as plt
-from mshr import *
-import mgis.behaviour as mgis_bv
-
-N = 50
-mesh = UnitSquareMesh(N, 2*N, "crossed")
-domain = Rectangle(Point(0, 0), Point(1, 1)) - Circle(Point(0.5, 0.5), 0.2, 40)
-mesh = generate_mesh(domain, N)
-plot(mesh, linewidth=0.2)
-plt.show()
-
-Vu = VectorFunctionSpace(mesh, "CG", 1)
-def top(x, on_boundary):
-    return near(x[1], 1) and on_boundary
-def internal(x, on_boundary):
-    return near((x[0]-0.5)**2+(x[1]-0.5)**2, 0.2**2, 0.05) and on_boundary
-Uimp = Expression(("0", "t"), t=0, degree=0)
-bcu = [DirichletBC(Vu, Constant((0, 0)), internal),
-       DirichletBC(Vu, Uimp, top)]
-
-Vd = FunctionSpace(mesh, "CG", 1)
-bcd = DirichletBC(Vd, Constant(0.), internal)
-
-u  = Function(Vu, name="Displacement")
-d = Function(Vd, name="Damage")
-dold = Function(Vd, name="Previous damage")
-
-
-material_u = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
-                                      "TensionCompressionSplit",
-                                      hypothesis="plane_strain",
-                                      material_properties={"YoungModulus": 200.,
-                                                            "PoissonRatio": 0.2})
-problem_u = mf.MFrontNonlinearProblem(u, material_u, quadrature_degree=0)
-problem_u.bc = bcu
-problem_u.integration_type = mgis_bv.IntegrationType.IntegrationWithSecantOperator
-d0 = problem_u.get_state_variable("Damage")
-psi =problem_u.get_state_variable("EnergyDensity")
-def solve_u(t):
-    Uimp.t = t
-    d0.interpolate(d)
-    problem_u.affect_state_variables()
-    problem_u.solve(u.vector())
-
-material_d = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
-                                      "PhaseFieldCoupled",
-                                      hypothesis="plane_strain",
-                                      material_properties={"RegularizationLength": 0.02,
-                                                            "FractureEnergy": 1.})
-
-problem_d = mf.MFrontNonlinearProblem(d, material_d, quadrature_degree=0)
-problem_d.register_gradient("Damage", d)
-problem_d.register_gradient("DamageGradient", grad(d))
-H = problem_d.get_state_variable("HistoryFunction")
-# prm = problem_d.solver.parameters
-# prm["report"] = False
-# prm["absolute_tolerance"] = 1e-2
-# prm["relative_tolerance"] = 1e-2
-# prm["maximum_iterations"] = 1
-
-
-
-def solve_d(bc=[]):
-    H.assign(local_project(Max(H, psi), H.function_space(), problem_d.dx))
-    problem_d.affect_state_variables()
-    problem_d.bc = bcd
-    dold.assign(d)
-    problem_d.solve(d.vector())
-    p=plot(d)
-    plt.colorbar(p)
-    plt.show()
-
-tol = 1e-2
-import numpy as np
-for (i, t) in enumerate(np.linspace(0, 2e-1,21)[1:]):
-    res = 1.
-    i = 1
-    print("Time step:", i)
-    while res > tol:
-        solve_d(bcd)
-        solve_u(t)
-        res = np.max(d.vector().get_local()-dold.vector().get_local())
-        print("   Iteration {}: {}".format(i, res))
-        i += 1
-
-

materials/PhaseField.mfront → materials/PhaseFieldDamage.mfront

 @DSL DefaultGenericBehaviour;
-@Behaviour PhaseField;
+@Behaviour PhaseFieldDamage;
 
 @Gradient TVector g;
 g.setEntryName("DamageGradient");
@@ -16,13 +16,18 @@ Y.setEntryName("EnergyRelease");
 
 @MaterialProperty real l0;
 l0.setEntryName("RegularizationLength");
+@MaterialProperty real Gc;
+Gc.setEntryName("FractureEnergy");
+
+@StateVariable real H;
+H.setEntryName("HistoryFunction");
 
 @ProvidesTangentOperator;
 @Integrator {
   // remove useless warnings, as we always compute the tangent operator
   static_cast<void>(computeTangentOperator_);
-  ∂q∕∂Δg = pow(l0, 2)* tmatrix<N, N, real>::Id();
-  ∂Y∕∂Δd = 1;
-  q = pow(l0, 2)*g;
-  Y = ∂Y∕∂Δd ⋅ d;
+  ∂q∕∂Δg = Gc*l0* tmatrix<N, N, real>::Id();
+  ∂Y∕∂Δd = Gc/l0+2*H;
+  q = Gc*l0*(g+Δg);
+  Y = ∂Y∕∂Δd ⋅ (d+Δd)-2*H;
 }

materials/PhaseFieldDisplacement.mfront

+@DSL Default;
+@Author Tran Thang DANG;
+@Date 05/02/2016;
+@Behaviour PhaseFieldDisplacement;
+// @ModellingHypothesis {PlaneStrain,Tridimensional};
+@Description{
+  "A damage model coupled with a phase "
+  "field approach"
+}
+
+@MaterialProperty stress yg;
+yg.setGlossaryName("YoungModulus") ;
+@MaterialProperty real nu ;
+nu.setGlossaryName("PoissonRatio") ;
+@Parameter real kk = 1e-6 ;
+
+@StateVariable real H;
+H.setEntryName("DamageDrivingForce");
+
+@StateVariable real de;
+de.setGlossaryName("Damage");
+
+@ProvidesSymmetricTangentOperator ;
+@Integrator{
+  constexpr const strain emin = 1.e-12;
+  // positive part
+  const auto f  = [](const real x){return x>0 ? x : 0;};
+  // derivative of the positive part
+  const auto df = [&emin](const real x)
+    {return std::abs(x)<emin ? 0.5 : ((x<0) ? 0 : 1);};
+  // update the damage
+  const auto d = de+dde;
+  // lame coefficients
+  const auto lambda = computeLambda(yg,nu);
+  const auto mu     = computeMu(yg,nu);
+  // computation of the stress, positive energy density and consistent
+  // tangent operator
+  const StrainStensor e_ = eto + deto;
+  const auto fdf     = e_.computeIsotropicFunctionAndDerivative(f,df,emin*0.1);
+  const auto& ep     = fdf.first;  // positive part of e_
+  const auto& dep_de = fdf.second; // derivative of the positive part of e_
+  const StrainStensor en = e_-ep;  // negative part of e_
+  // energy density
+  const strain tr  = trace(e_);
+  const strain tr_p = max(tr,real(0));
+  const strain tr_n = min(tr,real(0));
+  H = max(H,(((lambda/2)*(tr_p)*(tr_p))+(mu*(ep|ep))));
+  // stress
+  const auto deff = ((1-d)*(1-d))+kk;
+  sig = 2*mu*(deff*ep+en)+lambda*(deff*tr_p+tr_n)*StrainStensor::Id();
+  // consistent tangent operator (secant one here)
+  if(computeTangentOperator_){
+    if(tr>=0){
+      Dt=deff*(lambda*Stensor4::IxI()+2*mu*dep_de)+(2*mu*(Stensor4::Id()-dep_de));
+    } else {
+      Dt=deff*2*mu*dep_de+(lambda*Stensor4::IxI()+2*mu*(Stensor4::Id()-dep_de));
+    }
+  }
+  static_cast<void>(smt);
+} // end of @Integrator
+
+@APosterioriTimeStepScalingFactor {
+  // pas d'acroissement de l'endommagement sur le pas de temps
+  if (dde<1.e-4){
+    return {true,1.e-2/(max(dde,1e-4))};
+  }
+  return {true,1};
+}

mfront_wrapper/nonlinear_problem.py

@@ -135,31 +135,6 @@ class MFrontNonlinearProblem(NonlinearProblem):
         # tangent matrix quadrature space
         self.WCt = FunctionSpace(self.mesh, Wce)
 
-    # def strain_measure(self, v):
-    #     """ Strain measure associated with stress measure:
-    #         * small strain behaviour: linearized strain tensor epsilon = sym(grad(u))
-    #         * finite strain behaviour: transformation gradient F = Id + grad(u)
-    #     """
-    #     if self.axisymmetric:
-    #         r = abs(SpatialCoordinate(self.mesh)[0])
-    #         g = axi_grad(r, v)
-    #         E = symmetric_tensor_to_vector(sym(g))
-    #         if v.geometric_dimension()==2:
-    #             return as_vector([E[i] for i in range(4)])
-    #     else:
-    #         g = grad(v)
-    #         if self.finite_strain:
-    #             return transformation_gradient(g, dim=v.geometric_dimension())
-    #         else:
-    #             return symmetric_gradient(g)
-
-    # def strain_variation(self, v):
-    #     """ Variation of strain measure associated with stress measure:
-    #         * small strain behaviour: linearized strain tensor d_epsilon = sym(grad(du))
-    #         * finite strain behaviour: displacement gradient d_F = grad(du)
-    #     """
-    #     return derivative(self.strain_measure(v), v, v)
-
     def initialize_fields(self):
         for (g, size) in zip(self.material.get_gradient_names(), self.material.get_gradient_sizes()):
             try: