Modelling dielectric elastomer circuit networks for soft biomimetics

In order to obtain entirely soft bio-inspired robots, fully soft electronic circuits are needed. Dielectric elastomers (DEs) are electroactive polymers that have demonstrated multifunctionality. The same material can achieve different tasks like actuation, sensing, or energy harvesting. It has been shown that basic logic and memory functions can be realized with similar soft structures by combining multiple DE actuators and DE switches. Thus it would be possible to build, with the same materials and processes, a soft structure that mimics a biological being with all these capabilities. This contribution is focused on the modelling of the aforementioned soft electro-mechanical circuit networks. It is here reported the building process of a comprehensive SIMULINK model including the electro-mechanical behaviour of DE logic units and their interconnections. Conventional models deal with a single aspect of DEs, generating complex finite-element simulations. This contribution, based on a former model for an inverter-based DEO, shows how to integrate these various mathematical models and, with the help of direct measurements, create a software representation of DE circuit networks. This work is intended to demonstrate the validity of a recently introduced model and apply it to more complex circuit networks that have a higher number of components. Since, at the present state, the building processes are by hand, adding components generates more variability due to sample-to-sample variation and human error. Despite this, the model still shows a qualitatively good prediction of the devices’ behaviour. Furthermore, the introduction of new materials and automatic processes will help to reduce this variability, allowing the model to reach even better performance.


Introduction
Dielectric elastomers (DEs) are a class of electroactive polymers with mechanical properties very similar to biological muscles [1,2]. They use electricity, have a fast response time, comparable to the human muscle and can achieve manifold functionalities, such as actuation [1,2], sensing [3][4][5], energy harvesting [6,7], and basic signal processing [8,9], some of them even simultaneously [10,11]. For all these reasons they generated a lot of interest in the field of soft robotics, to be used as artificial muscles and to build totally autonomous soft structures. As an example, in [12] dielectric elastomer switches (DESs) are integrated with DE generators to obtain an entirely soft energy-harvester. In this case, DEs act both as power generators and driving circuit, reducing the need for external components (e.g. conventional diodes).
The field of DE electronics has not reached its full maturity yet, one instrument that it lacks is a comprehensive simulation tool. Due to DE's complex electro-mechanical behaviour and highly nonlinear material properties, modelling usually relies on time-consuming multi-physics finite-element simulations. In [13,14] Wissler and Mazza presented a finite-element analysis of the in-plane actuation of circular dielectric elastomer actuators (DEAs). For bending DE structures like the ones presented in [15,16], it is possible to predict the bending angle. Rizzello et al developed a model for out-of-plane actuating systems based on DEs, then used it to build a feedback controller [17]. In [18] a model for the leakage current through the elastomeric film of a DEA was developed and in [19,20] for the DEA's time-dependent viscoelastic response in actuation. The model proposed in this contribution has to take into account multiple aspects at the same time. The leakage current, viscoelasticity and in-plane actuation related to the DE material have to be considered, then, for circuit networks composed of multiple DE components, it also needs to consider the interconnections between the single elements. DE logic networks, in fact, are assembled from arrays of DESs and DEAs [8,9,21].
A previous work presented a first SIMULINK model of a DEO [22]. The model comprised submodel elements of the device components, DEA and DES, and the component interactions, overall representing the complicated electro-mechanical logic behaviour of the circuit. In this contribution, we show that a similar model can be applied to different and more complex DE structures with qualitatively good results. For a higher number of components, the pathway for the electrical signal propagation is more complex, and this is represented by additional model elements. With the handmade manufacturing process, the higher number of components increases errors due to human inaccuracies in the fabrication. Especially for DEOs, the oscillation's characteristics like frequency and the time to reach stable oscillation depend on some production details that will be discussed in sections 3 and 4, not only on the input values like it is for combinational circuits (e.g. logic gates). This contribution presents an updated model that can produce a qualitatively good representation of the overall electro-mechanical behaviour of a digital DEO [23]. The device has double the number of active electro-mechanical components in respect to the previously demonstrated sinusoidal DEO (the difference between the two will be analysed in section 2).
The biological role models, the central pattern generators (CPGs), are living neural circuits that can produce rhythmic signal outputs that drive muscle arrays to autonomously fulfil tasks [24]. A lot of biological activities that are performed in animals and humans every day are controlled by this kind of structures, without the need for centralized control by a central nervous system. Walking [25], respiration, swimming, flying and peristalsis in the gut are just some examples. Therefore, these types of structures are very important for building entirely autonomous soft robots that do not rely on external control or conventional circuits. In [21,26] has been demonstrated how a DEO can work as CPG for generating the crawling movement of a caterpillar-inspired robot, possessing soft actuators and soft circuits, based on DEs. Another important application of the digital DEO is a clock signal generator inside DE digital circuits. This functionality explains where its name comes from. To develop more complex circuits and functionalities of future bio-inspired soft robots, digital DEOs acting as clocks will become very important to control and synchronize processes.
With this work we want to confirm the validity of existing SIMULINK models, like the one presented in [22], for DE circuit networks. We demonstrate the correct functioning of the model with one of the most complex DE circuit network that can be obtained with the current processes, the digital DEO, and show its adaptability to different samples. Having such an easy-to-use simulation toolkit will be required to tailor the DEO characteristics to the desired application or to study new DE circuits, without the need to physically build them, saving materials and time.

