Grading Threads. Exploiting Viscous Thread Instability for the additive fabrication of Functionally Graded Structures via sensor-adaptive robotic control

We present a novel computational fabrication method for the production of Functionally Graded Structures (FGS) via robotic control of Viscous Thread Instability (VTI). Of interest in several fields and at different scales of application, the fabrication of FGS is often relying on offline fabrication workflows and on stable material conditions. By introducing partial control in the process of spatial deposition of an extruded clay thread in a state of instability, our method extends the design and fabrication possibilities of VTI to the production of FGS. Traditionally exploited for the industrial production of not-graded two-dimensional nonwoven textiles or for surface treatments in design-related 3d printing applications, we frame VTI as the main design and fabrication driver for the computational fabrication of functionally graded clay volumetric structures. Without relying on predictive physical simulation models, our method relies on feedback information provided by sensing equipment in combination with an industrial 6 axis robotic manipulator integrated with a numerically controlled clay extruder. The sensed information is used to retroactively update the inputs of a computational model programmed to guide the robotic additive fabrication of user-defined functional volumetric gradients. We illustrate the main design- and fabrication-related parameters and a set of material experiments designed to validate the accuracy of our model. We present a set of fabricated outputs to illustrate the flexibility of the model to accommodate a variety of design intentions and, finally, we discuss its potential for further research involving cross-scalar and trans-disciplinary applications.


Introduction
The fascination for complex processes induced by viscous material instabilities has well-known precedents in the broad field of the Visual Arts, both two-and three-dimensional. Connected to the use of flowing viscous paint as medium of choice, the research focus of artists such as Jackson Pollock and Perry Hall relates, at a deeper level, to the establishment of a form of soft control over the combined effects of motion and matter, constantly mediated by the artist's sensory supervision. Rather than forcing a pre-conceived shape on a passive and isotropic material substrate, both authors privilege the engagement in open-ended processes of guided material exploration, directly driven by non-linear and active material behaviors. Similar strategies have been developed by sculptors such as Richard Serra, Anish Kapoor and Eva Hesse. The behaviors of lead, hot wax, fresh concrete and free-falling ropes act as morphogenetic drivers informing a sculpting process proceeding "from the inside-out" [1].
In the last decade, similar yet more industry-oriented approaches to material investigation led to the development of several research projects in the field of Architectural Design and Construction, often situated at the intersection of disciplines such as Architectural Geometry, Computer Graphics, Robotics and Material Sciences [2,3]. In this context, the application of one-dimensional viscous-or elasticphase material relates to continuous, medium-to-large scale extrusion-based additive manufacturing processes exploiting the deposition of materials such as concrete, thermoplastic polymers, clay, water-based compounds and glass. During a 3d printing process involving viscous-phase materials, when the layer-to-layer vertical distance crosses a certain -variable, according to process and material parameters-threshold value, the extruded thread will snap out from its linear configuration, deforming in a variety of complex coiling and looping patterns. Causes and effects of this meandering behavior have been main subject of investigation of prior research in Physical Mechanics [4,5]. This complex material phenomenon involving a compound effect of gravity, viscosity, elasticity, surface tension and velocity is referred to as Viscous Thread Instability.
For understandable reasons, standard additive manufacturing applications tend to minimize any visible and/or structural effect of viscosity. A commonly undesired related effect is the progressive loss in printing height. This can be caused by a faulty calibration of the deposition process resulting in excessive downward pressure of the extruded material (over-feeding). Yet, the deformation under selfweight is an inherent property of 3d printing with viscous materials: as the height and the weight of the printed part increases, the relative layer-to-layer vertical distance also tends to increase, creating opportunities for the extruded coil to buckle under the effects of gravity. This behavior, if not counterbalanced, may lead to visible surface imperfections, large deformations or, in extreme cases, structural collapse of the printed part. In different contexts, more creative and experimental applications look at the gravitational effects produced by viscous materials as a design affordance that can enhance the aesthetic and/or functional qualities of the 3d printed object. These material effects are often the result of a linear (as opposed to circular, or feedback-based) computational process, specifically designed to parametrically generate one or more 3d printing toolpaths [6].
In recent years, research in Computer Graphics has been improving numerical simulation models meant for visual rendering and prediction of physical phenomena involving thread instability, from the flow of liquid jets to the kinematics of hair and fur [7,8]. The models underlying such simulation tools already proved successful integration in several architectural design-to-fabrication workflows involving elastic-phase construction materials [9]. Nevertheless, the models developed in the Computer Graphics area are usually not meant to promote meaningful design-related interactions, and their porting within 3d design and/or fabrication environments is still an expensive and specialized task [10]. Acknowledging their potential to play a fundamental role in the design process in future iterations of our research, we explicitly opted not to rely on simulation models in the context of this work. We instead decided to focus on a computational design strategy that relies on adaptation enabled by sensed information.
In order to cope with the complex behavior of materials in the context of continuous automated processes, traditionally the Manufacturing and, more recently the Construction sectors, have been developing Sensor-Adaptive strategies applied to processes of fabrication. Such applications focus on the design of sophisticated online production workflows which directly depend on the use of sensing equipment to provide real-time error-correction to meet production tolerances [11,12]. Other cases, usually situated in the Architecture and Construction area, focus on highly experimental, usually open-ended and retroactive computational fabrication processes to explore and exploit uncertainties of different material nature [13,14]. Following a similar research path, we frame Viscous Thread Instability as ideal testbed to adopt a methodological approach that still heavily depends on computational procedures, but which privileges sensing over predicting and adaptation over anticipation. Most importantly, our method actively exploits the aggregate effects of material behavior -namely, VTI-for the fabrication of volumetric Functionally Graded Structures.
In the following chapters we will reference our work to state-of-the-art research and applications within a broad multi-disciplinary context (Chapter 2). We will give an overview of the experimental hardware and software setups (Chapter 3) and we will expand on their inner functionalities by explaining and motivating their design-(Chapter 4), material-(Chapter 4), sensing-(Chapter 5) and fabricationrelated (Chapter 6) control functions and parameters. We will present a set of material outputs to validate our methods and assumptions (Chapter 7). We will conclude by discussing current limitations of our system and by envisioning future perspectives for our work (Chapter 8).

