Estimating Natural Frequencies of Cartesian 3D Printer Based on Kinematic Scheme

: Nowadays, 3D printers based on Cartesian kinematics are becoming extremely popular due to their reliability and inexpensiveness. In the early stages of the 3D printer design, once it is chosen to use the Cartesian kinematics, it is always necessary to select relative positions of axes and linear drives (prismatic joints), which would be optimal for the particular speciﬁcation. Within the class of Cartesian mechanics, many designs are possible. Using the Euler–Lagrange formalism, this paper introduces a method for estimating the natural frequencies of Cartesian 3D printers based on the kinematic scheme. Comparison with the ﬁnite element method and experimental validation of the proposed method are given. The method can help to develop preliminary designs of Cartesian 3D printers and is especially useful for emerging 3D-printing technologies.


Introduction
Robotic systems are becoming more and more popular in various applications since they perform fast and accurate operations while decreasing production costs and reducing tedious and potentially hazardous tasks. One rapidly developing field of application for robots is additive manufacturing, often referred to as 3D printing [1,2]. Having emerged as a tool for rapid prototyping, nowadays 3D printing is also widely used for various industrial applications [3,4].
Industrial robots are highly diverse in many parameters, and several classifications of them can be proposed [5]. Usually, researchers identify six significant kinematics of industrial robots [6]: articulated, SCARA, Cartesian, Parallel (or Delta), Cylindrical, Spherical.
Cartesian mechanical design became a predominant configuration in 3D printing [7,8]. The main advantage of Cartesian robots for additive manufacturing is the lowest cost compared to the other types of robots with the same accuracy and repeatability. The second popular design for 3D printers is the Delta robot, almost as inexpensive as Cartesian but more demanding on production quality and calibration accuracy. This issue sufficiently limits its popularity [9]; nevertheless, special implementation of this type of robot may yield unique advantages, e.g., energy efficiency [10]. Articulated robots gained the most popularity in such applications like welding, assembly, handling, and inspecting [11] and are used only in specific 3D printing applications [12]. SCARA robots, as well as cylindrical robots, are relatively rare in 3D printing but have some prospects in this field [13,14].
As can be expected for the rapidly developing area, the level of scientific comprehension of additive manufacturing sufficiently lags behind the practice. For example, a study of publications on the Design for additive manufacturing (DfAM) found that the main content of these works is guidelines, rules and certain practical aspects, papers are mainly concentrated in a small number of clusters and the level of international collaboration in them is weak [15]. There is a similar situation with publications on the theory of 3D printer design. This paper tries to introduce one more important theoretical aspect of additive manufacturing into scientific discourse.
The printing rate of the 3D printer is one of the most critical consumer characteristics, directly affecting its user experience and economic returns. A recent study by J. Go et al. proposed that three main factors limit the building rate in FDM/FFF technology: the rate of filament feed, the rate of filament heating, and velocity of the printhead positioning [16]. The first and second factors have been addressed in several recent studies. For example, a novel rotatable extruder with a slot-shaped nozzle was designed to increase the filament feed rate [17], and the heating rate was sufficiently increased by introducing a laser heating system [18]. In the article [19], the authors propose to use fuzzy adaptive control for balanced thermal distribution during extrusion, in order to achieve the highest quality rapid prototyping. However, increasing the speed of printhead positioning is probably the most challenging problem in mass-market printers. Despite the technical opportunity of using rigid and heavy CNC machine frames providing better printing quality at high speeds, the design of 3D printers always implies a trade-off between mechanical stiffness and the printer cost, which should be as minimal as possible for economic reasons. This is especially crucial in FDM/FFF 3D printers since this technology is the most inexpensive and the most demanded in the 3D printer mass market. Therefore, many various designs are present on the market. Unfortunately, almost each of them suffers from certain flaws limiting the building rate sufficiently below the theoretical maximum determined by the properties of the printing material.
The vibration characteristics of 3D printers and components directly affect the accuracy of the work, as reported in several studies, e.g., [20][21][22]. The main problem in FDM/FFF 3D printing originating from mechanical imperfectness is ghosting (ringing): the emergence of a repetitive pattern on the printed detail surface caused by mechanical vibration of the printer frame. This matter limits the maximal acceleration and speed available to a particular 3D printer [23,24]. Another defect also caused by mechanical imperfectness is corner swell, an outgrowth of the deposed material on a corner [24,25]. An example of these defects is given in Figure 1. They are negligible within a specific range of positioning settings, but printing errors become significant and practically unacceptable when they exceed particular values. The magnitude of the resulting vibrations directly depends on the rigidity of the structure, which is also affected by many other factors [26]. Nevertheless, reducing the vibration amplitude of the system and increasing the natural frequency of vibrations can improve the quality of 3D printing.
Several recent studies are dedicated to detecting and decreasing these effects [24,27,28]. The most effortless way to do that is to increase the natural frequency of the printer head vibration due to proper mechanical design, leading to decreasing mode shape displacements [29]. A more elaborate solution is using active vibration control, which is a common practice in CNC milling machines [30,31]. For 3D printing, some solutions based on feedforward control have been reported recently applied to vibration-prone 3D printers [32,33].
In addition, open-source Klipper firmware allows using one of the feed-forward resonance compensation algorithms for all compatible 3D printers [34].
However, how does an engineer understand which design is the most vibrationtolerant at the stage of preliminary design, where neither a complete CAD/CAE model exists nor any experimental specimens are available? This paper proposes using natural frequencies as the feature for distinguishing vibration-tolerant designs from vibrationprone ones. This feature depends on the relative positions and mobility of 3D printer parts. We propose a method for natural frequencies estimation from the 3D printer kinematic scheme, study in detail its bounds of applicability and give several illustrative examples including decision making on the 3D printer construction.
The paper is organized as follows. Section 2 describes the approach to calculating the natural frequency for a given Cartesian design. Section 3 presents examples of natural frequency analysis for Anycubic i3 Mega (Shenzhen Anycubic Technology Co., Ltd., Shenzhen, China) The natural frequencies have been obtained theoretically using the proposed approach and experimentally, confirming its applicability. Section 4 gives brief conclusions.

