#
=========================================================================
#
FEniCS
implementation:
Homogenized
model
of
a
single-edge
notched
plate
#
with
a
2D
honeycomb
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
subprocess
#
-------------------------------------------------------------------------
#
Create
mesh
#
-------------------------------------------------------------------------
dfact
=
12
#
Domain
scaling
hsize
=
0
.
25
#
Finite
element
size
[
mm]
subprocess
.
run
([
'
dolfin-convert
'
,
'
./mesh
/
sen_dfact_
'
+
str
(dfact)
\
+
'
_h_
'
+
str
(hsize)
+
'
.msh
'
,
'
./mesh
/
sen_dfact_
'
\
+
str
(dfact)
+
'
_h_
'
+
str
(hsize)
+
'
.xml
'
])
mesh
=
Mesh
(
'
./mesh
/
sen_dfact_
'
+
str
(dfact)
+
'
_h_
'
+
str
(hsize)
+
'
.xml
'
)
#
-------------------------------------------------------------------------
#
Input
parameters
#
-------------------------------------------------------------------------
#
Metamaterial
parameters
Lb
=
2
.
#
Bar
length
[
mm]
rb
=
0
.
1
*
Lb
#
Bar
width
[
mm]
Eb
=
430
.
#
Bar
Young's
modulus
[
N
/
mm^2]
Ab
=
1
*
rb
#
Bar
area
[
mm^2]
Ib
=
(
1
/
12
)
*
1
*
rb
**
3
#
Bar
MoM
[
mm^4]
#
Homogenized
properties
lmbda
=
(
1
/
(
2
*
sqrt
(
3
)))
*
Eb
*
Ab
*
(
-
12
*
Eb
*
Ib
+
Eb
*
Ab
*
Lb
**
2
)
/
(
12
*
Eb
*
Ib
*
Lb
\
+
Eb
*
Ab
*
Lb
**
3
)
mubar
=
4
*
sqrt
(
3
)
*
Eb
*
Ab
*
Eb
*
Ib
/
(
12
*
Eb
*
Ib
*
Lb
+
Eb
*
Ab
*
Lb
**
3
)
kappa
=
8
*
sqrt
(
3
)
*
Eb
*
Ib
/
(Lb
*
Lb
**
2
)
#
Applied
load
ubar
=
dfact
#
Displacement
[
mm]
#
-------------------------------------------------------------------------
#
Function
spaces
#
-------------------------------------------------------------------------
v_element
=
VectorElement
(
'
CG
'
,
mesh
.
ufl_cell
()
,
2
)
#
Deflections
t_element
=
FiniteElement
(
'
CG
'
,
mesh
.
ufl_cell
()
,
1
)
#
Rotations
U
=
FunctionSpace
(mesh
,
v_element
*
t_element)
V
,
T
=
U
.
split
()
#
1
DOF
for
rotations
V_0
,
V_1
=
V
.
split
()
#
2
DOFs
for
deflections
#
-------------------------------------------------------------------------
#
Lists,
DOFs
and
BCs
#
-------------------------------------------------------------------------
coords
=
mesh
.
coordinates
()
xmin
=
min
(coords[
:
,
0
])
ymin
=
min
(coords[
:
,
1
])
ymax
=
max
(coords[
:
,
1
])
top
=
CompiledSubDomain
(
"
near
(
x
[
1],
ymax)
&&
on_boundary
"
,
ymax
=
ymax)
bottom
=
CompiledSubDomain
(
"
near
(
x
[
1],
ymin)
&&
on_boundary
"
,
ymin
=
ymin)
left
=
CompiledSubDomain
(
"
near
(
x
[
0],
xmin)
&&
on_boundary
"
,
xmin
=
xmin)
#
Dirichlet
boundary
BC_v_left
=
DirichletBC
(V_0
,
Constant
(
0
.
0
)
,
left)
BC_v_bottom_0
=
DirichletBC
(V_0
,
Constant
(
0
.
0
)
,
bottom)
BC_v_bottom_1
=
DirichletBC
(V_1
,
Constant
(
0
.
0
)
,
bottom)
BC_v_top_0
=
DirichletBC
(V_0
,
Constant
(
0
.
0
)
,
top)
BC_t_left
=
DirichletBC
(T
,
Constant
(
0
.
0
)
,
left)
BC_t_bottom
=
DirichletBC
(T
,
Constant
(
0
.
0
)
,
bottom)
#
Non-homogeneous
Dirichlet
boundary
dispLoad
=
Expression
(
"
t
"
,
t
=
ubar
,
degree
=
1
)
BC_v_top_1
=
DirichletBC
(V_1
,
dispLoad
,
top)
BC
=
[BC_v_left
,
BC_v_bottom_0
,
BC_v_bottom_1
,
BC_t_left
,
BC_t_bottom
,
\
BC_v_top_0
,
BC_v_top_1]
#
-------------------------------------------------------------------------
#
Generalized
strains
and
stresses
#
-------------------------------------------------------------------------
def
epsilon_sym
(
v
):
strain_sym
=
sym
(
grad
(v))
return
strain_sym
def
rot_rel
(
v
,
theta
):
eps
=
skew
(
grad
(v))
star_eps_skw
=
0
.
5
*
(eps[
1
,
0
]
-
eps[
0
,
1
])
theta_rel
=
theta
-
star_eps_skw
return
theta_rel
def
sigma
(
v
):
eps_sym
=
epsilon_sym
(v)
stress
=
lmbda
*
tr
(eps_sym)
*
Identity
(
2
)
+
2
.
*
mubar
*
eps_sym
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
.
)))
*
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
/
(dfact
**
2
)
,
'
[
MPa]
'
)
#
Field
variables
vtkfile_v
=
File
(
'
./output
/
sen_v.pvd
'
)
vtkfile_v
<<
v_h
vtkfile_theta
=
File
(
'
./output
/
sen_theta.pvd
'
)
vtkfile_theta
<<
theta_h