- Email: [email protected]
Experimental characterization and modeling of the mechanical behavior of brittle 3D printed food
Journal Pre-proof Experimental characterization and modeling of the mechanical behavior of brittle 3D printed food N. Jonkers, J.A.W. van Dommelen, M.G.D. Geers
PII: DOI: Reference:
S0260-8774(20)30038-8 https://doi.org/10.1016/j.jfoodeng.2020.109941 JFOE 109941
To appear in:
Journal of Food Engineering
Received date : 2 July 2019 Revised date : 10 January 2020 Accepted date : 21 January 2020 Please cite this article as: N. Jonkers, J.A.W. van Dommelen and M.G.D. Geers, Experimental characterization and modeling of the mechanical behavior of brittle 3D printed food. Journal of Food Engineering (2020), doi: https://doi.org/10.1016/j.jfoodeng.2020.109941. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
Journal Pre-proof
*Highlights (for review)
Highlights The mechanical behavior and brittle failure of 3D printed samples is characterized A model to predict mechanical properties of brittle 3D printed food is presented Sample imperfections are taken into account to identify the model parameters
Jo
urn al P
re-
pro
of
• • •
Journal Pre-proof
of
*Conflict if Interest form
Jo
urn al P
re-
pro
The authors of this paper declare that they have no competing interests
Journal Pre-proof
*Credit Author Statement
Author statement Nicky Jonkers: Conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing – original draft.
of
Hans van Dommelen: conceptualization, funding acquisition, methodology, project administration, supervision, writing – review & editing.
Jo
urn al P
re-
pro
Marc Geers: conceptualization, funding acquisition, methodology, resources, supervision, writing – review & editing.
Journal Pre-proof
*Manuscript Click here to view linked References
3
N. Jonkersa , J.A.W. van Dommelena,∗, M.G.D. Geersa a Department
of Mechanical Engineering, Eindhoven University of Technology, The Netherlands.
pro
4
of
2
Experimental characterization and modeling of the mechanical behavior of brittle 3D printed food
1
Abstract
6
3D printing is a unique manufacturing method that enables food customization. The development of a
7
modeling framework to predict mechanical properties of food products is an invaluable tool in such a cus-
8
tomization process. To set up this framework, 3D printed samples are mechanically characterized by means
9
of compression testing. The observed phenomena are captured in a constitutive model that describes the
10
large deformation behavior and the brittle failure of the material. Due to the rough contact surface of 3D
11
printed samples, spatial homogeneity is lost and parameter identification is rendered not straightforward. To
12
incorporate this non-uniformity, the model is implemented in a finite element package. Simulations reveal
13
the influence of this geometrical effect, allowing to identify the model parameters by which the mechanical
14
behavior of the material is adequately described.
urn al P
re-
5
15
Keywords: 3D food printing, food design, material characterization, constitutive model, finite element
16
simulations
17
1. Introduction
18
Because the taste and appreciation of food is typically subjective and personal, there is a commercial driver
19
towards the customization of food (Sun et al., 2015; Godoi et al., 2016). Another important reason for food
20
customization is health (Lipton, 2017). 3D printing is a modern manufacturing method in which a product
21
is built layer by layer, which is especially suited for processing customized products. Besides creating
22
attractive, unique dishes, 3D food printing can be used to tailor properties of products. A good example
23
of this application is the texture-modified food for elderly people that require specific diets (Aguilera and
24
Park, 2016).
The technique that is currently most used for food printing is based on extrusion, as used to process,
26
for example, chocolate (Mantihal et al., 2017; Lanaro et al., 2017), breakfast spreads (Hamilton et al., 2018)
27
and cheese (Le Tohic et al., 2018). Although the presented work is applicable to all 3D printing techniques,
28
29
30
Jo
25
Selective Laser Sintering (SLS), which is based on the sintering of powder particles (Noort et al., 2017; Godoi et al., 2019), is chosen. This technique has the advantage that there is no need for support structures which are not part of the design but only needed for the process. In SLS, the unsintered powder provides ∗ Corresponding author Preprint submitted to Journal of Food Engineering Email address: [email protected] (J.A.W. van Dommelen)
January 10, 2020
Journal Pre-proof
the necessary support, which makes it suitable for complex geometries. Because the SLS process is more
32
complicated than the extrusion based method, SLS is not frequently used for food printing yet (Sun et al.,
33
2015).
µm
mm
pro
of
31
cm
Figure 1: Different length scales are important to understand the mechanics of food. Effective properties at the millimeter scale originate from the micrometer scale whereas the dimensions of the product are at the centimeter scale.
Le Tohic et al. (2018) showed that the mechanical properties of the printed material may be altered by
35
adjusting the printing process. To make use of the power of 3D printing to manufacture food with specific
36
designed or optimized properties, a modeling approach is beneficial. This requires in-depth understanding
37
of the mechanical behavior of the printed material and the ability to incorporate this knowledge into a
38
constitutive model. The mechanics of food is an inherently multi-scale problem (Liu and Scanlon, 2003;
39
Guessasma et al., 2011), as schematically illustrated in Figure 1. The mechanical properties result from the
40
microstructure, which is a result of the process conditions and the incorporated ingredients. The effect of the
41
microstructure on mechanical properties is studied by Attenburrow et al. (1989) for sponge cake, Afoakwa
42
et al. (2009) for chocolate, Sozer et al. (2011) for chocolate cake, Agbisit et al. (2007) and Warburton et al.
43
(1990) for cornstarch extrudates, Robin et al. (2011) for wheat flour based foams and Keetels et al. (1996),
44
Scanlon and Zghal (2001) and Zghal et al. (2002) for bread.
urn al P
re-
34
To test textural properties of food, the Texture Profile Analysis (TPA) test (Friedman et al., 1963) is
46
often used (see e.g. Rosenthal, 2010). A complication of the TPA test is that properties are extracted from
47
a force-displacement curve taken from an experiment with a probe that may have different geometries and a
48
sample with a non-standardized geometry. The experimental results are influenced by the chosen geometries
49
and settings of the TPA test, while the intrinsic mechanical properties are independent of the geometry
50
(Peleg, 2019). 3D printing is used to print specific structures, so here geometry matters. To be able to
51
predict the mechanical properties of a printed product (i.e. the effect of geometry, material, loading/contact
52
conditions), the starting point should be a constitutive model that captures the intrinsic mechanical behavior
53
of the material. Using a constitutive model in combination with the finite element method (FEM) can result
54
in quantitative predictions of properties for complex geometries and loading conditions (see e.g. Wismans
55
et al., 2009) and thereby the textural response. This is an essential step in enabling food customization.
56
57
58
Jo
45
Only a limited amount of papers address the constitutive modeling of food. Guessasma et al. (2011)
and Vancauwenberghe et al. (2018) used a linear elastic model in combination with FEM to predict effective properties of food products. Del Nobile et al. (2007) focused on the viscoelastic response and used a gen2
Journal Pre-proof
eralized Maxwell model to fit stress-relaxation experiments. Shirmohammadi et al. (2015) considered large
60
strains with a constitutive model of pumpkin peel and flesh that is based on an exponential strain energy
61
function. The predictive capability of the constitutive model was limited as different parameters are needed
62
to describe different loading conditions, so an accurate description of complex loading conditions does not
63
seem to be possible. Liu and Scanlon (2003) considered complex deformations, specifically indentation of
64
bread crumb, using a hyperelastic model based on the Ogden strain energy function.
pro
of
59
None of these models consider plastic deformation or have the ability to describe brittle failure of food
66
products. The aim of the current study is to obtain a constitutive model that describes the large-strain
67
material behavior of 3D printed starch-based foods based on careful experimental research. An important
68
aspect in modeling is the considered length scale, referring to the multi-scale problem in Figure 1. Here, the
69
focus is on the millimeter length scale at which the effective properties originate from the microstructure
70
which depends on the process parameters. As only the effective properties are considered, the microstructural
71
features are not explicitly taken into account. To quantify various aspects of the mechanical behavior, a set
72
of experiments is performed. The presented experimental methods are suitable for characterizing food in
73
general and are not restricted to the application of 3D printing. Based on these experiments, a 3D constitutive
74
model is formulated to describe the observed elasto-viscoplastic and damaging behavior of the material. An
75
important limitation, considering the layer-by-layer manufacturing process, of the presented model is that
76
material properties are assumed to be isotropic. The model is implemented in a finite element framework
77
to identify the material parameters, considering the experimental conditions. In this finite element study at
78
the sample scale, the effect of geometry on the mechanical response clearly emerges.
79
2. Experimental methods
80
2.1. Material and processing
81
The considered material is a mixture of 50% native wheat starch (Excelsior, Avebe), 40% maltodextrin
82
(Glucidex MD29, Roquette) and 10% palm oil powder (Admul PO 58, Kerry). This composition is similar
83
to the one that was used by Noort et al. (2017) and is chosen to have a structural component (native wheat
84
starch) in combination with a binder (maltodextrin and palm oil), which makes it suitable for SLS. The
85
powder is Selective Laser Sintered in an EOS P380 machine at room temperature. A pre-conditioning step
86
is required because the material is sensitive to environmental changes in moisture content. The moisture
87
content influences the binding between particles and the flowability of the powder. The powder is therefore
88
conditioned in a climate chamber at 50% RH and a temperature of 23 ◦ C for at least 16 hours. A layer
urn al P
Jo
89
90
re-
65
thickness of 0.3 mm, writing speed of 100 mm/s, line distance of 0.1 mm and laser power of 1.9 W are used. The scanning direction of the laser is alternated 90◦ in every subsequent layer to obtain the same
3
Journal Pre-proof
Hd
of
R
pro
H
D
Figure 2: Start of the compression test: the contact area Figure 3: Illustration of a sample with height H, diameter is small. D, roughness R and a damaged zone Hd .
91
mechanical properties in the two principal directions. The resulting printed cylinders have an average height
92
and standard deviation of H = 4.32 ± 0.07 mm and diameter D = 3.70 ± 0.05 mm. After sintering, the samples are baked in an oven at 180 ◦ C for 10 minutes.
94
2.2. Mechanical testing
95
Compression tests are performed with a rheometer (Discovery HR-3, TA instruments), having a flat bottom
96
plate and flat top plate (Petisco-Ferrero et al., 2019). The axial displacement is controlled, while measuring
97
the axial force, to obtain the force-displacement curves of the compressed samples. The printed samples are
98
again conditioned at 50% RH and a temperature of 23 ◦ C for at least 16 hours prior to testing. The samples
99
are not perfectly flat, see Figure 2, which results in imperfect contact between the non-flat sample and the
100
flat compression plate. The average height of the surface irregularities R, illustrated in Figure 3 and referred
101
to as roughness, is 0.31 mm and has a standard deviation of 0.09 mm. As a consequence, the contact area
102
varies during the initial part of the test.
urn al P
re-
93
103
To translate the force to a true stress and the displacement to true strain, the actual deformed cross-
104
sectional area is required because of the large applied strain. In this research however, the engineering
105
stress is used because it is not possible to correctly quantify the dynamic contact area between sample
106
and compression plate. Moreover, the samples are highly compressible and there is no visible deformation
107
of the cross-sectional area, suggesting that the engineering stress is close to the true stress. The initial
108
cross-sectional area A0 is therefore taken to calculate the engineering stress σeng , according to
σeng = F/A0 .
The true strain ε is calculated from the displacement u and sample height H (see e.g. Fenner, 1999) using
Jo
109
(1)
ε = ln(u/H + 1).
4
(2)
Journal Pre-proof
An exponentially decreasing velocity profile is prescribed such that the applied true strain rate ε˙ is constant
111
during the test. In addition, cyclic loading tests are performed. The sample is loaded and subsequently
112
unloaded with a constant velocity of 20 µm/s, increasing the maximum displacement for each loading-
113
unloading cycle, see Figure 4. In between cycles, the sample is not loaded during an interval that spans a
114
duration of 10 times the loading time, allowing the sample to relax.
pro
of
110
0.9 0.8 0.7 0.6 0.5 0.4 0.3
re-
0.2 0.1 0
500
1000
1500
Figure 4: Applied displacement in cyclic loading.
Finally, to measure linear viscoelastic effects, stress relaxation experiments are performed. A strain
116
between 1% and 2% is applied and kept constant, while measuring the stress as a function of time. The
117
resulting viscoelastic response is described by a series of stresses σi and relaxation times ti , according to
urn al P
115
σeng (t) =
M X i=1
σi · exp(−t/ti ),
(3)
in which M is the number of used modes.
119
3. Constitutive model
120
3.1. Material characterization
121
The engineering stress as a function of compressive true strain of six compression tests is shown in Figure
122
5, for a strain rate of approximately 10−3 s−1 . The effect of the non-flat surface is clearly visible in the
123
first part of the curves, where the sample surface is flattened and the contact area is increasing. This causes
124
a gradual start-up region, which is sample geometry dependent. The different measurements are therefore
125
126
127
Jo
118
shifted based on their maximum stiffness using the correction procedure from Pelletier et al. (2007). At the maximum stress, crack formation is observed and from this point onwards the structure of the sample is disintegrating. In the stress-strain response, this is manifested as a decrease of the stress as a function of
5
Journal Pre-proof
128
strain, also known as strain softening. From Figure 5, it is observed that the stiffness and fracture stress are quite reproducible. The start-up region due to the non-flatness of the contact surface prevents an accurate direct measurement of the Young’s modulus. The mean tangent stiffness ∂σ/∂ε between a strain of 0.025
131
and 0.030 for the tests with a strain rate of 10−3 s−1 is 50.5 MPa, and the standard deviation is 5.3 MPa.
132
The differences in sample roughness are expected to be the main cause for the variations in stiffness. For
133
the fracture stress, a mean value of 2.8 MPa and a standard deviation of 0.4 MPa are obtained.
pro
of
129
130
4
3
1
0 -0.1
-0.05
re-
2
0
0.05
0.1
0.15
0.2
urn al P
Figure 5: Compression testing of cylindrical samples: each curve represents a different sample and is horizontally shifted based on maximum stiffness.
134
The compression tests are repeated at a true strain rate of approximately 10−4 s−1 to investigate the
135
rate-dependence of the fracture stress. In Figure 5, four measurements with this smaller strain rate are
136
compared to the compression tests with a strain rate of 10−3 s−1 . No clear indication that the fracture stress
137
may depend on the applied deformation rate is observed. In the compression tests, the tangent stiffness is increasing as a function of the applied strain. Image
139
analysis of the deformation, taken during the tests, indicates that this is not only due to flattening of the
140
rough surface. The increase in stiffness originates from the porosity of the samples resulting from the laser
141
sintering process. With increasing applied deformation, the pores are collapsing, the material becomes denser
142
and thus the stiffness increases. This is confirmed in the cyclic loading tests. The result of such a test is
143
presented in Figure 6a, in which the loading curves are shown. In between cycles, the sample is not loaded
144
during an interval that spans a duration of 10 times the loading time, allowing the sample to relax. The
145
sample does not recover its original height indicating the presence of plastic deformation. The tangent
146
stiffness is determined for each cycle by the tangents shown in Figure 6a. The measurement is performed
147
three times resulting in the stiffness as a function of the applied strain, depicted in Figure 6b. The stiffness
148
149
150
Jo
138
is clearly increasing until the fracture stress is reached. After this, both stiffness and maximum stress of the damaged sample are decreasing. Finally, stress relaxation experiments are performed on three samples. Figure 7 shows the decay in stress 6
3
140
2.5
120
of
Journal Pre-proof
100
2
80 1.5 60 1
pro
40
0.5
20
0
0
0
0.05
0.1
0.15
0.2
0
(a)
0.05
0.1
0.15
0.2
(b)
re-
Figure 6: (a) Loading curves of a single sample (black) revealing tangent lines (blue) in cyclic loading and (b) the apparent stiffness for three different samples, determined from the tangents.
151
during such a measurement, indicating a viscoelastic material response. The response is described using four
152
modes (M = 4) in Equation 3 and this prediction is included in Figure 7. 0.3
urn al P
0.25 0.2
0.15
0.1
0.05
0 100
101
102
103
104
Figure 7: Stress relaxation measurement and model prediction using Equation 3.
To describe the material behavior, a constitutive model is needed that relates stress and strain as observed
154
in the experiments. The constitutive model is formulated in 3D in a tensorial form. The formulation allows
155
for high nonlinearity and large deformations. The model is represented by a single Maxwell element including
156
a damage model to take the brittle failure into account. The mechanical analogue of the Maxwell element
157
consists of a spring with shear modulus Gi and a dashpot with shear viscosity ηi = Gi ti as illustrated in
158
Figure 8.
Jo
153
7
Journal Pre-proof
Gi
of
ηi
Figure 8: A Maxwell element represented by a spring with shear modulus Gi and dashpot with viscosity ηi .
3.2. Tensor notations
160
In the derivation of the model, some notations are used which are briefly described here. A second-order
161
tensor A in the three-dimensional space is written as:
pro
159
A = Aij ~ei~ej ,
162
(4)
with Cartesian vector basis {~e1 , ~e2 , ~e3 } and summation over repeated indices. The index notation of this
tensor is Aij and the transpose of a tensor A, denoted as AT , equals Aji . The dot product or inner product
164
of two tensors A and B results in a second-order tensor C which is here defined as:
re-
163
C = A · B, with Cik = Aij Bjk ,
166
167
and the double dot product in a scalar c:
urn al P
165
(5)
c = A : B = Aij Bji .
(6)
tr(A) = Aii
(7)
1 Ad = A − tr(A)I, 3
(8)
The trace of a tensor A is
and the deviatoric part is obtained as
168
with the unit tensor I = ~ei~ei . The notation A˙ represents a differentiation with respect to time.
169
3.3. Kinematics
170
~ 0 ~u)T transforms the initial configuration to a deformed configThe deformation gradient tensor F = I + (∇
172
173
~ 0 denotes the gradient operator with respect to the initial configuration and ~u the displacement uration. ∇
Jo
171
vector. The starting point of the model is the multiplicative decomposition of the deformation gradient tensor in an elastic deformation Fe and a plastic deformation Fp (Lee, 1969), given by
8
Journal Pre-proof
Cc
Fp
C^
pro
Ci
of
F
Fe
Figure 9: Multiplicative decomposition of the deformation gradient tensor.
(9)
re-
F = Fe · Fp .
174
This is illustrated in Figure 9. The intermediate configuration Cˆ is obtained by deforming the initial
175
configuration Ci with Fp . The current configuration Cc is reached by applying Fe to this intermediate
176
configuration.
The velocity gradient tensor L is obtained as
urn al P
177
L = F˙ · F −1 = F˙e · Fe−1 + Fe · F˙p · Fp−1 · Fe−1 .
178
The plastic part of the velocity gradient tensor is defined in the intermediate configuration as ˆ p = F˙p · F −1 , L p
179
(10)
(11)
ˆ p in this intermediate configuration is and the plastic rate of deformation tensor D ˆ p = 1 (L ˆp + L ˆ T ). D p 2
(12)
180
Assuming that plastic deformation is isochoric, the volume change ratio J is calculated from the determinant
181
of Fe :
J = det(F ) = det(Fe ).
183
The assumption is made that the plastic part of the deformation does not have a rigid body rotation, i.e. Rp = I, with Rp the plastic rotation tensor. The plastic deformation gradient tensor is therefore given by
Jo
182
(13)
Fp = Rp · Up = Up , 9
(14)
Journal Pre-proof
in which Up is the plastic right stretch tensor.
185
3.4. Stress calculation
186
The Cauchy stress σ is the stress in the current configuration Cc . The total stress is split into two parts, a
187
hydrostatic part σ h and a deviatoric part σ d , according to
pro
of
184
σ = σh + σd .
(15)
188
The deviatoric part is obtained using a hyperelastic Ogden model (Ogden, 1972, 1984). This model is chosen
189
for its ability to model strongly nonlinear elastic behavior, in this case an increase in stiffness in compression,
190
see Figure 6b. The deviatoric stress is calculated from the isochoric elastic left stretch tensor V˜e : N 1X gi (V˜eai )d , J i=1
re-
σd =
(16)
191
in which gi and ai are the shear modulus and the exponent of Ogden term i. The total shear modulus G is
192
calculated from these parameters as
N X gi ai
urn al P
G=
i=1
2
.
(17)
193
The pair of parameters (gi , ai ) may also have negative values to influence the mechanical response in com-
194
pression, while positive values are determining the response in tension. By having multiple terms, it is
195
possible to address both tension and compression with a single set of model parameters. The left elastic
196
stretch tensor Ve is determined from Fe ,
Ve =
197
q
Fe · FeT ,
(18)
and the isochoric elastic left stretch tensor to the power ai using a
i V˜eai = J − 3 Veai .
(19)
The increase in shear modulus, due to the decreasing porosity during compression, is accompanied with
199
an increase in bulk modulus. The factor ai , which is numerically taking care for this stiffness increase, is
200
therefore included in the hydrostatic stress (Nedjar, 2016). The hydrostatic part of the stress is given by
Jo
198
σh =
N X Kgi i=1
2G
(J
2ai 3
10
−1
−J
−ai 3
−1
)I,
(20)
Journal Pre-proof
with K the bulk modulus.
202
3.5. Plasticity and damage
203
ˆ p is calculated from the stress-dependent viscosity function η(¯ The plastic rate of deformation tensor D τ)
204
and the deviatoric stress using ˆp = D
pro
of
201
1 ˆd S . 2η(¯ τ)
(21)
205
The deviatoric stress in the intermediate configuration is expressed in the elastic second Piola-Kirchhoff-like
206
stress tensor Sˆd which is calculated from the Cauchy stress:
207
The viscosity is calculated as
re-
Sˆd = (JFe−1 · σ · Fe−T )d .
τ¯ . γ˙ p
η(¯ τ) =
(22)
(23)
The plastic strain rate γ˙ p , parameterized by the reference shear rate γ˙ 0 , the power m and the characteristic
209
stress τc , is given by
urn al P
208
γ˙ p = γ˙ 0
210
τ¯ τc
m
,
(24)
1 ˆd ˆd S :S . 2
(25)
and the equivalent shear stress is calculated as
τ¯ =
r
211
The softening part of the stress-strain curve is described by a damage law. The evolution of the damage
212
parameter ξ, which is initially zero, is calculated from the plastic strain rate by ξ˙ =
ξ 1− dγ˙ p , ξ∞
(26)
in which d and ξ∞ are parameters describing the rate at which the damage develops and the limit value for the damage, respectively. The damage parameter ξ affects the current shear moduli gi and the characteristic
215
stress τc , such that
216
Jo
213
214
τc = τ0 (1 − ξ) and gi = g0,i (1 − ξ)
in which τ0 and g0,i are the initial, undamaged variables. 11
(27)
Journal Pre-proof
3.6. Time integration
218
The constitutive equations are integrated from time tn to tn+1 using an implicit integration scheme (Amiri-
of
217
Rad et al., 2019). The applied deformation is given by F and the elastic part of this deformation gradient
220
tensor is used to calculate the total stress. One determines the amount of plasticity, starting from the plastic
221
right Cauchy-Green deformation tensor Cp = FpT · Fp , and then looks to the remaining part of F , which
222
gives Fe . To this purpose, a Newton-Raphson scheme is used, in which first an initial guess for Cp,n+1 is
223
taken. From Cp,n+1 the plastic right stretch tensor Up,n+1 is determined:
Up,n+1 =
pro
219
p Cp,n+1 .
(28)
Using the polar decomposition of the deformation gradient, Equation 14, the elastic deformation gradient
225
tensor is recovered as
re-
224
−1 Fe,n+1 = Fn+1 · Up,n+1 .
(29)
d ˆ p,n+1 is found from the flow rule. The time The stress σn+1 is obtained with this Fe,n+1 and subsequently D
227
derivative of Cp,n+1 can now be calculated using
urn al P
226
T T C˙ p,n+1 = F˙p,n+1 · Fp,n+1 + Fp,n+1 · F˙p,n+1
=
228
T 2Up,n+1
(30)
ˆ p,n+1 · Up,n+1 . ·D
Using C˙ p,n+1 , the tensor Cp,n+1 is recomputed using the backward Euler method: Cp,n+1 = Cp,n + ∆tC˙ p,n+1 ,
(31)
in which ∆t is the time step size. The new value of Cp,n+1 enters the Newton-Raphson scheme again,
230
until convergence in the backward Euler scheme is obtained. The Frobenius norm of the residue (difference
231
between initial guess and recomputed tensor) is used as convergence criteria, repeating the time integration
232
until a norm smaller than 10−10 is obtained. The damage parameter ξ is updated at the end of each increment
233
by the explicit time integration of Equation 26.
234
3.7. Implementation
235
236
Jo
229
The constitutive model is implemented in the finite element package MSC.Marc using the user subroutine HYPELA2.
12
Journal Pre-proof
4. Results
238
4.1. Finite element model
239
The model performance is first assessed through the homogeneous deformation of a cube with sides of
240
1 mm. The element size is varied from a single element up to cubic elements with sides of 62.5 µm to
241
investigate mesh-sensitivity of the model. 3D linear hexahedral elements are used. The boundary conditions
242
are schematically presented in Figure 10a. Uniaxial compression simulations are performed by prescribing
243
the displacement at the top surface and fixing the bottom in z-direction. Moreover, the nodes on the left
244
face are fixed in x-direction and in the front face in y-direction, corresponding to symmetry conditions. An
245
initial time step size of 1.0 second is taken, which turned out to be sufficiently small (i.e. the final solution
246
is not dependent on it anymore), and is decreased during softening as a result of the high nonlinearity of
247
Equation 24.
re-
pro
of
237
248
Figure 2 reveals that the sample is not flat resulting in a small tangent stiffness at the start of the
249
compression test. To identify the effect of the geometry, resulting in inhomogeneous deformation, on the effective macroscopic response, finite element simulations are performed with a simplified geometry that represents the actual samples, depicted in Figure 10b. The average height H=4.32 mm, diameter D=3.70
252
mm and roughness R=0.31 mm, illustrated in Figure 3, of the samples are taken to model the geometry
253
with 3040 3D linear hexahedral elements. The bottom nodes are fixed in z-direction and the rotation of
254
the geometry is constrained. The compression is applied in the z-direction by contact with a plate with
255
a prescribed velocity such that the true strain rate is constant. The initial time step size is equal to 0.2
256
seconds. The effect of contact friction was tested with a Coulomb friction model, which turned out to be
257
small compared to the sample roughness and is therefore further neglected.
urn al P
250
251
y
z
z
x
(a)
y
x
(b)
Jo
Figure 10: Finite element model: (a) boundary conditions in homogeneous deformation and (b) representative geometry. Red arrows indicate fixed boundary condition and black arrows represent the applied compression.
13
Journal Pre-proof
4.2. Parameter identification
259
The model parameters and their identified values are listed in Table 1. From the experiments, a stiffness of
260
40 MPa is estimated at the moment the top surface was flattened. This value is assumed to be equal to the
261
Young’s modulus E. The flattening makes it hard to determine a Poisson’s ratio ν, which is estimated as
262
0.15. From these values, the shear modulus and bulk modulus are calculated using
G=
pro
of
258
E E , and K = , respectively. 2(1 + ν) 3(1 − 2ν)
(32)
The hyperelastic part of the model consists of two Ogden terms each with a shear modulus g0,i and
264
exponent ai . One term is dominant in tension (positive values) and the other is dominant in compression
265
(negative values). Because of the brittleness of the samples, it is not possible to perform tensile tests on this
266
particular material. A standard value corresponding to a1 = 2 is therefore adopted, which makes the term
267
equal to a Neo-Hookean model and g0,1 = 1, which has an arbitrary small value for the product a1 g0,1 as a
268
result.
re-
263
Table 1: Model parameters.
K [MPa] 19.0
g0,1 [MPa] 1.0
urn al P
G [MPa] 17.4 γ˙ 0 [s−1 ] 1.0·10−3
τ0 [MPa] 2.6
m [-] 15
5
a1 [-] 2.0
g0,2 [MPa] -0.73
d [-] 5.0
ξ∞ [-] 0.85
a2 [-] -45
ξ 0.85
4 3 2
0
0
0.05
0.1
0.15
y
m = 15
1
z
0.2
m = 50 x
0.0
Figure 11: Regularization of damage localization by decreasing m in a uniaxial compression test on a cube, showing (a) the macroscopic response and (b) the local damage after applying a strain of 0.13.
270
271
272
To uniquely determine the three parameters that describe the plastic strain rate as given in the specific
Jo
269
form in Equation 24, the value of one parameter can be chosen arbitrary. The reference shear rate γ˙ 0 is therefore taken equal to the experimentally applied strain rate. An effect of the damage is that it may localize, which can have mesh-dependent results as a consequence. This is regularized by introducing a 14
Journal Pre-proof
273
rate-dependent material behavior (Needleman, 1988). The rate-dependence of the material is captured by the model parameter m. Figure 11a shows the macroscopic response of a uniaxial compression test on the cube with sides of 62.5 µm and Figure 11b depicts the damage parameter after applying a compressive strain
276
of 0.13, for different values of m. For this simple test case, choosing m = 30 is sufficient to regularize the
of
274
275
damage. However, it is found that for a more complex geometry, a value of m = 15 is necessary to provide
278
a unique solution which is independent of the mesh size and time step size. The introduced rate-dependent
pro
277
279
behavior may be in contradiction with the rate-independent behavior observed in the experiments, as shown
280
in Figure 5. The chosen value m = 15 is sufficiently small to regularize the damage without introducing a
281
significant effect of the deformation rate on the fracture stress.
The five remaining parameters are g0,2 , a2 , τ0 , d and ξ∞ . Considering Equation 17, only a2 has to
283
be determined to also recover g0,2 , knowing the parameters of the first Ogden (i=1) term and the total
284
shear modulus. The number of parameters to determine is therefore reduced to four. The effect of these
285
parameters on the intrinsic material response, in terms of true stress (deformed area), is presented in Figure
286
12. As a reference, the parameter set given in Table 1 is used. Parameter a2 affects the stiffness increase as
287
a function of the applied compressive strain, and τ0 defines the fracture stress. The parameter d is identified
288
from the rate of damage, i.e. the slope of the stress-strain curve in the damaged region, and ξ∞ determines
289
the fully damaged state and corresponds to the ratio between the maximum stress and the plateau level at
290
large strains.
urn al P
re-
282
291
The main issue in identifying the model parameters is the local damage that is resulting from the non-
292
flat top surface. Figure 13 shows the damage parameter in a cross section of the sample after applying a
293
strain equal to 0.05. Locally, the material is already severely damaged, while a flat sample would be still
294
undamaged. The consequences are shown in Figure 14 for different surface topologies. The roughness R
295
is varied and different surface profiles are considered. In comparison with a flat sample, the stiffness and
296
fracture stress decrease and softening becomes more severe. Moreover, it is concluded that the influence of
297
R is more important than the surface topology. This confirms the choice for the definition of R and the
298
representative geometry.
To complete the identification, the influence of the different parameters on the stress-strain response of the
300
representative geometry is shown in Figure 15 and assessed through their influence on the intrinsic response,
301
depicted in Figure 12. The main difference is that the intrinsic response reveals a clear distinction between
302
elastic and plastic deformation. The parameters that correspond to plasticity do not affect the elastic part.
303
Adding roughness has the result that local plasticity is influencing the macroscopic stress-strain response,
304
including the deformation before the fracture stress is reached. Compared to Figure 12, the effect of a2 is
305
306
Jo
299
similar, but τ0 and d are clearly affecting the stiffness with increasing strain. Moreover, the fracture stress decreases with increasing d and the effect of d on the softening slope is fading. The softening is dominated
15
5
4
4
3
3
2
2
pro
5
of
Journal Pre-proof
1
1
0
0
0
0.05
0.1
0.15
0.2
0
5
0.05
0.1
0.15
0.2
0.1
0.15
0.2
5
4
4 3
2
re-
3
2
1
1
0
0
0
0.05
0.1
0.15
0.2
0
0.05
urn al P
Figure 12: Main characteristics of the model, showing the influence of a2 , τ0 , d and ξ∞ on the intrinsic mechanical behavior for a constant strain rate ε˙ = 10−3 s−1 . Note that for different values of d all curves will eventually reach the same plateau stress.
307
by geometrical effects and the material softening becomes less prominent. A new effect emerges, which are
308
the oscillations of the curve for larger values of d, induced by the many geometrical imperfections. The finally identified parameters are given in Table 1 and determined using the following characterization procedure. First, ξ∞ is determined because it is the only parameter that determines the ratio between the
311
fracture stress and the plateau stress after fracture. Moreover, ξ∞ is determining the depth of damage, Hd in
312
Figure 3. Next, d is chosen to approximate the jaggedness of the experimental results. Finally, τ0 is defined
313
to match the fracture stress, whereas a2 covers the necessary stiffness increase. Given the experimental
314
scatter, the obtained set of parameters adequately describes the overall experimentally obtained mechanical
315
behavior, as seen in Figure 16. Comparing the damaged zones in the experiments and simulation, Figure 17,
316
it is obvious that the simulation is also predicting this zone well. In comparison with the experiment, the
317
predicted damaged is more smooth. This difference is due to the fact that local heterogeneities and defects,
318
which may cause localization of damage, are not taken into account.
Jo
309
310
16
Journal Pre-proof
ξ
4.5
0.85
of
4 3.5 3 2.5 2
pro
1.5 1
0.5
0
0
0.0
0.05
0.1
0.15
0.2
re-
Figure 14: Comparison of the response (ε˙ = 10−3 s−1 ) Figure 13: Simulated damage in a cross section of the of samples with different roughness values. For each R sample after applying a strain of 0.05. value, three different surface profiles are simulated.
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
urn al P
0.5
0
0
0.05
0.1
0.15
0.5
0
0.2
3.5
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
0.05
0.1
0.15
0.2
0.05
0.1
0.15
0.2
0
0.05
0.1
0.15
0.2
3.5
3
0
0
0
Figure 15: Influence of model parameters on the macroscopic stress-strain response (ε˙ = 10−3 s−1 ) of the representative geometry.
5. Conclusion
320
An experimental procedure, based on compression testing, is presented and provides a solid basis for char-
321
322
Jo
319
acterizing the intrinsic mechanical properties and failure of 3D printed food. For the samples that are used in this research, the following observations are incorporated in a 3D constitutive model that allows for high
17
Journal Pre-proof
of
4
3
1
0 0
0.05
0.1
pro
2
0.15
0.2
Figure 16: Comparison between experiments and model prediction for ε˙ = 10−3 s−1 . ξ
re-
0.85
damaged
0.0
urn al P
Figure 17: Damaged zone at the surface in simulation and experiment after applying a strain of 0.2.
323
nonlinearity, inelastic deformation and damaging behavior. First, the stiffness increases as a function of the
324
applied strain caused by densification due to the porous microstructure. This is followed by the initiation of
325
brittle failure at the macroscopic fracture stress. During failure, softening is observed in which both stiffness
326
and stress decrease. After characterising some of the model parameters, a representative sample geometry,
327
taking the roughness of the samples into account, is modeled, allowing to identify the complete set of material
328
parameters. The presented model and characterization procedure are not only useful for the application of
329
3D food printing, but can be considered more generic and applicable to other food. A method to process or post-process flat samples without introducing any damage before testing would
331
enable a more accurate determination of the parameters. Characterization of the volumetric compressibility
332
of the material was not possible with the available experimental data and requires future investigation.
333
The value for the Poisson’s ratio is indicative only and the assumption that plasticity is volumetrically
334
incompressible may not be correct for this material. Although the constitutive model is able to distinct non-
335
linear behavior in both the compressive and tensile regime, the current characterization procedure focuses
336
on the more important compressive response only. An important future extension of the current model
337
338
Jo
330
could be the addition of anisotropy which results from the layer-by-layer manufacturing process. Despite these limitations, the model adequately predicts the mechanical response and damaged zone of the printed
18
Journal Pre-proof
samples. A powerful approach for predicting properties of specifically designed food products is therefore
340
obtained.
of
339
The mentioned aspects are important at the current length scale of modeling. This length scale only
342
considers effective properties originating from the microstructure. However, the microstructure of the 3D
343
printed food is highly heterogeneous. For future research, it would therefore be interesting to relate the
344
mechanical properties to the microstructure. The spatial distribution of material parameters can then be
345
incorporated by assigning different parameter values at different locations. Moreover, process parameters
346
(e.g. laser power) can be locally varied within a printed product. This enables tailoring the microstructure
347
(e.g. porosity, dimensions of voids and structural components) and, as a result, the mechanical properties
348
(e.g. Young’s modulus, fracture stress). Including process-microstructure-property relations should enrich
349
the current framework to become an invaluable tool in designing customized food.
350
Acknowledgements
351
The authors wish to thank TNO (Netherlands Organisation for Applied Scientific Research) for the financial
352
support, for the instructions on sample processing and for the access to the SLS machine. We also thank
353
Martijn Noort (Wageningen University and Research) for the preparation of the powder mixture, Ruth
354
Cardinaels (Eindhoven University of Technology) for the use of the rheometer and Susana Petisco Ferrero
355
(Eindhoven University of Technology) for the training on compression testing with the rheometer.
356
References
357
J. Sun, Z. Peng, W. Zhou, J. Y. H. Fuh, G. S. Hong, A. Chiu, A review on 3D printing for customized food
360
361
362
363
364
365
366
367
re-
urn al P
359
fabrication, Procedia Manufacturing 1 (2015) 308–319. F. C. Godoi, S. Prakash, B. R. Bhandari, 3d printing technologies applied for food design: Status and prospects, Journal of Food Engineering 179 (2016) 44–54. J. I. Lipton, Printable food: the technology and its application in human health, Current Opinion in Biotechnology 44 (2017) 198–201.
J. M. Aguilera, D. J. Park, Texture-modified foods for the elderly: Status, technology and opportunities, Trends in Food Science and Technology 57 (2016) 156–164.
Jo
358
pro
341
S. Mantihal, S. Prakash, F. C. Godoi, B. Bhandari, Optimization of chocolate 3D printing by correlating thermal and flow properties with 3D structure modeling, Innovative Food Science and Emerging Technologies 44 (2017) 21–29.
19
Journal Pre-proof
M. Lanaro, D. P. Forrestal, S. Scheurer, S. Liao, D. J. Slinger, S. K. Powell, M. A. Woodruff, 3D printing
369
complex chocolate objects: Platform design, optimization and evaluation, Journal of Food Engineering
370
215 (2017) 13–22.
372
C. A. Hamilton, G. Alici, M. in het Panhuis, 3D printing vegemite and marmite: Redefining breadboards, Journal of Food Engineering 220 (2018) 83–88.
pro
371
of
368
373
C. Le Tohic, J. J. O’Sullivan, K. P. Drapala, V. Chartrin, T. Chan, A. P. Morrison, J. P. Kerry, A. L.
374
Kelly, Effect of 3D printing on the structure and textural properties of processed cheese, Journal of Food
375
Engineering 220 (2018) 56–64.
378
379
380
381
382
383
384
385
M. W. Noort, J. V. Diaz, K. J. C. van Bommel, S. Renzetti, J. Henket, M. B. Hoppenbrouwers, Method for the production of an edible object using SLS. US2017/0266881 A1, 2017.
re-
377
F. C. Godoi, B. R. Bhandari, S. Prakash, M. Zhang, Fundamentals of 3D food printing and applications, Elsevier Ltd., London, 2019. doi:https://doi.org/10.1016/C2017-0-01591-4. Z. Liu, M. G. Scanlon, Predicting mechanical properties of bread crumb, Food and Bioproducts Processing 81 (2003) 224–238.
S. Guessasma, L. Chaunier, G. Della Valle, D. Lourdin, Mechanical modelling of cereal solid foods, Trends
urn al P
376
in Food Science and Technology 22 (2011) 142–153.
G. E. Attenburrow, R. M. Goodband, L. J. Taylor, P. J. Lillford, Structure, mechanics and texture of a food sponge, Journal of Cereal Science 9 (1989) 61–70.
386
E. O. Afoakwa, A. Paterson, M. Fowler, J. Vieira, Microstructure and mechanical properties related to
387
particle size distribution and composition in dark chocolate, International Journal of Food Science and
388
Technology 44 (2009) 111–119.
390
391
392
393
394
395
396
N. Sozer, H. Dogan, J. L. Kokini, Textural properties and their correlation to cell structure in porous food materials, Journal of Agricultural and Food Chemistry 59 (2011) 1498–1507. R. Agbisit, S. Alavi, E. Cheng, T. Herald, A. Trater, Relationships between microstructure and mechanical properties of cellular cornstarch extrudates, Journal of Texture Studies 38 (2007) 199–219. S. C. Warburton, A. M. Donald, A. C. Smith, The deformation of brittle starch foams, Journal of Materials Science 25 (1990) 4001–4007.
Jo
389
F. Robin, C. Dubois, D. Curti, H. P. Schuchmann, S. Palzer, Effect of wheat bran on the mechanical properties of extruded starchy foams, Food Research International 44 (2011) 2880–2888.
20
Journal Pre-proof
400
401
402
403
404
405
406
of
399
bread, Journal of Cereal Science 24 (1996) 15–26.
M. G. Scanlon, M. C. Zghal, Bread properties and crumb structure, Food Research International 34 (2001) 841–864.
M. C. Zghal, M. G. Scanlon, H. D. Sapirstein, Cellular structure of bread crumb and its influence on
pro
398
C. J. A. M. Keetels, K. A. Visser, T. van Vliet, A. Jurgens, P. Walstra, Structure and mechanics of starch
mechanical properties, Journal of Cereal Science 36 (2002) 167–176.
H. H. Friedman, J. E. Whitney, A. S. Szczesniak, The texturometer - a new instrument for objective texture measurement, Journal of Food Science 28 (1963) 390–396.
A. J. Rosenthal, Texture profile analysis - How important are the parameters?, Journal of Texture Studies 41 (2010) 672–684.
re-
397
407
M. Peleg, The instrumental texture profile analysis revisited, Journal of Texture Studies (2019) 1–7.
408
J. G. F. Wismans, J. A. W. van Dommelen, L. E. Govaert, H. E. H. Meijer, B. van Rietbergen, Computed
409
tomography-based modeling of structured polymers, Journal of Cellular Plastics 45 (2009) 157–179. V. Vancauwenberghe, M. A. Delele, J. Vanbiervliet, W. Aregawi, P. Verboven, J. Lammertyn, B. Nicola¨ı,
411
Model-based design and validation of food texture of 3D printed pectin-based food simulants, Journal of
412
Food Engineering 231 (2018) 72–82.
413
414
urn al P
410
M. A. Del Nobile, S. Chillo, A. Mentana, A. Baiano, Use of the generalized maxwell model for describing the stress relaxation behavior of solid-like foods, Journal of Food Engineering 78 (2007) 978–983.
415
M. Shirmohammadi, P. K. D. V. Yarlagadda, Y. T. Gu, A constitutive model for mechanical response
416
characterization of pumpkin peel and flesh tissues under tensile and compressive loadings, Journal of Food
417
Science and Technology 52 (2015) 4874–4884.
419
420
421
422
423
424
425
Z. Liu, M. G. Scanlon, Modelling indentation of bread crumb by finite element analysis, Biosystems Engineering 85 (2003) 477–484.
S. Petisco-Ferrero, R. Cardinaels, L. C. A. van Breemen, Miniaturized characterization of polymers: From synthesis to rheological and mechanical properties in 30 mg, Polymer 185 (2019) 121918. R. T. Fenner, Mechanics of Solids, CRC Press, Boca Raton, 1999.
Jo
418
C. G. N. Pelletier, E. C. A. Dekkers, L. E. Govaert, J. M. J. den Toonder, H. E. H. Meijer, The influence of indenter-surface misalignment on the results of instrumented indentation tests, Polymer Testing 26 (2007) 949–959.
21
Journal Pre-proof
E. H. Lee, Elastic-plastic deformation at finite strains, Journal of Applied Mechanics 36 (1969) 1–6.
427
R. W. Ogden, Large deformation isotropic elasticity: On the correlation of theory and experiment for
428
of
426
compressible rubberlike solids, Proceedings of the royal society A 328 (1972) 567–583.
R. W. Ogden, Non-linear elastic deformations, Chichester: E. Horwood; New York: Halsted Press, 1984.
430
B. Nedjar, On constitutive models of finite elasticity with possible zero apparent Poisson’s ratio, International
434
435
Viscoplastic Model for Short-Fiber Composites, Mechanics of Materials 137 (2019) 103141. A. Needleman, Material rate dependence and mesh sensitivity in localization problems, Computer Methods in Applied Mechanics and Engineering 67 (1988) 69–85.
re-
433
A. Amiri-Rad, L. V. Pastukhov, L. E. Govaert, J. A. W. van Dommelen, An Anisotropic Viscoelastic-
urn al P
432
Journal of Solids and Structures 91 (2016) 72–77.
Jo
431
pro
429
22
PII: DOI: Reference:
S0260-8774(20)30038-8 https://doi.org/10.1016/j.jfoodeng.2020.109941 JFOE 109941
To appear in:
Journal of Food Engineering
Received date : 2 July 2019 Revised date : 10 January 2020 Accepted date : 21 January 2020 Please cite this article as: N. Jonkers, J.A.W. van Dommelen and M.G.D. Geers, Experimental characterization and modeling of the mechanical behavior of brittle 3D printed food. Journal of Food Engineering (2020), doi: https://doi.org/10.1016/j.jfoodeng.2020.109941. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
Journal Pre-proof
*Highlights (for review)
Highlights The mechanical behavior and brittle failure of 3D printed samples is characterized A model to predict mechanical properties of brittle 3D printed food is presented Sample imperfections are taken into account to identify the model parameters
Jo
urn al P
re-
pro
of
• • •
Journal Pre-proof
of
*Conflict if Interest form
Jo
urn al P
re-
pro
The authors of this paper declare that they have no competing interests
Journal Pre-proof
*Credit Author Statement
Author statement Nicky Jonkers: Conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing – original draft.
of
Hans van Dommelen: conceptualization, funding acquisition, methodology, project administration, supervision, writing – review & editing.
Jo
urn al P
re-
pro
Marc Geers: conceptualization, funding acquisition, methodology, resources, supervision, writing – review & editing.
Journal Pre-proof
*Manuscript Click here to view linked References
3
N. Jonkersa , J.A.W. van Dommelena,∗, M.G.D. Geersa a Department
of Mechanical Engineering, Eindhoven University of Technology, The Netherlands.
pro
4
of
2
Experimental characterization and modeling of the mechanical behavior of brittle 3D printed food
1
Abstract
6
3D printing is a unique manufacturing method that enables food customization. The development of a
7
modeling framework to predict mechanical properties of food products is an invaluable tool in such a cus-
8
tomization process. To set up this framework, 3D printed samples are mechanically characterized by means
9
of compression testing. The observed phenomena are captured in a constitutive model that describes the
10
large deformation behavior and the brittle failure of the material. Due to the rough contact surface of 3D
11
printed samples, spatial homogeneity is lost and parameter identification is rendered not straightforward. To
12
incorporate this non-uniformity, the model is implemented in a finite element package. Simulations reveal
13
the influence of this geometrical effect, allowing to identify the model parameters by which the mechanical
14
behavior of the material is adequately described.
urn al P
re-
5
15
Keywords: 3D food printing, food design, material characterization, constitutive model, finite element
16
simulations
17
1. Introduction
18
Because the taste and appreciation of food is typically subjective and personal, there is a commercial driver
19
towards the customization of food (Sun et al., 2015; Godoi et al., 2016). Another important reason for food
20
customization is health (Lipton, 2017). 3D printing is a modern manufacturing method in which a product
21
is built layer by layer, which is especially suited for processing customized products. Besides creating
22
attractive, unique dishes, 3D food printing can be used to tailor properties of products. A good example
23
of this application is the texture-modified food for elderly people that require specific diets (Aguilera and
24
Park, 2016).
The technique that is currently most used for food printing is based on extrusion, as used to process,
26
for example, chocolate (Mantihal et al., 2017; Lanaro et al., 2017), breakfast spreads (Hamilton et al., 2018)
27
and cheese (Le Tohic et al., 2018). Although the presented work is applicable to all 3D printing techniques,
28
29
30
Jo
25
Selective Laser Sintering (SLS), which is based on the sintering of powder particles (Noort et al., 2017; Godoi et al., 2019), is chosen. This technique has the advantage that there is no need for support structures which are not part of the design but only needed for the process. In SLS, the unsintered powder provides ∗ Corresponding author Preprint submitted to Journal of Food Engineering Email address: [email protected] (J.A.W. van Dommelen)
January 10, 2020
Journal Pre-proof
the necessary support, which makes it suitable for complex geometries. Because the SLS process is more
32
complicated than the extrusion based method, SLS is not frequently used for food printing yet (Sun et al.,
33
2015).
µm
mm
pro
of
31
cm
Figure 1: Different length scales are important to understand the mechanics of food. Effective properties at the millimeter scale originate from the micrometer scale whereas the dimensions of the product are at the centimeter scale.
Le Tohic et al. (2018) showed that the mechanical properties of the printed material may be altered by
35
adjusting the printing process. To make use of the power of 3D printing to manufacture food with specific
36
designed or optimized properties, a modeling approach is beneficial. This requires in-depth understanding
37
of the mechanical behavior of the printed material and the ability to incorporate this knowledge into a
38
constitutive model. The mechanics of food is an inherently multi-scale problem (Liu and Scanlon, 2003;
39
Guessasma et al., 2011), as schematically illustrated in Figure 1. The mechanical properties result from the
40
microstructure, which is a result of the process conditions and the incorporated ingredients. The effect of the
41
microstructure on mechanical properties is studied by Attenburrow et al. (1989) for sponge cake, Afoakwa
42
et al. (2009) for chocolate, Sozer et al. (2011) for chocolate cake, Agbisit et al. (2007) and Warburton et al.
43
(1990) for cornstarch extrudates, Robin et al. (2011) for wheat flour based foams and Keetels et al. (1996),
44
Scanlon and Zghal (2001) and Zghal et al. (2002) for bread.
urn al P
re-
34
To test textural properties of food, the Texture Profile Analysis (TPA) test (Friedman et al., 1963) is
46
often used (see e.g. Rosenthal, 2010). A complication of the TPA test is that properties are extracted from
47
a force-displacement curve taken from an experiment with a probe that may have different geometries and a
48
sample with a non-standardized geometry. The experimental results are influenced by the chosen geometries
49
and settings of the TPA test, while the intrinsic mechanical properties are independent of the geometry
50
(Peleg, 2019). 3D printing is used to print specific structures, so here geometry matters. To be able to
51
predict the mechanical properties of a printed product (i.e. the effect of geometry, material, loading/contact
52
conditions), the starting point should be a constitutive model that captures the intrinsic mechanical behavior
53
of the material. Using a constitutive model in combination with the finite element method (FEM) can result
54
in quantitative predictions of properties for complex geometries and loading conditions (see e.g. Wismans
55
et al., 2009) and thereby the textural response. This is an essential step in enabling food customization.
56
57
58
Jo
45
Only a limited amount of papers address the constitutive modeling of food. Guessasma et al. (2011)
and Vancauwenberghe et al. (2018) used a linear elastic model in combination with FEM to predict effective properties of food products. Del Nobile et al. (2007) focused on the viscoelastic response and used a gen2
Journal Pre-proof
eralized Maxwell model to fit stress-relaxation experiments. Shirmohammadi et al. (2015) considered large
60
strains with a constitutive model of pumpkin peel and flesh that is based on an exponential strain energy
61
function. The predictive capability of the constitutive model was limited as different parameters are needed
62
to describe different loading conditions, so an accurate description of complex loading conditions does not
63
seem to be possible. Liu and Scanlon (2003) considered complex deformations, specifically indentation of
64
bread crumb, using a hyperelastic model based on the Ogden strain energy function.
pro
of
59
None of these models consider plastic deformation or have the ability to describe brittle failure of food
66
products. The aim of the current study is to obtain a constitutive model that describes the large-strain
67
material behavior of 3D printed starch-based foods based on careful experimental research. An important
68
aspect in modeling is the considered length scale, referring to the multi-scale problem in Figure 1. Here, the
69
focus is on the millimeter length scale at which the effective properties originate from the microstructure
70
which depends on the process parameters. As only the effective properties are considered, the microstructural
71
features are not explicitly taken into account. To quantify various aspects of the mechanical behavior, a set
72
of experiments is performed. The presented experimental methods are suitable for characterizing food in
73
general and are not restricted to the application of 3D printing. Based on these experiments, a 3D constitutive
74
model is formulated to describe the observed elasto-viscoplastic and damaging behavior of the material. An
75
important limitation, considering the layer-by-layer manufacturing process, of the presented model is that
76
material properties are assumed to be isotropic. The model is implemented in a finite element framework
77
to identify the material parameters, considering the experimental conditions. In this finite element study at
78
the sample scale, the effect of geometry on the mechanical response clearly emerges.
79
2. Experimental methods
80
2.1. Material and processing
81
The considered material is a mixture of 50% native wheat starch (Excelsior, Avebe), 40% maltodextrin
82
(Glucidex MD29, Roquette) and 10% palm oil powder (Admul PO 58, Kerry). This composition is similar
83
to the one that was used by Noort et al. (2017) and is chosen to have a structural component (native wheat
84
starch) in combination with a binder (maltodextrin and palm oil), which makes it suitable for SLS. The
85
powder is Selective Laser Sintered in an EOS P380 machine at room temperature. A pre-conditioning step
86
is required because the material is sensitive to environmental changes in moisture content. The moisture
87
content influences the binding between particles and the flowability of the powder. The powder is therefore
88
conditioned in a climate chamber at 50% RH and a temperature of 23 ◦ C for at least 16 hours. A layer
urn al P
Jo
89
90
re-
65
thickness of 0.3 mm, writing speed of 100 mm/s, line distance of 0.1 mm and laser power of 1.9 W are used. The scanning direction of the laser is alternated 90◦ in every subsequent layer to obtain the same
3
Journal Pre-proof
Hd
of
R
pro
H
D
Figure 2: Start of the compression test: the contact area Figure 3: Illustration of a sample with height H, diameter is small. D, roughness R and a damaged zone Hd .
91
mechanical properties in the two principal directions. The resulting printed cylinders have an average height
92
and standard deviation of H = 4.32 ± 0.07 mm and diameter D = 3.70 ± 0.05 mm. After sintering, the samples are baked in an oven at 180 ◦ C for 10 minutes.
94
2.2. Mechanical testing
95
Compression tests are performed with a rheometer (Discovery HR-3, TA instruments), having a flat bottom
96
plate and flat top plate (Petisco-Ferrero et al., 2019). The axial displacement is controlled, while measuring
97
the axial force, to obtain the force-displacement curves of the compressed samples. The printed samples are
98
again conditioned at 50% RH and a temperature of 23 ◦ C for at least 16 hours prior to testing. The samples
99
are not perfectly flat, see Figure 2, which results in imperfect contact between the non-flat sample and the
100
flat compression plate. The average height of the surface irregularities R, illustrated in Figure 3 and referred
101
to as roughness, is 0.31 mm and has a standard deviation of 0.09 mm. As a consequence, the contact area
102
varies during the initial part of the test.
urn al P
re-
93
103
To translate the force to a true stress and the displacement to true strain, the actual deformed cross-
104
sectional area is required because of the large applied strain. In this research however, the engineering
105
stress is used because it is not possible to correctly quantify the dynamic contact area between sample
106
and compression plate. Moreover, the samples are highly compressible and there is no visible deformation
107
of the cross-sectional area, suggesting that the engineering stress is close to the true stress. The initial
108
cross-sectional area A0 is therefore taken to calculate the engineering stress σeng , according to
σeng = F/A0 .
The true strain ε is calculated from the displacement u and sample height H (see e.g. Fenner, 1999) using
Jo
109
(1)
ε = ln(u/H + 1).
4
(2)
Journal Pre-proof
An exponentially decreasing velocity profile is prescribed such that the applied true strain rate ε˙ is constant
111
during the test. In addition, cyclic loading tests are performed. The sample is loaded and subsequently
112
unloaded with a constant velocity of 20 µm/s, increasing the maximum displacement for each loading-
113
unloading cycle, see Figure 4. In between cycles, the sample is not loaded during an interval that spans a
114
duration of 10 times the loading time, allowing the sample to relax.
pro
of
110
0.9 0.8 0.7 0.6 0.5 0.4 0.3
re-
0.2 0.1 0
500
1000
1500
Figure 4: Applied displacement in cyclic loading.
Finally, to measure linear viscoelastic effects, stress relaxation experiments are performed. A strain
116
between 1% and 2% is applied and kept constant, while measuring the stress as a function of time. The
117
resulting viscoelastic response is described by a series of stresses σi and relaxation times ti , according to
urn al P
115
σeng (t) =
M X i=1
σi · exp(−t/ti ),
(3)
in which M is the number of used modes.
119
3. Constitutive model
120
3.1. Material characterization
121
The engineering stress as a function of compressive true strain of six compression tests is shown in Figure
122
5, for a strain rate of approximately 10−3 s−1 . The effect of the non-flat surface is clearly visible in the
123
first part of the curves, where the sample surface is flattened and the contact area is increasing. This causes
124
a gradual start-up region, which is sample geometry dependent. The different measurements are therefore
125
126
127
Jo
118
shifted based on their maximum stiffness using the correction procedure from Pelletier et al. (2007). At the maximum stress, crack formation is observed and from this point onwards the structure of the sample is disintegrating. In the stress-strain response, this is manifested as a decrease of the stress as a function of
5
Journal Pre-proof
128
strain, also known as strain softening. From Figure 5, it is observed that the stiffness and fracture stress are quite reproducible. The start-up region due to the non-flatness of the contact surface prevents an accurate direct measurement of the Young’s modulus. The mean tangent stiffness ∂σ/∂ε between a strain of 0.025
131
and 0.030 for the tests with a strain rate of 10−3 s−1 is 50.5 MPa, and the standard deviation is 5.3 MPa.
132
The differences in sample roughness are expected to be the main cause for the variations in stiffness. For
133
the fracture stress, a mean value of 2.8 MPa and a standard deviation of 0.4 MPa are obtained.
pro
of
129
130
4
3
1
0 -0.1
-0.05
re-
2
0
0.05
0.1
0.15
0.2
urn al P
Figure 5: Compression testing of cylindrical samples: each curve represents a different sample and is horizontally shifted based on maximum stiffness.
134
The compression tests are repeated at a true strain rate of approximately 10−4 s−1 to investigate the
135
rate-dependence of the fracture stress. In Figure 5, four measurements with this smaller strain rate are
136
compared to the compression tests with a strain rate of 10−3 s−1 . No clear indication that the fracture stress
137
may depend on the applied deformation rate is observed. In the compression tests, the tangent stiffness is increasing as a function of the applied strain. Image
139
analysis of the deformation, taken during the tests, indicates that this is not only due to flattening of the
140
rough surface. The increase in stiffness originates from the porosity of the samples resulting from the laser
141
sintering process. With increasing applied deformation, the pores are collapsing, the material becomes denser
142
and thus the stiffness increases. This is confirmed in the cyclic loading tests. The result of such a test is
143
presented in Figure 6a, in which the loading curves are shown. In between cycles, the sample is not loaded
144
during an interval that spans a duration of 10 times the loading time, allowing the sample to relax. The
145
sample does not recover its original height indicating the presence of plastic deformation. The tangent
146
stiffness is determined for each cycle by the tangents shown in Figure 6a. The measurement is performed
147
three times resulting in the stiffness as a function of the applied strain, depicted in Figure 6b. The stiffness
148
149
150
Jo
138
is clearly increasing until the fracture stress is reached. After this, both stiffness and maximum stress of the damaged sample are decreasing. Finally, stress relaxation experiments are performed on three samples. Figure 7 shows the decay in stress 6
3
140
2.5
120
of
Journal Pre-proof
100
2
80 1.5 60 1
pro
40
0.5
20
0
0
0
0.05
0.1
0.15
0.2
0
(a)
0.05
0.1
0.15
0.2
(b)
re-
Figure 6: (a) Loading curves of a single sample (black) revealing tangent lines (blue) in cyclic loading and (b) the apparent stiffness for three different samples, determined from the tangents.
151
during such a measurement, indicating a viscoelastic material response. The response is described using four
152
modes (M = 4) in Equation 3 and this prediction is included in Figure 7. 0.3
urn al P
0.25 0.2
0.15
0.1
0.05
0 100
101
102
103
104
Figure 7: Stress relaxation measurement and model prediction using Equation 3.
To describe the material behavior, a constitutive model is needed that relates stress and strain as observed
154
in the experiments. The constitutive model is formulated in 3D in a tensorial form. The formulation allows
155
for high nonlinearity and large deformations. The model is represented by a single Maxwell element including
156
a damage model to take the brittle failure into account. The mechanical analogue of the Maxwell element
157
consists of a spring with shear modulus Gi and a dashpot with shear viscosity ηi = Gi ti as illustrated in
158
Figure 8.
Jo
153
7
Journal Pre-proof
Gi
of
ηi
Figure 8: A Maxwell element represented by a spring with shear modulus Gi and dashpot with viscosity ηi .
3.2. Tensor notations
160
In the derivation of the model, some notations are used which are briefly described here. A second-order
161
tensor A in the three-dimensional space is written as:
pro
159
A = Aij ~ei~ej ,
162
(4)
with Cartesian vector basis {~e1 , ~e2 , ~e3 } and summation over repeated indices. The index notation of this
tensor is Aij and the transpose of a tensor A, denoted as AT , equals Aji . The dot product or inner product
164
of two tensors A and B results in a second-order tensor C which is here defined as:
re-
163
C = A · B, with Cik = Aij Bjk ,
166
167
and the double dot product in a scalar c:
urn al P
165
(5)
c = A : B = Aij Bji .
(6)
tr(A) = Aii
(7)
1 Ad = A − tr(A)I, 3
(8)
The trace of a tensor A is
and the deviatoric part is obtained as
168
with the unit tensor I = ~ei~ei . The notation A˙ represents a differentiation with respect to time.
169
3.3. Kinematics
170
~ 0 ~u)T transforms the initial configuration to a deformed configThe deformation gradient tensor F = I + (∇
172
173
~ 0 denotes the gradient operator with respect to the initial configuration and ~u the displacement uration. ∇
Jo
171
vector. The starting point of the model is the multiplicative decomposition of the deformation gradient tensor in an elastic deformation Fe and a plastic deformation Fp (Lee, 1969), given by
8
Journal Pre-proof
Cc
Fp
C^
pro
Ci
of
F
Fe
Figure 9: Multiplicative decomposition of the deformation gradient tensor.
(9)
re-
F = Fe · Fp .
174
This is illustrated in Figure 9. The intermediate configuration Cˆ is obtained by deforming the initial
175
configuration Ci with Fp . The current configuration Cc is reached by applying Fe to this intermediate
176
configuration.
The velocity gradient tensor L is obtained as
urn al P
177
L = F˙ · F −1 = F˙e · Fe−1 + Fe · F˙p · Fp−1 · Fe−1 .
178
The plastic part of the velocity gradient tensor is defined in the intermediate configuration as ˆ p = F˙p · F −1 , L p
179
(10)
(11)
ˆ p in this intermediate configuration is and the plastic rate of deformation tensor D ˆ p = 1 (L ˆp + L ˆ T ). D p 2
(12)
180
Assuming that plastic deformation is isochoric, the volume change ratio J is calculated from the determinant
181
of Fe :
J = det(F ) = det(Fe ).
183
The assumption is made that the plastic part of the deformation does not have a rigid body rotation, i.e. Rp = I, with Rp the plastic rotation tensor. The plastic deformation gradient tensor is therefore given by
Jo
182
(13)
Fp = Rp · Up = Up , 9
(14)
Journal Pre-proof
in which Up is the plastic right stretch tensor.
185
3.4. Stress calculation
186
The Cauchy stress σ is the stress in the current configuration Cc . The total stress is split into two parts, a
187
hydrostatic part σ h and a deviatoric part σ d , according to
pro
of
184
σ = σh + σd .
(15)
188
The deviatoric part is obtained using a hyperelastic Ogden model (Ogden, 1972, 1984). This model is chosen
189
for its ability to model strongly nonlinear elastic behavior, in this case an increase in stiffness in compression,
190
see Figure 6b. The deviatoric stress is calculated from the isochoric elastic left stretch tensor V˜e : N 1X gi (V˜eai )d , J i=1
re-
σd =
(16)
191
in which gi and ai are the shear modulus and the exponent of Ogden term i. The total shear modulus G is
192
calculated from these parameters as
N X gi ai
urn al P
G=
i=1
2
.
(17)
193
The pair of parameters (gi , ai ) may also have negative values to influence the mechanical response in com-
194
pression, while positive values are determining the response in tension. By having multiple terms, it is
195
possible to address both tension and compression with a single set of model parameters. The left elastic
196
stretch tensor Ve is determined from Fe ,
Ve =
197
q
Fe · FeT ,
(18)
and the isochoric elastic left stretch tensor to the power ai using a
i V˜eai = J − 3 Veai .
(19)
The increase in shear modulus, due to the decreasing porosity during compression, is accompanied with
199
an increase in bulk modulus. The factor ai , which is numerically taking care for this stiffness increase, is
200
therefore included in the hydrostatic stress (Nedjar, 2016). The hydrostatic part of the stress is given by
Jo
198
σh =
N X Kgi i=1
2G
(J
2ai 3
10
−1
−J
−ai 3
−1
)I,
(20)
Journal Pre-proof
with K the bulk modulus.
202
3.5. Plasticity and damage
203
ˆ p is calculated from the stress-dependent viscosity function η(¯ The plastic rate of deformation tensor D τ)
204
and the deviatoric stress using ˆp = D
pro
of
201
1 ˆd S . 2η(¯ τ)
(21)
205
The deviatoric stress in the intermediate configuration is expressed in the elastic second Piola-Kirchhoff-like
206
stress tensor Sˆd which is calculated from the Cauchy stress:
207
The viscosity is calculated as
re-
Sˆd = (JFe−1 · σ · Fe−T )d .
τ¯ . γ˙ p
η(¯ τ) =
(22)
(23)
The plastic strain rate γ˙ p , parameterized by the reference shear rate γ˙ 0 , the power m and the characteristic
209
stress τc , is given by
urn al P
208
γ˙ p = γ˙ 0
210
τ¯ τc
m
,
(24)
1 ˆd ˆd S :S . 2
(25)
and the equivalent shear stress is calculated as
τ¯ =
r
211
The softening part of the stress-strain curve is described by a damage law. The evolution of the damage
212
parameter ξ, which is initially zero, is calculated from the plastic strain rate by ξ˙ =
ξ 1− dγ˙ p , ξ∞
(26)
in which d and ξ∞ are parameters describing the rate at which the damage develops and the limit value for the damage, respectively. The damage parameter ξ affects the current shear moduli gi and the characteristic
215
stress τc , such that
216
Jo
213
214
τc = τ0 (1 − ξ) and gi = g0,i (1 − ξ)
in which τ0 and g0,i are the initial, undamaged variables. 11
(27)
Journal Pre-proof
3.6. Time integration
218
The constitutive equations are integrated from time tn to tn+1 using an implicit integration scheme (Amiri-
of
217
Rad et al., 2019). The applied deformation is given by F and the elastic part of this deformation gradient
220
tensor is used to calculate the total stress. One determines the amount of plasticity, starting from the plastic
221
right Cauchy-Green deformation tensor Cp = FpT · Fp , and then looks to the remaining part of F , which
222
gives Fe . To this purpose, a Newton-Raphson scheme is used, in which first an initial guess for Cp,n+1 is
223
taken. From Cp,n+1 the plastic right stretch tensor Up,n+1 is determined:
Up,n+1 =
pro
219
p Cp,n+1 .
(28)
Using the polar decomposition of the deformation gradient, Equation 14, the elastic deformation gradient
225
tensor is recovered as
re-
224
−1 Fe,n+1 = Fn+1 · Up,n+1 .
(29)
d ˆ p,n+1 is found from the flow rule. The time The stress σn+1 is obtained with this Fe,n+1 and subsequently D
227
derivative of Cp,n+1 can now be calculated using
urn al P
226
T T C˙ p,n+1 = F˙p,n+1 · Fp,n+1 + Fp,n+1 · F˙p,n+1
=
228
T 2Up,n+1
(30)
ˆ p,n+1 · Up,n+1 . ·D
Using C˙ p,n+1 , the tensor Cp,n+1 is recomputed using the backward Euler method: Cp,n+1 = Cp,n + ∆tC˙ p,n+1 ,
(31)
in which ∆t is the time step size. The new value of Cp,n+1 enters the Newton-Raphson scheme again,
230
until convergence in the backward Euler scheme is obtained. The Frobenius norm of the residue (difference
231
between initial guess and recomputed tensor) is used as convergence criteria, repeating the time integration
232
until a norm smaller than 10−10 is obtained. The damage parameter ξ is updated at the end of each increment
233
by the explicit time integration of Equation 26.
234
3.7. Implementation
235
236
Jo
229
The constitutive model is implemented in the finite element package MSC.Marc using the user subroutine HYPELA2.
12
Journal Pre-proof
4. Results
238
4.1. Finite element model
239
The model performance is first assessed through the homogeneous deformation of a cube with sides of
240
1 mm. The element size is varied from a single element up to cubic elements with sides of 62.5 µm to
241
investigate mesh-sensitivity of the model. 3D linear hexahedral elements are used. The boundary conditions
242
are schematically presented in Figure 10a. Uniaxial compression simulations are performed by prescribing
243
the displacement at the top surface and fixing the bottom in z-direction. Moreover, the nodes on the left
244
face are fixed in x-direction and in the front face in y-direction, corresponding to symmetry conditions. An
245
initial time step size of 1.0 second is taken, which turned out to be sufficiently small (i.e. the final solution
246
is not dependent on it anymore), and is decreased during softening as a result of the high nonlinearity of
247
Equation 24.
re-
pro
of
237
248
Figure 2 reveals that the sample is not flat resulting in a small tangent stiffness at the start of the
249
compression test. To identify the effect of the geometry, resulting in inhomogeneous deformation, on the effective macroscopic response, finite element simulations are performed with a simplified geometry that represents the actual samples, depicted in Figure 10b. The average height H=4.32 mm, diameter D=3.70
252
mm and roughness R=0.31 mm, illustrated in Figure 3, of the samples are taken to model the geometry
253
with 3040 3D linear hexahedral elements. The bottom nodes are fixed in z-direction and the rotation of
254
the geometry is constrained. The compression is applied in the z-direction by contact with a plate with
255
a prescribed velocity such that the true strain rate is constant. The initial time step size is equal to 0.2
256
seconds. The effect of contact friction was tested with a Coulomb friction model, which turned out to be
257
small compared to the sample roughness and is therefore further neglected.
urn al P
250
251
y
z
z
x
(a)
y
x
(b)
Jo
Figure 10: Finite element model: (a) boundary conditions in homogeneous deformation and (b) representative geometry. Red arrows indicate fixed boundary condition and black arrows represent the applied compression.
13
Journal Pre-proof
4.2. Parameter identification
259
The model parameters and their identified values are listed in Table 1. From the experiments, a stiffness of
260
40 MPa is estimated at the moment the top surface was flattened. This value is assumed to be equal to the
261
Young’s modulus E. The flattening makes it hard to determine a Poisson’s ratio ν, which is estimated as
262
0.15. From these values, the shear modulus and bulk modulus are calculated using
G=
pro
of
258
E E , and K = , respectively. 2(1 + ν) 3(1 − 2ν)
(32)
The hyperelastic part of the model consists of two Ogden terms each with a shear modulus g0,i and
264
exponent ai . One term is dominant in tension (positive values) and the other is dominant in compression
265
(negative values). Because of the brittleness of the samples, it is not possible to perform tensile tests on this
266
particular material. A standard value corresponding to a1 = 2 is therefore adopted, which makes the term
267
equal to a Neo-Hookean model and g0,1 = 1, which has an arbitrary small value for the product a1 g0,1 as a
268
result.
re-
263
Table 1: Model parameters.
K [MPa] 19.0
g0,1 [MPa] 1.0
urn al P
G [MPa] 17.4 γ˙ 0 [s−1 ] 1.0·10−3
τ0 [MPa] 2.6
m [-] 15
5
a1 [-] 2.0
g0,2 [MPa] -0.73
d [-] 5.0
ξ∞ [-] 0.85
a2 [-] -45
ξ 0.85
4 3 2
0
0
0.05
0.1
0.15
y
m = 15
1
z
0.2
m = 50 x
0.0
Figure 11: Regularization of damage localization by decreasing m in a uniaxial compression test on a cube, showing (a) the macroscopic response and (b) the local damage after applying a strain of 0.13.
270
271
272
To uniquely determine the three parameters that describe the plastic strain rate as given in the specific
Jo
269
form in Equation 24, the value of one parameter can be chosen arbitrary. The reference shear rate γ˙ 0 is therefore taken equal to the experimentally applied strain rate. An effect of the damage is that it may localize, which can have mesh-dependent results as a consequence. This is regularized by introducing a 14
Journal Pre-proof
273
rate-dependent material behavior (Needleman, 1988). The rate-dependence of the material is captured by the model parameter m. Figure 11a shows the macroscopic response of a uniaxial compression test on the cube with sides of 62.5 µm and Figure 11b depicts the damage parameter after applying a compressive strain
276
of 0.13, for different values of m. For this simple test case, choosing m = 30 is sufficient to regularize the
of
274
275
damage. However, it is found that for a more complex geometry, a value of m = 15 is necessary to provide
278
a unique solution which is independent of the mesh size and time step size. The introduced rate-dependent
pro
277
279
behavior may be in contradiction with the rate-independent behavior observed in the experiments, as shown
280
in Figure 5. The chosen value m = 15 is sufficiently small to regularize the damage without introducing a
281
significant effect of the deformation rate on the fracture stress.
The five remaining parameters are g0,2 , a2 , τ0 , d and ξ∞ . Considering Equation 17, only a2 has to
283
be determined to also recover g0,2 , knowing the parameters of the first Ogden (i=1) term and the total
284
shear modulus. The number of parameters to determine is therefore reduced to four. The effect of these
285
parameters on the intrinsic material response, in terms of true stress (deformed area), is presented in Figure
286
12. As a reference, the parameter set given in Table 1 is used. Parameter a2 affects the stiffness increase as
287
a function of the applied compressive strain, and τ0 defines the fracture stress. The parameter d is identified
288
from the rate of damage, i.e. the slope of the stress-strain curve in the damaged region, and ξ∞ determines
289
the fully damaged state and corresponds to the ratio between the maximum stress and the plateau level at
290
large strains.
urn al P
re-
282
291
The main issue in identifying the model parameters is the local damage that is resulting from the non-
292
flat top surface. Figure 13 shows the damage parameter in a cross section of the sample after applying a
293
strain equal to 0.05. Locally, the material is already severely damaged, while a flat sample would be still
294
undamaged. The consequences are shown in Figure 14 for different surface topologies. The roughness R
295
is varied and different surface profiles are considered. In comparison with a flat sample, the stiffness and
296
fracture stress decrease and softening becomes more severe. Moreover, it is concluded that the influence of
297
R is more important than the surface topology. This confirms the choice for the definition of R and the
298
representative geometry.
To complete the identification, the influence of the different parameters on the stress-strain response of the
300
representative geometry is shown in Figure 15 and assessed through their influence on the intrinsic response,
301
depicted in Figure 12. The main difference is that the intrinsic response reveals a clear distinction between
302
elastic and plastic deformation. The parameters that correspond to plasticity do not affect the elastic part.
303
Adding roughness has the result that local plasticity is influencing the macroscopic stress-strain response,
304
including the deformation before the fracture stress is reached. Compared to Figure 12, the effect of a2 is
305
306
Jo
299
similar, but τ0 and d are clearly affecting the stiffness with increasing strain. Moreover, the fracture stress decreases with increasing d and the effect of d on the softening slope is fading. The softening is dominated
15
5
4
4
3
3
2
2
pro
5
of
Journal Pre-proof
1
1
0
0
0
0.05
0.1
0.15
0.2
0
5
0.05
0.1
0.15
0.2
0.1
0.15
0.2
5
4
4 3
2
re-
3
2
1
1
0
0
0
0.05
0.1
0.15
0.2
0
0.05
urn al P
Figure 12: Main characteristics of the model, showing the influence of a2 , τ0 , d and ξ∞ on the intrinsic mechanical behavior for a constant strain rate ε˙ = 10−3 s−1 . Note that for different values of d all curves will eventually reach the same plateau stress.
307
by geometrical effects and the material softening becomes less prominent. A new effect emerges, which are
308
the oscillations of the curve for larger values of d, induced by the many geometrical imperfections. The finally identified parameters are given in Table 1 and determined using the following characterization procedure. First, ξ∞ is determined because it is the only parameter that determines the ratio between the
311
fracture stress and the plateau stress after fracture. Moreover, ξ∞ is determining the depth of damage, Hd in
312
Figure 3. Next, d is chosen to approximate the jaggedness of the experimental results. Finally, τ0 is defined
313
to match the fracture stress, whereas a2 covers the necessary stiffness increase. Given the experimental
314
scatter, the obtained set of parameters adequately describes the overall experimentally obtained mechanical
315
behavior, as seen in Figure 16. Comparing the damaged zones in the experiments and simulation, Figure 17,
316
it is obvious that the simulation is also predicting this zone well. In comparison with the experiment, the
317
predicted damaged is more smooth. This difference is due to the fact that local heterogeneities and defects,
318
which may cause localization of damage, are not taken into account.
Jo
309
310
16
Journal Pre-proof
ξ
4.5
0.85
of
4 3.5 3 2.5 2
pro
1.5 1
0.5
0
0
0.0
0.05
0.1
0.15
0.2
re-
Figure 14: Comparison of the response (ε˙ = 10−3 s−1 ) Figure 13: Simulated damage in a cross section of the of samples with different roughness values. For each R sample after applying a strain of 0.05. value, three different surface profiles are simulated.
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
urn al P
0.5
0
0
0.05
0.1
0.15
0.5
0
0.2
3.5
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
0.05
0.1
0.15
0.2
0.05
0.1
0.15
0.2
0
0.05
0.1
0.15
0.2
3.5
3
0
0
0
Figure 15: Influence of model parameters on the macroscopic stress-strain response (ε˙ = 10−3 s−1 ) of the representative geometry.
5. Conclusion
320
An experimental procedure, based on compression testing, is presented and provides a solid basis for char-
321
322
Jo
319
acterizing the intrinsic mechanical properties and failure of 3D printed food. For the samples that are used in this research, the following observations are incorporated in a 3D constitutive model that allows for high
17
Journal Pre-proof
of
4
3
1
0 0
0.05
0.1
pro
2
0.15
0.2
Figure 16: Comparison between experiments and model prediction for ε˙ = 10−3 s−1 . ξ
re-
0.85
damaged
0.0
urn al P
Figure 17: Damaged zone at the surface in simulation and experiment after applying a strain of 0.2.
323
nonlinearity, inelastic deformation and damaging behavior. First, the stiffness increases as a function of the
324
applied strain caused by densification due to the porous microstructure. This is followed by the initiation of
325
brittle failure at the macroscopic fracture stress. During failure, softening is observed in which both stiffness
326
and stress decrease. After characterising some of the model parameters, a representative sample geometry,
327
taking the roughness of the samples into account, is modeled, allowing to identify the complete set of material
328
parameters. The presented model and characterization procedure are not only useful for the application of
329
3D food printing, but can be considered more generic and applicable to other food. A method to process or post-process flat samples without introducing any damage before testing would
331
enable a more accurate determination of the parameters. Characterization of the volumetric compressibility
332
of the material was not possible with the available experimental data and requires future investigation.
333
The value for the Poisson’s ratio is indicative only and the assumption that plasticity is volumetrically
334
incompressible may not be correct for this material. Although the constitutive model is able to distinct non-
335
linear behavior in both the compressive and tensile regime, the current characterization procedure focuses
336
on the more important compressive response only. An important future extension of the current model
337
338
Jo
330
could be the addition of anisotropy which results from the layer-by-layer manufacturing process. Despite these limitations, the model adequately predicts the mechanical response and damaged zone of the printed
18
Journal Pre-proof
samples. A powerful approach for predicting properties of specifically designed food products is therefore
340
obtained.
of
339
The mentioned aspects are important at the current length scale of modeling. This length scale only
342
considers effective properties originating from the microstructure. However, the microstructure of the 3D
343
printed food is highly heterogeneous. For future research, it would therefore be interesting to relate the
344
mechanical properties to the microstructure. The spatial distribution of material parameters can then be
345
incorporated by assigning different parameter values at different locations. Moreover, process parameters
346
(e.g. laser power) can be locally varied within a printed product. This enables tailoring the microstructure
347
(e.g. porosity, dimensions of voids and structural components) and, as a result, the mechanical properties
348
(e.g. Young’s modulus, fracture stress). Including process-microstructure-property relations should enrich
349
the current framework to become an invaluable tool in designing customized food.
350
Acknowledgements
351
The authors wish to thank TNO (Netherlands Organisation for Applied Scientific Research) for the financial
352
support, for the instructions on sample processing and for the access to the SLS machine. We also thank
353
Martijn Noort (Wageningen University and Research) for the preparation of the powder mixture, Ruth
354
Cardinaels (Eindhoven University of Technology) for the use of the rheometer and Susana Petisco Ferrero
355
(Eindhoven University of Technology) for the training on compression testing with the rheometer.
356
References
357
J. Sun, Z. Peng, W. Zhou, J. Y. H. Fuh, G. S. Hong, A. Chiu, A review on 3D printing for customized food
360
361
362
363
364
365
366
367
re-
urn al P
359
fabrication, Procedia Manufacturing 1 (2015) 308–319. F. C. Godoi, S. Prakash, B. R. Bhandari, 3d printing technologies applied for food design: Status and prospects, Journal of Food Engineering 179 (2016) 44–54. J. I. Lipton, Printable food: the technology and its application in human health, Current Opinion in Biotechnology 44 (2017) 198–201.
J. M. Aguilera, D. J. Park, Texture-modified foods for the elderly: Status, technology and opportunities, Trends in Food Science and Technology 57 (2016) 156–164.
Jo
358
pro
341
S. Mantihal, S. Prakash, F. C. Godoi, B. Bhandari, Optimization of chocolate 3D printing by correlating thermal and flow properties with 3D structure modeling, Innovative Food Science and Emerging Technologies 44 (2017) 21–29.
19
Journal Pre-proof
M. Lanaro, D. P. Forrestal, S. Scheurer, S. Liao, D. J. Slinger, S. K. Powell, M. A. Woodruff, 3D printing
369
complex chocolate objects: Platform design, optimization and evaluation, Journal of Food Engineering
370
215 (2017) 13–22.
372
C. A. Hamilton, G. Alici, M. in het Panhuis, 3D printing vegemite and marmite: Redefining breadboards, Journal of Food Engineering 220 (2018) 83–88.
pro
371
of
368
373
C. Le Tohic, J. J. O’Sullivan, K. P. Drapala, V. Chartrin, T. Chan, A. P. Morrison, J. P. Kerry, A. L.
374
Kelly, Effect of 3D printing on the structure and textural properties of processed cheese, Journal of Food
375
Engineering 220 (2018) 56–64.
378
379
380
381
382
383
384
385
M. W. Noort, J. V. Diaz, K. J. C. van Bommel, S. Renzetti, J. Henket, M. B. Hoppenbrouwers, Method for the production of an edible object using SLS. US2017/0266881 A1, 2017.
re-
377
F. C. Godoi, B. R. Bhandari, S. Prakash, M. Zhang, Fundamentals of 3D food printing and applications, Elsevier Ltd., London, 2019. doi:https://doi.org/10.1016/C2017-0-01591-4. Z. Liu, M. G. Scanlon, Predicting mechanical properties of bread crumb, Food and Bioproducts Processing 81 (2003) 224–238.
S. Guessasma, L. Chaunier, G. Della Valle, D. Lourdin, Mechanical modelling of cereal solid foods, Trends
urn al P
376
in Food Science and Technology 22 (2011) 142–153.
G. E. Attenburrow, R. M. Goodband, L. J. Taylor, P. J. Lillford, Structure, mechanics and texture of a food sponge, Journal of Cereal Science 9 (1989) 61–70.
386
E. O. Afoakwa, A. Paterson, M. Fowler, J. Vieira, Microstructure and mechanical properties related to
387
particle size distribution and composition in dark chocolate, International Journal of Food Science and
388
Technology 44 (2009) 111–119.
390
391
392
393
394
395
396
N. Sozer, H. Dogan, J. L. Kokini, Textural properties and their correlation to cell structure in porous food materials, Journal of Agricultural and Food Chemistry 59 (2011) 1498–1507. R. Agbisit, S. Alavi, E. Cheng, T. Herald, A. Trater, Relationships between microstructure and mechanical properties of cellular cornstarch extrudates, Journal of Texture Studies 38 (2007) 199–219. S. C. Warburton, A. M. Donald, A. C. Smith, The deformation of brittle starch foams, Journal of Materials Science 25 (1990) 4001–4007.
Jo
389
F. Robin, C. Dubois, D. Curti, H. P. Schuchmann, S. Palzer, Effect of wheat bran on the mechanical properties of extruded starchy foams, Food Research International 44 (2011) 2880–2888.
20
Journal Pre-proof
400
401
402
403
404
405
406
of
399
bread, Journal of Cereal Science 24 (1996) 15–26.
M. G. Scanlon, M. C. Zghal, Bread properties and crumb structure, Food Research International 34 (2001) 841–864.
M. C. Zghal, M. G. Scanlon, H. D. Sapirstein, Cellular structure of bread crumb and its influence on
pro
398
C. J. A. M. Keetels, K. A. Visser, T. van Vliet, A. Jurgens, P. Walstra, Structure and mechanics of starch
mechanical properties, Journal of Cereal Science 36 (2002) 167–176.
H. H. Friedman, J. E. Whitney, A. S. Szczesniak, The texturometer - a new instrument for objective texture measurement, Journal of Food Science 28 (1963) 390–396.
A. J. Rosenthal, Texture profile analysis - How important are the parameters?, Journal of Texture Studies 41 (2010) 672–684.
re-
397
407
M. Peleg, The instrumental texture profile analysis revisited, Journal of Texture Studies (2019) 1–7.
408
J. G. F. Wismans, J. A. W. van Dommelen, L. E. Govaert, H. E. H. Meijer, B. van Rietbergen, Computed
409
tomography-based modeling of structured polymers, Journal of Cellular Plastics 45 (2009) 157–179. V. Vancauwenberghe, M. A. Delele, J. Vanbiervliet, W. Aregawi, P. Verboven, J. Lammertyn, B. Nicola¨ı,
411
Model-based design and validation of food texture of 3D printed pectin-based food simulants, Journal of
412
Food Engineering 231 (2018) 72–82.
413
414
urn al P
410
M. A. Del Nobile, S. Chillo, A. Mentana, A. Baiano, Use of the generalized maxwell model for describing the stress relaxation behavior of solid-like foods, Journal of Food Engineering 78 (2007) 978–983.
415
M. Shirmohammadi, P. K. D. V. Yarlagadda, Y. T. Gu, A constitutive model for mechanical response
416
characterization of pumpkin peel and flesh tissues under tensile and compressive loadings, Journal of Food
417
Science and Technology 52 (2015) 4874–4884.
419
420
421
422
423
424
425
Z. Liu, M. G. Scanlon, Modelling indentation of bread crumb by finite element analysis, Biosystems Engineering 85 (2003) 477–484.
S. Petisco-Ferrero, R. Cardinaels, L. C. A. van Breemen, Miniaturized characterization of polymers: From synthesis to rheological and mechanical properties in 30 mg, Polymer 185 (2019) 121918. R. T. Fenner, Mechanics of Solids, CRC Press, Boca Raton, 1999.
Jo
418
C. G. N. Pelletier, E. C. A. Dekkers, L. E. Govaert, J. M. J. den Toonder, H. E. H. Meijer, The influence of indenter-surface misalignment on the results of instrumented indentation tests, Polymer Testing 26 (2007) 949–959.
21
Journal Pre-proof
E. H. Lee, Elastic-plastic deformation at finite strains, Journal of Applied Mechanics 36 (1969) 1–6.
427
R. W. Ogden, Large deformation isotropic elasticity: On the correlation of theory and experiment for
428
of
426
compressible rubberlike solids, Proceedings of the royal society A 328 (1972) 567–583.
R. W. Ogden, Non-linear elastic deformations, Chichester: E. Horwood; New York: Halsted Press, 1984.
430
B. Nedjar, On constitutive models of finite elasticity with possible zero apparent Poisson’s ratio, International
434
435
Viscoplastic Model for Short-Fiber Composites, Mechanics of Materials 137 (2019) 103141. A. Needleman, Material rate dependence and mesh sensitivity in localization problems, Computer Methods in Applied Mechanics and Engineering 67 (1988) 69–85.
re-
433
A. Amiri-Rad, L. V. Pastukhov, L. E. Govaert, J. A. W. van Dommelen, An Anisotropic Viscoelastic-
urn al P
432
Journal of Solids and Structures 91 (2016) 72–77.
Jo
431
pro
429
22