Output Force Density Saturation in COMSOL Simulations of Biomimetic Artificial Muscles

Many modern applications, such as undersea drones, exoskeletal suits, all-terrain walker drones, prosthetics, and medical augments, would greatly benefit from artificial muscles. Such may be built through 3D-printed microfluidic devices that mimic biological muscles and actuate electrostatically. Our preliminary results from COMSOL simulations of individual devices and small arrays (2 × 2 × 1) established the basic feasibility of this approach. Herein, we report on the extension of this work to N × N × 10 arrays where Nmax = 13. For each N, parameter sweeps were performed to determine the maximal output force density, which, when plotted vs. N, exhibited saturation behavior for N ≥ 10. This indicates that COMSOL simulations of a 10 × 10 × 10 array of this type are sufficient to predict the behavior of far larger arrays. Also, the saturation force density was ~9 kPa for the 100 μm scale. Both results are very important for the development of 3D-printable artificial muscles and their applications, as they indicate that computationally accessible simulation sizes would provide sufficiently accurate quantitative predictions of the force density output and overall performance of macro-scale arrays of artificial muscle fibers. Hence, simulations of new geometries can be done rapidly and with quantitative results that are directly extendable to full-scale prototypes, thereby accelerating the pace of research and development in the field of actuators.


Introduction
Traditional robotic actuators are electric motors and hydraulic systems. Electric step motors [1] provide precision of motion, utilize convenient power, and are scalable to a degree, making them the common choice in small robotic systems. As the actuation force is the Lorentz force, a significant magnetic field must be applied, usually by running large electric currents through copper solenoids. That generates significant Joule heating and leads to high power requirements. Electromagnetic motors have been found to be most effective for large displacements and less effective in small displacements [2].
The other major approach Is the use of hydraulic systems, e.g., in derricks, cranes, bulldozers, forklifts, industrial robots, and the US Army Mule walking robot [3,4]. However, hydraulic actuators scale disadvantageously with system size, because shrinking the system decreases the piston area, so the same pressure produces less force. Hence, hydraulic systems are unsuitable for miniaturization. In addition, pneumatic actuators output less force than alternative solutions at larger displacements [2]. Hydraulic actuators perform comparably to pneumatic artificial muscles [5]. Finally, biomimetic applications, e.g., exoskeletons and prosthetics, require fluidity and fine motion control, which are difficult to achieve through hydraulics [6]. 2 of 8 In view of the above limitations of traditional actuators, many modern applications, such as undersea drones, exoskeletal suits, all-terrain walker drones, prosthetics, and medical assists and augments, look for alternative solutions, such as artificial muscles [7,8]. For example, Dielectric Elastomer Actuators (DEAs) [9][10][11][12] utilize the deformation of a polymer slab due to the electrostatic force between the charges built-up on electrodes under externally applied voltage. This force scales as the inverse square of the distance between the electrodes. Therefore, it is small for macroscopic distances but increases disproportionally at microscopic scales. (Incidentally, the same principle can be used to build a capacitive high-sensitivity gauge of mechanical deformation [13,14], as well as a pressure sensor [15].) Such devices can be arrayed longitudinally, to increase elongation distance, and laterally, to increase force output. However, the unit cost would make the simple arraying and stacking of individually fabricated devices prohibitively expensive if traditional manufacturing was used. As a result, more recent attempts tend to be limited to small devices, such as single-chamber microrobots [16][17][18] and HASEL tentacle actuators [19,20]. Multi-material 3D printing is promising [21][22][23][24], but still in its infancy.
We pioneered an alternative solution [25]: artificial muscles built by 3D printing microfluidic channels, wherein the channels become embedded wires and electrodes for micro capacitors arranged to form muscle fibers [26]. Applying voltage to an array of micro capacitors makes them contract in unison, producing an electrostatic actuation force. That force is transmitted to the outside world via the surrounding polymer material, functioning as tendons. The resulting electrostatic actuators would be scalable, energy-efficient, and offer a high force-to-weight ratio [25]. The COMSOL simulations [25] indicated up to 33 MPa force density at current advertised limits of fabrication capability. Those simulations were limited to 1 × 1 × 1 and 2 × 2 × 1 devices [25]-devices containing one or two micro capacitors in each of the lateral directions and just one micro capacitor in the longitudinal direction of the artificial muscle.
The output force density scales are the inverse square of the separation distance between the plates in the micro capacitor, so best results are achieved when the individual capacitors are microscopic. Simultaneously, the total force of a muscle fiber bundle is the number of fibers times the generated force in each fiber, as the forces are additive when fibers are arranged in parallel [25,26]. As the force in each fiber is relatively small, a larger number of fibers need to be arrayed. This means that a practical macro-scale muscle would have at least tens of thousands of micro capacitors. Finally, COMSOL is very computationally intensive, particularly when multiple nodes are used simultaneously, and when the number of objects is large. Hence, directly simulating a macro-scale muscle made of micro capacitors is technically unfeasible due to insufficient computational power.
To solve this problem, we came up with the hypothesis that a saturation of output force density can be reached as the size of the array is increased, but before we run out of computational power. If this hypothesis is correct, then that saturation value can be used as a good estimate for the output force of much larger arrays, avoiding the need to simulate them directly. That would circumvent the above problem of insufficient computational power for macro-scale muscles.
To test this idea, we set up and ran simulations of artificial muscle structures at a 100 µm scale, built as N × N × 10 arrays of micro-capacitors, where N was varied from 1 to 13, while tendon thickness within each array was swept. The maximal force density of each sweep was plotted vs. N to produce an overall curve that indicated saturation around and beyond N = 10, thereby proving our starting hypothesis. Hence, a 10 × 10 × 10 array is sufficiently representative and predictive of far larger arrays of such artificial muscle devices. This result is of great practical significance to further research in this field and to building practical artificial muscles of this type. Furthermore, the saturation level was determined to be~9 kPa at the 100 µm scale, which suggests a force density of~1 MPa at the 10 µm scale, indicating the feasibility of strong, energy-efficient, low-density muscles for a wide range of practical applications. Figure 1 shows the basic geometry and parameter definitions used in the simulations. The overall bulk material (outer rectangular prism) and alternating pairs of electrode chambers (inner rectangular prisms) were situated in parallel within the material [26]. The structures were joined together using the union function. Tendon thicknesses e (at device edges) and h (between fibers) were locked to be equal and swept together between 100 and 200 µm, in steps of 10 µm. Top and bottom edge thickness c was set to 100 µm. The non-flexed lateral widths (plateX, plateY) of micro capacitor plates were set equal to 1000 µm. The non-flexed plate thickness (plateZ) and plate distance D were set equal to 100 µm. Silicone (polydimethylsiloxane (PDMS)) was assigned to the bulk dielectric. Liquid water was assigned to the electrode chambers [25].

