%
=========================================================================
%
Stabil
implementation:
DNS
of
a
perforated
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
%
=========================================================================
clear
close
all
clc
%%
Beam
parameters
Types
=
{
1
'
beam
'
}
;
Lb
=
2
;
%
Bar
length
[
mm]
rb
=
0
.
1
*
Lb
;
%
Bar
width
[
mm]
Eb
=
430
;
%
Bar
Young's
modulus
[
N
/
mm^2]
Ab
=
rb
*
1
;
%
Bar
area
[
mm^2]
Ib
=
rb
^
3
*
1
/
12
;
%
Bar
MoM
[
mm^4]
%
Sections
=
[
SecID
A
ky
kz
Ixx
Iyy
Izz]
(
Euler-Bernoulli
beam)
Sections
=
[
1
Ab
Inf
Inf
0
0
Ib]
;
Materials
=
[
1
Eb
0
.
0
]
;
%%
Load
mesh
dfact
=
12
;
%
Scale
the
domain
Nodes
=
load
([
'
mesh
/
holeinplate_dfact_
'
,
num2str
(dfact)
'
_Nodes.dat
'
])
;
Elements
=
load
([
'
mesh
/
holeinplate_dfact_
'
,
num2str
(dfact)
'
_Elems.dat
'
])
;
%%
Lists
and
DOFs
tol
=
Lb
/
50
;
TopNodeIDs
=
find
(
abs
(Nodes(
:
,
3
)
-
max
(Nodes(
1
:
end
-
1
,
3
)))
<
tol)
;
RightNodeIDs
=
find
(
abs
(Nodes(
:
,
2
)
-
max
(Nodes(
1
:
end
-
1
,
2
)))
<
tol)
;
LeftNodeIDs
=
find
(
abs
(Nodes(
:
,
2
)
-
min
(Nodes(
:
,
2
)))
<
tol)
;
BottomNodeIDs
=
find
(
abs
(Nodes(
:
,
3
)
-
min
(Nodes(
:
,
3
)))
<
tol)
;
DOF
=
getdof(Elements
,
Types)
;
seldof
=
[Nodes(
:
,
1
)
+
0
.
03
;
Nodes(
:
,
1
)
+
0
.
04
;
Nodes(
:
,
1
)
+
0
.
05
;
...
BottomNodeIDs
+
0
.
02
;
LeftNodeIDs
+
0
.
01
;
LeftNodeIDs
+
0
.
06
;
BottomNodeIDs
+
0
.
06
;
TopNodeIDs
+
0
.
06
]
;
DOF
=
removedof(DOF
,
seldof)
;
K
=
asmkm(Nodes
,
Elements
,
Types
,
Sections
,
Materials
,
DOF)
;
%%
Displacement
loads
(
use
constraints)
ubar
=
dfact
;
%
Applied
displacement
(
scaled)
[
mm]
nTop
=
numel
(TopNodeIDs)
;
%
Constr
=
[
Constant
Coef1
DOFs
...]
Constr
=
[ubar
*
ones
(nTop
,
1
)
ones
(nTop
,
1
) TopNodeIDs
+
0
.
02
]
;
P
=
zeros
(
size
(K
,
1
)
,
1
)
;
[K
,
P]
=
addconstr(Constr
,
DOF
,
K
,
P)
;
%%
Solver
U
=
K
\
P
;
%%
Post
[aDOF]
=
getdof(Elements
,
Types)
;
seldof
=
[Nodes(
:
,
1
)
+
0
.
03
;
Nodes(
:
,
1
)
+
0
.
04
;
Nodes(
:
,
1
)
+
0
.
05
;
Nodes(
end
,
1
)]
;
aDOF
=
removedof(aDOF
,
seldof)
;
Ldof
=
selectdof(aDOF
,
DOF)
;
Ua
=
Ldof
.'
*
U
;
ux
=
Ua(
1
:
3
:
end
,
1
)
;
uy
=
Ua(
2
:
3
:
end
,
1
)
;
omz
=
Ua(
3
:
3
:
end
,
1
)
;
caxux
=
[
-
0
.
9
,
0
]
;
caxuy
=
[
0
,
1
]
;
caxoz
=
[
-
0
.
05
,
0
]
;
colmap
=
turbo
;
h
=
figure
;
set
(h
,
'
color
'
,
[
1
1
1
])
;
pointsize
=
10
;
scatter
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
pointsize
,
ux
/
dfact
,
'
filled
'
)
;
axis
off
axis
equal
colormap
(colmap)
cb
=
colorbar
;
cb
.
FontSize
=
18
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
caxis
(caxux)
h
=
figure
;
set
(h
,
'
color
'
,
[
1
1
1
])
;
pointsize
=
10
;
scatter
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
pointsize
,
uy
/
dfact
,
'
filled
'
)
;
axis
off
axis
equal
colormap
(colmap)
cb
=
colorbar
;
cb
.
FontSize
=
18
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
caxis
(caxuy)
h
=
figure
;
set
(h
,
'
color
'
,
[
1
1
1
])
;
pointsize
=
10
;
scatter
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
pointsize
,
omz
,
'
filled
'
)
;
axis
off
axis
equal
colormap
(colmap)
cb
=
colorbar
;
cb
.
FontSize
=
18
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
caxis
(caxoz)
E_e
=
energy_bars(Nodes
,
Elements
,
Types
,
Sections
,
Materials
,
DOF
,
U)
;
fprintf
([
'
Energy:
'
,
num2str
(E_e
/
(dfact
^
2
))
,
'
[
MPa]
\n
'
])