LAUREA MAGISTRALEQuesta tesi di Laurea Magistrale indaga l'applicazione del Metodo degli Elementi Finiti con Particelle (PFEM) per affrontare le sfide nella stampa 3D del calcestruzzo (3DCP), concentrandosi in particolare sui modelli di instabilità dei fluidi di Bingham durante il processo di stampa. Sfruttando le capacità del PFEM per la simulazione della fluidodinamica non newtoniana, questo studio fornisce un'analisi dei fattori che influenzano la stampa 3D del calcestruzzo basato sull'estrusione, comprese le proprietà dei materiali e i parametri di estrusione, come velocità, posizione e dimensione del dispositivo di estrusione.
Attraverso il confronto tra simulazioni numeriche e risultati sperimentali presenti in letteratura, la tesi presenta riflessioni sull’ottimizzazione dei parametri di stampa 3D per una migliore integrità strutturale e accuratezza geometrica dei costrutti stampati.
I risultati ottenuti contribuiscono a prevedere correttamente il processo di deposizione del filamento conoscendo alcuni parametri geometrici e cinematici chiave.This Master's thesis investigates the application of the Particle Finite Element Method (PFEM) to address challenges in 3D Concrete Printing (3DCP), particularly focusing on the instability patterns of yield stress fluids during the printing process.
By exploiting the PFEM's capabilities for simulating non-Newtonian fluid dynamics, this study provides an analysis of the factors influencing the extrusion-based 3DCP, including material properties and extrusion parameters, such as velocity, position and dimension of the extrusion device.
Through the comparison between numerical simulations and experimental results from literature, the thesis presents novel insights into optimizing 3D printing parameters for improved structural integrity and geometric accuracy of printed constructs.
The findings contribute to predict correctly the filament deposition process knowing some key geometric and kinematic parameters

Minora, Adriano

English

POLITesi - Archivio digitale delle tesi di laurea e di dottorato (Politecnico di Milano)