Functionality of dielectric elastomer oscillators
DEs are composed of an elastic membrane, sandwiched between two compliant electrodes, forming a compliant capacitor. When a high input voltage V in (in the order of kilovolts) is applied to the electrodes, the opposite charges on different sides attract each other and the like charges on each side repel each other. This results in a constrictive force, the so-called Maxwell pressure, that makes the structure expand inplane and contract in thickness. This functionality of actuation is the basic concept of DEAs. It is possible to pair this structure with a DES, as firstly demonstrated in [27]. The working principle is shown in figures 1(a) and (b). The DES itself is an area on top of the membrane, covered with a piezo-resistive carbongrease composite electrode, that can change its electrical resistance, due to mechanical strain, over several orders of magnitude. When no voltage is applied to the input terminal of the DEA, there is no actuation in the DEA, the resistor is in its rest state and the resistance is high. This is due to the percolation effect: in the rest state, the conductive particles are dispersed with low density and do not form conductive paths in the matrix, thus they do not allow the current to flow through them [28]. When the DEA is actuated, the DES is compressed and, thus, the resistance lowers drastically as shown in figure 1(b). Such a combination of DEA and DES will be further addressed as unit or gate.
By combining multiple DEA-DES units and connecting them in different ways, it is possible to build electronic components, such as logic gates [8] and basic memory structures [9]. It has also been demonstrated that a DEO can be realized by linking together an odd number of DE-inverters [23,26,29]. In this contribution, a digital DEO composed of six DESs and six DEAs (figure 2) is investigated. This  structure has been firstly proposed in [23]. Based on the research, done over the last years, we reported here an optimised version of it, with different units, like the ones used in [8,26] and shown in figures 1 and 2. The digital DEO has an output voltage more similar to a square wave, compared to the version with only three DES-DEA units that performs a sinusoidal-like oscillation [22].
As described in [8], when using this kind of DEs for logic gates, there is significant time delay. This is because of the highly viscoelastic behaviour of the material used, (VHB 4905, by 3M). When a high voltage is applied to a DEA, the actuator will not reach its full expansion instantaneously, but it will slowly creep to the maximum. In a DES this translates to slower resistance changes. For the opposite transition (when the applied voltage is reduced, or absent), instead, there is almost zero delay because the elastic restoring force dominates over the viscoelasticity (the difference between the two delays has been shown in [8]). Now, if we consider the DEO with three DES, in each inverter the resistor of the DES (R DES ) is in a voltage divider with a fixed resistor R S (usually in the order of 10 7 Ω [22]) as shownin figure 3(a). If V CC is on and V in is zero, then V out is high because all the voltage drops across R DES . When V in is raised, R DES starts to decrease but it takes some time for it to reach the minimum. This means it will take some time for V out to reach the low logical output value (we need R DES R S to have a low output). In the DEO with six DES, instead of R S there is a second DES (figure 3(b)), and applying a high voltage means that while R DES2 is slowly decreasing, R DES1 will experience a very fast transition (because this is a relaxation transition) to its maximum high value, since this is way higher than R S was (see figure 1), less time is required in order to saturate the output value (meaning reaching logical '0' in this case). Considering now the whole DEO, this affects the shape of the output wave that will look more like a square wave in the six DES case, because of a shorter time required for the saturation of the output.
However, also because the building process is all manual, adapting the model to this DEO is more challenging: there is double the number of structures to simulate and more human-induced error.