Vibration Analysis of 3D Printers
This section introduces kinematic schemes for 3D printers. We show that a simple approach based on a flexible joint model is applicable in the case of composite structures with considerably stiff beams, such as constructions of the aluminum extruded profile, often used in custom 3D printing buildings.

Problems of Modal Analysis for Constructions with Bolted Joints
Knowing its dimensions and material, the vibration characteristics of any monolithic component can be calculated with one of the well-established approaches. In the case of beams, the basic tool for their normal mode estimation is the elastic beam theory, including the Euler-Bernoulli, Rayleigh, Timoshenko, and other models [35]. However, stiffness calculation for the composite structures with this approach is a nontrivial problem [36]. Usually, estimation of vibration frequency distinguishes several types of boundary constraints for beam couplings: hinged, fixed, sliding, and free. They imply an absolute coupling stiffness. In practice, components are usually coupled with different bracing (bolts, angles, screws, etc.). In such systems, joints, backlashes, surface friction, and nonlinear deformations play a huge role, making the problem of normal mode estimation extremely complex. Only very recent studies deal with efficient calculations of natural frequencies of bolted structures, given their nonlinear character. For example, the work [37] proposes a machine learning approach to nonlinear modal analysis. The paper [38] describes an artificial neural network for predicting parameters of normal contact stiffness, penetration limit, and contact. Recent work [39] proposed an efficient algorithm for estimating damping in bolted joints providing good accuracy compared to simulations in commercial finite element software. Nevertheless, this solution is still not implemented in commercially available engineering packages.
In several papers, the determination of natural frequencies for bolted structures using FEM analysis shows unsatisfactory results, with errors ranging from 9 to 46.9% in determining natural frequencies [40]. A technique of using equivalent material can be applied to improve the accuracy [41]. Equivalent lower-dimensional models can also be efficient for simulating bolted constructions [42]. In civil engineering, to avoid detailed analysis of each element joint, a concept of effective stiffness is introduced reflecting the overall stiffness of a joint [43].
In our previous study, we did not find low-frequency vibrations for the 3D printer Anycubic i3 Mega using finite-elements analysis [44] since we did not simulate joints, which played a leading role in the occurrence of this low-frequency mode.
Therefore, we will consider an approach based on the equivalent model of the entire construction. Since we consider the preliminary design stage, we will omit any implementation details and focus on the effective stiffness of the joint. We will use the term "stiffness" for simplicity, meaning exactly "effective stiffness". Nevertheless, the proposed approach is valid only if the beams can be treated as absolutely stiff. Therefore, we must ensure that at least their first normal modes are far above the normal modes of the entire setup, i.e., the condition must be satisfied where w b is the beam first normal mode, w is the setup lower frequency. Elastic beam theory provides a well-established tool for verifying this assumption via estimating normal beam modes. The following section briefly recalls the main concepts of the Euler-Bernoulli beam theory and its application to our problem.

