Data-Driven Anisotropic Biomembrane Simulation Based on the Laplace Stretch

: Data-driven simulations are gaining popularity in mechanics of biomaterials since they do not require explicit form of constitutive relations. Data-driven modeling based on neural networks lacks interpretability. In this study, we propose an interpretable data-driven finite element modeling for hyperelastic materials. This approach employs the Laplace stretch as the strain measure and utilizes response functions to define constitutive equations. To validate the proposed method, we apply it to inflation of anisotropic membranes on the basis of synthetic data for porcine skin represented by Holzapfel-Gasser-Ogden model. Our results demonstrate applicability of the method and show good agreement with reference displacements, although some discrepancies are observed in the stress calculations. Despite these discrepancies, the proposed method demonstrates its potential usefulness for simulation of hyperelastic biomaterials.


Introduction
Material behaviour is of great importance in many fields of engineering, especially for soft materials.Examples include biomedical engineering [1,2] and soft robotics [3].A hyperelastic material model is often used to describe rubber-like materials [4,5] and biological soft tissues [4,6].Hyperelasticity implies the existence of an elastic potential that defines the constitutive equations, i.e., the mechanical behaviour of the material [7].The elastic potential has to obey the objectivity and the symmetries of the material, as well as the polyconvexity of the elastic potential that provides a sufficient condition for the existence of solutions to boundary value problems in hyperelasticity [7].
Dozens of hyperelastic models have been proposed for soft materials [4].Most of these models meet the necessary requirements such as objectivity, material symmetry and polyconvexity through an invariant-based approach [6].This means that a set of invariants is chosen to define the strain measure and the elastic potential is the function of the invariants.It is the common practice to use the invariants of the right/left Cauchy-Green tensors [6].
One of the main advantages of hyperelastic models is that there is no real need to know the exact form of the elastic potential.To define the constitutive equations of a hyperelastic material, we only need to know the derivatives of the elastic potential with respect to the chosen strain measure, so-called response functions [8,9].The response functions may be constructed from experimental data on mechanical behaviour.To characterize nonlinear materials, one test them under a wide range of deformation modes.These include uniaxial, biaxial, shear and sometimes triaxial deformation [10].The resulting data from these tests are stress-strain curves.There are few approaches to constructing response functions from experimental curves.One of them, the most popular, is the phenomenological approach.In other words, the forms of the response functions are suggested a priori as analytical functions with model parameters, the latter are found by fitting to the experimental curves.
The main disadvantages of this approach are the non-uniqueness of the optimal set of model parameters [11] and the lack of any guidance on how to select the form of the function from a variety of expert-constructed models [12], since no universal model has yet been proposed.
Data-driven (model-free) approaches overcome such problems and are becoming more attractive with the development of full-field experimental techniques.Data-driven computational mechanics is an actively developing field (see e.g., [13][14][15][16][17]).One of the main trends in data-driven hyperelasticity is constitutive equations based on neural networks.Input data for neural networks are given strain measures and output data are elastic potential or response functions or some stress measures [13,18,19].There are several neural networks that consider physically based requirements such as polyconvexity [13,19,20].However, a constitutive equation based on neural networks is a black box and such approach raises the question of its non-interpretability [21,22], i.e., transparency of the mechanism by which the model works [23].
In the scope of data-driven hyperelasticity we want to overcome the problem of interpretability.Namely, we want to perform finite element forward simulations for datadriven constitutive equations.To this end, we choose as the strain measure the Laplace stretch based on the QR-decomposition of the deformation gradient [24,25].In this case the response function is equal to the corresponding component of the stress tensor [25].This means that we can reformulate our experimental data explicitly in terms of the response functions and the Laplace stretch, which are nothing else but data-driven constitutive equations.To perform data-driven finite element forward simulations, we use the hyperelastic nodal force version of the P 1 finite element method [26][27][28] with a simple interpolation method for response functions in the 3D space of the Laplace stretch variables that is requested in the iterative solution of the equilibrium problem.For the case of an isotropic material such approach was proposed in [29].However, the advantage of using the Laplace stretch as the strain measure and the response function based approach is its applicability to both isotropic and anisotropic materials.In this study we test our approach for datadriven simulation of anisotropic membranes using the similar steps proposed in [29] for the isotropic material, except for virtual experimental protocols applicable to real-world scenarios and material anisotropy.To this end, we use synthetic data based on virtual experiments for Holzapfel-Gasser-Ogden model [30] with parameters corresponding to porcine skin [31].For virtual experiments we use experimental protocols proposed for soft tissue [32] which are enriched for the sake of data completeness.The obtained experimental synthetic data is used explicitly for the finite element forward simulation of the inflation of a square membrane defined by the data-driven constitutive equation without any prescription about material anisotropy.The corresponding nonlinear algebraic system is solved by a simple relaxation method where an inverse weighted distance interpolation method finds requested values of the response functions at each iteration at any point of the 3D Laplace stretch space.The relative difference between computed and reference displacements is very small, whereas the relative difference between computed and reference stresses may achieve 10% at the membrane corners and diagonals.
It is important to highlight that the utilization of the actual protocols established for cruciform specimens [32] serves multiple purposes.First and foremost, the authors [32] have proposed a "series of mechanical tests designed to maximize inhomogeneous strain fields and in-plane shear forces", which means that it is not necessary to introduce holes in order to achieve a heterogeneous deformation field as it was done previously for isotropic case by [13,29].Furthermore, through the implementation of high resolution imaging and digital image correlation techniques, it becomes possible to capture accurately the complete displacement field.Secondly, the biaxial extensions of the specimens offer a wide range of options that can be employed in real experimental setups.At the same time virtual experiments act as an effective tool to provide additional insights that can aid in improving the actual test protocols in terms of sufficiency of data.Naturally, in real-life scenarios, our knowledge about the mechanics of materials is limited before conducting tests; however, if we can operate within the scope of hyperelastic behavior that is generally accepted for biomaterials [4,6], we believe that utilizing virtual experiments can also contribute to refining the design of experiments.
The main focus of this paper is to demonstrate the viability of our proposed forward data-driven simulation, which utilizes Laplace stretch and response functions, in the absence of any prior knowledge regarding the symmetry of the material.We aim to determine whether the anisotropic deformation of the membrane can be achieved using the constitutive equations provided in the tabulated data.
The outline of the paper is as follows.In Section 2 we briefly introduce the hyperelastic nodal force method, our approach to data-driven constitutive modeling, the virtual experiment for the 2D biaxial extension of a membrane sample.In Sections 3 and 4 we verify our approach for the data-driven forward simulation of square patch inflation.Section 5 summarizes the proposed approach, discusses its advantages and drawbacks, and present future directions of our research.