Executive Summary of the ThesisNumerical Simulation of Flow Instability Patternsin 3D Concrete Printing with PFEMDouble Degree Master Thesis in:Civil Engineering - Ingegneria CivileMathematical Engineering - Ingegeria MatematicaAuthor: Adriano MinoraAdvisor: Prof. Massimiliano CremonesiCo-advisor: Giacomo RizzieriAcademic year: 2022-20231. IntroductionThe adoption of 3D Concrete Printing (3DCP)represents a significant advancement in the fieldof construction. However, the inherent chal-lenges associated with managing instability pat-terns during the extrusion of yield stress fluidsis an obstacle to achieving optimal outcomes.The aim of this work is to show the validity ofthe Particle Finite Element Method (PFEM) asa tool for investigating the complex instabilitypatterns that can arise in additive manufactur-ing processes.By exploiting the PFEM’s capabilities in simu-lating non-Newtonian fluid dynamics, this studyexplores the interplay between material proper-ties and extrusion parameters in the context of3DCP.The starting point of this contribution is thework done by Geffrault et al. in [3]. This pa-per proposes a so-called universal mapping plot :a very simple yet powerful tool for determiningthe instability patterns in 3DCP based on fewgeometric and kinematic parameters. Numericalsimulations subsequently presented are intendedboth as an example of how the PFEM can be ef-fectively used in the context of 3DCP and as amean to further discuss the validity of such auniversal mapping plot.2. Extrusion-based 3DCPThe existing additive manufacturing approachesinvolving cement-based materials can be clus-tered in several groups. Among these, the meth-ods that have shown particular success are mate-rial extrusion, particle bed binding and materialjetting. In particular, extrusion-based 3DCPstands as a primary solution both in researchand first applications worldwide.2.1. Types of extrusionFocusing on the extrusion-based additive man-ufacturing, three situations can be idealized,based on the characteristics of the material tobe extruded:1Executive summary Adriano Minora(a) (b) (c)Figure 1: Comparison between different filamentcross sections based on the material properties.• Infinite brick extrusion of stiff material,where the extrude filament cross sectionmatches the extrusion device shape (see 1a);• Free-flow deposition of a flowable material,where the extruded material loses its shapesand it spreads over the substrate (see figure1c);• An intermediate situation, in which thefilament cross-section undergoes a certaindegree of deformation, without completelylosing its shape (see figure 1b).The letter case is the most frequent in practicalapplications and the one faced in this work. Inthis scenario the printed material is typically ayield stress fluid.2.2. The effects of GravityFollowing [4], the effects of gravity on the ex-truded filament can be quantified by the ratiobetween the gravitational stresses ρgh and theyield stress of the fluid τ0.Here, ρ is the material density, g the accelerationdue to gravity, and h the height of the depositedlayer.Based on these quantities, two limit scenarioscan be depicted, depending on their relativemagnitude:• Dominant yield stress τ0 ≳ ρgh:The layer of cementitious material remainsunaltered following deposition. This situa-tion is applicable when the material is care-fully and relatively slowly placed.• Dominant gravity stress τ0 << ρgh:The material undergoes deformation until itachieves a balance between the yield stressand gravity-induced stresses and the newlydeposited layer spreads over the previouslyprinted layer.2.3. Layer pressing and free layersFigure 2 shows a sketch of the different casesthat can arise based on the vertical position ofthe nozzle in the printing process.Figure 2: Example of free layer extrusion (left)and layer pressing extrusion (right)If the distance from the nozzle to the base planeis safely larger than the nozzle diameter, thenthe filament can freely bend and no flowingoccurs. As a consequence the filament cross-section remains essentially constant during theprocess. This scenario is referred to as free layer.On the other hand, if the distance from the noz-zle and the plane is less or equal to the nozzlediameter, a certain degree of flowing is neededby the filament in order to squeeze between thenozzle and the base plane. This case is referredto as layer pressing.3. Instability PatternsThe work by Geffraul Bessaies-Bey, Roussel andCoussot [3] establishes a general rule for pre-dicting instability patterns in 3D printing usingyield stress fluids.This rule is summarized in the universal map-ping plot, depending on two non-dimensionalvariables that condense the influence of nozzlediameter, nozzle height, printhead horizontal ve-locity and extrusion velocity. The plot is derivedby the analysis of several experiments involvingthe printing of a yield stress fluid filament on aplane.In this paragraph, the main concepts that leadto these findings are reported.3.1. Fundamental hypothesesFollowing the procedure depicted in [3], the sim-plest case is studied, so that the influence ofcomplex factors can be excluded. The main as-sumptions under which this analysis holds arethe following:1. Dominant yield stress flow,2. Negligible surface tension effects,3. No layer pressing,4. Single layer printing process,5. Circular nozzle cross-section2Executive summary Adriano Minora3.2. Experimental settingsThe yield stress model fluid used in the exper-imental campaign is a Kaolin suspension. Bychanging the Kaolin concentration it is possibleto obtain different values of yield stress, withoutvarying significantly the viscosity and the den-sity. The suspensions adopted has yield stressvalues ranging from 150 Pa to 1300 Pa, and den-sities from 1390 kg/m3 to 1590kg/m3.The experimental device is composed by a com-pressor, a nozzle with interchangeable outletshapes, and a moving plane. In a real appli-cation scenario it is usually the nozzle to moveabove the substrate and not vice-versa. This setup is justified by the hypothesis that for com-mon velocities used in 3DCP inertia effects arenegligible, so that the only horizontal velocity toplay a role is the relative velocity between nozzleand plain.In such a setting, the control variables are thefollowing:• The height of the nozzle from the plane H,• The diameter of the nozzle D,• The extrusion velocity Ve,• The translational relative velocity betweenthe plane and the nozzle Vt,• The rheological characteristics of the fluid,condensed in the yield stress τ0.Geometric and kinematic quantities can be con-densed in two non-dimensional variables:H∗ =HD, V ∗ =VtVe.3.3. Universal mapping of patternsThe patterns obtained by a various combinationof the control variables can be categorized in twomain distinct groups (see figure 3).If the height of the nozzle is higher than a certaincritical height Hc, then the filament breaks be-fore touching the plane and drops are formed.On the contrary, if the height of the nozzledoes not exceed Hc, then the filament can reachthe plane without braking and consequently fivepatterns can form:• Discontinuous line: due to the transla-tional velocity the filament shrinks at cer-tain sections, where eventually it brakes.• Straight line: the filament maintains astraight line and a constant cross section.• Meanders: the filament presents some pe-riodic waves in the horizontal plane withoutFigure 3: Universal mapping of patterns as afunction of the non-dimensional variables (usinghere 1− εc = 0.8).any overlapping of material.• Alternated loops: the filament undergoesa periodic pattern presenting loops that fallalternatively on one side or on the other.• Translated loops: the periodic loops ap-pear to be translated in space.From the analysis of the experimental resultsperformed in [3], the universal mapping plot offigure 3 is obtained. The frontier curves are ex-pressed as a function of the critical height Hcand the correlated critical elongational strainεc.4. The PFEM strategy for3DCPIn order to perform numerical simulations andassess the validity of the universal mapping plot,an appropriate numerical method is needed.The Particle Finite Element Method is a versa-tile method for solving fluid dynamics problems,particularly efficient for problems regarding free-surface flows, breaking waves, fluid–structure in-teractions, and fluid–object interactions [1, 2].4.1. Generalities on PFEMThe strength of this method stands in its ca-pability to combine the adaptability of a La-grangian approach with an effective remeshingalgorithm (more details in [1]). Essentially,each node of the computational mesh representsa physical particle, carrying its own inheritedproperties. When the method is applied to flu-ids, nodes of the mesh can exhibit large displace-ments with respect one another, rapidly degrad-ing the mesh quality.An efficient re-meshing algorithm is needed inorder to restore a valid computational mesh3Executive summary Adriano Minorawhen needed. Each time a re-meshing is per-formed a new boundary identification process iscarried out. This peculiarity, obtained by meansof the Alpha Shape method, makes the PFEMparticularly valid when dealing with changingtopology.4.2. Mathematical formulationThe fluid dynamics problem is governed by theNavier-Stokes equations. Being Ωt an evolvingdomain for t in [0 : T ], denoting by u = u(x, t)the velocity of the fluid, σ = σ(x, t) the Cauchystress tensor, b = b(x, t) the external bodyforces and ρ the density, Navier-Stokes equationsread: {ρDuDt = divσ + b in Ωt × [0 : T ]DρDt + ρdivu = 0 in Ωt × [0 : T ]where DDt is the material derivative.The Cauchy stress tensor can be written in termsof its isotropic and deviatoric parts:σ = −pI+ τwhere p = p(x, t) is the pressure field and τ isthe deviatoric stress tensor.The constitutive law usually has the form of:τ = µ(ε)ε,where µ is the viscosity, which in general can befunction of the symmetric part of the strain ratetensor ε.The system is closed by the following compati-bility equation:ε(u) =12(∇u+∇uT )− 13Idivu.The discretized nodal velocity and nodal pres-sure vectors can be defined as:ui(x, t) =∑aNua (x)Uia(t)p(x, t) =∑aNpa (x)Pa(t),where Nu and Np are linear shape functions anda stands for the a− th element in the computa-tional domain. The Navier-Stokes equations interms of the nodal vectors read:{MuU̇+KU+DTP = FMpṖ+ κDU = 0where Mu and Mp are the mass matrices of ve-locity and pressure, K is the fluid stiffness ma-trix, D is the discretized divergence operatorand κ is the bulk modulus of the fluid. VectorF accounts both for body forces and boundaryconditions. The definition of mass matrices andstiffness matrix are:Mu(i,j)(t) =∫ΩtρNui(x) ·Nuj(x) dΩMp(i,j)(t) =∫ΩtNpi (x)Npj (x) dΩK(i,j)(t) =∫Ωt∇Nui(x) : ∇Nuj(x) dΩSince the equal order (linear) finite elements donot satisfy the inf-sup condition, a proper stabi-lization method is needed.The time discretization follows an explicit Eu-ler’s scheme, so that the system of equation be-comes:{Muk+1Uk+1 = MukUk +∆tkRkMk+1p Pk+1 = MkpPk −∆tkκDkUkwhere Rk = Fk −KkUk −DTPk and the timestep ∆tk has to respect the Courant-Friedrichs-Lewy stability condition.Mass and stiffness matrices depend on the cur-rent domain configuration, so they are constantin time between one re-mashing and the other.Moreover, with appropriate element basis func-tion they are also lumped and an efficient par-allelization can be adopted (details in [2]). Themathematical formulation is complete with anappropriate set of boundary condition.4.3. Computational toolsEach time the mesh is regenerated the new do-main boundaries have to be identified, this isdone by the Alpha Shape method. After con-structing a Delaunay Triangulation a convex do-main is obtained. For retrieving non convexshapes, this method sets a threshold to the el-ements distortion and it excludes elements thatare beyond this threshold (see figure 4). The dis-tortion is generally defined as the ratio betweenthe radius of the circumsphere of the elementand a reference mean dimension of the elements:αe =Rehmean4Executive summary Adriano Minora(a) ᾱ −→ ∞ (b) ᾱ=20(c) ᾱ=10 (d) ᾱ=1Figure 4: Alpha Shape: decreasing values of thethreshold ᾱ lead to smaller outer angles.One peculiarity of the PFEM is the ability todeal with isolated nodes that can join or leavethe fluid mass based on their position. It is alsopossible to redefine the nodal positions in orderto remove nodes where their concentration is toohigh in favour of more sparse areas.The contact with the base plane, positioned atz = 0, is enforced by monitoring the points co-ordinates. When the z coordinate of a point isclose to zero, the node is set to be bounded. Thiseliminates the need to model the base plain withphysical points, procedure that would generatefictitious fluid volume near the contact area.In order to improve the computational efficiencyof the method for 3DCP problems, an ad hocprocedure for lowering the number of nodes tobe solved has been implemented. Practically,nodes that have negligible velocities are seen asbounded and so excluded from the solving phase.This is done by freezing nodes horizontally farfrom the nozzle (see figure 5).4.4. Validation of PFEMFor validating the PFEM for non-Newtonian flu-ids, five slump tests have been simulated. Re-sults in terms of radius and height of the finaldeformed shape are compared with the corre-sponding values computed with the formula pro-posed by Roussel and Coussot in [4]. Figure 6Figure 5: Frozen and active regions of the fluid(a) t = 0 s (b) t = 0.2 s(c) t = 1 s (d) t = 5 sFigure 6: Slump simulationshows an example of a slump simulation, whiletable 1 reports the comparison between simu-lated results (R and H) and theoretical results(R, H).5. Numerical SimulationsTwelve numerical simulations have been carriedout. The relative values of the control vari-ables and of the non-dimensional variables arereported in table 2. These simulations repro-duce the actual situation in which the nozzle ismoving over a fixed base plane. For this numer-ical campaign, only line patterns are analyzed,excluding the case of falling drops.G1 G2 G3 G4 G5 AverageR [cm] 26.9 26.4 24.9 24.5 24.4R [cm] 26.3 25.6 23.9 24.3 23.5Error 2.5% 3.0% 4.1% 1.0% 3.7% 2.8%H [cm] 5.49 5.63 5.71 5.39 5.32H [cm] 4.75 4.98 5.75 5.55 5.9Error 15.7% 13.1% -0.6% -3.0% -10.4% 2.9%Table 1: Slump test: comparison between simu-lated values and theoretical values.5Executive summary Adriano MinoraFigure 7: Simulation resultsThe positions of the simulations in the universalmapping plot is shown in figure 7, discerning thecases for which the simulated pattern is coher-ent with the position in the plot (circle mark-ers), from the cases in which different patternsare obtained ("x" markers). One did not showany repetitive pattern ("?" marker). This is dueto the particularly long transient phase, not yetvanished at the end of the numerical simulation.The obtained patterns are compared with theexpected patterns in table 3.6. ConclusionsIt can be stated that the universal mapping plotis a useful tool for predicting the instability pat-terns of a yield stress fluid extruded filament. Itis also evident that the PFEM proves to be ade-quate in simulating such a scenario. Still, somediscrepancies between numerical results and thereference experimental campaign exist and theyincrease when the scale of the problem rises.Numerical simulations showing different resultsτ0 [pa] H [cm] D [cm] Vt [cm/s] Ve [cm/s] H∗ V ∗C1 150 1.0 0.5 2.0 2.0 2.00 1.00C2 150 2.0 0.5 2.0 5.0 4.00 0.40C3 150 3.0 0.5 2.0 5.0 6.00 0.40C4 300 1.5 1.0 3.0 1.5 1.50 2.00C5 300 2.5 1.0 3.0 6.0 2.50 0.50C6 300 4.0 1.0 1.0 3.3 4.00 0.30C7 700 5.1 1.7 4.0 4.0 3.00 1.00C8 700 5.1 1.7 2.0 6.0 3.00 0.33C9 700 8.0 1.7 0.5 2.5 4.70 0.20C10 1300 6.0 3.0 5.0 3.0 2.00 1.67C11 1300 14.0 3.0 5.0 7.5 4.67 0.67C12 1300 14.0 3.0 5.0 6.5 4.67 0.77Table 2: Simulations summaryExpected pattern Simulated patternC1 Straight Line Straight LineC2 Alternated Loops Alternated LoopsC3 Alternated Loops UndefinedC4 Discontinuous Line Discontinuous LineC5 Meanders MeandersC6 Alternated Loops Translated LoopsC7 Straight Line Straight LineC8 Alternated Loops Alternated LoopsC9 Translated Loops Translated LoopsC10 Discontinuous Line Discontinuous LineC11 Meanders Alternated LoopsC12 Meanders MeandersTable 3: Expected and simulated patternsfrom what expected present high values of ve-locities or height. Probably, even in this lim-ited range of velocities and heights, some iner-tia effects could play a role. If so, it is soundto consider numerical simulations more reliable,as they intrinsically consider inertia effects andthey are able to well capture transient phases.Anyway, when velocities rise, the computationalcost of simulations becomes sensible, posing apractical limit to this approach exactly wherethey seem to be more accurate. Future develop-ments could find solutions for reducing the com-putational costs of simulations.References[1] M. Cremonesi, A. Franci, S. Idelsohn, andE Oñate. State of the Art Review of theParticle Finite Element Method (PFEM).Archives of Computational Methods in En-gineering, 27:1709–1735, 2020.[2] M. Cremonesi, S. Meduri, U. Perego, andA Frangi. An explicit Lagrangian finite el-ement method for free-surface weakly com-pressible flows. Comp. Part. Mech., 4:357–369, 2017.[3] A. Geffrault, H. Bessayes-Bay, N. Roussel,and P. Coussot. Printing by yield stress fluidshaping. Additive Manufacturing, 75:103–752, 2023.[4] N. Roussel and P. Coussot. “Fifty-centrheometer” for yield stress measurements:From slump to spreading flow. Journal ofRheology, 49:705–718, 2005.6

Numerical simulation of flow instability patterns in 3D concrete printing with PFEM

https://www.politesi.polimi.it/retrieve/6d93a90b-8d6e-46d1-95a6-76cbaa5c0d22/2024_04_Minora_Executive%20Summary_02.pdf

Numerical simulation of flow instability patterns in 3D concrete printing with PFEM

Abstract

Similar works

Full text

Available Versions

POLITesi - Archivio digitale delle tesi di laurea e di dottorato (Politecnico di Milano)