CC BY
18
DOI: 10.17516/1997-1397-2021-14-6-712-725 УДК 624
Numerical Simulation of the Fused Deposition Modeling for the Manufacturing of Parts with both High Geometric Fidelity and Mechanical Quality
Holm Altenbach* Gabor Janiga^
Otto-von-Guericke-Universitat Magdeburg Magdeburg, Germany
Mario Beiner* Rene Androsch Paul-Maximilian Runge^
Martin-Luther-Universitat Halle-Wittenberg Halle (Saale), Germany
Received 26.06.2021, received in revised form 21.08.2021, accepted 23.09.2021 Abstract. With increasing usage of additive manufacturing methods for mechanical parts the need for precise and reliable simulations of the manufacturing process increases as well. In this paper various computations suited for simulating the fused deposition modeling process are considered in two dimensions. In fused deposition modeling a molten polymer is laid down on a prescribed path before the cooling of the melt begins. The occuring flows are treated as multiphase flows. To model the deposition of the filament, methods of computational fluid dynamics are used in ANSYS-Fluent, namely the volume of fluid method (VOF). Different numerical experiments are simulated.
Keywords: volume of fluid method (VOF), non-Newtonian fluid, two-phase flow.
Citation: H. Altenbach, G. Janiga, M. Beiner, R. Androsch, P.-M. Runge, Numerical Simulation of the
Fused Deposition Modeling for the Manufacturing of Parts with both High Geometric Fidelity and
Mechanical Quality, J. Sib. Fed. Univ. Math. Phys., 2021, 14(6), 712-725.
DOI: 10.17516/1997-1397-2021-14-6-712-725.
1. Motivation
To model the geometric fidelity of a mechanical part manufactured by fused deposition modeling, the flow of the molten filament when it is extruded from the printing nozzle and laid down on the printing bed has to be considered. In particular, the printed strand of filament forms a free surface in contact with air. The form of that interface can be captured by means of a multiphase model of the problem, with two fluids, namely the molten filament and air. In order to compute the form of the interface between the fluids, in [3] the volume of fluid (VOF) method
§ [email protected] ^[email protected] © Siberian Federal University. All rights reserved
is employed in a three-dimensional isothermal simulation to find optimal process parameters for a sharp corner of ninety degrees in the first layer on the printing bed, as the polymer is modeled as a Newtonian fluid. In [12] a computational fluid dynamics simulation is carried out to simulate a printed part of multiple layers laid down above each-other, including a temperature gradient from the hot polymer to the surrounding air. In that work, the viscosity is considered to be temperature and shear rate dependent. In the continuing paper, [11] solidification and shrinkage of the laid-down material is considered, together with an improved nozzle geometry. In [2] the dependence of the cross-sectional area of an extruded strand in fused deposition modeling is shown to be influenced by temperature-dependent viscosity of the printing material and results are compared to experiments. In this article various numerical calculations are carried out in the computational fluid dynamics program ANSYS-Fluent. In Sect. 4 a simulation is carried out to validate the VOF method in ANSYS-Fluent based on the experiments in [5] as they are found in [6]. The results of the simulation will be compared to the experimental results from literature qualitatively. In Sect. 5, the isothermal flow of a arbitrarily chosen non-Newtonian fluid through a two-dimensional nozzle will be considered. In all simulations the VOF is used in ANSYS-Fluent. One critical part of that simulation technique is the transport of the interface between the two considered liquids. To advect that interface correctly in the finite volume method, the compressive interface capturing scheme for arbitrary meshes (CICSAM) firstly published in [8] is used in ANSYS-Fluent for all calculations. A short overview of this technique will be given according to [9] since it is suited for multiphase simulations with a high viscosity ratio [10]. In Sect. 6 a simplified two-dimensional model of the FDM-process without the need of a moving mesh is considered.
2. Volume of fluid method
The VOF is an interface-capturing method which does not define the interface between two immiscible fluids as a boundary but as a value of a boundary fraction a of cells of a fixed grid that are partially filled [4]. The computation is performed on a fixed grid, which covers the two fluids on both sides of the free surface. The shape of the free surface is determined by computing the fraction of each near-interface cell that is partially filled [4]. The computation in this paper is performed on such a fixed grid of computational cells, which cover both of the two fluids on each side of the fluid-fluid interface. The shape of that interface is determined by computing the volume fraction a of the traced fluid in each computational cell and by advercting that volume fraction by a transport equation, which is introduced in the VOF model [6]. The value of the volume fraction a is based on an indicator function 4 which can only take the values 0 and 1 according to the following formula [2]
I 0 for x e Dgas at time t, 4 (x,t) = { gas
[1 for x e Diiquid at time t,
where Dgas and Di;qu;d are domains of gas and liquid, respectively. The indicator function accounts for the instantly changing material properties in the domain at the liquid-gas interface. Averaging the indicator function given for a computational cell Q leads to the volume fraction a of the liquid in that cell [2]
1
a (x,t) = -1 i 4 (x,t) dV. vq Jq
The volume fraction a of the liquid phase is computed for each computational cell in the domain (Fig. 1). The value of the volume fraction of the liquid a varies between 0 and 1, 0 < a < 1 (cf. [9]). Does a computational cell Q haves a value of a = 1, it is completely filled by the liquid.
(a)
(b)
Fig. 1. Scheme for volume fraction [7]: (a) true interface, (b) volume fractions
The value of a = 0 indicates that the liquid is absent and the cell is filled with the second, the gaseous phase. However, if its value varies between zero and one 0 < a < 1 for a given cell Qj, that cell encloses the interface of the multiphase (cf. [9]). The VOF treats the multiphase system of liquid and gas as one single fluid with a density p and a dynamic viscosity With the volume fraction of the liquid phase, those parameters are varied in each location x and for all times t in the computational domain [9]
p (X,t) = piiquida (x,t)+[1 - a (X,t)] pgas, n (x,t) = ^liquida (x,t) + [1 - a (x,t)] figas.
The location of the interface can change in time and thereby the volume fraction a varies as well. The transport of the volume fraction is derived by the mass-transport equation for the averaged density p and holds in the following form [2]
da „_,
— + u -Va = 0. dt
Now the movement of the fluid with given averaged density p and averaged dynamic viscosity H can be given by one single Navier-Stokes-equation and one single continuity equation for an averaged fluid [6]
Vu
0,
dpu ^ _2 „_,
~J—+ pu -Vu + HV2U = -Vp + pg + fa, dt
where u is the velocity vector of the averaged fluid and p is the pressure field. Both fluids always travel with the same velocity u in the VOF model [1]. The vector g represents the gravitational acceleration. The last vector on the right-hand side of the momentum-transport equation accounts for surface tension effects on the interface [6]. In this work however, surface-tension-effects are neglected. The conservation transport equations of mass, momentum, and volume fraction hold only for an isothermal mixture. In case of a temperature-difference between the fluids like in the extrusion of a molten polymer in additive manufacturing, an additional
equation has to be considered to account for heat-transfer. In ANSYS-Fluent this equation is the scalar energy equation [1]
d
- (pE) + [u (pE + p)] = (keSVT) + Sh,
where the energy E is treated as a mass-averaged value [1]
apiEi + (1 - a) p2E2
E =
api + (1 - a) p2
In this equation the energies of each phase Ei and E2 are based on the specific heat capacity cp of that phase and the between the phases shared temperature T [1]. The density p is shared between the phases according to the equation above as well as the effective thermal conductivity keff. The source term, Sh contains contributions from radiation, as well as any other volumetric heat sources [1].
3. Scheme of numerical implementation
To solve the given system of partial differential equations, the software suite ANSYS-Fluent is used. The crucial issue for modelling of the multiphase flow is a proper solution of the transport-equation of the volume fraction a, that means the temporal discretization of the time derivative and convective term [9]. The chosen discretization has to fulfill two properties, that is avoiding a smearing of the step interface profile across the domain and assuring that the boundedness criterion is satisfied [9]. This is shown here for the Crank-Nicolson method for integrating in time and the finite volume method for spatial discretization as it is in [9]. Discretizing the advection equation for the volume fraction a
da „
— + u•Va = 0
dt
gives the algebraic equation [9]
VQ t+At , 1 'n t+At o t VQ t 1 to t
AtaQ +2 af Sfui,f = AtaQ - afSfUif,
f=i f=i
where VQ is the volume of one computational cell, Sf is the area of one face of this cell and At denotes the time step size. The convective flux af (value of the volume fraction at the face of the control volume) at the cell faces has to be interpolated [9]. In this work the compressive interface capturing scheme for arbitrary meshes (CICSAM), [8] is employed as an option in ANSYS fluent [1] to obtain the value of the volume fraction at the face of the control volume af, which represents the convective flux.
3.1. Normalized variable diagram
The normalized variable diagram (NVD) provides the foundation of the CICSAM scheme [9]. It is based on the convective boundedness criterion (CBC) that states that the variable distribution between the centers of the neighbouring computational volumes, e.g. Donor (D) and acceptor (A) should remain smooth aD < af < aA (Fig. 2). Using this constraint and information about the value of the volume fraction a.u in the upwind control volume a.u, normalized variables are introduced [9]:
af
af — au a a — au' ad — au a a — au1
which allow the CBC to be rewritten in the form aD ^ af ^ 1, which is represented in the right side of Fig. 2.
Fig. 2. Boundedness criterion with normalized variables [9]
The CBC is only satisfied in the grey area of the right-hand side diagram. Now it is possible to find the value of the volume fraction at the face of the control volume af [9]
af = (1 - /3f) a.D + f aA, af — aD
ßf
1 — aD
This expression for the interpolation of af is used in the algebraic form of the transport-equation for the volume fraction a.
3.2. The compressive interface capturing scheme for arbitrary meshes
In the case of the CICSAM scheme, additional assumption about the dependence where the CBC is satisfied on the CFL condition is used [9]. The local value of the Courant number, defined at the face of the control volume Cf = uf Sf At/Vn together with the CBC changes the following formula aD < a f < 1 to [9]
aD ^ af ^ min (1, aD/Cf),
which can be plotted in the following NVD (Fig. 5). In this figure, the grey area is again the region where the face-fluxes are bounded. It is visible, that the stable region becomes larger with a smaller Courant number. The CICSAM scheme is built of two components, firstly the HYPER-C scheme and secondly hte ULTIMATE-QUICKEST scheme [9]. The HYPER-C scheme
Fig. 3. CICSAM: dependence of the CBC on the local CFL condition [9]
is represented as [9]
af = aD, 0 <aD,a.D > 1,
afCBC
I, 0 < aD < 1.
The HYPER-C scheme satisfies the CBC criterion and is shown to be compressive, which means that it changes any gradient to a step profile due to the downwind differencing scheme employed (DDS) [9]. The compressive character of that scheme has a drawback in certain flow configurations, that is, if the interface is tangential to the flow direction, it tends to artificially deform its shape [9]. So there has to be employed a less compressive formulation. In the case of CICSAM, this second component is the ULTIMATE-QUICKEST scheme [9]
af = ad, 0 < aD,aD > 1,
aUQ = \ . ( 8Cf aD + (1 ~ Cf )(65p + 3) _ ) _
1 mm i —---^--, afCBC j , 0 < aD < 1.
Finally, to switch smoothely between both schemes, linear blending is used, with a blending factor of Yf: 0 < Yf ^ 1 [9]. The value of that blending factor depends on the angle Of between the unit vector normal to the interface n = VaD/\VaD | and the unit vector parallel to the line between centers of the donor cell D and the acceptor cell A d = DA/\DA\. When the interface position is normal to the direction of the flow, the blending factor yields Yf = 1 and the HYPER-C scheme is used. In the case of tangential orientation of the interface, the blending factor takes Yf =0 and the second scheme, ULTIMATE-QUICKEST is employed. So we get for the interpolated face-flux [9]:
af = Yf afCBC + (1 - Yf) afUQ, 0 < Yf ^ 1,
Of = arccos \ n • d \,
. (1 +cos(2Of) )
Yf = mm I -2- , 1 ) •
This derivation was carried out for one-dimension. To extend the formulation to multiple dimensions, not the Courant number at the face, Cf is used but the cell Courant number is introduced
n
CD = max(Cf, 0) [9], where n denotes the number of control volume faces.
f=1
4. Collapse of a water column in two dimensions
Firstly, the breaking of a dam is simulated in two dimensions (Fig. 4) and the results compared to an experiment (results taken from [6]) conducted by Koshizuka in 1995 [5]. The shapes of the
Fig. 4. Experimental set-up [6]
water-air interface were captured by a camera after flowtimes t = 0 s, t = 0.30 s, t = 0.40 s and t = 0.50 s. At those precise times the simulated shapes of the boundary between water and air is plotted. The form of the interface is qualitatively compared.
4.1. Set-up of the experiment
At the beginning of the experiment (see Fig. 4), a column of water is trapped behind a membrane and the remaining space of the domain is filled with air. At time t = 0 s that membrane is removed and the column collapses due to gravity. As the water hits the barrier a complex interface shape develops which is captured by experiment and simulation. As the grid-spacing (see Fig. 5) is very fine with an edge length of a finite volume of Ae = 6 mm the timestep size had to be chosen accurately small with a stepsize of At = 0.00025 s to satisfy the CLF condition [1]. The CICSAM scheme was used to treat the convective term of the transport-equation of the volume fraction a of the water-phase. It was of particular interest if that compressive interface-capturing scheme could resolve the interface accurately. The phases of the two phase mixture are water and air taken from the database in ANSYS-Fluent. The material parameters of the density p and the dynamic viscosity n were directly taken from the internal material-database in ANSYS-Fluent 19.
Fig. 5. Computational mesh
4.2. Results
In the following, the shapes of the water-air interface in the experiment are compared to simulation-results after flowtimes t = 0 s, t = 0.30 s, t = 0.40 s and t = 0.50 s in a qualitative manner. The material properties are given in Tab. 1. The results of the simulation and the
Table 1. Material properties (ANSYS-Fluent internal)
density p, [p] = kg m 3 dynamic viscosity p, [p] = Pas
air 1.225 1.7894e-05
water 998.2 0.001003
comparison with experimental data are shown ion Fig. 6.
As one can see, the simulation led to good agreement with the experiment. Especially the forming of the air pocket on the right-hand side of the domain after the collapsing water hits the right wall was captured well. However, the intense spray of water in the last frame was not captured. Intense spraying of a liquid in a liquid-gas flow will not be of interest since the focus will be on high-viscosity fluids which are unlikely to spray intensely.
5. Squeezing of an arbritary non-Newtonian fluid through a nozzle
A second numerical experiment was carried out to study the capture of non-Newtonian flow behaviour through a nozzle which is comparable to the printing-nozzle used in fused deposition modeling. That non-Newtonian behavior was of the power-law-type [1]
n (Y) = kYn-1,
Y =
with the rate-of-deformation tensor D and the shear rate Y. The investigation was isothermal.
a)
Fig. 6. Comparison of experimental data [5] (left) with simulation (right): a) after t = 0s, b) after t = 0.30 s, c) after t = 0.40 s, d) after t = 0.50 s
5.1. Expermimental set-up
It was desired to model a shear-thinning fluid (so called pseudo-plastic, Fig. 7), so the shear-coefficient was chosen to n = 0 [1] and the coefficient k to k = 4.6852 . It should be pointed out that it was not the goal of this numerical experiment to model a realistic fluid, but the qualitative behaviour of a shear-thinning fluid. The fluid is pushed into the cavity through a velocity-inlet
Fig. 7. Numerical setup
in ANSYS-Fluent with a steady flow velocity of w;nlet = 0.01113ms_1. At the bottom of the cavity, the narrowest opening has a diameter of d = 0.8 mm. All numerical parameters are of the calculation are unchanged from the previous simulation of the collapsing water column.
5.2. Results
At first, characteristic flow-stages are given shortly (Fig. 8). The results show the prescribed non-Newtonian behaviour which is visible from the last two figures on the bottom left and the bottom right. As the fluid passes through the nozzle, high shear rates occur which lad to a decease in viscosity of the fluid. In the last picture on the bottom right, the fluid touches the bottom wall of the domain. The shear rate is minimal at that bottom portion of the domain, which lets to a increase in the dynamic viscosity of the fluid.
6. Simple two dimensional model of the FDM-process
The flow of a non-Newtonian fluid through a nozzle in two dimensions was considered previously and in that simulation the source of the fluid-flow was held stationary, so the grid of computational cells was stationary as well, which is a foundation of the VOF [4]. But an intrinsic property of the FDM process is that the nozzle is moving with time. To model such a behavior in ANSYS-Fluent, an additional program has to be provided by the user which governs the exact movement of the mesh. The numerical procedure might be simplified and the need for an additional program avoided by defining a region of interest and multiple nozzles adjacent to each other.
a) b)
Fig. 8. Numerical results at various flow-times: a) nozzle-setup, b) pushing in, c) pushing through, d) pushing contact
6.1. Set-up of the experiment
The region of interest of this experiment is the green coloured transparent plane in the left-hand side figure (Fig. 9). The grey lines in that picture are representative of the printing path and are derived directly from the code, which itself was made by the slicing program CURA. When the nozzle moves according to that code, it passes through the green highlighted region of interest in its normal direction. On the right-hand side of the picture, five such nozzles are modeled in their initial configuration at the beginning of the simulation. They are already partially filled with a non-Newtonian fluid with a fluid-source at the top edge of each nozzle. The polymer has a temperature of $ = 210°C. The pressure outlet at the top of the domain has a temperature of $ = 26.85°C. Since it is the bottom layer of the printed part which is modeled here, the printing bed temperature at the bottom edge of the computational domain has to be considered as well and it is set to $ = 60°C. One major drawback of the described model of adjacent nozzles is that the distance between nozzles is far too wide compared to the distance of two adjacent printed
Fig. 9. Region of interest (green plane) and simplified numerical setup
paths. The function for temperature dependence of the dynamic viscosity was taken from the ANSYS-Fluent internal database [1].
6.2. Results
The described set-up is very simple and the most important result of it is the interaction of the individually printed strands after they spread on the printing bed and make contact with each other, as it can be seen that they hinder one another in a free flow (Fig. 10). Since the
Fig. 10. Numerical results
entire domain is modelled symmetrically, the resulting pattern of flow is symmetrical as well. However, on the bottom edge (the printing domain) numerical errors can be observed, which is possibly a reason of a numerical grid which is too coarse.
Conclusions
In this work, several two dimensional models of a two fluid mixture consisting of a liquid and a gas were considered. To solve this problem, ANSYS-Fluent was applied. To solve the advection equation of the volume fraction of the liquid, the CICSAM algorithm was chosen. A simplified two dimensional model of manufacturing a bottom layer by FDM was considered. Numerical
errors occured in the region where the liquid made contact with the heated bottom wall of the domain. In future works, those models will be extended to three dimensions with proper material constants suited for the materials handled in fused deposition modeling and refined meshes will be employed.
The present work is financially supported by the European Regional Development Fund (ERDF) and the Investitionsbank Sachsen-Anhalt (procedure number ZS/2020/01/105113). This support is gratefully acknowledged.
References
[1] ANSYS 13 Handbook. Ansys, Inc., 2007.
[2] B.Behdani, M.Senter, L.Mason, M.Leu, J.Park, Numerical study on the temperature-dependent viscosity effect on the strand shape in extrusion-based additive manufacturing, Journal of Manufacturing and Materials Processing, 4(2020), no. 2, 46.
DOI: 10.3390/jmmp4020046
[3] R.Comminal, M.P.Serdeczny, D.B.Pedersen, J.Spangenberg, Numerical modeling of the material deposition and contouring precision in fused deposition modeling, In Proceedings of the Annual International Solid Freeform Fabrication Symposium, US, Austin, 2018, 18551864.
[4] J.H.Ferziger, M.Peric, R.L.Street, Computational Methods or Fluid Dynamic, Springer, 2020.
[5] S.Koshizuka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dynamics J., 4(1995), 29-46.
[6] J.Marino-Salguero, M.Schafer, A modified normalized weighting factor method for improving the efficiency of the blended high-resolution advection schemes in the context of multiphase flows, Experimental and Computational Multiphase Flow, 3(2021), no. 3, 208-255. DOI: 10.1007/s42757-020-0074-2
[7] J.E.Pilliod, E.G.Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, Journal of Computational Physics, 199(2004), no. 2, 465-502.
[8] O.Ubbink, Numerical prediction of two fluid systems with sharp interfaces, Phd thesis, Imperial College, 1997.
[9] T.Waclawczyk, T.Koronowicz, Modeling of the wave breaking with CICSAM and HRIC high-resolution schemes, In P. Wesseling, E. Onate, J. Periaux, eds., European Conference on Computational Fluid Dynamics ECCOMAS CFD 2006, 2006.
[10] T.Waclawczyk, T.Koronowicz, Comparison of CICSAM and HRIC high-resolution schemes for interface capturing, Journal of Theoretical and Applied Mechanics, 46(2008), 325-345.
[11] Huanxiong Xia, J.Lu, G.Tryggvason, Fully resolved numerical simulations of fused deposition modeling. Part II - Solidification, residual stresses, and modeling of the nozzle, Journal of Manufacturing and Materials Processing, 24(2018), no. 6, 973-987.
DOI: 10.1108/RPJ-11-2017-0233
[12] H.Xia, J.Lu, S.Dabiri, G.Tryggvason, Fully resolved numerical simulations of fused deposition modeling. Part I, Fluid flow Rapid Prototyping Journal, 24(2018), no. 2, 463-476.
Численное моделирование наплавленного осаждения для производства деталей с высокой геометрической точностью и механическим качеством
Хольм Альтенбах Габор Янига
Университет им. Отто фон Герике Магдебург, Германия
Марио Байнер Рене Андрош Поль-Максимилиан Рунге
Галле-Виттенбергский университет им. Мартина Лютера
Галле (Заале), Германия
Аннотация. С увеличением использования аддитивных методов производства механических деталей также возрастает потребность в точном и надежном моделировании производственного процесса. В этой статье различные расчеты, пригодные для моделирования процесса наплавленного осаждения, рассматриваются в двух измерениях. При моделировании наплавленного осаждения расплавленный полимер укладывается по заданной траектории до начала охлаждения расплава. Возникающие потоки рассматриваются как многофазные. Для моделирования отложения нити в Ansys Fluent используются методы вычислительной гидродинамики, а именно метод объема жидкости (VOF). Проведено моделирование различных численных экспериментов.
Ключевые слова: объемный жидкостный метод (VOF), неньютоновская жидкость, двухфазный поток.