Equilibrium Equations and Finite Element Discretization
We consider deformation of an incompressible thin hyperelastic membrane with thickness H. Its mid-surface kinematics fully characterizes the deformation of the membrane.Let X and x denote positions of a material point in the reference (initial) Ω 0 and current (deformed) Ω t membrane mid-surface configurations, respectively.The deformation is defined by the mapping x = x(X), the displacement of the mid-surface is u = x − X, and the surface deformation gradient is F 2d = ∂x/∂X.
We use the Laplace stretch ξ = (ξ 1 , ξ 2 , ξ 3 ) as the strain measure for surface deformation [24].The Laplace stretch is based on the QR-decomposition of F 2d and is connected with the right Cauchy-Green deformation C 2d = F T 2d F 2d via Cholesky factorization [25].The form of the constitutive equations depends on the choice of the strain measure.In the considered case of hyperelastic membrane the constitutive equation reduces to the definition of a hyperelastic potential as a function of the Laplace stretch ψ = ψ ξ F T 2d F 2d .For a body at rest in the current configuration we impose mixed boundary conditions on where T is the Cauchy tension tensor, n t is the unit outward normal to ∂Ω t , ū and t are given displacement and tension on corresponding parts of the boundary.
The equilibrium equations for a hyperelastic body can be written in the following variational form [7]: where b is the density of body forces.
For the approximate solution of Equation ( 2) we use the hyperelastic nodal force version of the P 1 finite element method [26][27][28].For the reader's convenience, we briefly review its membrane variant in terms of the Laplace stretch [29].We consider the deformation of a consistent triangulation of the initial planar configuration Ω 0 .A triangle T P with vertices P 1 , P 2 , P 3 deforms into a triangle The Jacobian of the deformation is J = A Q /A P , where A P is the area of the undeformed triangle T P and A Q the area of the deformed triangle T Q .If λ 1 (X), λ 2 (X), λ 3 (X) are the barycentric coordinates of a material point X ∈ T P , then for X and x = x(X) ∈ T Q , it holds (3) Due to the linearity of P 1 finite element basis functions the deformation gradient F 2d and the elastic potential ψ ξ F T 2d F 2d are constant on each triangle T P , and the internal energy is for any point X ∈ T P .The static equilibrium (2) implies a nonlinear system for new positions Q i , i = 1, . . ., K of the mesh nodes where Σ i is a set of triangles sharing the i-th node, F i (T P ) and F i,ext (T P ) are the elastic force and the external force at the i-th node which are computed on triangle T P .
Using the chain rule we can obtain for the hyperelastic nodal force where ∂ψ/∂ξ s are so-called response functions [8,9].The other factors ∂ξ s /∂Q i are defined completely by concise formulas [29].Let the initially flat membrane in flat configuration Ω 0 be deformed so that its current configuration Ω t be defined in the three-dimensional space with basis vectors e 1 , e 2 , e 3 of global (fixed) Cartesian coordinates.Then for each triangle T P one has The shape vectors D i = ∂λ i /∂X are orthogonal to each triangle edge and completely defined by the reference mid-surface geometry [33] with notations P 4 := P 1 , P 5 := P 2 , R := n P × (P i+1 − P i+2 ), n P for the unit normal to the plane of triangle T P , and Due to (5), the response functions (the derivatives of the elastic potential for a given strain measure) fully define the constitutive relationships (mechanical behaviour) of a hyperelastic material.The response functions may be derived from experimental data [8].In case of the Laplace stretch, the derivation can be done explicitly from the deformation curves, since in the 2D case the response function is equal to the corresponding component of the stress tensor [25].It is important to note that Equations ( 6) and (7) do not distinguish between anisotropic and isotropic materials and are valid for any hyperelastic material.
The simplest relaxation method x n+1 = x n + δ • R n is used in our numerical experiments for the solution of the nonlinear system (4), where R n ∈ R 3K is the residual vector for system (4), x n is the vector of mesh node positions at the n-th iteration, δ = 4 × 10 −4 ≪ 1 is the iteration parameter.We note that the solution of nonlinear system (4) may be obtained by the Jacobian-free Newton-Krylov methods [34] as well.

Our Approach to Data-Driven Constitutive Modeling
Experimental data for the response functions ∂ψ/∂ξ s and explicit formulae for ∂ξ s /∂Q i (6) and (7) provide forward simulations of hyperelastic membrane deformation on the basis of data-driven constitutive models.The latter do not require assumptions on exact forms of elastic potentials or response functions.Namely, for forward simulations we use experimentally obtained points (ξ, ∂ψ/∂ξ) which define data-driven constitutive relations via ( 5) and (6).To solve Equation (4), we use the weighted k-nearest neighbors interpolation method [35] for calculating response functions ∂ψ/∂ξ at any requested point (ξ 1 , ξ 2 , ξ 3 ).In our numerical experiments we use k = 10 and weights 1/d 2 , where d is the distance from the requested point to its neighbor.For small deformations (Euclidean norm ||ξ|| < 0.03) we use a linear elasticity model instead of interpolating data from the experimentally obtained points [36].The stiffness tensor components for the linear elasticity model are obtained by the least squares method under the condition of positive definiteness of the matrix corresponding to the stiffness tensor [37].

Virtual Experiments for 2D Membrane Deformation and Their Relevance for Data-Driven Constitutive Modeling
We use synthetic experimental data to test our data-driven approach for anisotropic membrane inflation.Namely, we generate data by virtual experiments on planar extensions of a sample and use it as the input data for data-driven constitutive relations of a hyperelastic membrane without any extra knowledge on sample anisotropy and form of the potential.
Generation of rich data from mechanical testing experiments is one of the main directions in experimental mechanics [42].This is essential for data-driven modelling as well.To design our virtual experiments, we analyzed the protocols [32] proposed for the soft tissue study, in terms of the sufficiency of experimental data for data-driven modelling.The sample geometry and the sizes for the virtual sample correspond to realworld analogs [32]: we stretch the cruciform sample with even arms and collect data only in a small central region, rf. Figure 2. The schematic representation of the protocols is shown in Figure 3, where w i ∈ [0, 1], i ∈ 1, 4 is the percentage of the given maximum displacement u max for the i-th arm: w i = 0 yields the fixed arm and w i = 1 yields the maximum displacement.By varying w i , we can obtain a range of different experiments.In our virtual experiments, we gradually apply displacement with certain step ∆s until the maximum displacement is reached.Displacement w i • n • ∆s is applied to the i-th arm at the n-th step, where n = 1, . . ., N, N = u max /∆s is the number of steps.The triangular mesh for the sample is quasiuniform with mesh size h f it = 0.25 mm, the maximum displacement u max = 2 mm and ∆s = 0.025 mm.At each step we collect (ξ, ∂ψ/∂ξ) for all triangles belonging to the central region.Since we use the linear (P 1 ) finite elements, (ξ, ∂ψ/∂ξ) are constant on each triangle.To analyze sufficiency of experimental data, we studied protocols (Table 1 in [32]) represented by combinations of four values w1 : w2 : w3 : w4 = 1 : 1 : 1 : 1; 1 : 1 : 0 : 0; 1 : 0 : 0 : 0; 0 : 1 : 0 : 1; 1 : 0.67 : 0 : 1 with all possible permutations.In this case, we lack experimental points in the region ξ 2 < ξ 1 , rf. left plots in Figure 4. We varied the weights w i to obtain a more even distribution of points in the experimental data cloud.Our proposed testing protocol assumes five experiments: w1 : w2 : w3 : w4 = 1 : 1 : 1 : 1 w1 : w2 : w3 : w4 = 1 : 0.75 : 1 : 0.75 w1 : w2 : w3 : w4 = 0.75 : 1 : 0.75 : 1 (13) w1 : w2 : w3 : w4 = 1 : 0.5 : 1 : 0.5 w1 : w2 : w3 : w4 = 1 : 1/3 : 1 : 1/3 Figure 4 shows the virtual experimental data consisting of points ξ = (ξ 1 , ξ 2 , ξ 3 ), where we know the values of the response functions ∂ψ/∂ξ = (∂ψ/∂ξ 1 , ∂ψ/∂ξ 2 , ∂ψ/∂ξ 3 ).As one can see, our protocol (13) allows us to distribute the experimental data more evenly.Note that we can always add more data by varying w i .Virtual experiments based on the proposed protocol (13) form the data-driven constitutive equations or InputData for data-driven simulation.In total we obtain 36,472 points in the InputData.