Related work
Our work stands at the intersection of several fields of both fundamental and applied research, including Architecture and Construction, Computer Graphics, Industrial Manufacturing and Material Sciences. We reference state-of-the-art research from such disciplines, organized in four topic-based paragraphs: Functionally Graded Structures, Extrusion 3D Printing, Applied Viscous Thread Instability and Sensor-Adaptive Fabrication.
Functionally graded structures Single-material or composite structures that display a spatially heterogenous material distribution are usually referred to as Functionally Graded Structures (FGS). Numerous applications in fields ranging from mechanical, optical, bio-medical engineering to robotics and aerospace benefit from the use of engineered materials that exhibit the same graded properties at the microstructural scale. They are known as Functionally Graded Materials (FGM). The use of various high-resolution additive manufacturing technologies represents the ideal production method for the fabrication of structures with a complex and small internal geometric distribution of porosity [15]. Applications at the desktop scale, usually produced using FDM 3d printers, are developed in the Computer Graphics community to prove the efficiency of computational methods for the design of internal lattice structures to control the elasticity of free-form physical models [16]. In the field of Architecture and Construction, we reference the seminal work done at ILEK at the University of Stuttgart on large-scale graded concrete structures. Guidelines for applied material research and development [17] as well as full scale prototyping and characterization of several production methods have been published in the last five years [18]. This reference clearly indicates large research potential existing in the development of new production processes involving a deposition of fiber-like materials. Explicitly referring to extrusion-based additive processes, specific hardware solution for the control of production parameters such as fiber orientation, density and mass-volume fraction, should be developed and closely tied to coherent computational design strategies [17].
Extrusion 3D printing Acknowledging the positive impact that research around clay 3d printing has on larger scale concrete 3d printing [19], we restrict the field of investigation to applications involving these two materials in the context of Architecture and Construction. Here, applications involving one-dimensional threads of viscous materials often relate to extrusion-based additive processes. Traditionally, computational strategies for the generation of 3d printing toolpaths rely on the assumption that the digital input and physical output are, within negligible tolerances, identical. When this assumption does not apply, efforts are made to close such error gap, especially in large-scale extrusion processes involving viscous material mixes. As an example, the vertical deformation under self-weight is a known fabrication issue [20,21] that can lead to mild negative effects such as surface imperfections, to severe structural implications such as insufficient layer adhesion or structural collapse of the printed object. To tackle these critical effect, several offline and online strategies have been investigated, from fabrication-aware design environments [10] to realtime feedback [22]. Other research in the field of Design and Architecture aims at the production of medium sized parts, usually designed to withstand lower or null external stress. These conditions allow to look at the gravitational effects produced by viscous materials as a design opportunity to enhance functional and/or aesthetic qualities of the 3d printed objects. These material effects are usually obtained with linear (i.e., not involving feed-forward from simulations nor feed-back from sensing) parametric design strategies for the printing toolpath, such as uniform slicing [6], 3d discretization [23] or agent-based simulation [24]. Even if these examples leverage viscous instability, the resulting effects are always fully controlled in terms of material response and usually localized on the object's shell, affecting the visible surface with two-dimensional textured effects.

