Commit 5472d67f by Jeremy BLEYER

### Validated multiphase model

parent 32685cbd
 ... ... @@ -10,9 +10,10 @@ import mfront_wrapper as mf import numpy as np from ufl import diag width = 1. width = 0.5 height = 1. mesh = RectangleMesh(Point(0., 0.), Point(width, height), 10, 10) mesh = RectangleMesh(Point(0., 0.), Point(width, height), 100, 100) # mesh = BoxMesh(Point(0., 0., 0.), Point(width, height, height), 100, 2, 2) Ve = VectorElement("CG", mesh.ufl_cell(), 1) V = FunctionSpace(mesh, MixedElement([Ve, Ve])) ... ... @@ -37,13 +38,12 @@ bc = [DirichletBC(V.sub(0).sub(0), Constant(0), left), facets = MeshFunction("size_t", mesh, 1) ds = Measure("ds", subdomain_data=facets) mat_prop = {"MatrixYoungModulus": 70e3, "MatrixPoissonRatio": 0.2, "FiberYoungModulus": 200e3, mat_prop = {"MatrixYoungModulus": 10., "MatrixPoissonRatio": 0.45, "FiberYoungModulus": 10000, "FiberPoissonRatio": 0.3, "FiberVolumeFraction": 0.1, "Size": 0.1} "FiberVolumeFraction": 0.01, "Size": 1/16.} material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so", "MultiphaseModel", ... ... @@ -53,16 +53,40 @@ material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so", problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0) problem.bc = bc problem.register_gradient("MatrixStrain", sym(grad(u1))) problem.register_gradient("InclusionStrain", sym(grad(u2))) problem.register_gradient("FiberStrain", sym(grad(u2))) problem.register_gradient("RelativeDisplacement", diag(u2-u1)) # E = Function(problem.W0) # E.interpolate(Expression("20e3", degree=1)) # mat_prop["MatrixYoungModulus"] = E # problem.material.update_material_properties(mat_prop) problem.solve(u.vector()) x = np.linspace(width/2, width, 21) x = np.linspace(0, width, 21) import matplotlib.pyplot as plt plt.figure() 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"])) plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "ob", label=r"\$u_x^\textrm{matrix}\$ EF") plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "or", label=r"\$u_x^\textrm{fiber}\$ EF") E1 = mat_prop["MatrixYoungModulus"] nu1 = mat_prop["MatrixPoissonRatio"] E2 = mat_prop["FiberYoungModulus"] nu2 = mat_prop["FiberPoissonRatio"] eta = mat_prop["FiberVolumeFraction"] s = mat_prop["Size"] mu1 = E1/2/(1+nu1) mu2 = E2/2/(1+nu2) lmbda1 = E1*nu1/(1+nu1)/(1-2*nu1) lmbda2 = E2*nu2/(1+nu2)/(1-2*nu2) kappa = 12*mu1*mu2/((1-eta)*mu2+eta*mu1)/s**2 M1 = lmbda1+2*mu1 M2 = lmbda2+2*mu2 Er = eta*E2/(1-nu2**2) ahom = M1+Er l = sqrt(M1*Er/ahom/kappa) um = lmbda1/ahom*(x+l*(Er/M1)*np.sinh(x/l)/np.cosh(width/l)) ur = lmbda1/ahom*(x-l*np.sinh(x/l)/np.cosh(width/l)) plt.plot(x, um, '--b', label=r"\$u_x^\textrm{matrix}\$ ref.") plt.plot(x, ur, '--r', label=r"\$u_x^\textrm{fiber}\$ ref.") plt.legend() plt.show() \ No newline at end of file
 ... ... @@ -46,7 +46,6 @@ 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 ... ...
 @DSL DefaultGenericBehaviour; @Behaviour MultiphaseModel; @ModellingHypotheses {PlaneStrain, Tridimensional}; @ModellingHypotheses {PlaneStrain}; @Gradient StrainStensor e₁; e₁.setEntryName("MatrixStrain"); ... ... @@ -33,8 +33,8 @@ nu1.setEntryName("MatrixPoissonRatio"); Y2.setEntryName("FiberYoungModulus"); @MaterialProperty real nu2; nu2.setEntryName("FiberPoissonRatio"); @MaterialProperty real rho; rho.setEntryName("FiberVolumeFraction"); @MaterialProperty real ρ; ρ.setEntryName("FiberVolumeFraction"); @MaterialProperty real s; s.setEntryName("Size"); ... ... @@ -42,32 +42,42 @@ s.setEntryName("Size"); @Integrator { // remove useless warnings, as we always compute the tangent operator static_cast(computeTangentOperator_); const Stensor4 Cₕ = { 1, 2, 3, 4, 5, 6, // 1ère ligne 5, 6, 7, 8, 9, 10, // 2ème ligne 1, 2, 3, 4, 11, 12, // 3ème ligne 1, 2, 3, 4, 13, 14, // 4ème ligne 1, 2, 3, 4, 13, 14, // 4ème ligne 1, 2, 3, 4, 13, 14 // 4ème ligne }; const auto λ₁ = computeLambda(Y1,nu1); const auto μ₁ = computeMu(Y1,nu1); const auto λ₂ = rho*computeLambda(Y2,nu2); const auto μ₂ = rho*computeMu(Y2,nu2); constexpr const auto λ₁₂ = 0.; constexpr const auto μ₁₂ = 0.; const auto kappa = μ₁/12/(1-rho)/pow(s,2); const auto C1 = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄; const auto C2 = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄; const auto idC = invert(C2-C1); const auto avgC = rho*C1+(1-rho)*C2; const auto D11 = (1-rho)*C1-(C1*idC*(Cₕ-avgC)*idC*C1); ∂σ₁∕∂Δe₁ = C1; ∂σ₁∕∂Δe₂ = λ₁₂ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁₂ ⋅ I₄; ∂σ₂∕∂Δe₁ = ∂σ₁∕∂Δe₂; ∂σ₂∕∂Δe₂ = C2; const auto λ₂ = computeLambda(Y2,nu2); const auto μ₂ = computeMu(Y2,nu2); const auto Eₒₑ¹ = λ₁+2*μ₁; const auto Eₒₑ² = λ₂+2*μ₂; const auto Eₒₑ = (1-ρ)*Eₒₑ¹ + ρ*Eₒₑ²; const auto iEₒₑ = (1-ρ)/Eₒₑ¹ + ρ/Eₒₑ²; const auto λ = (1-ρ)*λ₁ + ρ*λ₂; const auto λiEₒₑ = (1-ρ)*λ₁/Eₒₑ¹ + ρ*λ₂/Eₒₑ²; const auto λ2iEₒₑ = (1-ρ)*pow(λ₁,2)/Eₒₑ¹ + ρ*pow(λ₂,2)/Eₒₑ²; const auto iμ = (1-ρ)/μ₁ + ρ/μ₂; const auto κ = 12/iμ/pow(s,2); // interaction stiffness const auto C₁₁₁₁ʰᵒᵐ = Eₒₑ - λ2iEₒₑ + pow(λiEₒₑ,2)/iEₒₑ; const auto C₁₁₂₂ʰᵒᵐ = λiEₒₑ/iEₒₑ; const Stensor4 Cʰᵒᵐ = { C₁₁₁₁ʰᵒᵐ , C₁₁₂₂ʰᵒᵐ, C₁₁₂₂ʰᵒᵐ, 0., // C₁₁₂₂ʰᵒᵐ, 1/iEₒₑ, λ, 0., // C₁₁₂₂ʰᵒᵐ, λ, Eₒₑ, 0., // 0., 0., 0., 2/iμ }; const auto C¹ = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄; const auto C² = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄; const auto iΔC = invert(C²-C¹); const auto Cᵃᵛᵍ = (1-ρ)⋅C¹ + ρ⋅C²; const auto H = iΔC ⋅ (Cᵃᵛᵍ-Cʰᵒᵐ) ⋅ iΔC; const auto D¹¹ = (1-ρ)⋅C¹-(C¹ ⋅ H ⋅ C¹); const auto D²² = ρ⋅C²-(C² ⋅ H ⋅ C²); const auto D¹² = C¹ ⋅ H ⋅ C²; ∂σ₁∕∂Δe₁ = D¹¹; ∂σ₁∕∂Δe₂ = D¹²; ∂σ₂∕∂Δe₁ = transpose(D¹²); ∂σ₂∕∂Δe₂ = D²²; σ₁ = ∂σ₁∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₁∕∂Δe₂ ⋅ (e₂ + Δe₂); σ₂ = ∂σ₂∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₂∕∂Δe₂ ⋅ (e₂ + Δe₂); I = kappa*(V + ΔV); ∂I∕∂ΔV = kappa ⋅ I₄; I = κ*(V + ΔV); ∂I∕∂ΔV = κ ⋅ I₄; }
 @DSL DefaultGenericBehaviour; @Behaviour MultiphaseModel; @ModellingHypotheses {Tridimensional}; @Gradient StrainStensor e₁; e₁.setEntryName("MatrixStrain"); @Flux StressStensor σ₁; σ₁.setEntryName("MatrixStress"); @Gradient StrainStensor e₂; e₂.setEntryName("FiberStrain"); @Flux StressStensor σ₂; σ₂.setEntryName("FiberStress"); @Gradient StrainStensor V; V.setEntryName("RelativeDisplacement"); @Flux StressStensor I; I.setEntryName("InteractionForce"); //@TangentOperatorBlocks {∂σ₁∕∂Δe₁,∂σ₁∕∂Δe₂,∂σ₂∕∂Δe₁,∂σ₂∕∂Δe₂,∂I∕∂ΔV}; @TangentOperatorBlock ∂σ₁∕∂Δe₁; @AdditionalTangentOperatorBlock ∂σ₁∕∂Δe₂; @AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₁; @AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₂; @AdditionalTangentOperatorBlock ∂I∕∂ΔV; @MaterialProperty stress Y1; Y1.setEntryName("MatrixYoungModulus"); @MaterialProperty real nu1; nu1.setEntryName("MatrixPoissonRatio"); @MaterialProperty stress Y2; Y2.setEntryName("FiberYoungModulus"); @MaterialProperty real nu2; nu2.setEntryName("FiberPoissonRatio"); @MaterialProperty real ρ; ρ.setEntryName("FiberVolumeFraction"); @MaterialProperty real s; s.setEntryName("Size"); @ProvidesTangentOperator; @Integrator { // remove useless warnings, as we always compute the tangent operator static_cast(computeTangentOperator_); const auto λ₁ = computeLambda(Y1,nu1); const auto μ₁ = computeMu(Y1,nu1); const auto λ₂ = computeLambda(Y2,nu2); const auto μ₂ = computeMu(Y2,nu2); const auto Eₒₑ¹ = λ₁+2*μ₁; const auto Eₒₑ² = λ₂+2*μ₂; const auto Eₒₑ = (1-ρ)*Eₒₑ¹ + ρ*Eₒₑ²; const auto iEₒₑ = (1-ρ)/Eₒₑ¹ + ρ/Eₒₑ²; const auto λ = (1-ρ)*λ₁ + ρ*λ₂; const auto λiEₒₑ = (1-ρ)*λ₁/Eₒₑ¹ + ρ*λ₂/Eₒₑ²; const auto λ2iEₒₑ = (1-ρ)*pow(λ₁,2)/Eₒₑ¹ + ρ*pow(λ₂,2)/Eₒₑ²; const auto iμ = (1-ρ)/μ₁ + ρ/μ₂; const auto μ = (1-ρ)*μ₁ + ρ*μ₂; const auto κ = 12/iμ/pow(s,2); // interaction stiffness const auto C₁₁₁₁ʰᵒᵐ = Eₒₑ - λ2iEₒₑ + pow(λiEₒₑ,2)/iEₒₑ; const auto C₁₁₂₂ʰᵒᵐ = λiEₒₑ/iEₒₑ; const Stensor4 Cʰᵒᵐ = { C₁₁₁₁ʰᵒᵐ , C₁₁₂₂ʰᵒᵐ, C₁₁₂₂ʰᵒᵐ, 0., 0., 0., // C₁₁₂₂ʰᵒᵐ, 1/iEₒₑ, λ, 0., 0., 0., // C₁₁₂₂ʰᵒᵐ, λ, Eₒₑ, 0., 0., 0., // 0., 0., 0., 2/iμ, 0., 0., // 0., 0., 0., 0., 2/iμ, 0., // 0., 0., 0., 0., 0., 2*μ // }; const auto C¹ = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄; const auto C² = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄; const auto iΔC = invert(C²-C¹); const auto Cᵃᵛᵍ = (1-ρ)*C¹ + ρ*C²; const auto H = iΔC ⋅ (Cᵃᵛᵍ-Cʰᵒᵐ) ⋅ iΔC; const auto D¹¹ = (1-ρ)*C¹-(C¹ ⋅ H ⋅ C¹); const auto D²² = ρ*C²-(C² ⋅ H ⋅ C²); const auto D¹² = C¹ ⋅ H ⋅ C²; const auto D¹²ᵀ = C² ⋅ H ⋅ C¹; ∂σ₁∕∂Δe₁ = D¹¹; ∂σ₁∕∂Δe₂ = D¹²; ∂σ₂∕∂Δe₁ = D¹²ᵀ; //transpose(D¹²); ∂σ₂∕∂Δe₂ = D²²; σ₁ = ∂σ₁∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₁∕∂Δe₂ ⋅ (e₂ + Δe₂); σ₂ = ∂σ₂∕∂Δe₁ ⋅ (e₁ + Δe₁) + ∂σ₂∕∂Δe₂ ⋅ (e₂ + Δe₂); I = κ*(V + ΔV); ∂I∕∂ΔV = κ ⋅ I₄; }
 ... ... @@ -12,7 +12,7 @@ INCLUDES := -I../include \ CXXFLAGS := -Wall -Wfatal-errors -ansi \$(shell tfel-config --oflags) -fPIC \$(INCLUDES) SRCCXX = IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx MultiphaseModel-generic.cxx MultiphaseModel.cxx PhaseField-generic.cxx PhaseField.cxx PhaseFieldCoupled-generic.cxx PhaseFieldCoupled.cxx StationaryHeatTransfer-generic.cxx StationaryHeatTransfer.cxx TensionCompressionSplit-generic.cxx TensionCompressionSplit.cxx SRCCXX = IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx MultiphaseModel-generic.cxx MultiphaseModel.cxx MultiphaseModel3D-generic.cxx MultiphaseModel3D.cxx PhaseField-generic.cxx PhaseField.cxx PhaseFieldCoupled-generic.cxx PhaseFieldCoupled.cxx StationaryHeatTransfer-generic.cxx StationaryHeatTransfer.cxx TensionCompressionSplit-generic.cxx TensionCompressionSplit.cxx makefiles1 = \$(SRCCXX:.cxx=.d) makefiles2 = \$(makefiles1:.cpp=.d) ... ... @@ -22,7 +22,7 @@ makefiles = \$(makefiles2) all : libBehaviour.so libBehaviour.so : MultiphaseModel-generic.o MultiphaseModel.o PhaseField-generic.o PhaseField.o PhaseFieldCoupled-generic.o PhaseFieldCoupled.o IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o StationaryHeatTransfer-generic.o StationaryHeatTransfer.o TensionCompressionSplit-generic.o TensionCompressionSplit.o libBehaviour.so : MultiphaseModel-generic.o MultiphaseModel.o PhaseField-generic.o PhaseField.o PhaseFieldCoupled-generic.o PhaseFieldCoupled.o IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o StationaryHeatTransfer-generic.o StationaryHeatTransfer.o TensionCompressionSplit-generic.o TensionCompressionSplit.o MultiphaseModel3D-generic.o MultiphaseModel3D.o @\$(CXX) -shared \$^ -o \$@ -L"\$(strip \$(shell tfel-config --library-path))" \$(patsubst %,-l%,\$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants)) install : ... ...
 ... ... @@ -19,7 +19,9 @@ sources : { "StationaryHeatTransfer-generic.cxx", "StationaryHeatTransfer.cxx", "TensionCompressionSplit-generic.cxx", "TensionCompressionSplit.cxx" "TensionCompressionSplit.cxx", "MultiphaseModel3D-generic.cxx", "MultiphaseModel3D.cxx" }; cppflags : { "\$(shell tfel-config --cppflags --compiler-flags)" ... ... @@ -65,7 +67,8 @@ epts : { "TensionCompressionSplit_Axisymmetrical", "TensionCompressionSplit_PlaneStrain", "TensionCompressionSplit_GeneralisedPlaneStrain", "TensionCompressionSplit_Tridimensional" "TensionCompressionSplit_Tridimensional", "MultiphaseModel3D_Tridimensional" }; }; headers : { ... ... @@ -96,6 +99,10 @@ headers : { "MFront/GenericBehaviour/TensionCompressionSplit-generic.hxx", "TFEL/Material/TensionCompressionSplit.hxx", "TFEL/Material/TensionCompressionSplitBehaviourData.hxx", "TFEL/Material/TensionCompressionSplitIntegrationData.hxx" "TFEL/Material/TensionCompressionSplitIntegrationData.hxx", "MFront/GenericBehaviour/MultiphaseModel3D-generic.hxx", "TFEL/Material/MultiphaseModel3D.hxx", "TFEL/Material/MultiphaseModel3DBehaviourData.hxx", "TFEL/Material/MultiphaseModel3DIntegrationData.hxx" }; };
 import mgis.behaviour as mgis_bv import dolfin mgis_hypothesis = {"plane_strain": mgis_bv.Hypothesis.PlaneStrain, "plane_stress": mgis_bv.Hypothesis.PlaneStress, ... ... @@ -13,15 +14,21 @@ class MFrontNonlinearMaterial: self.name = name # Defining the modelling hypothesis self.hypothesis = mgis_hypothesis[hypothesis] if self.hypothesis == mgis_bv.Hypothesis.Tridimensional: self._tensorsize = 9 else: self._tensorsize = 5 # if self.hypothesis == mgis_bv.Hypothesis.Tridimensional: # self._tensorsize = 9 # else: # self._tensorsize = 5 self.material_properties = material_properties self.external_state_variables = external_state_variables # finite strain options bopts = mgis_bv.FiniteStrainBehaviourOptions() bopts.stress_measure = mgis_bv.FiniteStrainBehaviourOptionsStressMeasure.PK1 bopts.tangent_operator = mgis_bv.FiniteStrainBehaviourOptionsTangentOperator.DPK1_DF # Loading the behaviour self.behaviour = mgis_bv.load(self.path, self.name, self.hypothesis) self.finite_strain = self.behaviour.getBehaviourType()=="StandardFiniteStrainBehaviour" if self.finite_strain: self.behaviour = mgis_bv.load(bopts, self.path, self.name, self.hypothesis) def set_data_manager(self, ngauss): # Setting the material data manager ... ... @@ -34,7 +41,13 @@ class MFrontNonlinearMaterial: self.material_properties = material_properties for s in [self.data_manager.s0, self.data_manager.s1]: for (key, value) in self.material_properties.items(): mgis_bv.setMaterialProperty(s, key, value) if type(value) in [int, float]: mgis_bv.setMaterialProperty(s, key, value) else: if isinstance(value, dolfin.Function): value = value.vector().get_local() mgis_bv.setMaterialProperty(s, key, value, mgis_bv.MaterialStateManagerStorageMode.LocalStorage) def update_external_state_variable(self, external_state_variables=None): if external_state_variables is not None: ... ... @@ -54,12 +67,9 @@ class MFrontNonlinearMaterial: def get_gradient_sizes(self): return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.gradients] def get_flux_sizes(self): return [self._tensorsize if (svar.name == "Stress" and self.finite_strain) else \ mgis_bv.getVariableSize(svar, self.hypothesis) \ for svar in self.behaviour.thermodynamic_forces] return [mgis_bv.getVariableSize(svar, self.hypothesis) for svar in self.behaviour.thermodynamic_forces] def get_tangent_block_names(self): return [(t[0].name, t[1].name) for t in self.behaviour.tangent_operator_blocks] def get_tangent_block_sizes(self): return [tuple([self._tensorsize if (tt.name == "Stress" and self.finite_strain) \ else mgis_bv.getVariableSize(tt, self.hypothesis) for tt in t]) \ return [tuple([mgis_bv.getVariableSize(tt, self.hypothesis) for tt in t]) \ for t in self.behaviour.tangent_operator_blocks] \ No newline at end of file
 ... ... @@ -200,13 +200,13 @@ class MFrontNonlinearProblem(NonlinearProblem): f, g = block t = self.fluxes[f].tangent_blocks[g] block_shape = self.flattened_block_shapes[i] if self.material.finite_strain and f == "Stress": Ctv = t.vector().get_local() mgis_bv.convertFiniteStrainTangentOperator(Ctv, self.material.data_manager, mgis_bv.FiniteStrainTangentOperator.DPK1_DF) t.vector().set_local(Ctv) else: t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten()) # if self.material.finite_strain and f == "Stress": # Ctv = t.vector().get_local() # # mgis_bv.convertFiniteStrainTangentOperator(Ctv, self.material.data_manager, # # mgis_bv.FiniteStrainTangentOperator.DPK1_DF) # t.vector().set_local(Ctv) # else: t.vector().set_local(self.material.data_manager.K[:,buff:buff+block_shape].flatten()) buff += block_shape def update_fluxes(self): ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment