An isogeometric scaled boundary plate formulation for the analysis of ionic electroactive paper

In recent years, electroactive paper emerged as a new alternative in the field of smart actuators. It is based on a cellulose material which is able to bend under the influence of an external electric field similarly as ionic polymer metal composites. The bending mechanism is mainly attributed to the migration of ionic charges over the thickness of a thin sheet of paper. The present contribution proposes a numerical framework for the simulation of electroactive paper. It is based on a scaled boundary plate formulation for isogeometric analysis. In contrast to the standard scaled boundary plate approach, the scaling direction is solved numerically by a B-Spline approximation. This allows to render nonlinear effects over the plate thickness as well as displacement fields of higher continuity. The model is applicable to very thin structures such as electroactive paper, and it also captures the nonlinear ionic charge distribution which is coupled to the bending mechanism of the actuator.


Introduction
The term electroactive paper (EAPap) firstly appears in [1] in which cellulose is proposed as a smart material. Previous research has already shown that wood, which is made out of cellulose, possesses piezoelectric properties [2]. More specifically, shear piezoelectricity is observed in cellulose crystallites [3]. Furthermore, it is also known that cellulose is a hygroscopic material which is capable to absorb water from ambient humidity [4].
EAPap can be fabricated by dissolving a cellulose film in sodium hydroxide. In this process, free sodium ions are attached to water molecules which are able to move freely under the influence of an external electric field [5]. Consequently, EAPap can be subdivided into ordered and unordered regions. The ordered regions are defined by the monoclinic crystal structures of cellulose with dipolar orientation which are responsible for the piezoelectric effect. This effect can be improved by the consideration of electric and magnetic fields as well as mechanical stretching [6]. The unordered regions are defined by water molecules which are attached to ions such as sodium. These regions are also affected by an external electric field by means of an ionic migration [8]. This means that the water molecules, which are attached to ions, are able to travel to the anode of the actuator creating a heterogeneous density and electric charge redistribution over the actuator thickness. Consequently, the bending effect of the EAPap actuator can be attributed to shear piezoelectricity as well as the ionic migration. Nevertheless, it is assumed that the main cause of the bending deformation is defined by ionic migration [8].
Regarding the mechanical properties of EAPap, it should be mentioned that cellulose is an anisotropic material due to the fibrous structure. The fiber orientation can be defined during the production process in order to control the Young's modulus in the longitudinal direction of the actuator. The range of the Young's modulus is given between 4 and 10 GPa in dependency of the fiber orientation [1,5]. Furthermore, EAPap can also be fabricated with hybrid materials such as multi-walled carbon nanotubes, conducting polymers and ionic liquids in order to enhance the ionic migration under dry ambient conditions [6]. The electrodes of the actuator are generally produced by coating the EAPap material with a thin gold film in the range of nanometers [8]. Also, here an actuation enhancement is possible by considering fishbone patterns for the electrodes structure [9].
EAPap can be regarded as an ionic electroactive polymer (EAP) whose actuation is driven by the displacement of ions inside of the polymer. A popular example of this technology is ionic polymer metal composites (IPMCs), see [15,16]. A second group of EAPs are defined by dielectric materials which are sandwiched between compliant electrodes. The actuation mechanism is defined by a transverse elongation of the film which is caused by the attraction of the electrodes [18]. The main drawback of this technology is given by the very high actuation voltages. Therefore, heterogeneous microstructures have been proposed in order to reduce the actuation voltage [19,20].
In this respect, EAPap has the advantage of a lower actuation voltage and consequently a reduced power consumption. The required electric fields are of one to two orders of magnitude lower than for electrostrictive polymers like silicone rubber or thermoplastic polyurethane [1]. Additional advantages of the material are that they are lightweight, inexpensive and biodegradable. In general, it can be used as a piezoelectric sensor, an artificial muscle or an energy harvester for microelectromechanical systems [10][11][12]. Some applications include micro-insect robots, micro-flying objects, biosensors, flexible electrical displays, etc.
In [1], first experimental results of an EAPap actuator under dynamic loading are presented. It is shown that the tip of a 30 mm cantilever strip of EAPap with a thickness of 30 µm can achieve an oscillation amplitude of 4.2 mm. The applied time-dependent voltage has an amplitude of 7 V at a frequency of 7 Hz. For the given sample thickness, the applied electric field corresponds to 0.23 V/µm. Furthermore, the ambient relative humidity in the experiment chamber is of 95 % and therefore enhances the bending deformation.
In [5,7], an algebraic model for the deformation of EAPap is proposed and compared with a quasi-static experiment. The model consists of a multilayered approach which is able to calculate the bending curvature of an EAPap cantilever actuator. Additionally, shear piezoelectricity is also taken into account. A good agreement of the model with the experimental data is presented, but the model constants are not provided. The experimental results show the dependency of the bending deformation on the ambient humidity. At an actuation voltage of 1.8 V, the bending tip displacement is of 0.15 mm at 50 % relative humidity and it increases to 0.35 mm at 60 % relative humidity.
To the knowledge of the authors, no numerical Galerkin-based methods can be found in literature for the simulation of EAPap. This motivates the attempt to propose a novel approach for this kind of material. It should be able to model a very thin structure with a slenderness of 1/1000. Additionally, the model must also be capable to render the heterogeneous ionic distribution through the thickness of the actuator.
The purpose of the present contribution is to provide a numerical framework for the computation of the deformation of EAPap in 3D. The proposed model neglects the piezoelectric effect and focuses on the heterogeneous distribution of ions. Since the ionic migration can be regarded as a 1D phenomenon which takes place through a very thin structure, a scaled boundary plate framework as offered in [13,14] is considered. The kinematics of this formulation allows for a scale separation of the thickness from the in-plane coordinates of the plate. Furthermore, the scaled boundary approach can also be applied to shell structures as presented in [17]. In this case, the method is considered for three dimensions by providing a scaling center for the shell. In contrast to the mentioned literature, the numerical approximation is performed by isogeometric analysis (IGA) instead of the finite element method (FEM). Furthermore, the thickness direction is approximated numerically by B-Spline functions. The electromechanic coupling is defined by the concept of a mechanic volume force density which depends directly on the ionic distribution. This concept is adapted from the model of IPMCs, which has a very similar actuation mechanism as EAPap by means of an ionic migration [15,16]. The proposed model is based on the assumption of linear isotropic elasticity for small deformations. It should be pointed out that this is a simplification, since EAPap can be produced with a defined fiber orientation [6] and experimental data show that EAPap exhibits viscous material properties [21,22]. Furthermore, the modeling of large deformations should also be considered in future works. A further assumption is the electromechanic coupling which is defined in a weak sense. Also the time dependent ionic migration in not discussed in this contribution. The ionic charge distribution is defined as static and represents the converged state of the migration. x(θ α , ζ) The purpose of this section is to provide the derivation details of the stiffness matrix of the proposed scaled boundary plate model. To this end, the mapping concept of the plate is introduced in first place. Then, the kinematics are derived and a linear elastic constitutive law is proposed. In order to solve the numerical problem by IGA, the approximation scheme with NURBS and B-Splines is presented. Since the formulation allows for a scale separation, the in-plane direction is approximated with NURBS and the plate thickness direction is approximated with B-Splines. Finally, the principle of virtual work is applied on the previously introduced concepts in order to obtain the stiffness matrix for the numerical implementation.

Scaled boundary plate mapping
The fundamental concept of the geometric description of the scaled boundary plate is given by scaling the plate reference plane in thickness direction. This means that in contrast to the standard scaled boundary method, a scaling center is not defined anymore since it is located in infinity. The reference plane of the plate is defined by a mapping from a parametric space spanned by the vectors θ α (α = 1, 2) into the 3D physical space spanned by an orthonormal basis system e i (i = 1, 2, 3) 1 .
In order to define the normal vector to the reference plane, the tangent vectors need to be calculated by partial derivation of the position vector w.r.t. the parameters, Then, the normal vector to the reference plane can be obtained by the normalized cross product of the tangent vectors, It should be noted that G 3 is a constant vector for a plate geometry which defines the scaling and thickness direction of the plate. Therefore, the plate description is obtained by introducing a scaling parameter ζ in the direction of G 3 , see Fig. 1. Thus, a position vector is defined by means of the scaled boundary plate kinematics by, The vectors G i define a covariant base system in the plate body. Moreover, these vectors also define the rows of the Jacobi matrix, In order to obtain the corresponding contravariant basis system, the matrix J just needs to be inverted. The columns of J −1 define the contravariant basis vectors, Furthermore, the determinant of the Jacobi matrix is introduced as J = det( J). It maps an infinitesimal volume element from the parametric to the physical space by

Scaled boundary plate kinematics and constitutive law
In order to describe the kinematics of the scaled boundary plate, the displacement vector for a material point is introduced as u = [u x , u y , u z ] T . Additionally, the second-order strain tensor is defined by means of Voigt notation as a vector ε = ε x x , ε yy , ε zz , ε xy , ε yz , ε xz T , in which the first three components describe the normal strains while the last three ones define the double shear strains. In the context of linear strain theory, the strain is defined as the symmetric part of the displacement gradient. Thus, Here, the L-operator is introduced as a matrix which contains the partial derivative operators. By means of the scaled boundary plate kinematics, the chain rule can be applied on the L-operator which can then be rewritten as in which the matrices b α and b 3 are defined by, Therefore, the strain which was introduced in (8) can now be rewritten by means of the scaled boundary plate as in which the displacement vector is derived in the in-plane direction α and in the thickness direction ζ . Analogously to the strain, the concept of the Voigt notation is also applicable to the Cauchy stress tensor, thus σ = σ x x , σ yy , σ zz , σ xy , σ yz , σ xz T . Here the first three components describe the normal stresses, while the last three components define the shear stresses. The constitutive law is considered in this work as a linear relation between the stress and the strain, Here C defines the 3D elasticity tensor written as a matrix in Voigt notation. The matrix C is a constant matrix which is defined by the Young's modulus E and the Poisson's ratio ν.

Isogeometric approximation
The approximation of the geometry and the solution field is performed by means of Isogeometric analysis using NURBS as shape functions. NURBS is based on B-Spline functions which are defined recursively, in which p represents the polynomial order of the functions, ξ is the parameter of the function and ξ i is defined as a knot which subdivides the parametric space. The complete set of subdivisions is defined by the knot vector = ξ 1 , . . . , ξ n+ p+1 in which n represents the total number of B-Spline functions. In order to obtain interpolatory functions at the boundaries, the first and the last knots are repeated p + 1 times which is known as an open knot vector. Furthermore, it should be noted that the knot vector defines the mesh of the geometry. In contrast to the FEM, the mesh is defined in the parametric space by subdividing it with knots. A NURBS function for surface geometries is based on weighted B-Splines according to (14) in which w i j are the weights of the functions and n, m represents the total number of functions in each parametric direction. In contrast to B-Splines, NURBS is rational polynomials which allow for the representation of conic sections. A NURBS surface is represented as a linear combination of shape functions. Therefore, the reference plane of the plate, as introduced in (1), can be expressed by In contrast to FEM, the n × m vectors a i j do not represent nodal coordinates but so-called control points. The difference is given by the fact that NURBS in general does not interpolate the control points. Only the boundary control points are interpolated by the application of open knot vectors on the definition of the NURBS functions. In addition, it should also be mentioned that B-Splines and NURBS are generally C p−1 continuous. Therefore, smooth surfaces can be represented by this method. In isogeometric analysis, the solution field is also approximated by B-Splines and NURBS. For this reason, solution fields of higher continuity can be obtained. In the context of the scaled boundary plate, the solution field is discretized by a multiplicative decomposition, in which N is the matrix which contains the NURBS functions used for the reference plane discretization, The approximation of the solution in the thickness direction is represented byû(ζ ). In the present work and by means of IGA, the thickness direction is approximated by B-Spline functions which introduce control points over the thickness of the plate. In this sense, the proposed method differs from the standard scaled boundary plate formulation introduced in [13] [14], in which the displacement in thickness direction is provided by an analytical approach. The main reason is given by the fact, that in this work, heterogeneous volume loads are to be applied over the thickness of the plate. Therefore, the solution field in thickness direction is approximated byû The solutions at the control points are contained in U and cp defines the number of control points in thickness or scaling direction. A graphical representation for the proposed scaled boundary plate formulation is depicted in Fig. 2.
In the scope of the multiplicative decomposition of the solution field, the local derivatives are obtained by