Applied viscous thread instability
We previously mentioned the interest of the Computer Graphics community for investigating the effects of viscous and elastic instability. Advancements in this field are of deep interest for future developments of our research, but, as previously stated, we excluded them from the scope of our current application. Prior research in the field of Experimental 3d printing addresses VTI for the design and production of patterned glass artifacts via a simplified geometric model [25]. Oriented to smaller scale applications, [26] and [27] set a relevant precedent as well as a context for the application of VTI in the production of open cellular foam structures by using extrusion-based 3d printing. A material characterization of thermoplastic polyurethane and composite siliconebased materials is provided, including measurement of the stress/strain ratio of different geometries and coiling patterns under both tension and compression, as well as a quantification of their densities and pore sizes. Our work expands the scale of application to an architecture-related context and follows up on some relevant aspects underlined by these references. Namely, the implementation of a process that can adaptively combine coiling densities also in, but not limited to, the vertical direction. In fact, by introducing a 6-axis robotic manipulator coupled with a depth sensor, we expand this possibility to further spatial configurations and density gradients for the production of volumetric FGS.
Applied research in the area of Automation and Industrial Manufacturing proves to be of particular interest for the purposes of our work, especially within the field of textile production of nonwoven fabrics. The direct focus on the application of the effects of viscous instability in a context of automated production, often relying on vision systems for performance assessment, is of great interest for the framing the methodological approach to our work. As illustrated in [28], nonwoven textile products show a wide range of applications (ranging from hygiene ad care items to medical implants to industrial air filters) and manufacturing techniques. Particularly interesting are processes of web-formation by spun-laying -such as electrospinning-which rely on a unified process of extrusion, collection and bonding of a single extruded thermoplastic polymer thread produced from the application of a strong electric field on a resin solution [29]. Crucial in the production of nonwovens are two main parameters informing their material performance: uniformity in fibers distribution and fibers orientation [30]. Several computational methods developed in this field leverage vision-based systems for the assessment of the spun fibers uniformity. Some of them rely on high-resolution scanning procedures to achieve a full three-dimensional description of the fibers [31]. Other faster but less accurate methods leverage statical evaluation of two-dimensional grayscale images based on reflected/transmitted light for the assessment of spatial mass uniformity of fibrous substrates [32]. Given all the design-and material-related differences existing between the abovementioned and our research, we can trace enough technical and methodological similarities to frame our work as a hybrid process enabling a numerically controlled transition between a fiber-spinning and a FDM 3d printing production process. The crossbreeding and the possibility of a modulation between these two fabrication methods enables the production of material structures that exhibit a variable and functionally tunable mass distribution and fiber orientation, enabling the spatial grading of properties depending on material density. The key for the actualization of such a process lies in the development of a common computational framework enabling the interoperability of robotic motion, material deposition and sensing in a continuous cycle of feedback.
Sensor-adaptive fabrication Recent developments in Architecture and Construction are, more and more often, taking advantage of the opportunities granted by a general broader access to affordable hardware for 1d, 2d and 3d environmental sensing as well as to established open-access software libraries for computer vision and image or point-cloud processing. Even if our current work does not rely on Machine Learning techniques, recent advances in the field are supporting the proliferation of sensor-adaptive robotic fabrication applications involving image and/or depth data streams as main design and fabrication inputs. Since our work shares analogous input data types and hardware setup, we reference an overview of such applications [33] to offer a perspective towards future scenarios for our work. Similar applications are conceived as processes either oriented to 1) online or offline error-correction, where design and production goals are either given a-priori within precise production and/or assembly constraints; or 2) open-ended and/or interactive design tasks, where the design space is unfolding during production time. The former often involves robotic setups where real-time online control over the fabrication [11] and/ or the assembly [12] of non-standard architectural scale components is usually the main target. The latter category privileges open-ended experimental design and production workflows, which are strictly enabled by retroactive computational fabrication strategies. Our research contributes to the advancement of sensor-adaptive robotic applications with strong connections to material instabilities, where lack of predictive accuracy in the digital design phase is often the necessary condition to involve sensing to track the spatial progression of the fabricated model. We share similar motivations and material constraints with project related to clay additive manufacturing. In [13] the authors rely on a 3d point cloud representation of the unpredictable behavioral effects provoked by the robotic throwing of soft loam projectiles deforming under the effects of stacking. In [34], a continuous process of spatial clay 3d printing is enabled by the online feedback provided by a low-range LRF sensor detecting local height deflections of the underlying lattice structure. The sensor information feeds back to a parametric model used for online recalibration of the robotic trajectory. The closed loop process is backed up by an ad-hoc software framework for real-time robotic control developed by the authors [35]. Fundamental precedent for our research is the project "Unspecified Clay" [14]. The project -inscribed within the framework of the inFORMed Matter Research initiated by Co-de-iT and the Responsive Matter program of EnsadLab in Paris-demonstrates the possibility of " […] designing a computational system able to negotiate and collaborate with this autonomous expression of clay behavior, with the objective of developing emergent material structures" (p. 235), in which "feedback processes would allow a robotic manufacturing unit to build artefacts through repeated cycles of deposition, scanning, and computation" (p.234).
Conceived as an open-ended process oriented to explore the emergent properties of a material system, "Unspecified Clay" leverages feedback strategies based on sensed point cloud data to address the adaptive and additive fabrication of objects characterized by variable "density, texture and articulation" (p. 238). By focusing on strategies to exploit material behavior rather than to explore it, we frame our work as a technical advancement of a comparable system. We introduce several upgrades in the experimental setup and in the vision processing unit, namely higher sensing resolution, full 6 axis robotic toolpath management, customized electro-pneumatic extrusion control and closed-loop data exchange pipeline between pc and robotic control. Such technical advancements allow us to narrow the research focus and to specifically shift it towards the production of FGS within a fully automated production workflow.

System overview
In order to guide the spatial grading of an unstable clay thread, we developed a computational fabrication system which enables the integrated computational design and digital production of customized volumetric Functionally Graded Structures (FGS). Our system provides a simple user interface within the Rhinoceros 3d modelling environment and allows the parametric design exploration of a series of input functions for the generation of FGS. The design instructions are seamlessly translated to machine instructions for the fabrication process, which relies on a robotic setup able to sense and dynamically react to the material conditions. The additive fabrication process unfolds in consecutive cycles of material deposition, the number and properties of which are described by the parametric input functions. Given the retroactive and parametric nature of the system, the properties of the deposition can range from an uncontrolled state of viscous instability to a fully controlled process of layer-based 3d printing. Therefore, by stacking and/or juxtaposing consecutive material deposition cycles with variable degrees of instability, the process can produce structures with varying degree of aesthetic-and performance-related features. In the following paragraphs we present an overview of the computational design input functions, as well as a description of the robotic fabrication setup.

