Mathematical Modeling of a Tactile Sensor with Applications in Minimally Invasive Surgery

This paper describes mathematical modeling of a novel supported tactile sensor. This system has only three or four sensing elements instead of an array configuration. By writing a code, together with using the piezoelectric capabilities of the finite element analysis software, it is possible to extract the stresses developed in the tactile membrane along the edges of the sensing elements, and convert them to output electrical charges. The comparison between the proposed mathematical modeling and the experimentally obtained data shows good agreement.


INTRODUCTION
In biomedical engineering applications of tactile sensors, defining the state of manipulating or gripping of a biological tissue is of great importance [1] . This is achieved by determination of two important physical parameters, i.e., force and position signatures [2] . To accomplish this goal, different types of tactile sensors can be utilized in detecting the presence or absence of a grasped object/tissue [3] . In some cases, an array of sensors might even provide the capability of measuring a number of tactile properties. By designing of a suitable sensor, it is possible to map a complete tactile image of the structure under investigation [4] . These systems have a considerable number of potential application areas including minimally invasive surgery (MIS) [5] . MIS is now being widely used as one of the most preferred choices for various types of operations. In MIS, any inhibitions on the surgeon's sensory abilities might lead to undesirable results. MIS has many advantages; however, it decreases the sensory perception of the surgeon and the surgeon might accidentally cut or incur damage to some of the tissues [6][7][8][9][10] . This effect is more pronounced when the surgeon approaches a target tissue while moving past other healthy tissues.
Currently, various kinds of MIS, such as laparoscopy and arthroscopy, are being routinely employed as the preferred choices, comparing to the conventional surgeries [11][12][13] . Because of the nature of MIS, which makes use of both visual and tactile capabilities of the surgeons, any impairment on these abilities should be dealt with properly. The success of a MIS procedure depends highly on the surgeon's ability in feeling the tissues and detecting the presence of blood vessels and ducts while performing the surgery. Therefore, for these delicate biological organs to be handled more safely, the measurement of the magnitude and location of the applied loads can convey a useful feedback to the surgeon [14][15][16][17] .
Most investigators have attempted to design a tactile sensor using an array of sensing elements arranged in matrix form [18][19][20] . The output voltages from individual sensing elements represent the magnitude of the applied force and the number of sensing elements pressed indicates the position of the object in relation to the tactile sensor. There are several problems associated with this type of tactile sensors [21,22] . These include cross-talk between sensing elements, fragility, and complexity.
Following the above research activities and our previous work on the design and manufacturing of a tactile sensor with only four sensing elements, we are now proposing a novel approach for mathematical analysis of the designed system [7] . A geometric mapping process is introduced, thereby the loci of the isocharge contours for the sensing elements are determined by applying force on various points of the sensor surface.

Triangulation approach:
The identification of a point on a two-dimensional plane requires the value of its coordinates to be known. The position of a point on the plane can be found from knowledge of its distance from three reference points, the coordinates of which are known. This approach used for position detection is similar to the concept behind global positioning system (GPS). The need for the distances from three reference

SCI-PUBLICATIONS Author Manuscript
Open Access

