The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Simultaneous Optimization of Topology and Orientationof Anisotropic Material Using
Isoparametric Projection Method
Tsuyoshi Nomura1
1

Toyota Central R&D Labs., Inc. Aichi, Japan, nomu2@mosk.tytlabs.co.jp

November 27, 2015

1/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Acknowledgments
Collaborators• Dr. Ercan M. Dede (Toyota Research Institute of North America)• Prof. Shintaro Yamasaki (Osaka University)• Dr. Atsushi Kawamoto (Toyota Central R&D Labs., Inc.)• Dr. Tadayoshi Matsumori (Toyota Central R&D Labs., Inc.)• Prof. Kazuhiro Izui (Kyoto University)• Prof. Shinji Nishiwaki (Kyoto University)

Funding
This project is partially supported by Strategic Innovation Promotion (SIP)
program, Innovative Design and Production Technology Project

2/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

4/36

Recent fabrication technologies and Orientation design
Tailored Fiber Placement⇒Super Robotics Sewing
MachineResin transfer molding

https://www.behance.net/gallery/11063817/L1-lightweight-stool
http://www.ipfdd.de/en/ departments/institute- of-polymer-materials/
composite-materials/ fields-of-research/ complex-structuralcomponents-tailored-fibre- placement-tfp/
aaa

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Recent fabrication technologies and Orientation design
Continuous fiber printing

Matsuzaki Lab / Tokyo University of Science1
http://www.rs.tus.ac.jp/rmatsuza/

5/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Recent fabrication technologies and Orientation design
Continuous fiber printing

Markforged MarkOne 1
http://www.markforged.com/

5/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Recent fabrication technologies and Orientation design
Toolpath depended anisotropy (Even without fiber)aa

http://www.re3d.org/Topology and Orientation design of 2D / shell is becoming practical andmore important

6/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Challenges of orientation design problems
Initial
Early stage
Final configuration
Variation in space∂Duy∂Dttx• Large local change in orientation• Has to be flexible for local change• Complex structure due to topological change• Spatial regularization (filtering) is necessary

7/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Parameterization schemes for Orientation design 1/2
Continuous Fiber Angle Optimization

Homogenization Design Method

Pedersen P.1989

Bendsøe M. P and Kikuchi N, 1988

Directly use angle as design variableθ

Local minima issue
cos(θ) − sin(θ)R(θ) =,sin(θ) cos(θ)

CGijkl =2
X
I,J,K ,L=1

CHIJKL RIi (θ)RJj (θ)RKk (θ)RLl (θ)

8/36

9/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, JapanParameterization schemes for Orientation design 1/2
More restrictive
Discrete Material Optimization (DMO)
Stegmann J and Lund E. 2005
• Based on multimaterial topology optimization

Cijkl =!
ρp1+ ρp2+ ρp3+ · · · + ρpnρ 1 + ρ 2 + ρ3 + · · · + ρn + 6 1• Each orientation corresponds to a material

More permissive
Free Material Optimization• Design all independent tensor element under loose constraints(1)

10/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Parametarization schemes for Orientation design 2/2• Vector element parameterization

Orientation vector
" #x
B = Br ry

−1 6 {x, y} 6 1 ,
kxk + kyk 6 1

JS Choi, J Yoo, 2008

yBoxconstraintx(2)Linearconstraint

Unit circle Expensive (Pointwise nonlinear constraint)
Box constraint Low cost (free) but insufficient (infeasible properties)
Linear constraint Acceptable cost but too tight
Converges into discrete design but depends on global
coordinate axis

11/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Advantage of vector parameterization

Vector element

Orientation angle

Distribution

After filtering
θinit

θinit
Blurring
Alias angle

Discontinuity−π, π−π, π0

θoptxy0θopt

No discontinuitiesθinitatan(y/x)θopt

No alias angles
xyatan(y/x)

14/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, JapanPDE density filter / Vector case
Actual design variable   ξ¯ξξ˜• Precursor variable,, Filtered variable,, Projected variable,η̄ηη̃• Helmholtz filtering2FilterIradius   ξ¯ξ˜ξ˜− Rυ ∇+=,η̄η̃η̃
2• Heaviside projection
Differentiable Heaviside function

1-10-1

1

ξ(x)
η(x)



2 H̃

=
2 H̃

 ˜ξ(x)−1 .
η̃(x) −1(5)(6)