Computational setup
In the following pages we provide an overview of our computational system by offering a description of the main input functions and by schematically outlining the information flow of our pipeline.
Our computational system relies on three main parametric input functions and on one main numerical input: As illustrated in Fig. 1, each iteration of our workflow is structured as follows: 1. A depth frame of the current working volume is captured as full resolution 3d point cloud. 2. According to VF, the point cloud is clipped and converted to quadrilateral mesh. 3. According to GF, the quadrilateral mesh is downsampled to a specific resolution. 4. According to MDF, each face of the mesh is assigned a scalar value relative to the material density. 5. The robotic toolpath and the deposition instructions are computed and the process executed. 6. The process is repeated from point 1 with updated variables according to VF, GF and MDF. 7. The process is stopped if conditions set by VF, GF or IN are met.
For diagnostic and analytical purposes, at each iteration our system outputs and archives the following information within an organized directory structure: 1. A log file with the main design and fabrication process parameters. 2. A 3d file of the last-captured depth frame. 3. A 3d file of the last-input mass distribution map. 4. A file of the last-executed robot code. 5. A top view of the last deposition round as a raster image.
As visible in Fig. 2, at each iteration (i.e., material deposition cycle) the system allows the tunability of the three input functions to enable the fabrication of a controllable spatial material gradation. Each input function can be tuned by means of dependent and independent input parameters, which we will later discuss in detail by illustrating their operative ranges and material effects, both individually (Chapter 4) and in combination (Chapter 7). In order to visually grasp how the functions and their main parameters are expressed in the computational system, we refer to Fig. 2. The Volume Function is a 3d volume defining the physical boundaries of the graded structure. The Grading Function is an interactive 2d curve which defines the unfolding of Printing Height and Gradient Resolution of each deposition. The Gradient Resolution is expressed as a quadrilateral mesh with a specified face density. The Mass Distribution integrates the quadrilateral mesh by adding a scalar value to each of its faces to indicate the amount of material at each location. The Printing Height defines the distance of the material extrusion from the previous deposition and it is represented as the distance of the robot toolpath from the quadrilateral mesh. In Fig. 2, the three depicted iterations have three different GR/PH values and three uniform MDFs scalar values. A uniform MDF is represented with a uniform mesh face color (Fig. 2c) and translates into an approximately flat material deposition (Fig. 2g). A non-unform MDF is represented with a grayscale mesh color map an translates into an articulated relief-like 2.5d deposition (See Fig. 9).
Our computational workflow is developed within the McNeel Rhinoceros 3d modelling software. It leverages the high degree of customization and automation afforded by the Rhino Common API, as well as its compatibility with the Python and C# scripting languages. The user interface is implemented within the Grasshopper parametric environment. The input parameters for the design, sensing, processing and actuating functions are exposed inside an interactive GUI. Robot kinematics and extruder deposition instructions are encoded in a.pdl file (proprietary Comau Robotics programming language), which is written in the Rhinoceros environment, sent to the Comau Control Unit and automatically translated and executed by a set of.pdl online routines running on the Robot Control Unit. A File Transfer Protocol (FTP) is responsible for the information exchange between the PC (client) and the Robot Control Unit (server) in both directions. Concurrently, we exploit a User Datagram Protocol (UDP) to bridge the deposition instructions from the server/Robot Control Unit to the client/Clay Extruder, activated by a script loaded on an Arduino Uno board.

Fabrication setup
As visible in Fig. 3

Material gradient parameters
Guiding the unfolding of an unstable clay thread under the influence of gravity, motion and pressure, while aiming at the production of a volumetric graded structure, requires a prior basic understanding of the relationship occurring between causes and effects of such instability. Within the operational boundaries dictated by the material conditions at our disposal, we identified three main input functions -and their related parametric domains -able to inform the design and fabrication of volumetric FGS while exploiting VTI. We present a description of the input functions and we support the definition of their operative ranges by presenting empirical data collected from specific material experiments. Additionally, we provide a concise description of how the function parameters are encoded in our computational system.

Printing Height < > Gradient Resolution
The fundamental role that height plays in triggering the instability of a viscous thread is clear and, within this context, it is exploited as the main independent fabrication parameter. Assuming a constant flow rate, a viscous thread falling freely on a planar surface will tend to coil with a curvature radius that proportionally increases with the distance of the thread from the surface. It follows that, if the initial position of the coil is fixed, the coil will tend to cover an increasing, approximately circular, surface area on the said planar surface, until some upper or lower boundaries is reached. On the lower side of the boundary, if the distance drops under a certain threshold -closer to value of the thread cross section diameter-the instabilities will cease, and the coil will degenerate into a point. On the upper side of the boundary, if a certain height threshold is passed, the coil diameter will tend to stabilize around a maximum value and, eventually, tear due to excessive tensional strain. In the context of the current work, we refer to the measure of the approximately circular footprint of the viscous coil as Gradient Resolution (GR). The Gradient Resolution is directly controlled by the Printing Height (PH), defined as the distance that separates the tip of the extrusion nozzle from the closest point laying on the uppermost surface of the previously deposited material layer. Our work relies on the assumption that repeated cycles of material deposition at variable GR/PH over the same spatial region will have a direct impact on the 1) visual, 2) spatial and 3) performative aspects of the produced objects. 1) refers to the size of the material trace left by the coil, consequently affecting the level of resolution that such trace can reach in the visual description of the produced object. 2) refers to the gain in printing height of the produced object. 3) refers to the physical and mechanical properties of the nonwoven object, such as stiffness, elasticity, porosity, density or fiber orientation. Therefore, the control over and the adaptation to the effects produced by the Printing Height / Gradient Resolution parameter stand at the core of our fabrication method. We illustrate an experiment to empirically correlate Printing Height and Gradient Resolution using a curve fitting method. The experimental setup consists of a stationary 2 mm diameter nozzle oriented perpendicularly to the printing surface, through which a 0% chamotte clay mix with 25% water content is extruded at the fixed rate of 1 rpm at 21 point-locations. We release a different amount of material at each location, manually interrupting the deposition after the stacking of 4 stable coils. The nozzle is progressively lowered from a height of 300 mm to a height of 0 mm at discrete steps of 15 mm. As shown in Fig. 4, we measured the resulting footprints by means of a 3d point cloud representation with ground resolution of 0.75 mm (later converted to grayscale heightfield mesh), extracted from an Intel RealSense L515 LiDAR camera module at the maximum depth resolution of 1024 × 768 pixels. Five consecutive frames have been captured and averaged to minimize noise-related artifacts. After clustering the point cloud in 21 individual groups, we chose to isolate a section of the point cloud comprised approximately between 4 and 6 mm from the lower printing surface in order to measure the X and Y dimensions of the bounding box in the centermost section of each deposited volume. By calculating the square root of their product, we reduce each X and Y measurement pair to a single data point, thus associating the approximately circular footprint of each coil to a square footprint. This choice aligns to the design decision of using quadrilateral meshes as the base geometric type of our computational system.
Finally, we correlate Printing Height and Gradient Resolution using a curve fitting process with a seconddegree polynomial function of the type y = ax 2 . The fitting returned the following coefficients: As already mentioned, GR and PH represent Gradient Resolution and Printing Height, both expressed in millimeters. The fitting process is graphically represented in Fig. 5.
It is worth remarking that this equation is calibrated to our specific physical setup. It follows that different extrusion hardware and/or material conditions would require  a new function to be fitted through a new set of measured data points. Given the unstable nature of the material conditions, we performed the described material experiment three times. The measurements of Figs. 3 and 4 are expressing the average of the mentioned three data sets of GR values. The recorded average deviations for each PH measurement are in the order of the tenth of a millimeter. To summarize, within the context of our work, the resulting equation stands at the core of our computational fabrication system, representing the key input for the design of a Grading Function which exploits the effects of Viscous Thread Instability.
As mentioned, the Printing Height is a numeric value expressed in millimeters. At each deposition cycle, it defines the properties of the extrusion process via the Gradient Resolution input values. Our system encodes the dual value of PH/GR within one 2d quadrilateral mesh per each process iteration.  The total surface area of each quadrilateral mesh represents the XY domain defined by the Volume Function and, most importantly, the edge length of each mesh face directly maps to the GR value and it indirectly maps to the PH value. In the context of this work, when referring to a quadrilateral mesh, we always refer to ordered UV Quadrilateral Meshes. Such mesh typology is generated by a regular UV subdivision of a mesh plane and is characterized by an ordered 1D array of faces which all have same surface area and aspect ratio. Therefore, we always assume a unique GR value within a single process iteration. We also assume that, in the following experiments, each individual material deposition cycle is fabricated at the same PH distance from the previously deposited layer. As visible in Fig. 6, the Grading Function is also responsible for the resolution at which the sensed geometry is processed, as it directly defines the downsampling of the original point cloud produced by the depth camera. We will discuss this specific aspect in the paragraph "Effects of the Grading Function".