Principle of virtual work
In the present section, the stiffness matrix for the scaled boundary plate element is derived for the global system. In order to obtain the element stiffness matrix, the integration domains just need to be subdivided by means of Isogeometric analysis with the use of the previously introduced knot vectors. According to the principle of virtual work, the virtual strain energy δU is balanced by the virtual work originated by external forces δW ; thus, The external volume load as well as the traction forces can be obtained from the term δW , while the stiffness matrix is contained in δU . The virtual strain energy of the plate is defined by the volume integral of the tensorial product of the virtual strain and the Cauchy stress. In matrix notation it reads, In analogy to (11), the virtual strain is obtained by replacing the displacement with the virtual displacement, The integration domain can be reparametrized by means of the Jacobian, see (7), and the virtual strain (25) can be substituted in (24), Now the numerical approximations of the derivatives of the virtual displacements given in (21) can be substituted in the previous equation and also the integration domain is separated in the scaling and in-plane direction, Next, the constitutive law for the scaled boundary plate given in (12) is inserted and also the derivatives of the displacement are replaced by the numerical approximation given in (20), In order to simplify the previous equation, the B-operators for scaled boundary plate kinematics are introduced. These are given by B 12 = b α N ,α and B 3 = b 3 N and are substituted in the previous equation leading to Furthermore, the in-plane integrals can be rewritten as in-plane stiffness matrices which are defined by the following expressions: It should be noted that these matrices already include the 3D constitutive law of the plate material. Substituting them in (29) simplifies the expression of the virtual strain energy to As observed, the strain energy is now reduced to an integral in scaling direction. At this point, independent approaches for the thickness direction can be chosen. In the present work, a thickness discretization by means of Isogeometric analysis is proposed. This allows for a continuous strain field along the thickness domain if the interpolation is at least of order r = 2. The thickness approximations are given in (22) and can be substituted in the previous equation. Thus, Finally, the global stiffness matrix is given in this expression by the integral in scaling direction, The algorithmic implementation of the scaled boundary stiffness matrix is performed by two subsequent numerical integration procedures. In the first step, the matrices E 0 , E 1 and E 2 are integrated over the reference surface of the plate by considering the NURBS functions which are also used for the discretization of the reference surface geometry. In the second step, the full stiffness matrix is obtained by the thickness integration with an independent approximation approach. In the present work, a B-Spline approximation of order r = 2 is proposed in order to obtain a C 1 -continuous displacement field in thickness direction. Additionally, knot insertion can be performed in order to increase the numerical resolution. This feature is very useful for the consideration of heterogeneous volume force densities along the thickness of the plate which is caused by ionic migration in electroactive paper and IPMC. Furthermore, material heterogeneities can also be modeled by these means. As a last remark, it should also be mentioned that a resorting of the stiffness matrix may be required in dependency of the convention which assigns the degrees of freedom on the mesh. Due to the separate integration scheme, the degrees of freedom are sorted in increasing order along the thickness of the plate.

Numerical examples
In this section, three benchmarks are provided to demonstrate the capabilities of the proposed isogeometric scaled boundary formulation under heterogeneous volume loads. The first example is of academic nature and consists of a cantilever beam which is loaded by a linear volume force density over the thickness. In this example, the numerical result is compared to an analytical one. In the second benchmark, the coupling of the ionic charge distribution to a volume force density is introduced. This concept is taken from the bending mechanism theory of IPMC, and therefore, the formulation is compared against experimental data of an IPMC actuator. In the last example, the coupling mechanism of IPMC is transferred to EAPap. It demonstrates that the proposed concepts are applicable to the modeling of the bending behavior of an EAPap actuator by comparing the numerical results with experimental data.

Cantilever beam under volume moment load
The purpose of the first numerical example is to provide a first check of the proposed plate formulation. To this end, an Euler cantilever beam as depicted in Fig. 3 is considered. The length of the beam in x-direction is given by l, the width in y-direction is b, and the thickness in z-direction is t. Since the cross section of the beam is rectangular, the area moment of inertia is defined by I = bt 3 /12. Additionally, the Young's modulus E is given, while the Poisson's ratio is neglected.
In order to test the implementation for non-constant volume forces across the thickness, a linear force density is defined according to, in which a is the magnitude of the force density at the top and bottom surface of the beam (z = ±t/2). Consequently, a bending moment per length results over the cross section of the beam which is obtained by the integration over the cross section as The bending behavior of the Euler beam under a bending moment load can then be described by a second order differential equation as follows, in which w describes the beam deflection in dependency of the x-coordinate. The equation is solved for w by integrating the equation two times, and the cantilever boundary conditions w(x = 0) = 0 and w (x = 0) = 0 are applied. Additionally, the area moment of inertia I is substituted. The maximal deflection at the right-hand side of the beam (x = l) is then obtained as, Furthermore, the Young's modulus can be prescribed in dependency of the beam thickness by E = 10/t N/m 2 in order to obtain a maximal deflection of w max = 1/10 m independently of the beam thickness. For the numerical simulation, the following thicknesses are considered: t 1 = 1/10 m, t 2 = 1/100 m, t 3 = 1/1000 m. The beam is then discretized by one single isogeometric scaled boundary plate element as introduced in the previous sections. The order of the approximation functions is adjusted to the nature of the beam deflection. This means that the order in x-direction is p = 3, in y-direction q = 1 and in z-direction r = 2. The total amount of degrees of freedom for this discretization is 72. The applied force density is integrated over the thickness of the element by Gaußpoint quadrature in which continuous B-Splines of order r = 2 are used as approximation functions. Furthermore, the Poisson's ratio is neglected, ν = 0.
The results of the numerical examples yield a maximal deflection of w num = 0.1000 m for the three proposed thicknesses. Hence, it can be stated that the implementation is capable to capture a non-linear volume force over the beam thickness up to t 3 = 1/1000 m. The contour plot of the deformation of a beam with thickness t 1 = 1/10 m is depicted in Fig. 4. Additionally, the displacement components in x-and z-direction are presented in Fig. 5.

