Hi everyone, Hi interested,
I trying to build some example to experiment and learn how to approach problems with multiple domains.
Therefore, I would like to build a stackup of different materials with different behaviours.
This is done with the following code:
Define the material stackup
from pymaterial.materials import IsotropicMaterial, TransverselyIsotropicMaterial, OrthotropicMaterial
mat0 = {
"name": "steal",
"thickness": 0.3,
"stiff_tens": IsotropicMaterial(Em=2e11,nu=0.3, density=7850).get_stiffness() # from ANSYS
}
mat2 = {
"name": "E-glass",
"thickness": 0.3,
"stiff_tens": OrthotropicMaterial(E_x=4.5e10,
E_y=1e10,
E_z=1e10,
nu_xy=0.3,
nu_yz=0.4,
nu_xz=0.3,
G_xy=5e9,
G_yz=3.8462e9,
G_xz=5e9,
density=2000
).get_stiffness() # from ANSYS Epoxy E-Glass UD
}
mat3 = {
"name": "carbon fibre",
"thickness": 0.3,
"stiff_tens": OrthotropicMaterial(E_x=2.3e11,
E_y=2.3e10,
E_z=2.3e10,
nu_xy=0.2,
nu_yz=0.4,
nu_xz=0.2,
G_xy=9e9,
G_yz=8.2143e9,
G_xz=9e9,
density=1800
).get_stiffness() # from ANSYS Carbon Fiber (230 GPa)
}
stack = [
mat1,
mat2,
mat3
]
Build the model and load it into fenics
import numpy as np
import gmsh
from mpi4py import MPI
gmsh.initialize()
model = gmsh.model()
name = "Material Stack"
model.add(name)
model.setCurrent(name)
width, length = 1, 1
thickness = 0
gdim = 3 # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:
boxes = []
coms = []
for i in range(len(stack)):
elem = stack[i]
elem_thickness = elem["thickness"]
start_height = thickness
end_height = thickness + elem_thickness
boxes.append(model.occ.addBox(0, 0, start_height, width, length, elem_thickness))
thickness = end_height
coms.append([width/2, length/2, start_height + elem_thickness/2])
model.occ.synchronize()
for i in range(len(boxes)-1):
model.occ.fragment([(3,boxes[i])], [(3,boxes[i+1])])
for dim_,id_ in model.getEntities(3):
com = gmsh.model.occ.getCenterOfMass(dim_, id_)
# print(com)
for i in range(len(coms)):
if not np.allclose(com, coms[i]):
continue
tag = i
model.addPhysicalGroup(3, [id_], tag=tag)
model.setPhysicalName(3, tag, stack[i]["name"])
break
model.mesh.generate(dim=3)
# gmsh.model.mesh.optimize("Netgen")
# Create DOLFINx distributed mesh
from dolfinx.io.gmshio import model_to_mesh
mesh, cell_markers, facet_markers = model_to_mesh(model, mesh_comm, model_rank)
gmsh.finalize()
Defining the problem
In the particular case one approach would be to define a Tensorfunction which represents the stiffness tensor C.
Like:
import ufl
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def strain2voigt(e):
"""e is a 2nd-order tensor, returns its Voigt vectorial representation"""
return ufl.as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s):
"""
s is a stress-like vector (no 2 factor on last component)
returns its tensorial representation
"""
return ufl.as_tensor([[s[0], s[2]],
[s[2], s[1]]])
def sigma(u):
return voigt2stress(ufl.dot(stiff_tens, strain2voigt(epsilon(u))))
and
from dolfinx.fem import VectorFunctionSpace, TrialFunction, TestFunction
import ufl
V = VectorFunctionSpace(mesh, ('Lagrange', 1))
# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
u = Function(V, name='Displacement')
a = ufl.inner(sigma(u), epsilon(v))* ufl.dx
The problem: generating the spacial dependent tensor function
But whatever I tried to build such a tensor function it did not work:
from dolfinx.fem import TensorFunctionSpace, Function, Constant
import ufl
from petsc4py.PETSc import ScalarType
Q = TensorFunctionSpace(mesh, ("CG", 2)) # CG-Lagrange element, 2-degree of freedom (2nd order element)
stiff_tens=Function(Q)# often symbolized by "C"
material_tags = np.unique(cell_markers.values)
for tag in material_tags:
cells = cell_markers.find(tag)
if tag < len(stack):
material = stack[i]
stiff_tens.x.array[cells]=np.array(list([material["stiff_tens"]])*len(cells), dtype=ScalarType)
which throws the following error:
ValueError Traceback (most recent call last)
Cell In[61], line 12
10 if tag < len(stack):
11 material = stack[i]
---> 12 stiff_tens.x.array[cells]=np.array(list([material["stiff_tens"]])*len(cells), dtype=ScalarType)
ValueError: shape mismatch: value array of shape (402,6,6) could not be broadcast to indexing result of shape (402,)
Another question:
dolphin had something called Subdomains, where I could also define other partial differential equations on. For me this means that I could perform the integration with one stiffness tensor over each material domain separetly and then add all matrices together. How would this look in dolphinx or is this not possible?
Thanks in advance!