Grading Function
The Grading Function input is a 2d function that represents the unfolding of the Printing Height/Gradient Resolution values across consecutive iterations of the fabrication process. By defining 1) the number of fabrication iterations as a numerical sequence and 2) the configuration of the said function as a 2d curve, the user can guide the spatial development of the material gradient in the Z dimension by manipulating an interactive 2d curve. Figure 7 illustrates the 2d curves (and the corresponding material effects) describing two simple descending Grading Functions: Linear (left) and Quadratic (right). Both samples have been produced across 5 printing iterations, within a Printing Height / Gradient Resolution domain ranging between a maximum of 100 mm and a minimum of 2 mm. The total mass value per iteration is equal and uniformly distributed within the XY plane at each printing iteration.  Established in the previous paragraph that a single Gradient Resolution map encodes the Printing Height values at each mesh face location during a single iteration of the fabrication process, it is trivial to visualize the Grading Function as a "stack" of consecutive Gradient Resolution maps, represented as an ordered series of quadrilateral meshes. In Fig. 8, we graphically illustrate the unfolding of a Linear GF across 5 iterations by juxtaposing the quadrilateral mesh representation to its material instantiation.

Mass Distribution Function
The Mass Distribution Function (MDF) controls the material distribution across consecutive iterations of the fabrication process in the horizontal plane. It is numerically expressed as a 2d scalar field of integer values populating the spatial regions of the XY fabrication area. The number of spatial regions at each iteration is the direct result of the subdivision of the X and Y dimensions of the fabrication area by the Gradient Resolution value. The values of the scalar field are defined by the user as a list of 3 numbers at each iteration: 1) Total Mass within the entire fabrication area; 2) Maximum Mass percentage; 3) Minimum Mass percentage. We treat mass as a dimensionally unspecified quantity, associating the units of the scalar field to an integer counter representing the robot end-effector visits to the spatial regions. In Fig. 9, we illustrate the material effects of 2 simple MDFs, both calculated over an arbitrarily defined XY fabrication area of 50 × 50 mm 2 at a GR of 5 mm, resulting in 25 square spatial regions measuring 10 × 10 mm 2 . In both examples the Total Mass is set to 175 units, Maximum Mass is 80% and Minimum Mass is 20%. The Printing Height Value is constant and equals 13 mm. The produced objects are the result of a single printing iteration, executed with a 0.84 mm nozzle.
With these experiments we demonstrate how the Mass Distribution Function at low Gradient Resolution values is able to guide Viscous Thread Instability towards the design definition of simple volumetric objects, materialized by coarse material formations. We also demonstrate the possibility of designing the MDF to guide the production of Synclastic and Anticlastic 2.5d relief-like configurations.
As shown in Fig. 9, the color of each mesh face encodes a local mass value. We thus represent the mass values as a scalar field, visually rendered as grayscalecolored map, where white mesh faces represent higher mass concentrations. A grayscale quadrilateral mesh is thus encoding the Mass Distribution Function at one deposition cycle. It follows that, for multiple deposition cycles, the MDFs are represented by multiple, stacked, 2d quadrilateral meshes.
The use of structured data as the main input data type-either represented as spatially organized stacks of quad meshes or bitmaps-facilitates the integration of our system with established computational workflows which inherently handle similar data types. Referring to the Architecture and Construction field, often the output of energy-related analysis is 3d geometry with associated per-face or per-pixel information related to structural, acoustic, or fluid dynamic performance. Furthermore, the use of Image Stacks to represent volumetric data is common practice in the Biological and Medical analysis.