Euler-Bernoulli First Normal Modes
The natural frequency of the beam using the Euler Bernoulli theory is calculated by the formula [45]: where E is elastic modulus, I is the second moment of the beam's cross-section area, ρ is beam material density, A is the cross-sectional area, and a is a wavenumber of Euler-Bernoulli mode. Wavenumber a for a beam fixed at one end and corresponding to the first mode equals 1.875 [45].
In our previous study, we found the experimental natural frequency of the aluminum profile beam to test the theoretical results (see Figure 2). The beam was attached to a fixed base using hidden aluminum profile corner brackets.
We investigated two options for fastening the beam to a fixed base: using one or two hidden brackets. The first mode frequency for the case with one hidden bracket was 22.1 Hz, and for the case with two hidden brackets, it was 45.2 Hz. As a result, the beam's natural frequency value calculated with the Euler-Bernoulli approach was significantly higher than the values obtained experimentally (the theoretical value is 135.86 Hz).
It means that almost all energy received by the beam from the impact hammer was turned into vibration caused by the joint flexibility rather than the beam elasticity. Furthermore, various beam mounting options affect the obtained values of natural frequencies, which verifies this assumption.
The flexible mounting option is not considered in the Euler-Bernoulli theory; however, it can be extended to such a case. Meanwhile, theoretical calculation confirms that the condition (1) is satisfied. Experimental results show that vibration on Euler-Bernoulli's first normal mode of the beam is almost indistinguishable compared to the vibration caused by a finite joint stiffness k. So we can conclude that the condition (1) allows considering the aluminum profile beams used in the 3D printer construction as a rigid body.