Author Manuscript
Am. J. Applied Sci., 4 (10): [779][780][781][782][783][784][785]2007 780 points can be easily seen by considering the case where the distances from only two reference points are known. These distances define a circle about each of the reference points, and their intersection gives the position of the unknown point. There are, however, two intersections so that these distances do not unambiguously define the position of the unknown point. This is shown in Fig. 1, where points a, b, and c represent three arbitrary points on a plane. This means the triangulation approach for position identification requires a minimum of three data. The equation used for position determination is: where x n and y n are known positions of the three sensing elements, r n is the distance between applied force's position and sensing element, and finally both x p and y p represent the unknown positions of the applied force. If the radii are known for the three circles and the position of the sensing elements are known, it will be possible to determine the position of the applied force, since we will have three equations (with n = a, b, c) and only two unknowns (x p ,y p ). With only three sensing elements, obviously, it will be difficult to obtain the shape as in the matrix arranged sensing elements, but this is just to illustrate the fundamentals of the triangulation approach and its potential to simplify the sensor's structure. From the triangulation equation, it is possible to know the distance between the position of applied force and sensing element. Obviously, if the applied force changes position, so does the radius. Therefore, the key idea of this sensor design is using triangulation method to find these distances, which will be discussed in the following section. The mapping scheme is shown in Fig. 2. First, one of the electrodes is selected. Then, a number of lines and circles are drawn with respect to the electrode center. The radial and angular increments are, respectively, 4 mm and 10°. Calibration of the sensor: The layout of the sensor is shown in Fig. 3. The PVDF film is considered as a circle which is clamped all around its perimeter. The center of each of the three sensing elements, A, B, and C, is 10 mm away from the center of the membrane, or the origin of the coordinate system. If the centers of all three electrodes are connected to each other, an equilateral triangle will be created. Lines connecting the center of the membrane to the centers of electrodes B and C make angles of 60° and 120° with the x-axis (or the drawn direction of the biaxial PVDF film), respectively. The center of electrode A is located exactly on the y-axis (or the transverse direction of biaxial PVDF film).   Finite Element Analysis: For this purpose, ANSYS software (version 10.0) was extensively used. Since ANSYS is capable of running multidisciplinary simulations like piezoelectric problems, at first an attempt was made to run a piezoelectric analysis. This analysis enables the user to input load or deflection and directly obtain the generated output charges (or voltages) developed on the surfaces of electrodes. There are a number of elements available in ANSYS which are ideally suited for piezoelectric analysis of volumes such as SOLID 5, SOILD 98, SOLID 226, and, SOLID 227. However, regardless as to the type of piezoelectric element used, all of the meshing attempts failed due to the extremely large aspect ratio of the circular membrane, i.e., (radius/thickness) = (45 mm/25µm)= 1800. We used an alternative solution which was to run a structural analysis, extract the resultant stresses, and finally convert them to charges using piezoelectric formulae. The most suitable element for a structural analysis of shells and membranes is SHELL 63 element. This element has both bending and membrane capabilities in which both in-plane and normal loads are permitted. After creation of the PDVF film and its sensing elements, the shape of the probe, either triangular or rectangular, was added. For the sake of brevity, most of results presented in this work belong to the triangular probe. Once we generated all the areas, a meshing process was performed. To achieve this, an element size of 0.75 mm was selected for meshing the probe and four square sensing elements. For the remainder of the membrane area, a size of 1 mm was chosen. Because of the nonlinearity of the static problem, all of the areas were meshed using triangular elements. For a planar element such as SHELL 63, it is often possible to choose two geometries: quadrilateral and triangular, but triangular element proved to be more appropriate since it resulted in a considerably more uniform meshing as depicted in Fig. 5. Although the number of generated elements varied from one model to another, depending upon the shape and location of the probe, the number typically fell in the range of 13,000 to 14,000. This ensures that the produced mesh is sufficiently uniform and captures the entire geometry well enough.
Structural nonlinearity: Geometrically speaking, the problem at hand can be categorized as a nonlinear type. In a broader sense, structural nonlinearity or nonlinear structural behavior arises from a number of causes including changing status, material nonlinearities, and geometric nonlinearities. In the case of changing status, many common structural features exhibit nonlinear behavior that is status-dependent. For material nonlinearity, this was not applicable because PVDF exhibits good linearity. The relationship between the applied force and the output charge is linear for stresses up to 40 MPa. In our case, however, it is the geometric nonlinearity issue, which is the main reason for nonlinear behavior of the model. It is known that if a structure experiences large deformations, its changing geometric configuration can make it respond nonlinearly. For the PVDF model of this project, geometric nonlinearity can be explained on the grounds of dimensional properties of the model and also the magnitude of the applied load. Although a load of 1 N might seem too small to cause large deflections and any sort of nonlinearity, it must be realized that this load is applied only to the area of the probe, which is very small, hence, generating noticeable pressures. This coupled with the dimensional properties of the PVDF film, leads to a geometrically nonlinear system. For the PVDF film and SHELL 63 element, due to the extremely small thickness of the membrane, it cannot be taken as being a large strain but rather as one of large-deflection.
The stress stiffness matrix was computed based on the stress state of the previous equilibrium iteration. Thus, to generate a valid stress-stiffened problem, at least two iterations were normally required of which the first was used to determine the stress state that will be used to generate the stress stiffness matrix of the second iteration. If this additional stiffness affected the stresses, further iterations were required to reach a solution. For the PVDF film, the film was modeled as a circular membrane which is clamped all around its perimeter. Using the "solid modeling" feature of ANSYS, boundary conditions were applied as constraining the outer perimeter of the circle. Instead of counting the number of nodes generated on the probe area and distributing the lateral pressure among them uniformly, the pressure was applied to the whole probe area as one entity.