Sensing parameters
To validate the efficiency and the accuracy of our sensing setup -which relies on an Intel RealSense L515 LiDAR camera module -we isolated two crucial parameters affecting, respectively, the pre-and the post-processing side of the sensing output, namely 1) Depth Sensor Resolution and 2) Median Filter Settings. As shown in Fig. 10, we established a ground truth by producing a high-resolution point cloud of a sample object of approximate size 250 × 250x50mm by means of a 7 For processing efficiency, we used the CloudCompare software to decimate the original high resolution point cloud from 3.1 M to 100 K vertices, making sure to keep a maximum deviation from the original below 0.1 mm. we clipped it at a height of 1 mm parallel to the zero plane to minimize the effects of false positives produced by the high correlation of the vertices laying on the planar working area. In order to compare the outputs produced by our sensing setup to the ground truth model, after orienting the point cloud to a shared reference system, we evaluated the smallest point-to-mesh distance at each vertex of the ground truth scan against the polygon meshes produced by different values of the two selected sensing parameters. We used the point-to-mesh mean deviation as comparison metric. To preserve constant ground resolution, all depth frames captured by our setup are performed at the same, fixed distance of 550 mm from the zero plane.

Effects of Sensor Resolution
We tested the three default depth sensor resolution settings-1024 × 768, 640 × 480 and 320 × 240 pixels-in both single capture mode and averaged mode (we captured three consecutive full-frame point clouds and averaged the spatial coordinates of each vertex), for a total of six test cases. Figure 11 shows a grayscale mapping of the three "single-capture mode" mesh outputs and the corresponding color visualization of the deviations rendered on the point cloud vertices. As visible in Fig. 12, the chosen camera settings (resolution 1024 × 768 pixel, single capture) produced a mean deviation from the high resolution point cloud of 1.64 mm. Given a nozzle diameter-and a consequent coil cross section-of 2 mm, this value is considered satisfactory. Furthermore, the chosen settings confine the 10% of the measured points in the upper 60% of the detected distance range (4.76 to 10.2 mm).

Effects of Median Filter
To reduce noise artifacts, our system applies a Convolutional Median Filter to the sensed point cloud before converting it to a regular quadrilateral mesh. The filter  function exposes two main parameters: Number of Iterations, and Kernel Size. By iterating through each vertex of the point cloud, the filter replaces each point Z coordinate with the average of the Z coordinates of a Kernel-Sizewide square region of neighboring vertices. The larger the number of iterations, the more averaged, or smooth, the new point cloud will result. The wider the Kernel, the lesser local shape features will be preserved. We performed several experiments to deduce optimal tradeoff settings between noise reduction, shape preservation and robot pose continuity. We post-processed the same initial point cloud (captured at maximum depth resolution) with four Kernel Size values (3 × 3, 5 × 5, 7 × 7, 9 × 9), applying the filter for a variable number of Iterations (1,2,4,8 times). As in the previous experiment, we computed the mean deviation of the resulting post-processed point cloud (converted to quad mesh with maximum edge length of 1 mm) from the ground truth point cloud. Figure 13 shows that 4 out of the 16 total tests have a mean deviation value larger than 2 mm.
By running some printing tests in a full 6 axis robot configuration and by evaluating Fig. 14, we opted for the 5 × 5 Kernel Size and 2 Iterations setting combination. Since our setup relies on an industrial robotic manipulator involved in a continuous process of material deposition -inherently sensitive to rapid pose changes in the attempt of minimizing deviations from the target velocity-we deliberately sacrificed sensing accuracy for a less stressful machine kinematic response due to a smoother mesh-face normal vector transition ensured by the 2 filter passes.

Effects of the Grading Function
The Grading Function does not only affect the actuation phase, but it also directly impacts the sensing process. As a result, the captured point cloud and the resulting quad mesh are downsampled to a number of spatial regions that matches the one produced by the Mass Distribution Function (See paragraph 4.3). This design decision is coherent to the nature of our adaptive system: for large values of Printing Height/Gradient Resolution, a highly resolute scanned geometry can be considered superfluous, whereas for small PH/GR values, a high resolution mesh representation of the sensed point cloud is necessary for an accurate deposition.

Actuation parameters
We present our fabrication process by briefly introducing the computational strategy for the design of the robotic toolpath and the consequences of its main parameters on the spatial and material grading of the produced Functionally Graded Structures.

Toolpath design
The only geometric input needed for the generation of the robotic toolpath is the 3d quadrilateral mesh produced by the sensing process, downsampled by the Grading Function and which embeds per-face scalar values defined by the Mass Distribution Function. We developed a simple agent-based strategy that produces a single-line toolpath expressing the necessary mass information to guide the material grading process. The algorithm triggers the search process of a single point/agent which, by iteratively visiting all the faces of the input mesh, reduces their mass count by one unit until all faces reach the zero value. We designed the algorithm so that the agent decision privileges the neighboring face with the highest mass value when performing its next step. If all neighboring faces have equal mass value-as depicted in Fig. 15a, b -a user-defined input influences the choice towards either the X or the Y direction. We discuss Fiber Orientation in Paragraph 6.2. For non-uniform MDFs -as in Fig. 15c-the output toolpath is not guaranteed to be Hamiltonian for each possible Mass Distribution Function. Nevertheless, the algorithm can flag the excessive per-face visits of the agent by performing a simple segment length check, and act during fabrication phase. When such exception occurs, in order not to exceed the planned MDF per-face values, our system can either 1) stop the extrusion at that indexed toolpath segment, or 2) increase the speed of the robotic arm at the flagged robot move in order to minimize the mass-related impact of the locally excessive material.
As overall computational design strategy, our algorithm directly acts on the 3d mesh and, therefore, can extract the normal vector for each visited face and use it to construct a 3d plane. Its origin point and Z vector represent sufficient information to compute the robot position and orientation along the printing toolpath in the necessary tool frame, thus expanding to 6 degrees of freedom the spatial possibilities of the robotic fabrication process (see Fig. 15a'').