How Cartesian 3D Printer Kinematic Scheme Affects Its Dynamics
Cartesian 3D printers consist of two independent parts: the frame with the printhead mounted on it and the printing surface (the printbed). These parts move relative to each other. The system requires the mobility of 3D printer parts along three orthogonal axes to be functional. Using the standard graphical notation of kinematic schemes [47], we can depict a kinematic structure of any Cartesian design as a very simple object with three translational degrees of freedom provided by three orthogonal groups of prismatic joints (redundant systems will not be considered). Examples of such schemes are given in Figure 3, left column. Figure 3a refers to the simplest structure used in very cheap 3D printers, Figure 3b refers to popular gantry kinematics with one closed kinematic chain, Figure 3c presents a variant of more stiff kinematics, containing three closed kinematic chains. These kinematic schemes figure out only how the structure is intended to move.
It is not obvious from each kinematic scheme in the left column, which one is better in terms of vibration tolerance unless we consider an additional factor: a finite stiffness of joints between beams, linear guides and other elements of construction. As the simplest approximation, these joints can be considered spherical joints with internal stiffeners. Such a point of view is reflected in Figure 3d-f. From these extended schemes, the advantage of closed kinematic chains becomes obvious: the more spherical joints are involved in possible angular movement the higher the overall stiffness is.
For better illustration, we introduce a simple formalism describing differences in relative axis mobility and axes structure for 3D printers and present all possible kinematic designs within this formalism in Appendix A. Examples of kinematic schemes for Cartesian 3D printers. Rectangles denote prismatic joints, circles denote spherical joints, a black triangle denotes the printhead, and a black parallelepiped denotes the printbed. Panels (a-c) refer to functional kinematic schemes, while the panels (d-f) additionally show frame and guide joints that cannot be considered absolutely stiff. Now, let us consider how the kinematics of the 3D printer affects its vibration characteristics and therefore printing quality. Denote the vector of angular rotations in spherical joints as θ = (θ 1 , θ 2 , . . . , θ n ) , and a vector of linear translations in prismatic joints in local coordinate systems associated with each direction of prismatic joints as q = (x L , y L , z L ) , where N is the number of spherical joints and n ≤ N is the number of rotational degrees of freedom, and the number of translational degrees of freedom is 3. Denote the global coordinates of the 3D printer working point as x = (x, y, z) . The overall dynamics of the 3D printer correspond to a linearized ODE: where J(x) is a diagonal matrix of moments of inertia, M is a diagonal matrix of masses, C θ , C x are diagonal matrices of damping coefficients, K θ , K x are matrices of linearized stiffness, f θ is a vector of disturbance torques, q is a vector of relative translations in local coordinate systems attached to the prismatic joints, f x is a vector of disturbance forces, x is an absolute position of the 3D printer working point, e is the positioning error. The kinematic scheme determines the number of degrees of freedom and the function J(x).
Let us show more exactly how dynamical error is related to mechanical movement. Consider a case where stiffness matrices are diagonal, so each equation of the system (3) is independent. Then, one rotational degree of freedom of a linearized system is described as follows. Let the system be affected by the torque f θ = Ae jΩt , where Ω is the frequency of vibrations induced by the 3D printer parts' translational movements. In case of small angles the function L is linear, so we imply e = Lθ, K θ = k, C θ = 0. The Equation (3) then reads: Substituting a known solution θ(t) = θ 1 e jωt + θ 2 e jΩt into the Equation (4), where ω = √ k/J is a natural frequency of the printer vibration, we obtain: The error amplitude maximal value is e * = L(|θ 1 | + |θ 2 |). The low-frequency vibrations caused by massive 3D printer parts movements satisfy the condition ω 2 > Ω 2 . The Equation (5) explains the profit from the natural frequency enhancement: increasing the natural frequency lowers down the amplitude of the dynamical error.
The way how real Cartesian systems are implemented suggests using more complicated structures than the ones shown in Figure 3. A typical Cartesian system can be illustrated by the gantry-type kinematics often used in architectural 3D printers or CNC milling machines. Its functional kinematic scheme is given in Figure 4a following literature on CNC machine design [48,49]. Meanwhile, a more detailed kinematic scheme given in Figure 4b shows some issues with practical implementations. In addition to actuated prismatic joints, each axis is supplied with passive linear translation mechanisms for several reasons: to decrease the load on actuated joints, to prevent skew and jamming in actuated joints, and to improve the overall stiffness. The latter is important in the context of the current paper. A study for the CNC machine [49] confirms the importance of passive joint stiffness. gives a detailed kinematic scheme of its implementation with actuated prismatic joints (shown with white rectangles) and passive rail guides (shown with grey rectangles). Circles denote passive spherical joints with internal stiffeners, a black triangle denotes the printhead, and a black parallelepiped denotes the printbed The engineering practice of 3D printer designing implies using multiple additional passive joints as well. Several types of passive linear guides are common including round rods with ball bearing carriages, profiled rails with ball bearing carriages and aluminum profiles with wheeled carriages. Their properties such as backlashes, precision and price are highly variable and should be considered individually for each manufacturer and each design, and their application is also highly variable. For example, in inexpensive gantry 3D printers a brintbed with round rod guides is used in (Prusa Research, Prague, Czech Republic) and Anycubic i3 Mega (Shenzhen Anycubic Technology Co., Ltd., Shenzhen, China) 3D printers [50,51], while Sovol SVO3 (Sovol 3D, Shenzhen, China) has a wheeled carriage on an aluminum profile [52]. None of these designs should be preferred without an additional investigation.

