%
=========================================================================
%
Stabil
implementation:
DNS
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
%
=========================================================================
clear
close
all
clc
%%
Beam
parameters
Types
=
{
1
'
beam
'
}
;
Lc
=
0
.
75
;
%
Cell
length
[
mm]
Eb
=
430
;
%
Bar
Young's
modulus
[
N
/
mm^2]
rb
=
0
.
065
/
2
;
%
Bar
radius
[
mm]
Ab
=
pi
*
rb
^
2
;
%
Bar
area
[
mm^2]
Ib
=
pi
*
rb
^
4
/
4
;
%
Bar
MoI
[
mm^4]
Lb
=
Lc
/
sqrt
(
2
)
;
%
Bar
length
[
mm]
%
Sections
=
[
SecID
A
ky
kz
Ixx
Iyy
Izz]
(
Euler-Bernoulli
beam)
Sections
=
[
1
Ab
Inf
Inf
2
*
Ib Ib Ib]
;
Materials
=
[
1
Eb
0
.
0
]
;
%%
Load
mesh
dfact
=
6
;
Nodes
=
load
([
'
mesh
/
rotcube_dfact_
'
,
num2str
(dfact)
'
_Nodes.dat
'
])
;
Elements
=
load
([
'
mesh
/
rotcube_dfact_
'
,
num2str
(dfact)
'
_Elems.dat
'
])
;
%%
Lists
and
DOFs
tol
=
Lc
/
50
;
TopNodeIDs
=
find
(
abs
(Nodes(
1
:
end
-
1
,
4
)
-
max
(Nodes(
1
:
end
-
1
,
4
)))
<
tol)
;
BottomNodeIDs
=
find
(
abs
(Nodes(
1
:
end
-
1
,
4
)
-
min
(Nodes(
1
:
end
-
1
,
4
)))
<
tol)
;
RightNodeIDs
=
find
(
abs
(Nodes(
1
:
end
-
1
,
2
)
-
max
(Nodes(
1
:
end
-
1
,
2
)))
<
tol)
;
LeftNodeIDs
=
find
(
abs
(Nodes(
1
:
end
-
1
,
2
)
-
min
(Nodes(
1
:
end
-
1
,
2
)))
<
tol)
;
BackNodeIDs
=
find
(
abs
(Nodes(
1
:
end
-
1
,
3
)
-
max
(Nodes(
1
:
end
-
1
,
3
)))
<
tol)
;
FrontNodeIDs
=
find
(
abs
(Nodes(
1
:
end
-
1
,
3
)
-
min
(Nodes(
1
:
end
-
1
,
3
)))
<
tol)
;
DOF
=
getdof(Elements
,
Types)
;
seldof
=
[BottomNodeIDs
+
0
.
01
;
BottomNodeIDs
+
0
.
02
;
BottomNodeIDs
+
0
.
03
;
TopNodeIDs
+
0
.
03
]
;
DOF
=
removedof(DOF
,
seldof)
;
K
=
asmkm(Nodes
,
Elements
,
Types
,
Sections
,
Materials
,
DOF)
;
%%
Displacement
loads
(
use
constraints)
vartheta
=
2
*
pi
/
180
;
CoR
=
[
0
.
5
*
(
max
(Nodes(
1
:
end
-
1
,
2
))
+
min
(Nodes(
1
:
end
-
1
,
2
)))
...
0
.
5
*
(
max
(Nodes(
1
:
end
-
1
,
3
))
+
min
(Nodes(
1
:
end
-
1
,
3
)))]
.'
;
appDispX
=
zeros
(
length
(TopNodeIDs)
,
1
)
;
appDispY
=
zeros
(
length
(TopNodeIDs)
,
1
)
;
for
iN
=
1
:
length
(TopNodeIDs)
nCoords
=
Nodes(TopNodeIDs(iN)
,
2
:
3
)
.'
;
shiftedCoords
=
nCoords
-
CoR
;
R
=
[
cos
(vartheta)
,
-
sin
(vartheta)
;
sin
(vartheta)
,
cos
(vartheta)]
;
appDisp
=
R
*
shiftedCoords
-
shiftedCoords
;
appDispX(iN)
=
appDisp(
1
)
;
appDispY(iN)
=
appDisp(
2
)
;
end
Constr
=
[appDispX
ones
(
numel
(TopNodeIDs)
,
1
) TopNodeIDs
+
0
.
01
;
appDispY
ones
(
numel
(TopNodeIDs)
,
1
) TopNodeIDs
+
0
.
02
]
;
P
=
zeros
(
size
(K
,
1
)
,
1
)
;
[T
,
Q0]
=
tconstr(Constr
,
DOF)
;
%%
Solver
[Ur
,
FLAG
,
RELRES]
=
gmres
((T
.'
*
K
*
T)
,
(T
.'
*
(P
-
K
*
Q0))
,
100
,
1e-8
,
100
)
;
U
=
T
*
Ur
+
Q0
;
%%
Post
[aDOF]
=
getdof(Elements
,
Types)
;
seldof
=
[Nodes(
end
,
1
)
+
0
.
0
]
;
aDOF
=
removedof(aDOF
,
seldof)
;
Ldof
=
selectdof(aDOF
,
DOF)
;
Ua
=
Ldof
.'
*
U
;
ux
=
Ua(
1
:
6
:
end
,
1
)
;
uy
=
Ua(
2
:
6
:
end
,
1
)
;
uz
=
Ua(
3
:
6
:
end
,
1
)
;
ox
=
Ua(
4
:
6
:
end
,
1
)
;
oy
=
Ua(
5
:
6
:
end
,
1
)
;
oz
=
Ua(
6
:
6
:
end
,
1
)
;
kinmap
=
turbo
;
minX
=
min
(Nodes(
1
:
end
-
1
,
2
))
;
minY
=
min
(Nodes(
1
:
end
-
1
,
3
))
;
minZ
=
min
(Nodes(
1
:
end
-
1
,
4
))
;
maxX
=
max
(Nodes(
1
:
end
-
1
,
2
))
;
maxY
=
max
(Nodes(
1
:
end
-
1
,
3
))
;
maxZ
=
max
(Nodes(
1
:
end
-
1
,
4
))
;
gap
=
[
0
.
16
0
.
03
]
;
margh
=
[
0
.
12
0
.
01
]
;
margw
=
[
0
.
012
0
.
012
]
;
csize
=
18
;
fsize
=
20
;
pointsize
=
25
;
caxux
=
[
-
0
.
05
,
0
.
05
]
;
caxuy
=
[
-
0
.
05
,
0
.
05
]
;
caxuz
=
[
-
0
.
004
,
0
.
004
]
;
caxox
=
[
-
0
.
018
,
0
.
018
]
;
caxoy
=
[
-
0
.
018
,
0
.
018
]
;
caxoz
=
[
0
,
0
.
035
]
;
h
=
figure
;
set
(h
,
'
Position
'
,
[
680
558
900
700
])
;
set
(h
,
'
color
'
,
[
1
1
1
])
;
ax(
1
)
=
subtightplot(
2
,
3
,
1
,
gap
,
margh
,
margw)
;
scatter3
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
Nodes(
1
:
end
-
1
,
4
)
,
...
pointsize
,
ux
/
dfact
,
'
filled
'
)
;
colorbar
off
colormap
(kinmap)
caxis
(caxux)
axis
equal
;
axis
off
;
view
(
37
.
5
-
90
,
30
)
cb
=
colorbar
(
'
location
'
,
'
south
'
,
'
fontsize
'
,
csize)
;
cb
.
AxisLocation
=
'
out
'
;
cb
.
Ruler
.
Exponent
=
-
0
;
cb
.
Position
=
cb
.
Position
-
[
-
0
.
02
0
.
05
0
.
04
0
]
;
cb
.
Title
.
Interpreter
=
'
latex
'
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
cb
.
Title
.
String
=
'
Deflection
$v_{1}$
[
mm]
'
;
cb
.
Title
.
Position
=
cb
.
Title
.
Position
-
[
10
75
0
]
;
cb
.
Title
.
FontSize
=
fsize
;
ax(
2
)
=
subtightplot(
2
,
3
,
2
,
gap
,
margh
,
margw)
;
scatter3
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
Nodes(
1
:
end
-
1
,
4
)
,
...
pointsize
,
uy
/
dfact
,
'
filled
'
)
;
colorbar
off
colormap
(kinmap)
caxis
(caxuy)
axis
equal
;
axis
off
;
view
(
37
.
5
-
90
,
30
)
cb
=
colorbar
(
'
location
'
,
'
south
'
,
'
fontsize
'
,
csize)
;
cb
.
AxisLocation
=
'
out
'
;
cb
.
Ruler
.
Exponent
=
-
0
;
cb
.
Position
=
cb
.
Position
-
[
-
0
.
02
0
.
05
0
.
04
0
]
;
cb
.
Title
.
Interpreter
=
'
latex
'
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
cb
.
Title
.
String
=
'
Deflection
$v_{2}$
[
mm]
'
;
cb
.
Title
.
Position
=
cb
.
Title
.
Position
-
[
10
75
0
]
;
cb
.
Title
.
FontSize
=
fsize
;
ax(
3
)
=
subtightplot(
2
,
3
,
3
,
gap
,
margh
,
margw)
;
scatter3
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
Nodes(
1
:
end
-
1
,
4
)
,
...
pointsize
,
uz
/
dfact
,
'
filled
'
)
;
colorbar
off
colormap
(kinmap)
caxis
(caxuz)
axis
equal
;
axis
off
;
view
(
37
.
5
-
90
,
30
)
cb
=
colorbar
(
'
location
'
,
'
south
'
,
'
fontsize
'
,
csize)
;
cb
.
AxisLocation
=
'
out
'
;
cb
.
Ruler
.
Exponent
=
-
0
;
cb
.
Position
=
cb
.
Position
-
[
-
0
.
02
0
.
05
0
.
04
0
]
;
cb
.
Title
.
Interpreter
=
'
latex
'
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
cb
.
Title
.
String
=
'
Deflection
$v_{3}$
[
mm]
'
;
cb
.
Title
.
Position
=
cb
.
Title
.
Position
-
[
10
75
0
]
;
cb
.
Title
.
FontSize
=
fsize
;
ax(
4
)
=
subtightplot(
2
,
3
,
4
,
gap
,
margh
,
margw)
;
scatter3
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
Nodes(
1
:
end
-
1
,
4
)
,
...
pointsize
,
ox
,
'
filled
'
)
;
colorbar
off
colormap
(kinmap)
caxis
(caxox)
axis
equal
;
axis
off
;
view
(
37
.
5
-
90
,
30
)
cb
=
colorbar
(
'
location
'
,
'
south
'
,
'
fontsize
'
,
csize)
;
cb
.
AxisLocation
=
'
out
'
;
cb
.
Ruler
.
Exponent
=
-
0
;
cb
.
Ticks
=
[
-
0
.
018
0
0
.
018
]
;
cb
.
Position
=
cb
.
Position
-
[
-
0
.
02
0
.
05
0
.
04
0
]
;
cb
.
Title
.
Interpreter
=
'
latex
'
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
cb
.
Title
.
String
=
'
Rotation
$
\t
heta_{1}$
[-]
'
;
cb
.
Title
.
Position
=
cb
.
Title
.
Position
-
[
10
75
0
]
;
cb
.
Title
.
FontSize
=
fsize
;
ax(
5
)
=
subtightplot(
2
,
3
,
5
,
gap
,
margh
,
margw)
;
scatter3
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
Nodes(
1
:
end
-
1
,
4
)
,
...
pointsize
,
oy
,
'
filled
'
)
;
colorbar
off
colormap
(kinmap)
caxis
(caxoy)
axis
equal
;
axis
off
;
view
(
37
.
5
-
90
,
30
)
cb
=
colorbar
(
'
location
'
,
'
south
'
,
'
fontsize
'
,
csize)
;
cb
.
AxisLocation
=
'
out
'
;
cb
.
Ruler
.
Exponent
=
-
0
;
cb
.
Ticks
=
[
-
0
.
018
0
0
.
018
]
;
cb
.
Position
=
cb
.
Position
-
[
-
0
.
02
0
.
05
0
.
04
0
]
;
cb
.
Title
.
Interpreter
=
'
latex
'
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
cb
.
Title
.
String
=
'
Rotation
$
\t
heta_{2}$
[-]
'
;
cb
.
Title
.
Position
=
cb
.
Title
.
Position
-
[
10
75
0
]
;
cb
.
Title
.
FontSize
=
fsize
;
ax(
6
)
=
subtightplot(
2
,
3
,
6
,
gap
,
margh
,
margw)
;
scatter3
(Nodes(
1
:
end
-
1
,
2
)
,
Nodes(
1
:
end
-
1
,
3
)
,
Nodes(
1
:
end
-
1
,
4
)
,
...
pointsize
,
oz
,
'
filled
'
)
;
colorbar
off
colormap
(kinmap)
caxis
(caxoz)
axis
equal
;
axis
off
;
view
(
37
.
5
-
90
,
30
)
cb
=
colorbar
(
'
location
'
,
'
south
'
,
'
fontsize
'
,
csize)
;
cb
.
AxisLocation
=
'
out
'
;
cb
.
Ruler
.
Exponent
=
-
0
;
cb
.
Ticks
=
[
0
0
.
035
/
2
0
.
035
]
;
cb
.
Position
=
cb
.
Position
-
[
-
0
.
02
0
.
05
0
.
04
0
]
;
cb
.
Title
.
Interpreter
=
'
latex
'
;
cb
.
TickLabelInterpreter
=
'
latex
'
;
cb
.
Title
.
String
=
'
Rotation
$
\t
heta_{3}$
[-]
'
;
cb
.
Title
.
Position
=
cb
.
Title
.
Position
-
[
10
75
0
]
;
cb
.
Title
.
FontSize
=
fsize
;
E_e
=
energy_bars(Nodes
,
Elements
,
Types
,
Sections
,
Materials
,
DOF
,
U)
;
fprintf
([
'
Energy:
'
,
num2str
(E_e
/
dfact
^
3
)
,
'
[
MPa]
\n
'
])