Ionic polymer metal composite bending actuator
The second example consists in the reproduction of the bending behavior of an IPMC actuator as provided in [15]. The structure of the actuator is defined by different regions as shown in Fig. 6. The main region consists of Nafion, which is an EAP which is able to bend under the influence of an external electric field. The top and bottom regions define the platinum electrodes of the actuator on which an electric potential difference is applied in order to produce an electric field on the Nafion material. Additionally, a diffusion region is included between the Nafion material and the electrodes due to the effects on the material properties of the ionic diffusion. The material properties of the three regions are specified in [15] as follows: E Naf = 200 MPa, ν Naf = 0.49, E Pt = 168 GPa, ν Pt = 0.38, E dif = 84 GPa, ν dif = 0.42.
If an electric potential difference is applied on the electrodes, an ionic diffusion takes place in which counter ions are redistributed over the thickness of the actuator leading to an ionic imbalance. The diffusion process of the counter ions can be modeled by the so called Nernst-Planck equation. The steady-state distribution of counter ions under the application of a potential difference of 2V can be reviewed in [15]. It is observed that a charge imbalance only occurs in the diffusion regions while in the Nafion region the charges remain balanced. Similar to EAPap, the bending actuation is attributed to the ionic migration. In order to couple the ionic redistribution to the mechanical bending, the concept of an in-plane volume force density is introduced, see [15] [16]. The volume force density depends on the charge density according to in which ρ describes the net charge density. The coupling constants are defined in [15] as G = 110 and H = 10 for the IMPC actuator. It is mentioned that these parameters are obtained by fitting simulations to experimental data.
In the present contribution, the net charge density distribution is approximated to the one proposed in [15] by considering a piecewise linear function as shown in Fig. 7. Furthermore, it is also assumed that the integral of the net charge density over the thickness vanishes. This implies that negative and  and q = 1 since the bending occurs along the x-axis. Therefore, the total amount of in-plane control points is 8. In the thickness direction, the order of the B-Spline functions is r = 2. The platinum electrodes as well as the Nafion material are each discretized with one element. The top diffusion region is also discretized with one element since the linear charge imbalance distribution in this region coincides with the thickness of the diffusion region. h − h 2 = 1 10 h = 0.01 mm. The bottom diffusion region is discretized with 8 elements since the linear charge imbalance only covers 1/8 of the diffusion region. h 1 − h = 1 80 h = 1.25 × 10 −3 mm. This means that the total amount of elements along the thickness direction is 12. The discretization in thickness direction is performed according to the isogeometric paradigm by means of knot insertion. Therefore, the shape functions are of higher continuity. Dirichlet boundary conditions are applied by means of the cantilever beam by restricting the displacements for the numerical points at x = 0. Furthermore, a volume force density as described in (41) is applied on the IPMC body. To this end, the charge density distribution as described in Fig. 7 is considered.
The solution of the numerical implementation by means of the proposed isogeometric scaled boundary plate element yields a tip displacement of w num = 11.4337 mm. As a reference, the solution of [15] is considered in which the tip displacement is of w ref ≈ 11 mm. Under this perspective, it can be stated that the proposed isogeometric scaled boundary plate formulation is capable to render the bending behavior of IPMC under the concept of a volume force density which couples the ionic distribution to the bending behavior. Additionally, the contour plot of the deformed IPMC is shown in Fig. 8.