RESULTS AND DISCUSSION
Since the type of analysis was selected as "structural", evidently the generated results are stresses and/or deflections. By writing our own code, it was possible to obtain directly the output electrical charges that surfaced on the sensing elements. For the meshing process, the element edge size for each electrode was 0.75 mm. The formula for computing the electrical charge is: where A is the area of the electrode surface and Q is the charge. If the average stress components along the outermost edges are found and plugged into the piezoelectric formulae, the corresponding obtained charge is precise. In accordance with these discussions, considering σ x for electrode A as an example, first line integrations of σ x along its right and left edges were done. Then, the result for each edge was divided by its length and the average of two sides was considered as the final value for σ x . This is better illustrated in the following equation: In the above equation, s is the total length of the path (which is the edge of the electrode) and d s is the incremental length of each edge. So, the value for s is 3 mm for all the four edges of all the four electrodes A, B, C, and D. Using the coordinates of nodes created on the right and left edges of electrode A, each of the edges is defined as a path.
The triangular probe dislocates between intersections of lines and arcs. Here is how it works: for the first set of simulations the probe is moved along the right-most line from one intersection to another with radial increments of 4 mm. For each line, which is drawn at a certain angle with respect to the x-axis, the resulting output electrical charges on the electrode A are obtained. This is followed by plotting a number of curves for each angle as shown in Fig. 6. In this figure, it is shown how the output charge changes against distance while moving on lines drawn with different angles, from 35º to 90º, with respect to the drawn direction. Another mapping is carried out, with respect to electrode B, which is shown in Fig. 7. Using the same procedure leads to the following curves shown in Fig. 8 for rectangular probe. In the above figure, the lines along which the force moves, make angles with the drawn direction which vary in the range of 140º to 240º. Throughout the entire simulations, the load transferred through the probe was assumed as 1 N. It is worth mentioning that the uppermost curves correspond to the smallest angles; i.e., 140° while the lower-most curves represent the largest angles; i.e., 240°. The next step is utilizing the triangulation technique. As mentioned earlier, this is a method for determining the location of a point in a plane. The earlier description of triangulation approach was based on the distances of an arbitrary point within a plane from three reference points. Here, it is dealt with in a different way. That is, if a force is applied onto the surface of a PVDF film via an arbitrary shape, stresses and consequently electrical charges will be developed on the electrodes. For a given shape of probe coupled with its fixed dimensions and a known magnitude of applied force, it is possible to find out the loci of points which all share something in common. If a particular force is applied at each of these points via the predetermined probe, a certain amount of charge will appear on the electrode. The connection of all those points leads into what is called an "isocharge contour". Therefore, the procedure for detecting the location of a force is redefined as follows: given the output electrical charges on the three electrodes are, namely, X, Y, and Z, for each electrode, the related isocharge contour is drawn. The intersection point of them yields the desired data.
Creation of the isocharge contours can be done by using the charge versus distance curves depicted earlier.
As an example, it is assumed that the charge on electrode A is 650 pC. A horizontal line corresponding to 650 pC is drawn so that it crosses all the graphs. The intersection points are what we need to construct the isocharge contour. Assuming this was done in Fig. 6, the following pairs of data could be extracted: 35°, 7 . Each of these points should be identified in a polar system, the origin of which is the center of the referred electrode. Eventually, they are connected to each other and form an isocharge contour. A polynomial was fit to each set of data and then their corresponding coefficients were found. Following this, the obtained equation was set equal to the desired output charge, and its root was found. The answer is one point of an isocharge contour. Later, the points were joined in AutoCAD using the SPLINE function. A triangulation is presented in Fig. 9 for a triangular probe. Fig. 9: Simulation with a triangular probe.
Since electrodes B and C are symmetrical with respect to the y-axis or transverse direction, it is possible to create the isocharge contour relative to electrode B and then use the MIRROR command of AutoCAD and mirror it relative to transverse direction. In Fig. 10, all three isocharge contours are drawn. The three contours are not concurrent, i.e., they do not intersect in one single point. Hence, it was decided to take the coordinates of all three intersecting points and compute their average as a unique yet imaginary intersection point. Using the DIST command of AutoCAD, the distance between the exact location of the center of the probe and the calculated average was found equal to only 0.122 mm, which is acceptable.
Before proceeding to the next simulations, two points should be cleared. First, in some simulations, the probe had a small overlap with electrode B. This may raise a concern whether the third piezoelectric coefficient, d 33 , should be taken into consideration or not because the electrode area is directly touched and under pressure of a force along the 3-or thickness direction. The answer is negative. That is so because, although electrode B undergoes a force along a direction perpendicular to its face, it is part of a very thin membrane (just 25 µm). As a result, there is no stress or strain in the 3-or thickness direction, which means no charge is produced because of d 33 effect. The second point to consider is the pyroelectric effect. When a foreign object or probe touches the electroded area of a PVDF film, a sudden change of temperature can occur. That leads to a secondary source of charge generation. At times, this effect can be noticeably large and even dominate other sources like stresses. Throughout simulations, it was assumed that, even if the probe comes into contact with an electrode area, no temperature alteration happens. This means, in effect, that the pyroelectric effect in the entire set of simulations has been ignored.
For all the cases that we studied in four sensing element layout, the results were in accord with those obtained in our experiments that were presented previously.
In Fig. 11, numerous isocharge contours for all four sensing elements are put on the same picture. These contours were created by moving a force of 1 N across different points of the central area of the sensor. The isocharge contours shown here were obtained by moving a circular probe of 1 mm radius. For probes with different shapes, the isocharge contours may change according to the shape. Nonetheless, this figure gives a clear view of how isocharge contours are created and located on the sensor area. Fig. 11: Isocharge contours for all four sensing elements under a force of 1 N. Figure 12 shows the simulations results for charge versus distance curves for electrode A, when a triangular probe used. As noticed in the figure, the charges vs. distance curves are more distinct and separated. This is a proof of how effective the piezoelectric coefficients can be in defining the behavior of the system. In these simulations, the load was taken as 1 N. The mentioned angles are the angles between drawn direction and the lines along which the probe moves.
Effect of thickness of the PVDF film: Considering the results reported, the main reason for the nonlinear variation of output charges with respect to the force can be stated as "nonlinear geometrical deformation" of the membrane. This problem becomes an issue when a film or membrane is extremely thin. To answer this question, first, the Web site of Good Fellow Co. (USA), a PVDF manufacturer, was consulted to see if there existed other commercially available, thicker PVDF films. It was noticed that PVDF films are produced in thicknesses of up to 1mm. As two arbitrary choices, two PVDF films with thicknesses of 110 µm and 1 mm were chosen. A three-electrode layout was assumed and various loads were applied through the same triangular probe. Our simulations show that an increase in thickness drives an enhancement in the linearity of the sensor response for all electrodes. However, the nonlinearity still exists. To verify that, the following equations are presented which relate charge to force when the thickness of PVDF film is 1 mm: The coefficients of 2nd and 3rd degree terms are not negligible at all so a linear approximation is prone to great inaccuracy. Furthermore, with the increase of film thickness higher strains and stresses will develop along the thickness direction of the film. This leads to a bigger role for the d 33 coefficient and a third source of charge generation. With all the above results in mind, it can be concluded that the triangulation technique can be used for even larger loads, at least up to 2 N. The only setback is that the related electrical charges on the electrodes do not vary linearly. Therefore, in order to find both magnitude and position of an applied load, the following solutions are recommended: (a) The sensor can be used in a specific band of forces, for instance, up to 0.5 N. In this range, the variation of charge against force is linear and interpolations are easier and more accurate. (b) A large database, or actually a lookup table, for different forces and positions should be created and stored in a computer. Then, for each unknown set of variables (force and position), the results should be compared or matched with the lookup table. Of course, in this method, the shape and size of an object in contact with the sensor must be known in advance.

CONCLUSION
Because of the complexities involved in closedform solution of such problems, Finite Element Analysis was the only viable theoretical approach. It was not possible to use piezoelectric elements. Instead, a structural analysis was performed. By writing a code, combined with using the basic piezoelectric formula, it was possible to extract the stresses developed in the membrane along the edges of the sensing elements, and convert them to output electrical charges. The triangulation method was used for detecting the position of the applied force. A mapping scheme was proposed in which, force was applied at different selected points within the central area of the sensor, confined to the edges of the sensing elements. Using the calibration method, and by applying the force through two shapes of probe, triangular and rectangular, a number of charge against distance curves were obtained. These curves were used for creating isocharge contours. By doing the triangulation four times for a four-sensing element layout, and taking the averages of the coordinates of the intersection points, a more accurate knowledge of the application point of the force was obtained. As for the variation of the output charge against force, the sensor was linear up to 0.5 N. For magnitudes of forces beyond this value, the variation of the output electrical charge versus the magnitude of the applied force was found to be a cubic polynomial.
Work is currently underway in our lab to detect the shape of the touched object using the proposed method.