Flexible Joint Model
The approach described in the previous section based on the flexible joint model is applicable to any construction made of rigid metal beams, in which natural frequencies are much higher than the natural frequencies emerging due to joint imperfectness. Unless other is claimed, we assume the system is strictly linear and symmetric and stiffness matrices in (3) are diagonal, so vibrations, which would occur in at least two orthogonal planes (denote them θ x (t) and θ y (t)) are independent and similar in case of similar initial conditions. Therefore, we first consider only a planar motion in one selected plane (denote the angle θ(t)) and only one corresponding frequency, and then extend this approach to more complex cases.
Vibration model of the beam. A joint model based on a linear flexible joint is proposed to find the natural frequency of vibrations, a. When the length of the beam is relatively small, it can be assumed that the structural profile is rigid. Then, according to Newton's second law, the motion of the beam is expressed in the following linear differential equation: where J = 1 3 M T L 2 is the beam moment of inertia, M T is the total mass of the beam, L is the beam length, and k is the torsional stiffness of the joint, θ is the angular displacement (see Figure 5). The analytical solution of Equation (6) has the form: Substituting (7) into (6), we obtain: hence, the vibration frequency of the extruded profile beam is expressed as: Equation (6) can be used to determine the stiffness of the joint by empirically determining the first mode frequency of the system and expressing the stiffness as a ratio: In the case of linear stiffness, (9) provides a simple way to determine k from experimentally estimated ω. In the case of nonlinear stiffness, the identification procedure should involve the waveform shape and would be more complicated [46].
Vibration model of the gantry structure. To find the natural vibration frequency of the gantry structure shown in Figure 6. Let us determine two independent mechanical stiffness of the joints. The stiffness of the lower joint is denoted by k 1 , and the stiffness of the upper joint is denoted by k 2 . The mass of the vertical beam is indicated M 1 , and the mass of the horizontal beam is represented by M 2 . The moment of inertia of the beam is equal to L 1 2 : The kinetic energy of motion of the gantry structure in this case is equal to: the potential energy: The Lagrange function, where L is the Langrangian, is equal to the difference of kinetic and potential energy: the Euler-Lagrange equation is written as: Using (7), we find the equation of the frequency of natural vibrations of the gantry structure: It should be noted that the length of the horizontal beam does not affect the frequencies of the gantry structure.
Vibration model of a parallelepiped. For simplicity, in the structure shown in Figure 7, the stiffness of the joints can be assumed and equal to k 1 and k 2 . The length of all vertical beams is equal to L 1 . The moments of inertia of the vertical and horizontal beams are equal, respectively, to: The kinetic energy of motion of the parallelepiped structure is equal to: potential energy: The Lagrange equation after transformations will be: The frequency of the parallelepipedal structure, in this case, is equal to: Taking into account several degrees of freedom. Usually, 3D printers have two, three or more degrees of freedom, including oscillations in directions of X and Y axes, oscillations of the printbed independently from the frame, moreover, some printer designs provide horizontal and vertical elements of the frame to oscillate independently as well, like the open kinematic chain design shown in Figure A3b. If we use linearized models and also assume that the stiffness matrices in the Equation (3) are diagonal, treatment of multiple degrees of freedom is fairly simple: each rotational degree of freedom is considered as if the other degrees of freedom are fixed. For example, let us calculate the gantry structure vibration in another plane, orthogonal to the plane considered before. The kinetic energy of motion of the gantry structure in this case is exactly the same but the potential energy, due to the fact that upper joints are not involved in the oscillations, is This results in the frequency: To simplify this approach, we present an Algorithm 1 for calculating Cartesian structures. For example, for the gantry and parallelepiped structures considered before longitudinal beams are vertical, transverse beams are horizontal, so Formulas (10)-(12) could be easily derived without Lagrangians.

Algorithm 1:
The algorithm for calculating Cartesian structure natural frequency.
Select the plane in which the vibration of the structure will be investigated. Select an appropriate degree of freedom.

2.
Decompose the construction into primitives:

3.
Calculate the stiffeners k j , j = 1 . . . N (N is the number of joints).

4.
Estimate the frequency using the formula:

Analysis of Sample Designs
The current section describes the application of the proposed method for estimating natural frequencies of 3D printers to some practically valuable examples.