Electroactive paper actuator
The purpose of this last example consists in transferring the volume force density distribution known from IPMC to the bending theory of EAPap. Furthermore, humidity and anisotropy are also factors which have a significant influence on the bending behavior but are neglected in this first modeling attempt. Also the time dependency of the ionic migration is not considered here but a steady-state distribution is assumed.
The experimental data provided in [5] are taken as a reference solution for the proposed formulation. In this work, an EAPap strip with the dimensions of 48 × 14 mm 2 and a thickness of 0.033 mm is considered. It should be noted that in comparison to IPMC, the thickness of EAPap is approx. 10 times smaller. Also in contrast to IPCM, the electrodes consist of a thin gold layer with a thickness of 200 nm which can be neglected for the bending stiffness. The electrodes cover the EAPap material only partially with a dimension of 40 × 10 mm 2 . Regarding the Young's modulus, it is mentioned that it is heterogeneous across the thickness of the actuator and can be averaged with 10 GPa.
For the numerical simulation, only the region which is sandwiched between the electrodes is considered. Therefore, the in-plane discretization consists of one element with the dimensions of 40 × 10 mm 2 and the order of the approximation functions are p = 3 and q = 1. Regarding the thickness direction, the dimension is set to 0.033 mm. Furthermore, quadratic B-Splines are employed in this direction, r = 2. The ionic charge density distribution is applied as for IPMC, see Fig. 7. In contrast to IPMC, a symmetric distribution is assumed. The value of h = 0.0165 mm corresponds to half of the thickness. The domain of the unbalanced region is set to a width of h 1 − h = h − h 2 = 0.08h. The values of the charge density limits are f a = −1.2 × 10 3 Fig. 9 Deformation of EAPap actuator magnified by a factor of 50 mol/m 3 , f b = 0 mol/m 3 , f c = 1.2 × 10 3 mol/m 3 . Since the charge density is subdivided into three piecewise linear functions, the thickness discretization is performed with three elements. As in the previous examples, the isogeometric paradigm is taken into account and thus the discretization is performed by knot insertion which leads to a continuous strain distribution over the thickness.
Regarding the material parameters, the Young's modulus is prescribed with E = 10 GPa and the Poisson's ratio is set to ν = 0.3. The coupling parameters for the volume force density, see (41), are estimated with G = 11 and H = 1.
The numerical result of the cantilever bending simulation yields a tip displacement of w num = 0.1786 mm, see Fig. 9. According to the experimental data provided in [5], a tip displacement of 0.1 mm is obtained for the application of 1 V on the electrodes while the application of 2 V causes a deflection of 0.2 mm. In comparison with experimental results, it can be observed that the numerical simulation predicts a bending tip displacement for an application of a voltage between 1 V and 2 V. It should also be noted that the ionic charge distribution as well as the electromechanic coupling constant is estimated in this example. Nevertheless, the capability of the model shows that it can be applied for the simulation of electroactive paper.

Conclusion
The present contribution proposes a novel approach for the prediction of the bending behavior of EAPap. The numerical framework is based on a FEM scaled boundary plate formulation as given in [13] [14]. The formulation is modified by substituting the FEM formulation with IGA which allows for a higher continuity of the deformation field. In addition, the scaling direction is treated numerically by means of a weak formulation. This allows the consideration of nonlinear effects across the plate thickness as observed in IPMC and EAPap. The electromechanic coupling is performed under the concept of a volume force density which is used for IPMC, see [15] [16]. It is a phenomenological model in which unbalanced ionic charge distributions are coupled to volume forces in longitudinal beam direction in order to create a bending moment.
The most valuable aspect of the proposed formulation lies in the separation of the thickness scale from the in-plane scale. This allows for a separate treatment of the solution in thickness direction, which can be adapted to the phenomenological aspects of a given physical problem. Some possibilities are given by an analytical solution, a collocation method or a weak formulation as proposed in this contribution. Furthermore, nonlinear phenomena can be well captured by providing appropriate approximation schemes. In this work, only piecewise linear functions are considered for the volume forces, but in general, functions of higher order can be considered as well in a straightforward manner.
The numerical results show that the proposed approach is capable to predict the bending behavior of IPMC. It is also shown that the current framework can transfer the IPMC coupling concept to EAPap which have a slenderness in the range of 1/1000.
It should be pointed out that the obtained results are promising but only provide an initial attempt to model the behavior of EAPap by an isogeometric scaled boundary plate formulation. Further enhancement is required to increase the modeling precision. Firstly, the ionic migration needs to be modeled properly by a diffusion equation in which the parameters are provided. It is observed that the bending behavior is very sensitive to the distribution of ionic charges. Also the consideration of humidity represents an important aspect in the actuator bending. Experimental data show that an increase in humidity enhances the ionic migration which leads to larger bending displacements. Furthermore, piezoelectricity and material anisotropy are also aspects which can be considered in the future research of the modeling of EAPap.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.