Materials and Methods
Appl. Sci. 2023, 13, x FOR PEER REVIEW 3 of 9 10 µm scale, indicating the feasibility of strong, energy-efficient, low-density muscles for a wide range of practical applications. Figure 1 shows the basic geometry and parameter definitions used in the simulations. The overall bulk material (outer rectangular prism) and alternating pairs of electrode chambers (inner rectangular prisms) were situated in parallel within the material [26]. The structures were joined together using the union function. Tendon thicknesses e (at device edges) and h (between fibers) were locked to be equal and swept together between 100 and 200 µm, in steps of 10 µm. Top and bottom edge thickness c was set to 100 µm. The non-flexed lateral widths (plateX, plateY) of micro capacitor plates were set equal to 1000 µm. The non-flexed plate thickness (plateZ) and plate distance D were set equal to 100 µm. Silicone (polydimethylsiloxane (PDMS)) was assigned to the bulk dielectric. Liquid water was assigned to the electrode chambers [25]. To set up the simulations, three physics nodes were used: electrostatics, solid mechanics, and moving mesh set to "normal". The solid mechanics node was applied to the bulk dielectric. The moving mesh node was applied to the electrode chambers for efficient computation in finite element analysis applications [27]. The electrostatics node was applied to all domains. Within the electrostatic node, a potential of 3000 V was applied to all surfaces of one set of electrodes and a ground potential was applied to all surfaces of the other set of electrodes. This choice was made for consistency with the previous work [25]. This voltage also ensures that even if the electrostatically induced deformation would shrink the distance between the plates of the micro capacitor, the resulting field would still be well-below the typical dielectric breach value of several hundred volts per micron.