14/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, JapanPDE density filter / Vector case
Actual design variable   ξ¯ξξ˜• Precursor variable,, Filtered variable,, Projected variable,η̄ηη̃• Helmholtz filtering2FilterIradius   ξ¯ξ˜ξ˜− Rυ ∇+=,η̄η̃η̃
2• Heaviside projection
Differentiable Heaviside function

1-10-1

1

ξ(x)
η(x)



2 H̃

=
2 H̃

 ˜
ξ(x)
−1
 .
η̃(x) −1(5)(6)

16/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Isoparametric projection method ξ(Projected variable)
ηSerendipity element−1 6 ξ 6 1, −1 6 η 6 1,Second order• Isoparametric projection,eight node quadrilateralc
ϑ== N ξ, η .
(7)elements• Vectorial design variable,• Orientation vector, c 2 + s 2 6 1yη
Isoparaξ metricprojection−1 6 ξ 6 1,−1 6 η 6 1yx

Shape function: N

Transform box constraints to circular constraints⊂x

c2 + s2 6 1

17/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Isoparametric projection method for discrete angle setyηyIsoparaξ metricprojection−1 6 ξ 6 1,−1 6 η 6 1

6 options / 2 elementsxx

Shape function: Nc2 + s2 6 18 options / 4 elements
η2211ξ
34

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Basic idea

• Cijkl has 4th order pronominal of cos(θ) and sin(θ) .

cos(θ)• These are the element of the orientation vector, i.e. e =.
sin(θ)
 c• We also have ϑ =, we can replacescos(θ) and sin(θ) by c and s , respectively

19/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

20/36

Material interpolation
Solid part (Anisotropic)
Void part 0−1C( ρ , ϑ ) = Cv + ρp (8)
 Ci + T̂ (ϑ̃) · Cu − Ci · T̂ (ϑ̃) − Cv  .
|{z}




Rotated Tensor
Isotropic component of Cu
Void partĈ11 ( c, s ) Ĉ12 ( c, s ) Ĉ16 ( c, s )T̂−1 (ϑ)CT̂0 (ϑ) = Ĉ(ϑ) =Ĉ12 ( c, s ) Ĉ22 ( c, s ) Ĉ26 ( c, s )Ĉ16 ( c, s ) Ĉ26 ( c, s ) Ĉ66 ( c, s )|{z}
4th order Polynomial functions of c and s
E depends onp 2(c + s 2 )(9)

21/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Material interpolation
Rotation tensor can be written as 2

cs22csT(ϑ) =  s 2 c 2 −2cs  ,−cs cs c 2 − s 2c2
s2csT0 (ϑ) =  s 2c2−cs  .
2−2cs 2cs c − s 2
(10)therefore, elements of Ĉ(ϑ) areĈ11 = C11 c 4 − 4C16 c 3 s + (2C12 + 4C66 )c 2 s 2 − 4C26 cs 3 + C22 s 4Ĉ12 = C12 c 4 + 2(C16 − C26 )c 3 s + (C11 + C22 − 4C66 )c 2 s 2 − 2(C16 − C26 )cs 3 + C12 s 4
Ĉ16 = C16 c 4 + (C11 − C12 − 2C66 )c 3 s − 3(C16 − C26 )c 2 s 2 + (C12 − C22 + 2C66 )cs 3 − C26 s 4Ĉ22 = C22 c 4 + 4C26 c 3 s + (2C12 + 4C66 )c 2 s 2 + 4C16 cs 3 + C11 s 4
Ĉ26 = C26 c 4 + (C12 − C22 + 2C66 )c 3 s + 3(C16 − C26 )c 2 s 2 + (C11 − C12 − 2C66 )cs 3 − C16 s 4Ĉ66 = 2(C16 − C26 )c 3 s + (C11 − 2C12 + C22 − 2C66 )c 2 s 2 − 2(C16 − C26 )cs 3 + C66 (c 4 + s 4 ).
(11)

4th order polynomials, no trigonometric functions

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

22/36

Mapping to −π/2
< θ < π/2 ccovers whole circle, −π < θ < π, but since Cijkl iss
non-directional, we want to map −π < θ < π into −π/2 < θ < π/2
Variable ϑ =2.
cos3 (θ/2) sin(θ/2) = sin(θ) (1 + cos(θ)) 44.
.
cos2 (θ/2) sin2 (θ/2) = (1 − cos(θ))(1 + cos(θ)) 4cos(θ/2) sin3 (θ/2) = sin(θ) (1 + cos(θ)) 4.
2sin4 (θ/2) = (1 − cos(θ)) 4 ,(Half angle formulas)cos4 (θ/2) = (1 + cos(θ))