Fiber configuration and orientation
Two production parameters relate to the possibility of controlling configuration and orientation of the viscous thread by means of the toolpath design. The effects of these two parameters directly inform the spatial fiber arrangement of the nonwoven/graded object, and therefore, its mechanical performance.
The Fiber Configuration is controlled via a single parameter affecting the linearity of the toolpath across the spatial regions/faces of the quad mesh. It is a single normalized floating-point. A value of 0 implies the toolpath to rigorously interpolate the centroids of the spatial regions. The value of 1 implies a maximum displacement of each toolpath waypoint from the centroid of the mesh face towards its boundaries. The maximum displacement value depends on the size of the mesh face (originally depending on the Gradient Resolution values) and equals half of its projected edge length, whereas the direction of the displacement vector is currently based on a random function.
The Fiber Orientation is a simple Boolean value which affects the choice of the next waypoint by privileging either the X or the Y direction. In Fig. 16 we illustrate the combined effect of the two mentioned parameters on the toolpath design (Fig. 16a and c) and on the fabricated objects ( Fig. 16b and d). The depicted material samples are the result of two consecutive iterations of the same fabrication process, executed over a printing area of 80 × 160 mm, using a 0.84 mm diameter nozzle at 12 mm Printing Height and 6 mm Gradient Resolution, and with a uniform Mass Distribution Function. One sample depicts a Fiber Configuration  (Fig. 16a), while the other represents a FC value of 0.75 and a X directed FO (Fig. 16c).

Results
We present a detailed documentation of the computational fabrication process of four volumetric Functionally Graded Structures produced by exploiting clay Viscous Thread Instability. With the first three samples (50 × 50x50 mm 3 produced with the 2 mm nozzle) we demonstrate the possibilities afforded by the manipulation of the Grading Function in the materialization of spatial gradients. In order to obtain a clear comparison, the three samples have been produced with the same design and fabrication parameters, except for the nature of the Grading Function itself. We kept Mass Distribution Function (tot density 100, min 0.5, max 0.5), Fiber Configuration (value 0), Fiber Orientation (X direction), Robot Speed (0.01 m/s) and Robot Kinematic Range (3 axis) constant across the 5 production iterations. In Fig. 17 we illustrate the effects of three different Grading Functions-Linear, Quadratic and Periodic -by capturing the material samples at progressive stages of fabrication and from different angles.
Observing the progressive top views in the three case studies we clearly notice the sequential variation in both Printing Height and Gradient Resolution, as well as a general coherence relative to the designed input functions. The top views-captured at each consecutive printing iteration-show that the curve fitting method used to define our PH/GR equation described in the paragraph "Printing Height < > Gradient Resolution" delivers consistent results. By looking at the iterations with PH lower than 75 mm we notice a regular spacing between the coils in the XY plane and, if we compare the Linear and the Quadratic gradients, we notice a clear distinction between the X and the Y fiber orientations. We assume that the deposition rounds with PH larger than 75 mm do not fully comply to the measured GR because of the combined effect of excessive robot speed and insufficient extrusion pressure.
By relying on experimental results produced with a stationary nozzle, we disregarded the time that the thread needs to fill the height gap between the nozzle tip and the desired spatial region. When the height difference is too large for the current motion/extrusion parameters, the nozzle moves to the next spatial region without having produced a coil in the previous one. This behavior results in a decreased material density and can be noticed in the Iterations 0 of the three samples, where the coiling pattern does not fully express the regularity of the designed toolpath. Observing the front views of the 3 samples we can notice a certain adherence to both the Grading Functions (with noticeable aesthetic manifestations of the material grading in terms of contrasting surface qualities) and the Mass Distribution Function (with a visible vertical progression of the consecutive depositions and a relatively linear horizontal subdivision of the printing volume). We also provide a cross section perpendicular to the fiber orientation of the 3 models as an initial attempt to evaluate the internal fibrous structure and pores distribution. We designed a fourth sample with the aim of producing a small-scale specimen (80 × 40x50 mm 3 using the 0.84 mm nozzle) expressing the full possibilities offered by the parameter space of our system. Differently from the previous three samples, in this experiment we kept most of the input functions and process parameters variable across the fabrication process, such as Mass Distribution (tot. mass variable, min. mass 0.1, max. mass 0.9), Grading Function (Printing Height from 15 to 1 mm), Robot Speed (from 0.005 to 0.02 m/s) and Kinematic Range (full 6 axis motion at Iteration 6), Fiber Configuration (1.0 to 0.0), Fiber Orientation (X and Y swap at Iteration 4).  In order to reinforce the digital-physical relationships enabled by our system, in Fig. 18 we present a more complete documentation. We juxtapose photographic material of the physical artifacts to the input Mass Distribution Function maps, to the sensed and processed geometry, as well as to the progressive steps of the robotic fabrication. We can observe a clear correlation between the top row of Fig. 18, representing the unfolding of the 2D input Mass Distribution Function and the bottom row, depicting the 3D output height gain of the fabricated structure with a side view of the sensed mesh geometry. Such distribution function embeds a material distribution in the XY plane characterized by a decreasing density in the positive X direction (left to right in the image space) and with a constant density in the positive Y direction (bottom to top). A material deposition guided by such a gradient will result in a three-dimensional volume with a decreasing height in the positive X direction and with a constant height in the Y direction.
We grayscale-painted the faces of the scanned mesh geometry to highlight the correlation between the 2D digital input (first row from the top) and the 3D physical object (last row from the top). In the plan views (third row from the top) we illustrate the decreasing linear progression of the Grading Function corresponding to an increase in the 3d printing process resolution and in surface treatment from bottom to top of our graded structure. Additionally, we can notice a change in the Fiber Orientation from Y to X from Iteration 3 to Iteration 4 (see closeup in Fig. 19). The robotic fabrication snapshot at Iteration 5 (second row from the top) offers evidence of a full 6 axis 3d printing process.