Data-Driven Forward Simulation of Square Patch Inflation
For data-driven simulation we consider inflation of a square patch of a hyperelastic membrane with size R = 10 cm.The constitutive equation is defined by InputData without any a priori knowledge of isotropy or anisotropy.The patch is clamped around its perimeter and inflated by varying pressure p defined via the dimensionless pressure p * = 3pR/(µH), where p * ∈ [0, 12] with increment in 1.
We compare the calculated membrane displacements and the second Piola-Kirchhoff tensors with the displacements and tensors of a reference solution.The reference solution for the HGO membrane was obtained by the hyperelastic nodal force version of the P 1 finite element method described in Section 2.1.The description of the method and its validation are presented in detail in [28].We consider two quasiuniform meshes: a coarse mesh with mesh size h = 4 mm (1448 elements, 776 nodes) and a fine with mesh size h = 1 mm (23,020 elements, 11,711 nodes).
We performed 600,000 iterations of the relaxation method for the fine grid and 150,000 iterations for the coarse grid.Further iterations do not improve the solution accuracy compared to the reference finite element solution.
We estimate the integral relative errors for the displacement u and the second Piola-Kirchhoff stress tensor S by integration over the reference configuration Ω 0 : where || • || 2 denotes Euclidean norm of a vector, || • || F denotes Frobenius norm of a matrix, subscript f denotes the finite-element reference solution, and subscript e denotes the data-driven forward simulation solution.