Materials and Methods
Within the solid mechanics node, a fixed-boundary node was applied to one side of the bulk dielectric that was parallel to the electrode plates (the bottom outer surface in Figures 1 and 2), while the opposite side (top outer surface in Figures 1 and 2) was allowed to deform. This method [25] significantly simplified the calculations for COMSOL while still producing valid results. Further, a boundary load from the electrostatic force was applied to all surfaces of the model with the outputs for each component force being equated to their respective Maxwell Surface Tensor Equations [28]. The boundary node allowed the electromagnetic forces to be included in the physical model. A quadrilateral linear mesh of uniform seed size was applied to the model. To set up the simulations, three physics nodes were used: electrostatics, solid mechanics, and moving mesh set to "normal". The solid mechanics node was applied to the bulk dielectric. The moving mesh node was applied to the electrode chambers for efficient computation in finite element analysis applications [27]. The electrostatics node was applied to all domains. Within the electrostatic node, a potential of 3000 V was applied to all surfaces of one set of electrodes and a ground potential was applied to all surfaces of the other set of electrodes. This choice was made for consistency with the previous work [25]. This voltage also ensures that even if the electrostatically induced deformation would shrink the distance between the plates of the micro capacitor, the resulting field would still be well-below the typical dielectric breach value of several hundred volts per micron.
Within the solid mechanics node, a fixed-boundary node was applied to one side of the bulk dielectric that was parallel to the electrode plates (the bottom outer surface in Figures 1 and 2), while the opposite side (top outer surface in Figures 1 and 2) was allowed to deform. This method [25] significantly simplified the calculations for COMSOL while still producing valid results. Further, a boundary load from the electrostatic force was applied to all surfaces of the model with the outputs for each component force being equated to their respective Maxwell Surface Tensor Equations [28]. The boundary node allowed the electromagnetic forces to be included in the physical model. A quadrilateral linear mesh of uniform seed size was applied to the model.
A table of parameters was created to avoid manually building each individual array and instead make proportional relationships between the parameters, so that increasing the array size would expand the overall structure to scale. The parameter N was utilized to define the lateral size of the model during our study in the number of micro capacitors deployed in each of the two lateral dimensions. Therefore, an N × N × 10 device contained N × N parallel columns, wherein each column contained 10 micro capacitors. The model was iterated from N = 1 to N = 13 at a step of 1, while the size of the model in the z direction was kept constant (10 micro capacitors tall). A table of parameters was created to avoid manually building each individual array and instead make proportional relationships between the parameters, so that increasing the array size would expand the overall structure to scale. The parameter N was utilized to define the lateral size of the model during our study in the number of micro capacitors deployed in each of the two lateral dimensions. Therefore, an N × N × 10 device contained N × N parallel columns, wherein each column contained 10 micro capacitors. The model was iterated from N = 1 to N = 13 at a step of 1, while the size of the model in the z direction was kept constant (10 micro capacitors tall).
For each N value, the tendon width between adjacent muscle fibers (h = u) was set equal to the tendon width e at the edges of the overall array (Figure 1), and the resulting parameter e was swept from 100 to 200 µm, in steps of 10 µm. Heat maps (as in Figure 2) of the resulting deformations were generated to visually and qualitatively confirm the correct behavior of the simulation.
For each swept value e, the total output force of the device was calculated as the surface integral of the Von Mises Stress along the unconstrained outer top surface of the device, opposite to the boundary-constrained bottom surface. The total deformed surface area was calculated as the surface integral of the area of the unconstrained outer top surface of the device. Dividing the total output force by the total deformed top surface area produced the average output force density.
For each value of N, the average force density was plotted as a function of the swept tendon thickness e. The result was a response curve for that N. These response curves were organized in a single plot vs. tendon thickness e (Figure 3). The maximum force density of each response curve determined the best performance for that N and the corresponding optimal value of the tendon thickness e. These maximal force densities were then plotted as a function of N, producing the force density saturation curve (Figure 4).
The overall procedure is outlined in a flowchart in Figure 5. For each N value, the tendon width between adjacent muscle fibers (h = u) was set equal to the tendon width e at the edges of the overall array (Figure 1), and the resulting parameter e was swept from 100 to 200 µm, in steps of 10 µm. Heat maps (as in Figure 2) of the resulting deformations were generated to visually and qualitatively confirm the correct behavior of the simulation.
For each swept value e, the total output force of the device was calculated as the surface integral of the Von Mises Stress along the unconstrained outer top surface of the device, opposite to the boundary-constrained bottom surface. The total deformed surface area was calculated as the surface integral of the area of the unconstrained outer top surface of the device. Dividing the total output force by the total deformed top surface area produced the average output force density.
For each value of N, the average force density was plotted as a function of the swept tendon thickness e. The result was a response curve for that N. These response curves were organized in a single plot vs. tendon thickness e (Figure 3). The maximum force density of each response curve determined the best performance for that N and the corresponding optimal value of the tendon thickness e. These maximal force densities were then plotted as a function of N, producing the force density saturation curve (Figure 4).

Results
The simulations were set up and run as described in Section 2, in accordance with the following considerations.
The COMSOL models were designed to satisfy a set of practical requirements. While specific 3D printers may promise resolution around a couple tens of microns with hard The overall procedure is outlined in a flowchart in Figure 5.

Results
The simulations were set up and run as described in Section 2, in accordance with the following considerations.
The COMSOL models were designed to satisfy a set of practical requirements. While specific 3D printers may promise resolution around a couple tens of microns with hard