.

(12)

and thereforeno(1+c)2 C11 − (1−c)2 C22 − 4 (1+c) sC16 − 4 (1−c) sC26 + 2 1−c 2 (2C66 +C12 )4no2221−c (C11 +C22 ) + 2 1+c C12 + 4cs (C16 −C26 ) + 4 c −1 C66
Ĉ12 =4no
22Ĉ16 = s (1+c) C11 − s (1−c) C22 − 2−2c−4c C16 + 2+2c−4c C26 − 2cs (C12 +2C66 )
4n
o
(13)
222
Ĉ22 = −(1−c) C11 + (1+c) C22 + 4s (1−c) C16 + 4s (1+c) C26 + 2 1−c (2C66 +C12 )
4n

o
22
Ĉ26 = s (1−c) C11 − s (1+c) C22 + 2csC12 + 2+2c−4c C16 − 2−2c−4c C26 − 4sC66
4no2221−c (C22 +C11 ) − 2 1−c C12 + 4cs (C16 −C26 ) + 4c C66
Ĉ66 =4.
Ĉ11 =

2nd order polynomial, no trigonometric functions!!

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Numerical examples / Short cantilever

hc

∂Du

Fixed design domain D∂Dttyxwc

Height hc = 1, Width hc = 3, Mesh size d = 0.02, Filter size d = 0.1d,24/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Short cantilever / Result 1/2
Volume fraction = 0.5
Continuous, C=4.17

Discrete 4 options C=4.98(0, π/4, π/2, 3π/4)Discrete 4 options C=5.52(0, π/4, π/2, 3π/4) + π/625/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Short cantilever / Result 2/2
Volume fraction = 0.75
Continuous, C=3.09

Discrete 4 options C=3.78(0, π/4, π/2, 3π/4)Discrete 4 options C=4.12(0, π/4, π/2, 3π/4) + π/626/36

27/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Numerical examples / MBB beam

∂DtthmFixed design domain D∂Dr∂Duyxwm

Height hc = 1, Width hc = 5, Mesh size d = 0.02, Filter size d = 0.1d,

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

MBB beam / Result
Isotropic, Volume fraction 0.25

Anisotropic, Volume fraction 0.25

Isotropic, Volume fraction 0.5

Anisotropic, Volume fraction 0.5

Isotropic, Volume fraction 0.75

Anisotropic, Volume fraction 0.75

28/36

29/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Numerical examples / Multiload beam

∂DtthmFixed design domain D∂Dt∂Duy∂Drtx

wm

Height hc = 1, Width hc = 5, Mesh size d = 0.02, Filter size d = 0.1d,

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Multiload beam / Result, Multiload caseIsotropic, Volume fraction 0.25

Anisotropic, Volume fraction 0.25

<4>Isotropic,Volume fraction 0.5

Anisotropic, Volume fraction 0.5

Isotropic, Volume fraction 0.75

Anisotropic, Volume fraction 0.75

30/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Numerical examples / Holed panel problemofofrf

Hole radius rf = 0.4

∂Dulf

rf

∂Dt
Fixed design
domain Dofyx

Side length lf = 2oflftHole offset of = 0.6
Mesh size d = 0.02
π π 3πAngle set (0, , , π) +4 2 4831/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

Holed panel problem / Results
Radius: 0.2d

Radius: 0.4d

Radius:0.1d

Radius: 0.2d

Radius:0.4d

Discrete

Continuous

Radius:0.1d

32/36

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, Japan

33/36

Future work
Toolpath design• Toolpath generation algorithm based on PDE is being developed

3D design
1

3D design of density + Layer by layer 2D orientation design

2

Orientation design on 3D surfaces• For current composite printers•34Full 3D orientation design• For Full 3D composite printerOpen platform composite printer is desired!!
• Open driver development, open G-code specification (or at least SDK)
• Open hardware based research platform

The First Workshop on Composite 2.0 / Composite Printing, Novbember 27, 2015 Tokyo, JapanConclusion• A method for simultaneous optimization of topology and orientation isproposed.
• Capable for both continuous design and discrete design• Comparison to orientation angle design• Compatible with density filtering• Less local optima issue• Comparison to DMO• Continuous design• Design variable is 2 per design point in 2D problem regardless to numberof options• DMO offers more degree of freedom in material choice• Comparison to FMO• Physically feasible design
General topology optimization method with continuous and discrete orientation
design using isoparametric projection
DOI: 10.1002/nme.4799

34/36