Results
The distribution of the relative differences in displacements and stresses for the fine mesh and the dimensionless pressure p * = 8 are shown in Figure 5.The Euclidean norm is used for the evaluation of displacements, the Frobenius norm is used for the evaluation of stress tensors.The integral relative errors for both meshes and for various dimensionless pressures are shown in Figure 6.According to Figure 5, the relative error in displacements is very small, less than 0.02% almost everywhere.The error in stresses is less than 5% for the most part of the patch, it can be higher than 10% in the corners of the patch and on its diagonals.This is due to the fact that for corner boundary triangles there is a need to obtain stress values at points ξ outside the convex hull of the data cloud, and data extrapolation results in incorrect stresses.The integral relative error of the data-driven approach is almost independent of the mesh size of the grid providing the calculations.As p * approaches its maximum, the integral relative errors begin to grow both for displacements and stresses.

Discussion
In the paper we test the applicability of the data-driven forward simulation approach based on the response functions and the Laplace stretch for anisotropic membranes.To the best of our knowledge, this approach has been used for anisotropic materials for the first time.To build the data-driven constitutive equation, we performed virtual experiments based on real protocols for cruciform sample extensions of biomaterials [32].The datadriven constitutive equation is provided by a set of experimental points (Laplace stretch, response functions) of some constitutive manifold.The sample material for the virtual experiments corresponds to porcine skin characterized by the Holzapfel-Gasser-Ogden model.
Since data-driven approaches need feasible data distribution, we analyzed experimental protocols in terms of data sufficiency and proposed new experimental protocols that give more uniform coverage of the Laplace stretch space (ξ 1 , ξ 2 , ξ 3 ) than the protocols proposed in [32].In particular, the use of a non-equally armed cruciform specimen should result in richer experimental data due to heterogeneous strain field.It is important to take note that the proposed protocols for virtual experiments are essentially shortened versions of the protocols used for studying soft tissue [32], with the addition of non-equibiaxial extension protocols.This means that the feasibility, potential limitations, and applicability of the proposed protocols are exactly the same as those for real experimental protocols used in studying biomaterials.When it comes to real experimental protocols for biomaterials, the main challenge lies in estimating the corresponding stress field.The majority of studies have relied on the assumption of a homogeneous stress field in regions of homogeneous deformation [43].However, biomaterials are typically heterogeneous, and it is more preferable to incorporate as much experimental information as possible when constructing constitutive equations.Hence, it is important to develop a reliable technique for evaluating the full stress fields, taking into account their heterogeneity.One possible method to tackle this challenge is by assuming a hyperelastic potential and utilizing inverse techniques, such as finite element model updating or the virtual field method [42].However, as mentioned earlier, there is currently no clear guidance on selecting the appropriate form of the hyperelastic potential from a range of established models proposed by experts, and there is also no universal model that has been put forth.Furthermore, there exist inverse techniques that allow for the evaluation of the heterogeneous stress field based on given experimental protocols [44].These techniques are based on finding an approximate solution to the weak form of the Cauchy stress balance equation for a given set of displacements derived from experiments.
The important feature of our approach is the use of a linear elastic model for small deformations ||ξ|| < 0.03.The elastic model provides us correct extrapolation of experimental data to the region of compression deformations which are not available in our virtual extension experiment [36].
We simulate membrane inflation with the data-driven constitutive relation without any knowledge of material symmetry (isotropy/anisotropy) or the form of the elastic potential.The data-driven constitutive equation in our framework is a set of known response functions (∂ψ/∂ξ 1 , ∂ψ/∂ξ 2 , ∂ψ/∂ξ 3 ) at points (ξ 1 , ξ 2 , ξ 3 ) and an interpolation scheme.Note that we do not introduce any additional invariants to describe the anisotropy, as this is done according to the popular and elegant invariant-based approach [6].The three response functions can provide a description of the mechanical behaviour of any hyperelastic material in two-dimensional setting.To interpolate the response functions when solving finite-element equilibrium equations, we use the simple weighted k-nearest neighbors interpolation.Any other interpolation method can be utilized as well.
According to our results, data-driven constitutive modeling is successful in finding displacements for the square membrane inflation.The relative error in displacements is very small, however, the relative error in stresses can be higher than 10% in the corners of the patch and on its diagonals because of using data extrapolation instead of data interpolation: we do not have enough information on ξ 3 in the vicinity of the corner points.Here we faced with the well-known issue of data-driven constitutive modeling that insufficiency of experimental data require extrapolation with large errors.Possible remedies are expanding experimental data and/or exploiting additional relationships for response functions such as material symmetry.Also, the errors may be caused by our simple interpolation/extrapolation scheme.More elaborated method should be used.
In the future, we will delve into the topic of piecewise linear interpolation utilizing Delaunay triangulation as a means of analysis.However, the primary concern that persists is the issue of extrapolation.To tackle this problem, there are a couple of strategies that can be employed.One such strategy involves supplementing the existing experimental data with novel protocols in order to broaden the dataset.Another approach entails imposing additional constraints on the response functions, thereby enhancing the accuracy of the results.As an illustration, when dealing with isotropic materials, it is crucial to ensure that the isotropic condition is satisfied in terms of both the left Cauchy-Green deformation tensor and the stress tensor [45].This isotropic condition can be rewritten in terms of response functions and the Laplace stretch.In the case of anisotropic materials, similar conditions on response functions will be explored and discussed further in the future.
Nevertheless, the results are promising for the application of the proposed data-driven approach to simulation of anisotropic biomaterials.The future work should be focused on better interpolation scheme; avoiding data extrapolation via virtual experiments that help to select appropriate protocols for data-driven simulation; account the polyconvexity constraint on the second derivative of the elastic potential; account the material symmetry in terms of response functions.
In addition to ensuring that our experiments are conducted under realistic conditions, our research will involve an investigation into the technique of heterogeneous stress field estimation put forth by [44].To gain a more comprehensive understanding, we plan to utilize various benchmark materials to estimate the errors that may arise from our interpolation and extrapolation methods.Furthermore, in order to validate and verify the effectiveness of our methodology, a wide range of deformation modes will be examined through actual experiments.It is worth mentioning that our proposed approach will not be limited to two-dimensional scenarios; rather, we will also extend it to encompass three-dimensional cases.