Experimental Analysis of a Gantry 3D Printer
It is possible to find the first lower mode for an existing 3D printer using the proposed approach. A gantry design Anycubic i3 Mega was chosen as an example. A distinctive feature of this 3D printer is the presence of a relatively heavy printhead mounted on two horizontal cylindrical guides working as a whole due to the rigid coupling of their ends (see Figure 8). The length of the vertical beams is labeled L 1 , and the height at which the printhead is located is labeled L 2 . This design has three types of joints with different stiffness.
The kinetic energy of motion of the gantry 3D printer is equal to: potential energy of motion: The Lagrange equation after transformations will look like this: The natural frequency of oscillation of the 3D gantry printer can be expressed using the formula: This formula can be also obtained using the Algorithm 1. Using the Formula (13), we can find the relationship between the printhead position on the Z-axis and the natural frequency of oscillations. A testbench for vibration data acquisition was assembled to measure the vibrations of the head, bed, and frame of the 3D printer. The testbench includes the following components: NI PXI 1042 controller with the PXI-4461 module for collecting vibroacoustic signals, piezoelectric accelerometer IMV VP-4200, accelerometer power amplifier, and custom mount for accelerometer printed on a 3D printer, which allows placing the sensor parallel to any of three axes X, Y, Z (see Figure 9).
A series of experimental data were recorded with a sampling rate of 1000 Hz. The program for recording data into a file was implemented in the LabVIEW environment. The data received from the accelerometer was processed in the MATLAB environment using the built-in FFT function.
We investigated the dependence of the natural frequency of oscillation on the height of the printhead. Since the printhead and related components involved in Z-axis movement (guides, belts, mounts) have a total mass of 1.7 kg, which is greater than the mass of the vertical beams (0.7 kg), the location of the printhead on the Z-axis affects the 3D printer natural frequency, much. In the experiment, we used the g-code to set a short-term extruder movement at each of the positions along the Z axis with a step of 10 mm, and the maximum height at which the printhead was positioned was 200 mm. Printer movement gave an impulse to the printer frame impact with the spectrum as close to that experienced by a 3D printer during printing as possible. The accelerometer was mounted on the X axis of the 3D printer in the right corner, see Figure 9. This position allows determining the frequency of frame vibrations affecting the printhead directly during operation.
The total stiffness of the 3D printer at zero position was found using Equation (13) to validate the theoretical approach based on the flexible joint model. Further, the theoretical natural frequencies of vibrations were found. Figure 10 shows the dependence of the experimental and theoretical natural frequencies on the position of the printhead along the Z axis.  Figure 10 shows that the proposed approach successfully predicts the first mode of the Anycubic i3 Mega. The relative error of the theoretical frequency was 2-6 percent (see Table 1).

Preliminary Analysis of a New Design
This subsection shows how natural frequency can be calculated from the classification scheme for an arbitrary 3D printer.
Suppose, natural frequencies of cubic printers depicted in Figure A2a,b are predicted. From Algorithm 1 it follows that the natural frequency equals to where k i is i-th joint stiffness, and J i is i-th moment of inertia. The formula for the frequency of open Z-axis design is: where M G is the guide mass and M P is the profile mass. Let the mass of a horizontal beam equal M h = ρ P L h , where L h is the length of the horizontal beam, and M v = ρ P L 1 is the mass of a vertical beam. From this and Equation (15) we can obtain: Substituting values from Table 2 yields: The frequency of the kinematics with closed Z-axes is calculated using the following formula (supposing that one motor is used for actuating both linear drives of the closed Z-axis on a printbed): Substituting values from Table 2, another value is obtained: The resulting difference is relatively small, so we can conclude that using the open Z-axis is reasonable for the design from the point of view of natural frequency in XY plane, and a closed Z will not yield much more stiffness.

Comparison of Common Designs
To compare the seven most common 3D printer designs, their CAD/CAE models have been developed. A more detailed description of these designs is given in Appendix A.
First, we used the finite-element method (FEM) to calculate the natural frequencies of the printers assuming that all joints are absolutely stiff. Then, we used the proposed method based on the kinematic schemes of the same structures and calculated the approximate lowest frequencies of these printers using Algorithm 1 with the parameters given in Table 2. Figure 11 shows small images of the printers in one panel with their ordinal numbers, and Figure 12 shows the corresponding lowest natural frequencies found with two methods: the FEM and the proposed one.  The experiment found that for the cuboid structures 1 and 6 the results of our method and FEM approximately coincide (see Table 3). This shows that profile elasticity has an almost similar impact on the vibration as the flexibility of the joints. For the other designs, the FEM overestimates the vibration frequency in comparison with the proposed method. This can be explained by the fact that these constructions are more light-weighted, and in the case of absolutely stiff joints the overall stiffness is high. In practice, flexibility in joints make such structures more vibration-prone, and this is confirmed by engineering practice: only a few cheapest 3D printers have an open kinematic chain structure 3, and structure 2 is also becoming less common in recent years. Of course, these results should be extrapolated on real designs with care. However, they clearly show the limitation of the FEM when joint flexibility is not taken into account. This is especially important when novel additive manufacturing machines are developed with heavy printheads, such as painting robots [53].

