#
=========================================================================
#
FEniCS
implementation:
Homogenized
model
of
a
cube
torsion
test
with
a
3D
#
octet-truss
microstructure
#
J.
Ulloa,
M.P.
Ariza,
J.E.
Andrade
&
M.
Ortiz
#
Homogenized
models
of
mechanical
metamaterials
#
Computer
Methods
in
Applied
Mechanics
and
Engineering,
August
2024
#
=========================================================================
from
dolfin
import
*
from
mshr
import
*
import
numpy
as
np
from
aux
import
rotate_elast
import
mshr
from
mshr
import
*
#
-------------------------------------------------------------------------
#
Create
mesh
#
-------------------------------------------------------------------------
hsize
=
0
.
12
#
Finite
element
size
[
mm]
cube
=
mshr
.
Box
(
Point
(
0
,
0
,
0
)
,
Point
(
3
,
3
,
3
))
mesh
=
mshr
.
generate_mesh
(cube
,
3
/
hsize)
#
-------------------------------------------------------------------------
#
Input
parameters
#
-------------------------------------------------------------------------
#
Metamaterial
parameters
Lc
=
0
.
75
#
Cell
length
[
mm]
Eb
=
430
.
#
Bar
Young's
modulus
[
N
/
mm^2]
rb
=
0
.
5
*
0
.
065
#
Bar
radius
[
mm]
Ab
=
pi
*
rb
**
2
#
Bar
area
[
mm^2]
Ib
=
0
.
25
*
pi
*
rb
**
4
#
Bar
MoM
[
mm^4]
#
Homogenized
properties
C11
=
(
4
.
*
Eb
*
Ab
*
Lc
**
2
+
24
.
*
Eb
*
Ib)
/
(
sqrt
(
2
.
)
*
Lc
**
4
)
C12
=
sqrt
(
2
.
)
*
(Eb
*
Ab
*
Lc
**
2
-
6
*
Eb
*
Ib)
/
(Lc
**
4
)
C44
=
(
2
.
*
Eb
*
Ab
*
Lc
**
2
+
12
.
*
Eb
*
Ib)
/
(
sqrt
(
2
.
)
*
Lc
**
4
)
kappa
=
48
.
*
sqrt
(
2
.
)
*
Eb
*
Ib
/
(Lc
**
4
)
#
Cubic
lattice
rotation
w.r.t.
(
x1,x2,x3)
ang
=
-
pi
/
4
o1
=
0
;
o2
=
0
;
o3
=
1
#
Applied
load
ubar
=
2
.
*
pi
/
180
.
#
Imposed
rotation
angle
at
the
top
[
rad]
#
-------------------------------------------------------------------------
#
Function
spaces
#
-------------------------------------------------------------------------
v_element
=
VectorElement
(
'
CG
'
,
mesh
.
ufl_cell
()
,
2
)
#
Deflections
t_element
=
VectorElement
(
'
CG
'
,
mesh
.
ufl_cell
()
,
1
)
#
Rotations
U
=
FunctionSpace
(mesh
,
v_element
*
t_element)
V
,
T
=
U
.
split
()
V_0
,
V_1
,
V_2
=
V
.
split
()
#
3
DOFs
for
deflections
T_0
,
T_1
,
T_2
=
T
.
split
()
#
3
DOFs
for
rotations
#
-------------------------------------------------------------------------
#
Lists,
DOFs
and
BCs
#
-------------------------------------------------------------------------
coords
=
mesh
.
coordinates
()
xmin
=
min
(coords[
:
,
0
])
xmax
=
max
(coords[
:
,
0
])
ymin
=
min
(coords[
:
,
1
])
ymax
=
max
(coords[
:
,
1
])
zmin
=
min
(coords[
:
,
2
])
zmax
=
max
(coords[
:
,
2
])
top
=
CompiledSubDomain
(
"
near
(
x
[
2],
zmax)
&&
on_boundary
"
,
zmax
=
zmax)
bottom
=
CompiledSubDomain
(
"
near
(
x
[
2],
zmin)
&&
on_boundary
"
,
zmin
=
zmin)
#
Dirichlet
Boundary
BC_bottom_0
=
DirichletBC
(V_0
,
Constant
(
0
.
0
)
,
bottom)
BC_bottom_1
=
DirichletBC
(V_1
,
Constant
(
0
.
0
)
,
bottom)
BC_bottom_2
=
DirichletBC
(V_2
,
Constant
(
0
.
0
)
,
bottom)
#
Non-homogeneous
Dirichlet
boundary
CoRX
=
0
.
5
*
(xmax
+
xmin)
CoRY
=
0
.
5
*
(ymax
+
ymin)
dispLoadX
=
Expression
(
"
(
x
[
0]
-
CoR_x)*cos
(
t)
-
(
x
[
1]
-
CoR_y)*sin
(
t)
-
\
(
x
[
0]
-
CoR_x)
"
,
t
=
ubar
,
CoR_x
=
CoRX
,
CoR_y
=
CoRY
,
degree
=
1
)
dispLoadY
=
Expression
(
"
(
x
[
0]
-
CoR_x)*sin
(
t)
+
(
x
[
1]
-
CoR_y)*cos
(
t)
-
\
(
x
[
1]
-
CoR_y)
"
,
t
=
ubar
,
CoR_x
=
CoRX
,
CoR_y
=
CoRY
,
degree
=
1
)
BC_top_0
=
DirichletBC
(V_0
,
dispLoadX
,
top)
BC_top_1
=
DirichletBC
(V_1
,
dispLoadY
,
top)
BC_top_2
=
DirichletBC
(V_2
,
Constant
(
0
.
0
)
,
top)
BC
=
[BC_bottom_0
,
BC_bottom_1
,
BC_bottom_2
,
BC_top_0
,
BC_top_1
,
BC_top_2]
#
-------------------------------------------------------------------------
#
Generalized
strains
and
stresses
#
-------------------------------------------------------------------------
def
epsilon_sym
(
v
):
strain_sym
=
as_tensor
([[ v[
0
]
.
dx
(
0
)
,
0
.
5
*
(v[
0
]
.
dx
(
1
)
+
v[
1
]
.
dx
(
0
))
,
0
.
5
*
(v[
0
]
.
dx
(
2
)
+
v[
2
]
.
dx
(
0
)) ]
,
[
0
.
5
*
(v[
0
]
.
dx
(
1
)
+
v[
1
]
.
dx
(
0
))
,
v[
1
]
.
dx
(
1
)
,
0
.
5
*
(v[
1
]
.
dx
(
2
)
+
v[
2
]
.
dx
(
1
)) ]
,
[
0
.
5
*
(v[
0
]
.
dx
(
2
)
+
v[
2
]
.
dx
(
0
))
,
0
.
5
*
(v[
1
]
.
dx
(
2
)
+
v[
2
]
.
dx
(
1
))
,
v[
2
]
.
dx
(
2
) ]])
return
strain_sym
def
rot_rel
(
v
,
theta
):
star_eps_skw
=
0
.
5
*
as_vector
((v[
2
]
.
dx
(
1
)
-
v[
1
]
.
dx
(
2
)
,
v[
0
]
.
dx
(
2
)
-
v[
2
]
.
dx
(
0
)
,
v[
1
]
.
dx
(
0
)
-
v[
0
]
.
dx
(
1
)))
theta_rel
=
theta
-
star_eps_skw
return
theta_rel
C_rot
=
rotate_elast
(C11
,
C12
,
C44
,
o1
,
o2
,
o3
,
ang)
def
sigma
(
v
):
eps_sym
=
epsilon_sym
(v)
epsilon_11
=
eps_sym[
0
,
0
]
epsilon_22
=
eps_sym[
1
,
1
]
epsilon_33
=
eps_sym[
2
,
2
]
epsilon_12
=
eps_sym[
0
,
1
]
epsilon_13
=
eps_sym[
0
,
2
]
epsilon_23
=
eps_sym[
1
,
2
]
#
Voigt
notation
(
considering
engineering
shear
strain)
sig_vec
=
np
.
array
([epsilon_11
,
epsilon_22
,
epsilon_33
,
2
*
epsilon_23
,
2
*
epsilon_13
,
2
*
epsilon_12])
sig_vec
=
np
.
dot
(C_rot
,
sig_vec)
stress
=
as_tensor
([[ sig_vec[
0
]
,
sig_vec[
5
]
,
sig_vec[
4
] ]
,
[ sig_vec[
5
]
,
sig_vec[
1
]
,
sig_vec[
3
] ]
,
[ sig_vec[
4
]
,
sig_vec[
3
]
,
sig_vec[
2
] ]])
return
stress
def
strain_energy
(
v
,
theta
):
eps_sym
=
epsilon_sym
(v)
sig
=
sigma
(v)
theta_rel
=
rot_rel
(v
,
theta)
energy
=
0
.
5
*
(
inner
(eps_sym
,
sig)
+
kappa
*
inner
(theta_rel
,
theta_rel))
return
energy
#
-------------------------------------------------------------------------
#
Variational
form
#
-------------------------------------------------------------------------
v
,
theta
=
TrialFunctions
(U)
vv
,
vtheta
=
TestFunctions
(U)
#
Bilinear
form
a
=
inner
(
epsilon_sym
(vv)
,
sigma
(v))
*
dx
\
+
kappa
*
inner
(
rot_rel
(vv
,
vtheta)
,
rot_rel
(v
,
theta))
*
dx
#
-------------------------------------------------------------------------
#
Solution
#
-------------------------------------------------------------------------
u_h
=
Function
(U)
L
=
dot
(vv
,
Constant
((
0
.
,
0
.
,
0
.
)))
*
dx
#
Zero
external
force
solve
(a
==
L
,
u_h
,
BC
,
solver_parameters
=
{
"
linear_solver
"
:
"
mumps
"
})
v_h
,
theta_h
=
u_h
.
split
()
#
-------------------------------------------------------------------------
#
Postprocessing
#
-------------------------------------------------------------------------
E_e
=
assemble
(
strain_energy
(v_h
,
theta_h)
*
dx)
print
(
'
Energy:
'
,
E_e
,
'
[
MPa]
'
)
#
Field
variables
vtkfile_v
=
File
(
'
./output
/
rotcube_v.pvd
'
)
vtkfile_v
<<
v_h
vtkfile_theta
=
File
(
'
./output
/
rotcube_theta.pvd
'
)
vtkfile_theta
<<
theta_h