Figure 2 .
Figure 2. The dimensions of the cruciform biomaterial sample.The radius of the cutout is the same for all the cutouts.The central 5 × 5 mm 2 region for data collection is marked by gray lines.

Figure 3 .
Figure 3.The experimental protocols of cruciform sample extensions.For each arm with index 1, 2, 3, 4 we introduce displacement ratio w i , i = 1, . . ., 4. The displacement ratio is a percentage of the given maximum displacement.The blue line indicates the direction of the material anisotropy.

Figure 4 .
Figure 4.The left plots demonstrate experimental data based on the protocol [32].The right plots demonstrate experimental data based on our protocol (13).The top plots show the distribution of points in (ξ 1 , ξ 2 , ξ 3 ) space.The bottom plots show the corresponding projection of the experimental data onto the plane ξ 3 = 0.

Figure 5 .
Figure 5. Distribution of relative differences in displacements and stresses for the fine mesh and the dimensionless pressure p * = 8.The Euclidean norm is used for the evaluation of displacements, the Frobenius norm is used for the evaluation of stress tensors.The subscript f denotes the finite element reference solution and the subscript e denotes the data-driven forward simulation solution.

Figure 6 .
Figure 6.Integral relative errors (15) for two meshes and various dimensionless pressures.The subscript f denotes the finite element reference solution and the subscript e denotes the data-driven forward simulation solution.