-
Notifications
You must be signed in to change notification settings - Fork 14
3. Creating material models in Florence
In this section, we will have a look at how to define material models in Florence. One of the core features of Florence is in its large amount of predefined material models for different types of variational formulations, as all analyses in Florence require defining a material model, even solving a Poisson equation. These set of material models are extremely high performant as they are implemented on the fly using the SIMD based tensor contraction framework Fastor. Fastor can generate the Hessian/tangent operators and gradients of different types of materials by employing FLOP reduction algorithm to determine the least amount of FLOP count for a given constitutive law and then produce data parallel SIMD code for them targeting all types of modern heterogeneous architectures.
To define an anisotropic material for Poisson type problems, i.e. anisotropic Poisson equation we can for instance do
material = AnisotropicIdealDielectric(ndim, e=np.eye(ndim))All material models need to know about the dimension of the problem that you are solving as they need to generate the right tangent operator/Hessian and work/energy conjugates. The parameter e in this case is the anisotropic matrix defining the the material constants for instance the anisotropic permittivity tensor in case of electrostatics. For standard Possion type problems the material model can be defined as
material = IdealDielectric(ndim, eps=1.)The linear elastic material model can be defined as
material = LinearElastic(ndim, youngs_modulus=1e9, poisson_ratio=0.32)or alternatively, using Lame constants
material = LinearElastic(ndim, mu=5e8, lamb=1e9)Note that, Florence strictly disallows using LinearElastic model with multiple load increments for quasi-static analyses, as solving a linear elastic problem in multiple step can just be done using linear superposition. In cases where there are prestresses in the material or the residual changes due to other external reasons, the equivalent IncremrentalLinearElastic model can be used instead, which behaves just like LinearElastic model but allows prestressed/non-zero residual state and also supports load increments.
material = IncremrentalLinearElastic(ndim, youngs_modulus=1e9, poisson_ratio=0.32)Transversely isotropic and anisotropic linear elastic models are also available
material = TranservselyIsotropicLinearElastic(ndim,
youngs_modulus=1e9,
poisson_ratio=0.32,
E_A=1e5, # out of plane Young's modulus
G_A=5e4, # out of plane shear modulus
)A suite of nonlinear hyperelastic material models are available
material = NeoHookean(ndim, mu=5e8, lamb=1e9) # or
material = NeoHookean(ndim, youngs_modulus=1e9, poisson_ratio=0.4)
material = RegularisedNeoHookean(ndim, lamb=1e9, mu=5e8)
material = CoerciveNeoHookean(ndim, lamb=1e9, mu=5e8)
material = MooneyRivlin(ndim, mu1=1e8, mu2=1e8, lamb=1e9)Generic anisotropic hyperelastic material models are available as well
material = TranservselyIsotropicHyperElastic(ndim,
youngs_modulus=1e9,
poisson_ratio=0.32,
E_A=1e5, # out of plane Young's modulus
anisotropic_orientations=numpy.random.rand(mesh.nelem,ndim)
)where anisotropic_orientations is the vector of anisotropy that is defined for every element of the computational mesh
A suite of nonlinear electro-hyperelastic material models are available
material = IsotropicElectroMechanics_101(ndim, mu=5e8, lamb=1e9, eps_1=2e-11) # NeoHookean based
material = IsotropicElectroMechanics_108(ndim, mu1=5e8, mu2=5e8, lamb=1e9, eps_2=2e-11) # MooneyRivlin based
material = IsotropicElectroMechanics_107(ndim, mu1=5e8, mu2=5e8, mue=1e7, lamb=1e9,
eps_1=1e-11, eps_2=2e-11, eps_e=1e-13) Generic anisotropic/piezoelectric hyperelastic material models are available as well
material = Piezoelectric_100(ndim, mu1=5e8, mu2=5e8, mu3=1e7, lamb=1e9,
eps_1=1e-11, eps_2=2e-11, eps_3=1e-13,
anisotropic_orientations=numpy.random.rand(mesh.nelem,ndim))where anisotropic_orientations is the vector of anisotropy that is defined for every element of the computational mesh
To check the corresponding energy of any material or constitutive law, in the ipython shell you can do:
In [0]: IsotropicElectroMechanics_107?
Out[0]:
"""Electromechanical model in terms of internal energy
W(C,D0) = W_mn(C) + 1/2/eps_1 (D0*D0) + 1/2/eps_2/J (FD0*FD0) + W_reg(C,D0)
W_mn(C) = u1*C:I+u2*G:I - 2*(u1+2*u2+6*mue)*lnJ + lamb/2*(J-1)**2 + mue*(C:I)**2
W_reg(C,D_0) = 2/eps_e/mue (F:F)(FD0*FD0) + 1/mu/eps_e**2*(FD0*FD0)**2
"""which explains all the material constants used and their energy form.
Support for user defined materials is limited at the moment, but will be incorporated into Florence soon, as this is a high priority. In essence, Florence, using Fastor backend should be able to generate optimised low level code for the tangent operator of any given energy and consistently linearise any complex expression on the fly. The user only has to supply the energy form and related parameters as a dictionary, like so:
user_material_context = Material()
user_material_context.Linearise(
dict("name":"MooneyRivlinType",
"energy_expression":"u_1*(F:F-3) + u_2*(H:H-3) - 2*(u_1+2*u_2)*ln(J)+lamb/2*(J-1)**2",
"fields":"mechanics",
"nature":"nonlinear",
"constants":"u_1,u_2,lamb",
"variables":dict("F":"deformation_gradient","H":"cofactor(F)","J":"determinant(F)")
)
)
user_material_context.Generate()This will generate a low level code for material model called MooneyRivlinType which can be used like any other material model in Florence
user_material = MooneyRivlinType(ndim, u_1=1, u_2=1, lamb=10)