Example of Cost Calculation
As an example of calculating the cost of an FDM 3D printer, a gantry printer was chosen, shown in Figure A3a, with a frame made of 20 × 20 mm aluminum profiles. The length of horizontal and vertical aluminum profiles is 400 mm. The length of lead screws was chosen as 350 mm. Table 4 shows the cost (minimum and maximum) of the main components of a 3D printer. The difference in the cost of the components is conditioned by their different quality, brands and other features. Among the components, there are obligatory ones used in any configuration of a 3D printer and not depending on the kinematic scheme: the power supply unit (PSU), the printbed, the printhead and the controller board. Table 4. The cost of a 3D printer in Figure A3a. The economic effect of creating a 3D printer with the specific kinematic scheme increases when cheap obligatory components are used, because in this case a rail guide system, motors, etc. will affect the total cost of the printer much more significantly. In the case of using top-level obligatory components, a more stiff kinematic scheme will not result in relatively much higher expenses, so it should be preferred. The exception here is the rail guide system, which is a rather costly component, but it can be replaced with cheaper cylindrical guides in many cases.

Conclusions
The paper proposes a simple method for calculating 3D printer natural frequencies based on the kinematic scheme. It is applicable in the case of high stiffness of elongated elements such as profiles, guides, and beams. This method allows replacing elaborate finite-element analysis of the detailed model on a preliminary design stage. The accuracy of this approach is relatively high and feasible for engineering purposes; the provided examples confirm its practical applicability. Several most common designs have been compared via their lowest natural frequency, and the variant with the greater value of vibration frequency is preferred using this technique. Additionally, we give an example of calculating the FDM printer cost. We conclude that using expensive obligatory parts such as the power unit, the printhead and the controller, may reduce the relative cost of mechanical components and therefore make using a more stiff kinematic scheme more reasonable. An example of decision-making concerning the type of 3D printer is given.
We also propose a classification scheme for the Cartesian mechanical systems used for additive manufacturing which considers variants of already existing structures and examples not yet applied to 3D printing. It can be used by developers of various additive machines, such as desktop 3D printers, architectural 3D printers, and others, to simplify the engineering process in a preliminary design stage. The algorithm for natural frequency estimation can be used for computer search of proper designs using optimization. A table of stiffness of various joint components obtained from the experiment can be developed which would help engineers to find optimal design solutions.
To summarize, the paper makes the following contributions. First, the paper provides theoretical and practical confirmation of the applicability of the flexible joint model and the necessity to take bolted joint flexibility into account during finite element analysis, which is a non-trivial task in most moderns software packages such as Fusion 360. Second, the paper gives a simplified and handy algorithm for natural frequency estimation of the 3D printer which does not require the application of the Euler-Lagrange formula directly and is suitable for both manual and machine usage. Third, a classification scheme for Cartesian 3D printers is proposed using two criteria: the type of kinematic chain associated with each axis (open/closed) and the mobility of each axis' elements (actuated/fixed) connected to the print head and the print bed.
Further work will be dedicated to testing some non-existing designs and considering other decision-making criteria, such as the system tolerance to the inaccuracy of the 3D printer components, in addition to the natural frequency criterion.

Acknowledgments:
The authors are thankful to anonymous Referees for their insightful reviews and useful recommendations.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Classification of Cartesian Designs: Principles and Examples
The usefulness of a specific classification scheme for practice depends on whether it is accurate enough to systematize all the existing knowledge and is capable of accommodating new subjects [54]. Basics of kinematic scheme topology classification and methods for mobility investigation are reported in papers [55][56][57]. For 3D printing, a detailed classification of possible designs may be helpful for engineers at the early design stage to choose the most feasible kinematics, for sellers to organize their catalogs, for purchasers to identify their requirements more efficiently, for researchers to place their novel designs among existing ones.
Cartesian 3D printers consist of two independent parts: the frame with the printhead mounted on it and the printbed.To distinguish between the coordinate axes of the printhead and the printbed, we denote them as X H , Y H , Z H for the printhead and X B , Y B , Z B for the printbed. These parts are separated in space by the printed part and move relative to each other. The system requires the mobility of 3D printer parts along three orthogonal axes to be functional. For example, if the coordinate system as a whole is stationarily attached to the printbed, then the printhead must move along the X H , Y H and Z H axes. Hereafter, for brevity, we will use the term "axis" not only for the direction in space but also to denote a set of mechanical elements (beams, guides) designed to move the printer's working body along a given coordinate axis.
Kinematics of modern Cartesian printers implies variation of mobility and immobility of the axes. In our analysis, the actuated axis is designated (A) and the fixed axis is designated (F). Table A1 shows all possible combinations of the printhead and the printbed axis mobility. Each row in the table represents a kinematics variant with three actuated axes. Because creating each layer of the printed part takes place in the XY plane, the Ox and Oy axes of the printhead and the printbed are interchangeable. Thus, the third row is synonymous with the fourth row, the fifth row is interchangeable with the seventh row, and only one version of each pair will be considered further.
Also, two mounting options for the actuated axes were considered: • Open axes (A O )-actuated axes that have one attachment point to a fixed base or elements of other axes. This corresponds to the open kinematic chain. • Closed axes (A C )-actuated axes that have at least two attachment points to a fixed base or elements of other axes. This corresponds to a closed kinematic chain.
The closed axis provides more rigidity to the construction but needs a more complicated design with additional guides, motors, etc. Table A1. Variability in axis mobility of the Cartesian 3D printer. Cells for the printhead are filled gray. Rows 4 and 7 are not investigated further.