Conclusion
We presented the physical outputs of an integrated computational fabrication system that exploits Viscous Thread Instability for the additive robotic production of singlematerial Functionally Graded Structures. Such outputs demonstrated the feasibility of an applied methodology that relies on the acquisition and processing of sensing information to guide the design-to-production process of volumetric structures with tunable aesthetic, spatial and mechanical properties. By synchronizing sensing and actuating-by means of a bespoke feedback-based computational system which includes a depth sensor, a 6 axis manipulator and a numerically controlled extruder-we gained design control over an additive fabrication process able to transition between a controlled, standard layer-based clay 3d printing process and a locally uncontrolled process of deposition under the influence of Viscous Thread Instabilities.
We provided a multi-disciplinary framework for possible applications of similar Functionally Graded Structures and addressed the focus towards large scale 3d printing in the Architecture and Construction field. Even if the documented output structures do not yet engage the architectural scale, we presented a computational framework that enables a potential integration with state-of-the-art architectural design workflows. We supported such statement with two major computational design decisions: 1) the choice of Rhinoceros as a main software environment; 2) the use of structured data (in our case, embedded in ordered UV quadrilateral meshes) as main workflow inputs. We consider both decisions as key aspects for further research developments in the architectural field, particularly when looking at upstream integration of energy-related simulation tools (structural, acoustic, fluid dynamic etc.), which often couple mesh geometries and per-pixel/per-voxel scalar quantities as main output data type.
Finally, we physically demonstrated the suitability of our fabrication setup for the production of small to medium scaled foam-like structures with tunable properties using clay materials. Yet, by indulging on methodological aspects, we addressed the potential of our fabrication system to address material applications not limited to clay: we are indeed confident that the presented experimental framework, if adapted to other viscous material systems, could be well suited for a variety of design and production scenarios in several fields of application.

Limitations
Currently our system shows various limitations that need to be addressed before moving to further stages of development. In particular, in order to address the Architecture and Construction sector, the issue of scale poses a set of unresolved challenges.
The production of larger and heavier structures will certainly limit the physical possibilities of gaining height by exploiting VTI. In order to successfully achieve the goal without sacrificing the integrity of the volumetric structure due to deformation and slumping, core implementations would have to be addressed, such as a controlled process of drying/curing during fabrication runtime. Similar technical implementations would have to be combined with a more focused assessment on the actual design space afforded by such material system. Assuming the availability of a larger robotic arm, the computational core of the pipeline would theoretically be neutral to an upscaling process, whereas the feeding and extrusion system would certainly require upgraded equipment to implement larger material volumes and extruder torque values. The fabrication tolerances afforded by the current sensing setup (oscillating between 1.5 and 2 mm) would intuitively react positively to a larger clay thread cross section and to wider surface areas. Nevertheless, the possibility of upgrading to a higher-class sensing technology would certainly increase production reliability and overall efficiency.
A characterization of the mechanical performance of the presented FGS has not been covered at this stage of development. Currently, the structural integrity (defined by the mechanical bonding between consecutive deposition iterations) of the graded structures produced with our method has still to be understood. Furthermore, mechanical testing would validate the hypothesis of a graded mechanical behavior. Also related to a thorough material characterization, is the still open question of developing a reliable process of assessment of the inner volumetric porous structure and its related properties.
On the computational side, the main bottleneck of our system is represented by the inefficiencies derived from interfacing the sensing device with the Rhino/Grasshopper environment. In particular, the Python server requires a large processing time for capturing a single depth frame at full resolution (approximately 10 s), whereas the rest of the design and fabrication pipeline requires one order of magnitude less time to compute (approximately 1 s). Another limitation is represented by the lack of sophistication of the current robotic toolpath-planner. The absence of advanced collision avoidance or path-planning strategies in a similar fabrication process, combined with the current computational design strategy, confines the existing production of FGS to simple box-like structures. To tackle similar issues, we are considering the integration of our computational system to other software platforms such as ROS. ROS would provide a framework which is specifically designed for developing experimental robotic processes, and that natively supports state-of-the-art pathplanning libraries and integrated communication protocols pipelines for sensors and actuators integration.

Future work
Within a shorter timeframe, we prioritize the attempt in closing the gap on the mentioned performance-related limitations such as the porosity assessment and the mechanical characterization of such category of FGS. Given the immediate availability of a larger-size Comau industrial robot (1200 mm horizontal reach and 10 kg wrist payload), we will also aim at the fabrication of a significantly larger clay graded structure. Considering possible future implementations on the material side, testing other viscous material systems such as concrete or different viscoelastic polymers could open fruitful applied research scenarios. Considering a longer research term, we believe that focusing on the integration of Reinforcement Learning strategies in a context of material instability could offer valuable insight on a category of automated fabrication methods based on dynamic adaptation rather than accurate prediction.