Results
The simulations were set up and run as described in Section 2, in accordance with the following considerations.
The COMSOL models were designed to satisfy a set of practical requirements. While specific 3D printers may promise resolution around a couple tens of microns with hard resin, the requirement to use soft resin [25] would inevitably degrade the achievable resolution to some extent. Furthermore, consistent quality of printing and surface roughness, as well as the need to clear out the sacrificial material or an uncured monomer, would decrease the practical resolution further. Finally, in preparation for experimental demonstration of the proof of principle, it made sense to conduct simulations at scales that were reasonably easily achievable with current state-of-the-art printers in soft resin.
Accordingly, while the Stratasys Object500 Connex3 3D printer has a professed resolution of~10 µm, the manufacturer recommends at least 100 µm to be used to consistently 3D-print internal channels in soft materials. This recommendation and the above considerations prompted the choice of the smallest dimension in our COMSOL simulations to be 100 µm. This number was used for the thickness of the capacitor plates (plateZ), the distance between the plates (D), the thickness of the top and bottom borders of the device (c), the minimal thickness of the tendon regions between neighboring fibers (h, or u in COMSOL), and the minimal thickness (e) of the tendon regions on the lateral edges of the device. Also, in accordance with empirical evidence from PDMS microfluidics [29][30][31][32][33][34][35][36], an aspect ratio of 10:1 between the lateral widths (plateX, plateY) and the longitudinal thickness of each capacitor plate was chosen to prevent the collapse of the hollow structures that would become the fluidic electrodes.
For every N, the tendon thickness e was swept ( Figure 3) to determine its optimal value that maximized the output force density (total force per unit area). These maxima were then plotted against N to produce a saturation curve (Figure 4). This curve suggests that the force density approaches saturation at and beyond N = 10, and that the saturation level for the output force density is~9 kPa at the 100 µm scale.

Discussion
The simulated numerical results in Figure 4 were validated by comparing them with the results of the previous simulations and first-principle back-of-the envelope calculations [25]. More specifically, the previous 1 × 1 × 1 simulation [25] predicted a 1.44 kPa force density, which is essentially in agreement with the 1 × 10 × 10 result here of~1.75 kPa. The~20% increase can be easily attributed to overall smaller edge effects. Simultaneously, back-of-the envelope calculations adjusted for effective muscle area predicted up to a 2.4 kPa force density [25], which is within a factor of two from both simulations' results.
The observed saturation has several important consequences. It proves the hypothesis that as the size of the array increases, the relative influence of edge effects decreases, so the output force density approaches an asymptotic level at array sizes that are still small enough to be successfully simulated in COMSOL. This is an important conclusion both for these artificial muscle devices and as a general method of simulating macro-actuators made of arrays of individual micro-actuators. It significantly alleviates limitations imposed by computational capacity and allows important determinations to be achieved with relatively modest resources.
For example, the presented N = 13 simulations took 4 days to complete on the Hamming cluster at NPS, while the presented N = 10 simulations took just 1 day. Knowing that N = 10 produces a sufficiently good estimate of the output force density already decreases computational time by at least 4×.
Furthermore, for this type of artificial muscle [25,26], the results show that the output force density at N = 10 is already approximately equal to the saturation value. Hence, simulations of laterally larger arrays should not be necessary to predict their behavior in terms of output force density. This array size can now be used as the optimal standard for simulations of such actuators, while the specifics of the individual micro capacitor can be varied.
Next, the presented results show that the saturation level of the output force density appears to be~9 kPa for the 100 µm scale of artificial muscles. Hence, a 10 cm × 10 cm artificial bicep printed at the 100 µm scale should output~90 Newtons of actuation force.
Furthermore, as the COMSOL works with linear homogeneous materials, the calculations scale with physical size. Hence, the same architectures printed at the 10 µm scale would have to produce 100× the force produced at the 100 µm scale, because the electrostatic force inside the capacitor scales as the inverse square of the plate distance. Consequently, the same artificial bicep from the above example should output~9 kN of actuation force. Such forces would be sufficient for all intended applications, including swimming robots [37].
The ability to quantitatively estimate the saturation force density using the shown method can be used to quantitatively estimate the performance of craft propelled by such artificial muscles. Hence, the presented results are an important step forward towards the practical implementation of this type of artificial muscle.
More generally, the presented method of stepwise increases in matrix size N with nested parameter sweeps to maximize the force density for each N and searching for asymptotic or saturation behavior is a good general method to use for other types of devices, and other geometries as well. So long as the asymptotic behavior can be achieved before the simulations run out of available computational power, the established performance and asymptotic quantitative output can be used to predict the performance and behavior of macro-scale arrays of the same type. Hence, what has been demonstrated is a general method for a much larger range of devices and simulations, as well as a specific tool for artificial muscles, and quantitative predictions for these types of artificial muscles.