. Classification of Cartesian Designs
In this subsection, we will overview all possible combinations of fixed and actuated axes and make a complete classification of Cartesian 3D printers. An example of the visual representation of each kinematics is presented in Appendix A in Figures A5-A10. It should be noted that the kinematic schemes can be depicted differently depending on where the entire structure is mounted (e.g., standing on the table or attached to the wall). All illustrated constructions shown in Appendix A are fixed on the floor. Tables A2-A7 show kinematic schemes with different variations. Since the X and Y axes are interchangeable, designs numbered 2 and 8 (marked yellow (*)) as well as 6 and 7 (marked red (**)) in some tables are the same, and one solution from each pair can be omitted. In Tables A4 and A5, there are no interchangeable designs because the axes X B and Y H (or X H and Y B ) belong to different 3D printer components.

. Survey of Existing Mechanical Designs
A review of existing 3D printer solutions showed that many of the variants presented in Table A2 of the kinematics are not used in practice, except for numbers three and four. Examples of a additive systems with closed X H , Y H and open Z H axes (see Figure A1a) are the machine painting systems [53,58] [62] and a custom printer with Core-XYZ kinamtics [63].
Creating 3D printers with bed fixation is sometimes not a convenient technical solution but rather a forced one. The choice of this configuration for architectural printers can be explained by the fact that we can not move the platform on which the building stands.   Figure A2a) is Anycubic 4max Pro [51] or Ultimaker 3. An example of a 3D printer with closed X H , Y H , Z B axes (see Figure A2b) is the Total Z Anyform 3D printer [64].These 3D printers differ only in the closed Z axis of the table, which affects the stiffness of the table. As a rule, 3D printers with a closed Z B axis are often significantly more expensive, which is due to using more precise ball screw drives instead of a trapezoidal screw and is not fully determined by the printer kinematics.  Table A4 shows two variants used in practice. The first variant is the one with closed axes X H , Y B and open Z H . This variant is trendy in the 3D-printers market due to its cheapness and simplicity of construction. Examples of a this type of 3D printers are Wanhao Duplicator i3 [65] and Anycubic i3 Mega [51] (see Figure A3a). The second variant with open X H , Z H axes and open Y B (see Figure A3b) is a simplified version of the previous model, and one of the cheapest variations of Cartesian 3D printers. An example of such a 3D printer is the Wanhao Duplicator i3 Mini 3D printer [66]. Examining the designs from Table A5, we found only one variant of the existing 3D printer with two actuated printbed axes. It is Felix 3.0 [67], a 3D printer with open Y B , Z B and closed X H (see Figure A4). The manufacturer claims that the 3D printer can develop a fairly high printing speed thanks to this design compared to its counterparts. Three-dimensional printers with the designs shown in Tables A6 and A7 were not found. Probably, such designs are not quite suitable for 3D printing. Moving the printing surface along all three axes can affect print quality since the printbed with the printed detail is much heavier than the printhead. However, it is possible that designs of this type are or will be used for other additive technologies, where complete fixation of the printhead is required due to its size and mass. A recent study published the design of a prototype 3D printer with a fully fixed extruder. However, this printer is not Cartesian and is not commercially available [68]. Now, make a complete classification scheme of the existing Cartesian 3D printers. It is given in Table A8. The classification scheme shows each case's fixed and actuated axes and the closed and open axes for each actuated axis.
Only seven of the 40 possible kinematics designs (8 interchangeable) were identified among existing designs. Table A8 shows that the actuated axes are often closed for greater structural stiffness. Most of the actuated axes relate to the printhead.