Materials and methods
A DES-DEA unit is built from VHB 4905 tape by 3M ( figure 4), this is a commercially available acrylicbased tape that can be readily used as elastomer for DEAs. A piece of this tape is mechanically prestretched equi-biaxially to a value of λ pre = 2. The piezo-resistive DES resistor pattern is then applied with a PDMS stamp using a mixture of 5:1 Molykote 44 medium grease to Cabot Vulcan XC72 carbon black. Afterwards, the membrane is further prestretched until the final pre-stretch ratio λ pre = 3.6 and stuck to PMMA frames. The remaining electrodes and conductive tracks are manually painted with carbon grease (Nye's Nyogel 756G). A more complete discussion about the methods is in [30,31].

Equivalent circuit
The system of a digital DEO can be represented by the equivalent circuit in figure 5, each unit is composed  The goal of the model is to calculate the output voltages of the three inverters: V out1 , V out2 and V out3 . To do this, it is possible to divide the DEO into three electrically independent circuits as highlighted in figure 5 (green, orange and yellow rectangles). In figure 6 a single stage is considered, where the contribution from the resistor for voltage measurement (R meas ) and the electrical behaviour of the viscoelastic membranes, through which there will inevitably be some leakage current, are added. The leakage is represented by a variable resistance in parallel to each capacitor (R mem1 and R mem2 ).
Therefore the DEO structure in the model is divided into three stages like figure 6 representing the three electrically independent and identical circuits present inside the DEO. However, these stages are mechanically dependent to one another, hence they have to be simulated altogether. The capacitors and resistors of a single inverter are mechanically dependent to one another (when a C DEA has a voltage applied, the actuator will actuate and compress the DES lowering the correspondent R DES ) but are part of different electrical circuits.

Model
The principles of modelling such structures were firstly presented in [22]. This chapter briefly describes its fundamentals and the modifications and  improvements that have been made, in this contribution, to adapt it to the digital DEO. The SIMULINK model will be explored in a top-down approach. Figure 7 describes the interaction between the three main blocks that the DEO is divided in, and the signal paths. Each stage outputs the values of actuation (u) of its DEAs and passes them to the next stage. This is needed because the actuation strains influence the resistors of the following stage as explained in the previous section. The individual stages contain circuits as depicted in figure 6. Figure 8 represents the inside the three stages, all three are identical. They are each composed of two identical sets of blocks ('DES', 'actuation', 'DEA', 'membrane') corresponding to the two gates that compose an inverter. Depending on their input, they calculate the resulting parameters and pass them on, to the electrical circuit where the calculations for the V out are executed.
The 'actuation' block is shown in figure 9(a). Its purpose is to calculate the actuation of the DEA that corresponds to the applied voltage, taking into account the viscoelastic behaviour of the membrane. As mentioned before, the viscosity of the membrane will resist the mechanical forces that try to deform it. This results in a general delay of actuation. In the case shown here, it is possible to see two main effects ( figure 9(b)): an incomplete actuation when the high voltage is applied, and a slow drift of the actuation to higher values over time, after the voltage is applied continuously for several minutes. Finally, when the voltage is removed, the actuation rapidly jumps back (typically to u ≈ 0.1 mm) due to the elastic restoring force. However, it takes around 5 min for the membrane to completely relax to its initial state (u = 0 mm).
The calculation in the 'actuation' block ( figure 9(a)) is based on a transfer function, which is   derived from a generalized viscoelastic model of the fifth-order [32]. The basic equation of the model is: where a i and b i are parameters to be determined.
To integrate this equation into SIMULINK, it is necessary, to transform equation (1) in the Laplace domain: (2) and then use the access transfer function data (tfdata) to obtain the transfer function:

U(s) =
As 5 + Bs 4 + Cs 3 + Ds 2 + Es + F Gs 5 + Hs 4 + Is 3 + Js 2 + Ks + L , where the coefficients from A to L have to be determined. This is done by measuring the step response of the DEA actuation, as shown in figure 9(b), with a high precision laser displacement sensor (optoN-CDT 1402 by Micro-Epsilon). When a high voltage step is applied, the DEA instantly moves to an initial displacement and then creeps slowly to a maximum actuation value. It is possible to extract the required parameters from the measured actuation values using Matlab. The actuation can be expressed as the sum of 5 exponential terms as shown in equation (1) (the fifth-order has been chosen to have a high enough precision, lower orders did not produce an accurate fit). The measured actuation curve is divided in 5 parts and for each part a fit with the equation a i 1 − e (−tb i ) is performed, until the final sum of the fits closely resembles the measurement like it is shown in figure 9(b). The Matlab fits calculate the coefficients a i and b i ; finally the coefficients A − L can be extracted following equations (2) and (3). A more complete discussion about this process can be found in [22]. The 'DES' block, shown in figure 10, calculates the R DES (u in ) parameter, as result of the DEA's actuation u in . This is done by using a lookup table: for a specified value of the actuation it gives back a value for the resistance of the DES, R DES .
To generate the data required for this lookup table, a measurement system composed of a single gate was set up. A high voltage generator specifically designed for DEs [33] was used to power the DEA and a precision laser displacement sensor optoNCDT 1402 from Micro-Epsilon to measure u(t). The data acquisition was done with the NI-9205 module (DAQ) by NI and data was recorded with a LabView program. In order to be able to measure high voltage signals, with a low voltage DAQ, the measured voltage was reduced by a custom-built pre-amplifier, possessing a voltage divider. The resistance of the DES and the actuation of the DEA with a varying input voltage were measured and used to build a lookup table in Matlab. Figures 11  and 12 depict the measured values, and comparison to simulated values. As input voltage, a triangular signal with a period T = 60 s was used. The long period was chosen to measure a quasi-static relation between actuation u in and R DES and to observe longterm behaviours due to the viscoelasticity. Figure 11 shows the correct prediction of the device's parameters in time. Figure 12 illustrates that the model reproduces the hysteresis typical of carbon composites inks [34], this is originated by the peculiar electrical conduction in such materials. The conductivity depends on the applied electric field and it needs some time to settle when there are variations of electric field. The next block is the 'DEA' in figure 13, it is the implementation of the formula for the variable  capacitance of the electric model of the DEA, a parallel plate capacitor: where 0 is the permittivity of vacuum, r the relative permittivity of the elastomer (based on [35]), A the area and d the distance between the two plates. This component monitors the variations on the actuation value and modifies accordingly the value of the capacitance.
In figure 14 there is the component for the estimation of the membrane resistance value, this is based on the formula discussed in [18,36]: where σ 0 is the value of the initial conductivity of the VHB, E B is its breakdown field value, E is the applied electric field, and A is the area. This block takes the geometric variations on the DEA and the electric field inside the capacitor into account and  obtains the resistance value by inverting equation (5) through Ohm's law. Finally, the block 'electrical circuit' in figure 8 represents the implementation of the circuit in figure 6 using the SIMULINK electrical package, here the voltage V out is calculated. In this block, all the components are linked together. The values of the varying electrical components are given by all the other blocks described in this chapter. The DEAs' capacitances are multiplied with a factor that accounts for their variations in area due to manual production. This is needed because, at the present state, the DEAs are hand-painted, thus their shape is irregular and also the thickness of the ink cannot be precisely determined; both these parameters are important for a precise calculation of C DEA . We aim to improve this aspect by introducing new and automatic processes (e.g. spray painting, doctor blading, ink-jet printing) for the device fabrication that will allow an easier estimation of DE properties. Figure 15 shows the performance of the model for a single NOT logic gate, this is obtained with the circuit in figure 3(b). The input voltage is correctly inverted at the output. The delay caused by the viscoelasticity of VHB is reproduced by the model and the simulation shows a good overall accuracy. As mentioned before, building a model for a combinational circuit (e.g. logic gates) is easier and gives better results than for DEOs. Combinational circuits, in fact, are governed by the input signal, any variation in the system is a consequence of a variation of the input. Whereas for a DEO, the output voltage is a consequence of the interactions between the internal components, the input is a DC voltage that does not change in time, thus it is critical to correctly predict the parameters of these components. This results in more complex analyses. Furthermore, looking at figure 3(b), the electrical circuit is composed only by a voltage divider, the DEAs are electrically independent from the DESs and do not need to be modelled (ground and V CC decouple the two parts of the electrical circuit). Thus, the model does not require subsystems for the leakage current and for the capacitance of the DEA. These aspects contribute to reduce the imperfections of the model that will be analysed in the remaining part of the section.

Results
Then, to deeper test the model, we present the results obtained from applying the model to three different digital DEOs (figures [16][17][18]. The DESs that are used in DEOs, as depicted in figure 2, can be made with different numbers of switching tracks and, thus, a different length of piezo-resistive membrane. Here we present the results for DEOs with five and nine tracks. In general, it was recognized that ninetracks DESs have a more stable performance over time because the DES resistance degrades more slowly in time (the ink gets dry with time and the overall  resistance increases). Figure 16 compares the measurement and simulation results for a five-track DEO with an input voltage V in = 2850 V applied. This DEO possesses pre-resistors R S = 50 MΩ in series to the upper DESs resistors to limit the maximum current. After a settling time, which is also shown by the model, the simulated signals fit the experiment. In particular, we can see how in this case neither the model nor the real DEO arrives at the maximum voltage available, this is due to the current limiting resistor R S . Furthermore, there is also a difference between the maximum of the measured and simulated voltage. This is most likely due to the imprecise leakage current prediction. For this parameter, equation (5) is used, this is an equation that works for standalone DEAs similar to the ones here presented (same materials) [18], but it is difficult to study how they behave inside a circuit network while working, thus we expect this parameter not to be completely precise. Also, the DES degradation over time contributes to the differences: over time, the minimum resistance starts to increase until the DES stops working as a switch (if the low and high resistances become of similar values there is no more  switching behaviour). We also calculated the frequency of the real DEO and the modelled DEO, the results are very similar as it can be seen in figure 16(a). For this measurement, we did not take into account the initial few seconds where the oscillation is not stable, and we used a Matlab script to perform the task, finding a relative error E r = |f meas −f model | f meas ≈ 0.037 (in percentage 3.7%) based on the frequency calculation with Matlab.
In figure 17 we present the same comparison for a different structure, the DESs now have nine tracks each and there are no pre-resistors R S . A V in = 3070 V is applied. A photograph of the DEO is shown in figure 2. For this DEO the output wave-forms look more similar to a square wave than the output waveforms of the DEO with five-tracks DESs. In particular, the measured V 1 reaches the maximum input voltage. This is because there was no R S , which was also reproduced by the model. V 2 and V 3 do not reach the maximum voltage output like V 1 , probably because their DESs were not of the same quality as the first ones. If the difference between low and high resistance of a DES is not high enough (three orders of magnitude or more) the voltage output of an inverter is not saturated to the maximum and minimum values (as can be seen in figure 3, the output voltage is given by a voltage division between two switches, one in the high-resistance state and one in the lowresistance state). However, there is good accordance between the model and the experimental values, in particular, the frequency in the two cases is very similar, the percentage error is around 1.2%.
The last test has been done with a second DEO with nine-tracks DESs. Because of the handmade process, there are small differences (e.g. the capacitors' shapes and thicknesses, the DESs quality) that will result in a different overall behaviour as it is shown in figure 18. The model is capable of adapting if its parameters are updated with measurements that fit better to the current structure (e.g. how the resistance varies with the input voltage and with actuation); indeed, also in this case, the frequency estimations are very close (error 2.1%).
Despite all the aspects that could be improved (e.g. production, deeper studies on the DEA-DES units), we have been able to fit the model to different structures with a maximum error of 3.7% for the most complicated structure (with also current limiting resistors R S ).

Discussion
We showed how to build a model that takes into account the main electro-mechanical properties of DEAs and DESs, puts them together, and simulates the behaviour of the resulting DE circuit network. Starting from an existing SIMULINK model [22], developed for sine-wave output DEOs, we adapted it to different circuit networks (digital DEOs) that are more complicated to simulate because composed of double the number of electro-mechanical elements. Modelling a single aspect of DEs is already a challenge [13,14,[18][19][20]; to simulate DE circuit networks it is required to combine multiple individual models for each DEA-DES unit and interface the units with each other. As shown in figures 11 and 12 the behaviour of DEAs and DESs is correctly taken into account by the model, then it is possible to interconnect these elements to obtain a working model of the final DE circuit that predicts its voltage output in time (figures [16][17][18]. At the present state, there are some limitations in the model. To obtain the best fitting results (the ones showed in figures [16][17][18] it is necessary to correct some internal parameters that cannot be known before the physical realization. The DEA capacitance C DEA and the DES resistance R DES are affected by variations induced by the manual production process. The DEAs' shape is not regular because it is handpainted, also its thickness cannot be precisely defined. The DESs do not always have constant values for the high-resistance state and low-resistance state, thus influencing the voltage output of the DE inverters ( figure 3); furthermore, the DESs degrade with time due to the current flowing through them that dries the ink and ruins the switching behaviour (this aspect is currently not implemented in the model due to the difficulty of estimating it). As an example, comparing figures 17 and 18, these two devices are made in the same way but the small differences we just mentioned change the overall behaviour. Thus, the model parameters need to be updated with measurements on the DEA-DES units (reproducing the same experiment shown in figures 11 and 12) each time a new DE circuit is built. Another problem is the maximum resistance of DES that, we think, our measurement setup does not allow to correctly measure. This is because it is difficult to precisely measure a resistor with such high variability in values and we accepted a compromise to have a good measurement on the overall curve. However, when modelling combinational circuits, the model performs better (as shown in figure 15). In this case, some aspects of DEs do not need to be modelled (e.g. leakage current, C DEA ), the only parameters to take into account are the variations of the switch resistances (R DES ) with the voltage applied. This reduces the number of components that need to be simulated and the related problems.
We showed that the model can reproduce the behaviour of complex DE circuit networks like digital DEOs and we believe that it will benefit from new materials and production processes for DEAs and DESs. With a reduced sample-to-sample variation, a more complete model of the single units could be made. This will lead to a more robust and precise model for DE circuit networks that could allow their design and investigation without the need to build them beforehand.