Finite strain porohyperelasticity: An asymptotic multiscale ALE-FSI approach supported by ANNs

The governing equations and numerical solution strategy to solve porohyperelastic problems as multiscale multiphysics media are provided in this contribution. The problem starts from formulating and non-dimensionalising a Fluid-Solid Interaction (FSI) problem using Arbitrary Lagrangian-Eulerian (ALE) technique at the pore level. The resultant ALE-FSI coupled systems of PDEs are expanded and analysed using the asymptotic homogenisation technique which yields three partially novel systems of PDEs, one governing the macroscopic/effective problem supplied by two microscale problems (fluid and solid). The latter two provide the microscopic response fields whose average value is required in real-time/online form to determine the macroscale response. This is possible efficiently by training an Artificial Neural Network (ANN) as a surrogate for the Direct Numerical Solution (DNS) of the microscale solid problem. The present methodology allows for solving finite strain (multiscale) porohyperelastic problems accurately using the direct derivative of the strain energy, for the first time. Furthermore, a simple real-time output density check is introduced to achieve an optimal and reliable training dataset from DNS. A Representative Volume Element (RVE) is adopted which is followed by performing a microscale (RVE) sensitivity analysis and a multiscale confined consolidation simulation showing the importance of employing the present method when dealing with finite strain poroelastic/porohyperelastic problems.


Introduction
Porohyperelastic materials present multi-scale biphasic problems in which a hyperelastic porous medium interacts with viscous fluid filling and percolating through its pores.
For decades, the nonlinear multi-scale phenomena including heterogeneous elasticity and poroelastic/porohyperelastic problems have been under investigation (especially, from the theoretical point of view) [1,2,3,4,5]. However, a robust multi-scale multi-physics methodology for solving finite strain poro-hyperelasticity (from theory to general computational mechanics solutions at both scales) was missing since the algorithms are multiplex requiring multiple domains of knowledge/expertise and great computational power. In [6] the ALE description is used to achieve a linear poroelastic model with infinitesimal porescale deformation assumption. This was not solved at the homogenised level although it was simplified by neglecting the effects of macroscopic deformation on RVE mechanical response. In [7] this problem is approached theoretically using asymptotic homogenisation of the Lagrangian description of the Fluid-Solid Interaction (FSI) at the microscale. However, since the fluid phase can introduce extremely large "deformations", this approach might be suffering from inaccuracy and numerical instability. Although the FSI problems in the Lagrangian setting have been addressed using advanced techniques (including mesh-moving and mixed-hybrid velocity-based formulations [8,9,10]) the complexity imposed by the multi-scale nature of porohyperelastic problems does not allow for embracing Lagrangian techniques.
In general, the structure of the multi-scale ALE-FSI approach includes the formulation of the FSI problem at the physical scale (before scales decoupling) using the ALE description shown, schematically, in Figure 1. This is followed by the multi-scale analysis of the achieved equations (ALE-FSI) using asymptotic (twoscales) homogenisation which results in three systems of PDEs for the homogenised (macroscale) problem, solid cell (RVE) problem, and fluid cell problem. The microscale solid response is required at every quadrature point and increment to determine the macroscale response, which is very time-consuming. Thus, we construct an Artificial Neural Network (ANN) as a surrogate for microscale problems delivering the required results in real-time. The mentioned techniques are explained further in the corresponding sections.
In this study, we employ the ALE method, which is widely embraced for fluid flow with moving boundaries [11,12], multidimensional fluid dynamics problems [13], FSI problems [14], etc.
The ALE method uses an arbitrary computational mesh that could be moved in any prescribed manner (including Lagrangian) or could be held fixed (Eulerian). We employ this approach since the fluid is flowing within a domain that is deforming/moving due to solid deformation. The arbitrary mesh is chosen to move with material (Lagrangian) in the solid skeleton and on the interface while it is held fixed within the fluid domain (except for the fluid-solid interface). The balance of mass and momentum equations of the fluid are also mapped to the reference domain (e.g. the fluid domain with undeformed solid phase) using the stress continuity on the interface (using Piola transformation), which is employed to establish the effective balance of linear momentum. This approach provides an economical, accurate, and numerically stable approximation of the model response as it does not directly include fluid displacement or fluid deformation. Since the problem has multiple length scales the ALE formulation of the FSI problem is non-dimensionalised to avoid dimension-related ambiguities. Furthermore, for the first time, we formulate this type of mathematical multi-scale problem directly using the firstorder derivative of a hyperelastic strain energy density function (neo-Hookean) with respect to the deformation gradient tensor (no explicit linearisation of the system) which again results in higher accuracy and numerical stability.
Assuming a sharp length scale separation (between micro and macro scales) and initial local periodicity, the problem can be regularised, thus standard homogenisation approaches can be adopted [15]. Here, the resulting dimensionless ALE-FSI formulation is analysed using two-scale asymptotic homogenisation techniques [6,7,16,17]. This includes differential operator decoupling and power series representation of the relevant fields which leads to the expanded form of the equations. Analysing the latter, for the macroscale, the effective stress, balance of linear momentum, mass conservation, and transformed form of Darcy's law are accompanied by the corresponding microscale hyperelastic solid and fluid RVE/cell systems of PDEs. This is followed by developing the weak formulation and the incremental analysis of the problem via the Finite Element (FE) method.
The general numerical strategy for solving the problem is similar to the FE2 method for composites (see e.g. [18,19]) in the sense that for each time increment and each macroscale numerical quadrature the responses of microscopic problems are required. In FE2, the volume average of stress is determined at the microscale (without any constitutive law at the macro level) while, here, we calculate the average microscale displacement gradient tensor using RVE problems at the microscale. The latter is supplied to the macroscale constitutive law to determine the effective stress. This has the advantage of constructing the deformation gradient and Piola transformation tensors directly at the macroscale, allowing us to solve the mass conservation equation and calculate the transformed pore pressure for the determination of the effective stress. This type of online calculation could be computationally very expensive and, in case of larger problems, infeasible. However, if we can obtain the average displacement tensor differently there is no need for solving for the full displacement field as part of the effective system of PDEs. This can be achieved with considerable speed-up by exploiting the predictive power of ANNs in the context of computational mechanics [20,21,22,23].
When constructing a suitable ANN to serve as a surrogate for the microscale system of PDEs (a continuous function which is relatively time-consuming to be solved) we consider the macroscale displacement gradient and pore pressure as the inputs while the average microscale displacement gradient tensor is the output. The ANN needs to be trained to deliver accurate and reliable outputs. The training procedure consists of tuning the ANNs parameters (weights and biases) by minimising a cost function. The latter represents the distance between the ANN outputs and the "exact" values of a sufficient number of samples that are provided in the training dataset. In fact, the final value of the minimised cost function (here, obtained using Adam optimiser [24]), shows how accurate the ANN represents the features introduced by the training dataset. However, in order to reach the desired fidelity, one should ensure that the training dataset also reflects the features of the original function. To this end, in this study, a simple real-time output density check algorithm is presented providing an optimal reliable density of the training dataset. Furthermore, to demonstrate the microscale response under a variety of possible macroscale displacement gradient components and pore pressure, we perform a sensitivity analysis on RVE problems and study the local phenomena (local concentration of strain energy density) in detail.
The above approach is implemented using the open-source package FEniCS [25] and the opensource Machine Learning package Pytorch [26] for ANNs. A confined consolidation/compression test is used to verify the numerical implementation and to gain a deeper understanding of the importance and robustness of the present method for finite strain porohyperelastic problems. The results of the numerical examples show that the maximum strain energy density and deformation are observed at the intersection of the pores with values considerably above the average. The deformation predicted by the present method is smaller than the one obtained using the linear poroelastic method (due to the strain stiffening feature of the neo-Hookean model). The steady-state is achieved in a considerably shorter time, and considerable variations in tangent properties in time and space can be observed. Finally, the non-dimensional variables are dimensionalised using the characteristic values of two scenarios of interest, namely, brain tissue and soil mechanics.
The ALE formulation together with asymptotic/two-scale homogenisation and the final macroscale and microscale governing PDEs are provided in Section 2. The weak formulation, incremental analysis and data-driven approach together with sensitivity analysis of the RVE response are carried out in Section 3. Section 4 presents a benchmark finite strain porohyperelastic problem (a consolidation test) using the present method. Finally, Section 5 provides the concluding remarks of the present work followed by some potential future directions.

Mathematical description
Let us assume a poroelastic body in physical (single) scale (domain Ω) which consists of a subdomain of hyperelastic porous skeleton Ω s (reference solid domain) and the complementary viscous fluid Ω f (reference fluid domain) filling and percolating the pores such that Ω = Ω s ∪ Ω f . In this section, the mathematical description of the ALE-FSI problem at the pore level is followed by asymptotic homogenisation/multi-scale analysis resulting in the effective governing system of PDEs and the corresponding solid and fluid RVE/cell problems.

ALE-FSI governing equations
The Finite strain solid problem is described using the Lagrangian method in which the observer follows the material movement. The first Piola-Kirchhoff stress (PK1) is employed which is a two-point tensor consistent with the Lagrangian displacement gradient. Furthermore, PK1 is calculated using the deformation gradient tensor which does not involve higher order displacement gradients. This renders it ideal for asymptotic homogenisation.
Since the fluid domain obeys the solid deformation (moving boundary/domain problem) the Eulerian method with the stationary coordinate system (which is usually employed in fluid mechanics) leads to inaccuracies. Due to the latter condition, (fluid flow with moving boundaries) we employ the ALE method which includes a mesh that follows the moving boundaries [11]. The ALE description is chosen to follow the material/Lagrangian description on the interface and within the solid domain while it uses the Eulerian/spatial description elsewhere in the fluid domain. Thus, one can formulate the fluid problem in the current configuration but within the reference fluid domain by pulling back the interface from the deformed configuration to the reference one using the ALE method. The fluid-solid coupling is provided by the interface transformation from the deformed configuration to the reference one using Nanson's formula.

Solid
The spatial variablesX andx refer to reference/undeformed and current/deformed configurations in physical scale (before multi-scale considerations), respectively. The solid problem is described using the balance of linear momentum where P = ∂Ψ ∂F in Ω s (2) and are the first Piola-Kirchhoff (PK1) stress and the deformation gradient tensor, respectively. Furthermore, Ψ and u indicate solid strain energy density function and displacement vector, respectively.

Interface and Fluid
Neglecting the inertial forces and convective term due to the standard fluid velocity scaling (explained in Section 2.2) the ALE form of the fluid problem description reads [11,12] The Cauchy fluid stress σ is defined as where µ f , v, and p are the fluid dynamic viscosity, velocity and pressure. In order to provide accurate interface conditions, the transformation of the interface area from deformed/current configuration to reference configuration is necessary. The latter is possible using Nanson's formula as follows where n and N indicate the normal vectors to the interface surface in deformed and reference configurations, respectively. The second order Piola transformation tensor G is Thus, the stress continuity on the interface is established as Furthermore, Nanson's formula together with divergence theorem provide a coordinate transformation for arbitrary vector V and tensor T Thus, Equations (4) and (5) can be pulled back to a fixed reference fluid domain as to be used in establishing the effective balance of linear momentum.

Non-dimensionalisation
The problem at hand involves two or more length scales that vary in an application-specific range (e.g. biomechanics, soil and rock mechanics). Moreover, due to the multi-physics nature of the matter, there is a complex interdependency between the units of measurement. Thus, the dimensions can become a potential source of ambiguity leading to a misunderstanding of the model parameters and response [27]. To address this issue, we consider the non-dimensionalisation procedure before proceeding to the multi-scale formulation. The non-dimensionalisation is similar to the ones in [22,28,29] but extended for ALE formulation and provided here for clarification and the reader's convenience. Let us define X and Y denoting, respectively, the formally independent macroscale and microscale reference spatial variables, defined as and x and y being the current counterparts. 0 < = d L 1 is the scale separation factor and d and L are, respectively, microscale (RVE) and macroscale representative lengths. We define Springer Nature 2021 L A T E X template two formally independent representative values for fluid dynamic viscosity µ f c and force f c such that where f is (any) force, and µ f is the interstitial fluid dynamic viscosity. All other parameters in this study are non-dimensionalised with respect to the mentioned independent characteristic values e.g.
where µ s and λ s are the solid Lamé constants (to be used in the solid constitutive law), ν is the solid Poisson's ratio, and r indicates the non-dimensional parameter. Substituting the non-dimensional fields in the governing equations and dropping the prime symbol, which indicates non-dimensionalised variables (for the sake of simplicity), all the equations will have the same form except for Equation (6) which, factorising and removing fc L 2 , takes the form [28,16] The fluid convective term and inertial effects are also scaled by ) and neglected from the governing equations.

Multi-scale expansion
Having achieved the dimensionless ALE-FSI governing equations, we can now apply asymptotic two-scale analysis which starts from the multiscale expansion of the relevant fields using the following steps. The differential operator is decoupled via 19) and the fields (ψ) are represented via power series We use Equation (19) and Equation (20) to achieve the expanded form of all relevant fields, which is provided in Appendix A. The integral average is defined as while the average can be obtained viā The integral average and average are applied to the parameters and fields with y-dependency to be employed at the homogenised level. We carry out the multi-scale expansion of the governing equations in the following.

Expansion of governing equations
Multi-scale expansion of the governing equations is provided here and will be served as the basis for upscaling and to achieve the RVE problems. Equation (1) could be written as Throughout this study, the coefficients of −1 , 0 , and 1 are referred to as C( −1 ), C(1), and C( 1 ), respectively. Equation (23) in C(1), C( ) and C( −1 ) reads, respectively, Expanding Equation (18) results in The expansion of Equation (4) renders Similarly, the mass conservation Equation (5) takes the form Substituting Equation (27) into (30) highlights that the pore pressure is constant with respect to the microscale space ( From multi-scale expansion of Equation (9) in C(1), C( ) we reach Having the expanded form of the governing equations, we can proceed to the upscaling process.

Upscaling
In this section, the PDEs governing the homogenised domain (Ω h ) are derived from the expanded governing equations. This is followed by establishing suitable cell problems whose solutions provide the required parameters of the homogenised system of PDEs.

Balance of linear momentum
From the sum of Equations (1) and (13) we have which, after the expansion, takes the form Applying the divergence theorem (considering the direction of the normal vector) we reach Considering Equations (27) and (36) we conclude where V s and V f indicate, respectively, volume of solid and fluid phases (|Ω s | and |Ω f |) in the reference configuration. Considering that p (0) and ∇ X u (0) are spatially constant in microscale space (y) and averaging the micro-dependent parts of the equations i.e.
the macroscale balance of linear momentum is where Ω h indicates the homogenised domain and We take the same approach as in [7], where the zeroth order of expansion of Equation (1) takes the form thus the constitutive law for solid stress reads where, using Equations (21) and (22),

Mass conservation and Darcy's law
Expanding Equation (14) reads which, applying the divergence theorem and considering local periodicity in the reference configuration, could be rewritten as Considering the compatibility condition on the interface and assuming no-slip interface condition (v =u on Γ) we can rewrite Equation (49) as which, again, using the divergence theorem takes the form Using Equation (21), Equation (52) can be written as which can be rewritten as where yet, this equation should be further processed to be employed in the macroscale system of equations.
Using the divergence identity ∇ · (Av) = v · ∇ · A T + A T :∇v for arbitrary tensor A and vector v, the last two terms of Equation (54) can be written as considering ∇ Yu (0) = 0 and Equation (109), the first and the last terms of the right hand side are zero. Furthermore, applying the divergence theorem twice, we Finally, substituting Equation (57) into Equation (54) and using the averaged ydependent quantities, the effective mass conservation takes the form which, at infinitesimal strains, holds the analogy with the infinitesimal poroelasticity in [30]. Darcy's law [31] is adopted as the standard Ansatz/constitutive law for the effective relative fluid velocity Thus, The term ∇ Yu (1) and hydraulic conductivity K are to be determined, respectively, using the RVE problems in solid and fluid domains (solid cell problem and fluid cell problem), which are provided in Section 2.5. We highlight that the term ∇ Yu (1) , implicitly, provides the fluid-solid coupling term.

RVE problems
In this section, the systems of PDEs to be solved in RVE domains to determine the micro-driven parameters (i.e. ∇ Y u (1) and hydraulic conductivity K) are provided.

Fluid phase
Since the fluid is assumed to be incompressible and Newtonian flowing at sufficiently small velocities making the inertial convective effects negligible, it is a linear problem. Thus, instead of solving the fluid problem at each pressure gradient, we calculate its tangent which provides the hydraulic conductivity for Darcy's law as in [16,28,30,32]. Substituting Equations (27) and (28) into (29) Considering the RVE geometry, local periodicity, no-slip interface condition, Equations (33) and (55) We consider the following Ansatz to determine the hydraulic conductivity which results in where the hydraulic conductivity of Darcy's law could be calculated via with the uniqueness conditions The solution of the fluid RVE problem can be obtained by concatenating the results of a componentwise analysis as in [30].

Solid phase
Equation (26) and (35), together with periodic boundary condition, construct a solid RVE problem as follows with the constitutive law Equation (45), repeated here for the reader's convenience, and At the microscale level (RVE problems), the macroscopic variables p (0) and ∇ X u (0) are knowns to impose a condition under which the microscopic response ∇ Y u (1) is to be determined by solving the RVE problem in the solid phase.

Incremental weak form supported by ANNs
In this section, we provide the formulation and details required for solving the sets of governing equations numerically. Basically, the macroscale problem is solved numerically via FEM which requires the weak form of the corresponding equations. Furthermore, due to fluid-structure interactions (multi-physics), the response is pathdependent which should be solved incrementally. The microscopic response is required to solve the macroscopic equations due to the multi-scale nature of the problem. This is included via the term ∇ Y u (1) and hydraulic conductivity in the effective governing equations (see Equation (47), (58) and (60)). The latter is a one-time calculation per initial properties (since the fluid model is linear) while the former (∇ Y u (1) ) is to be determined per initial properties and within each increment for each finite element since the response is nonlinearly correlated to the macroscopic inputs (∇ X u (0) and p (0) ). Such a high number of direct numerical simulations (DNS) of the microscopic RVE problems for one macroscopic problem renders it computationally unaffordable. This is compounded by the need for input-output tangents in the macroscopic iterative solver (e.g. Newton-Raphson). This problem is addressed by employing ANNs trained with numerical results of microscopic RVE problems in fluid and solid phases via DNS.
In this section, firstly, the procedures of solving RVE problems via FEM are explained and some RVE numerical tests are provided, secondly, the data-driven approach (as a surrogate for the RVE problems) is introduced and, finally, the solving strategy for the macroscale problem is presented.

RVE fluid problems
Following the componentwise analysis of the system of Equations (67)-(69), the final system of equations to be solved in the fluid domain reads for every i = 1, 2, 3, constructing three Stoke's problems each with a unit body force in the direction of i-th unit vector of the chosen system of coordinates (i.e. e i ). Consequently, the tensor K = K ji f can be constructed by concatenating the vectors K i f .
Finally, the weak form of Equations (76)-(78) reads where, δv and δq are arbitrary test functions. We highlight that, since the RVE fluid problem is linear, incremental analysis is not required. To accurately solve the fluid problems, the corresponding domain is discretised using 57273 tetrahedral elements.

RVE solid problem
The RVE solid problem described by Equations (72) is a non-symmetric problem with a relatively complex geometry of RVE. Solving the problem incrementally (applying BCs in small increments) improves the numerical stability, allowing us to consider larger macroscale deformations and porepressure and providing a dense dataset for ANN training. The incremental weak formulation reads Note that ∇ X u (0) and p (0) are y-constants given by macroscale problem that are imposed incrementally. δu (1) and u (1) are, respectively, the test function and the function (solution) of this problem. The output to be provided for macroscale problem is ∇ Y u (1) s . Furthermore, we assume the compressible neo-Hookean model as the leading order strain energy density function as follows where J (0) is defined via Equation (103) and I (0) 1 = Tr(F (0)TF(0) ). The constants µ s and λ s are the material parameters (the Lamé constants). The leading order strain energy density function (Equation (82)) can be non-dimensionalised using fields in (17) as which, dropping the prime of variables (for simplicity and consistency with the rest of the formulation), takes the form We highlight that although the unit of strain energy is f c L the strain energy density (strain energy per unit volume) has the same unit of stress ( fc L 2 ).

Local periodicity
Local periodicity is required within each finite element of the macroscale model at the Reference Lagrangian configuration. It allows us to exploit the potential of asymptotic homogenisation at each element. Since the RVEs in a neighbourhood within each finite element receive the same yconstant macroscale inputs (∇ X u (0) and p (0) ) the initial local periodicity is maintained throughout the problem.

RVE numerical tests
Localisation, which is one of the important points provided by the present methodology, allows us to compute the mechanical response field of the microstructure under arbitrary macroscopic deformation and porepressure. This can be served to study phenomena such as local strain energy density (or deformation/stress) distribution and concentration (e.g. in the context of local fracture/damage mechanism). Here, we consider a unit cube as microstructure RVE. The solid phase (a) Displacement (magnitude of u (1) ) field of the RVE due to the pore pressure (p (0) ) in the absence of any macroscopic deformation (∇ X u (0) ). The higher displacement magnitude at the intersection of the three cylinders is notable.
(b) The strain energy density field of the RVE corresponds to the same BCs as in Figure 2a. The energy concentration at the intersection of the fluid channels is of utmost importance, particularly, in the context of fracture/damage analysis.     is the cube from which three cylinders (the pores or fluid channels) with radius r f = 0.2 are subtracted resulting in a porosity of φ = 0.29 (for graphical representation of the RVE geometry see Figure 2b). The material model is the compressible neo-Hookean as in Equation (84) 12 ) macroscopic displacement gradients, as well as pore pressure (p (0) ). Since the latter conditions have strong interdependencies, we embrace the one-factor-at-a-time (OAT) sensitivity analysis method (varying only one macroscopic condition while the others are absent/zero) to study the role of each case.
As expected, the microscopic response is spatially heterogeneous with regions of high concentration. The colormaps of displacement magnitude and strain energy density shown in Figures 2,  3, 4 visualises the local concentration zones of each field. In general, there are nine displacement gradient elements yet we only provide the profile of the ones with considerable effects on the overall response of the medium. The latter is divided into "Corresponding elements" and "Non-corresponding elements" which are shown in Figure 5 (compare the indices of I with O). Figures 5 and 6 show, respectively, the variations of the average and maximum displacement gradient elements. Figure 7 provides the average and maximum strain energy density under different macroscopic conditions with further explanations provided in close captions.
Hydraulic conductivity is a critical parameter in determining the hydraulic response of a poroelastic medium which varies in cases under large porepressure/deformation. In the present methodology, this variation is taken into consideration using a transformation from the deformed configuration to the reference configuration. Considering the term 1 J in Equation (10), which is factorised and removed from mass conservation (in Equation (49) to simplify the mathematical procedure), and neglecting the effects of J (1) in the expanded term 1 J (0) + J (1) , the full form of the transformed hydraulic conductivity takes the form The variations of the main elements of the transformed hydraulic conductivity (normalised with respect to its initial value) due to the prescribed macroscale conditions and the resultant microscale effects are provided in Figure 8. Finally, for the sake of comparison with linear poroelasticity, the profiles of the forth order tensor M, the second order tensor Q, and Biot modulus (which are constants in linear poroelasticity) under ∇ X u (0) and p (0) are provided in Figure 9.

Data-driven approach using ANNs
In this section, we introduce ANNs as a surrogate for the solid cell problems that provide micro-macro scales link. We choose ANNs for this task because of the complexities enforced by the strong interrelation between the pore pressure and macroscale deformation gradient. This interrelation is due to, firstly, the pore pressure Neumann B.C. (introduced via Equation (73)) which depends on the leading order Piola transformation, secondly, the macroscale displacement gradient which results in a body force-type load and, thirdly, a neo-Hookean hyperelastic model which is nonlinear. Thus, a simple superposition of the effects is not valid. The RVE problems solved via DNS are too time-consuming to be included directly in the calculations of the provided concurrent multi-scale multi-physics approach. However, they are useful to provide dataset for training ANNs. Basically, the aim is to replace a continuous problem represented by a system of PDEs (RVE problem to be solved via DNS with a geometry described in Section 3.3) with a surrogate model (here, ANNs) that provides the response faster than the original approach. In general, the inputs of the ANN are the parameters affecting the microscopic response ∇ Y u (1) including geometrical features, material parameters, and macroscale inputs (∇ X u (0) and p (0) ). However, since providing a dataset for the general case is computationally expensive we adopt a case-specific approach in which only the varying parameters are introduced as the inputs of ANNs (see Figure 10 for a schematic graphical representation of the ANN). In this study, for verification/benchmark and comparison purposes, we adopt the relevant varying inputs of a confined consolidation test (e.g. Terzaghi's test), namely, ∂u (0) 2 ∂X2 and p (0) with the outputs being More complex examples will be considered in future studies. The suitable complexity of the model is chosen via ANN Hyperparameters tuning (so-called, Grid Search). In fact, we choose the smallest possible ANN architecture that delivers accurate results to 1-calculate the outputs efficiently and 2-to avoid overfitting due to unnecessarily large ANN capacity. The chosen ANN architecture has three hidden layers with 40 Neurones with Sigmoid activation functions in each layer. Furthermore, it is trained using Adam optimiser and MSE loss function.
Despite their efficiency, the accuracy/fidelity of deep learning approaches as the surrogate models are directly linked to the density (or sample rate) of the discrete sequence (training dataset) so that it reflects the features of the original function. This is provided by solving a sufficient number of cell problems described via Equation (72) and (73). The "sufficient" number of the problems to be solved (samples of training dataset) to provide the required density of the training dataset largely depends on the original problem. In case the latter shows a varying polynomial degree of nonlinearity (which is the case, here) the optimum density will be heterogeneous. Furthermore, covering the space of the high number of input elements at a given interval requires taking all the possible combinations into account. This renders an extreme computational cost if blind sampling techniques (such as equidistant discretisation at the finest required steps) are adopted.
Here, we develop a simple real-time output density check to adaptively refine the output density (adaptive sampling-step refinement). This ensures that the provided dataset reflects the features of the continuous problem at an optimum density so that the well-trained ANN will deliver the desired accuracy. Although the latter is not the focus of this study we briefly explain it. Let us adopt a simple equidistantly sequenced sampling through incremental analysis of the RVE (a) "Corresponding elements" of average microscale displacement gradient tensor.
ii ) under shear macroscale deformation (I = ∇ X u (1) 12 ) show one of the major benefits of using the introduced methodology (robustness and accuracy). We highlight that the latter is neglected using Eulerian/linear poroelastic formulations (e.g. by enforcing zero M ii12 in the calculation of effective stiffness tensor), (a) "Corresponding elements" of maximum microscale displacement gradient tensor.
(b) "Non-corresponding elements" of maximum microscale displacement gradient tensor. solid problem with one moving input element at a time (so-called one-factor-at-a-time OAT). With a default step/increment size of ∆ d i ∇ X u (0) and ∆ p i p (0) we define the residual of sample n as   where || r || F = [ i,j r2 i,j ] 1/2 indicates the Frobenius norm. The term (∇ Y u (1) ) n L is the result of a linear extrapolation which can be described by where ∆I n is the variation in the moving element of input array (∆I n = I n − I n−1 ) and J n−1 is the tangent acquired based on the results of increments n − 2 and n − 1 using Finally, if the residual is more than a pre-defined tolerance the moving input element will adopt the value I n = I n + I n−1 2 .
A training dataset suitable for the numerical reconstruction of a consolidation test (introduced in the following section) with two input variables is to be provided. To this end, we choose relatively coarse initial sampling steps (0.03 for provides a part of FEM Jacobian matrix of the macroscale problem. In linear infinitesimal and remodelling-based finite strain poroelasticity (see e.g. [30,22], respectively), this tensor determines the poroelastic parameters including effective elasticity tensor and Biot coefficient.
The integral average in the deformed configuration is calculated by multiplying the RVE volumetric strain (Ĵ = det(∇ Y u (1) + I)) to the integral average of M. The larger variations in this plot show that solving a poroelastic problem in reference configuration results in considerably smaller changes in the Jacobian matrix rendering the problem numerically more stable.

(c)
The second order tensor also provides a part of FEM Jacobian matrix of the macroscale problem. In this plot, similar to Figures 9a and 9b, it is shown that using ALE formulation the problem will be numerically more stable.  total of 1189, which, given the required time for each problem in Section 3.3, takes 406,44 minutes. The required time to train the ANN can vary, but it is usually around 10 minutes.
The results of the provided ANN is valid for many cases including any characteristic RVE size, solid material constant, and fluid viscosity (since they are non-dimensionalised). However, in cases with different parameters whose variations are not considered in this model (e.g. porosity in [21]) a more comprehensive dataset and network architecture is required. ) as a surrogate model for microscale solid RVE problem. We highlight that this ANN is only valid for a specific geometry of RVE, material model, and material parameters.

Macroscale problem
Equations (43), (44), (46), (58), and (60) construct a system of PDEs describing the macroscopic/effective problem. The weak formulation of the latter reads where δu (0) and δp (0) are the required test functions. Since the problem is history and timedependent we need to perform an incremental analysis which reads where N h is the normal vector to the surface of interest at macroscopic level in reference Lagrangian configuration. In fact, the problem reduces to finding p At this stage, it is necessary to calculate the increments of each field (stress, fluid velocity etc.) based on the increments of u (0) and p (0) , where u t )) employing the trained ANN as described in Section 3.4. As the ANN inputoutput tangent is required for the macroscopic iterative solver (Newton-Raphson) we consider the following transformation for a small increment of microscopic deformation where the fourth rank tensor M and the second rank one Q could be estimated numerically at the end of the last increment using where δ r (except for the test functions) in Equations (93) and (94) indicate a very small variation in r so that taking the partial derivatives in a numerical way. Finally, the term ∆ Ḡ (0) w f could be calculated via

Numerical example
We perform a confined compression (consolidation) test to, firstly, verify the implementation of the methodology and, secondly, to demonstrate to what extent employing the present model will improve the accuracy of the results when dealing with poroelastic/porohyperelastic problems under finite strain. Throughout this section, for the sake of abbreviation, we call the present method and the linear poroelastic equations the ALE and linear cases, respectively. We solve the problem in a non-dimensional form and, at the end of this study, we dimensionalise the variables using the representative/unit values corresponding to soil mechanics and brain tissue applications. Let us assume that the column of poroelastic/poro-hyperelastic material shown in Figure 12 is under mechanical pressure P = P N on Γ top (the external load) where the fluid drainage is allowed. The column is impermeable on all other surfaces (namely, Γ bottom and Γ sides ) with zero displacement B.C. on Γ bottom and zero displacement in abscissa and applicate coordinate axes but free displacement in ordinate axis.
We solve the mentioned problems using the present method with the neo-Hookean material model. We also solve the same problem using the linear poroelastic approach with the effective elasticity tensorC, Biot modulus M , Biot coefficient α, and hydraulic conductivity K identified using the same RVE properties as the ALE model (for more details on multi-scale linear poroelasticity see e.g. [33]). We expect the mechanical responses under smaller deformations to be close while they diverge at higher deformations. The dimensionless fluid dynamic viscosity µ f = 1e−3 is adopted resulting in initial hydraulic conductivity (K i ) jj = 1.943 j = 1, 2, 3. We discretise the homogenised domain with 106 elements (chosen through sensitivity analysis) in the longitudinal direction and one element in the other direction (with 3 nodal degrees of freedom) since the microscopic response ∇ Y u (1) s is 3D although the homogenised response ∇ X u (0) is 1D (see Figure  5b).
We define settlement as the negative displacement in ordinate direction (−u (0) 2 ). In order to ensure that the linearisation of the ANN (Equation (92)) is valid/accurate, we apply the external load in 10 increments at the beginning of the consolidation process with very small changes in time (one-tenth of the usual time increments) which results in a sheer (but incremental) increase in settlement and porepressure (see Figures 13 and 14).
We highlight that the difference between the final/steady-state deformation of the ALE and the linear case reflects the deviation of the neo-Hookean hyperelastic material model from the linear elastic case while the Piola transformation seems to have a strong role in determining the quality of reaching there (with respect to time). Since at small deformations the Piola transformation tensorḠ (0) approaches the identity tensor and the neo-Hookean constitutive equation is approximately the same as the linear elastic one, the solid mechanical response of the ALE and linear cases are expected to be very close which is shown in Figure 13b.  Figure 14a where the pore pressure vanishes considerably faster than the linear case. Concluding from Figures 13, 14, and 15, although the nonlinear material law has a considerable effect in determining the final solid deformation, the effects of the microscale and macroscale solid deformation (which constitutē G (0) andF (0) ) on the fluid part is even more critical. This is also reflected in the spatial profile of hydraulic conductivity at different times, shown in Figure 16, where the transformed hydraulic conductivity varies according to the microscale and macroscale deformations with the initial value being (K i ) jj = 1.943 j = 1, 2, 3. The spatial profiles of the pore pressure, settlement, macroscopic and microscopic displacement (a) As expected, the ALE case with neo-Hookean strain energy density function under finite strain exhibits stiffening behaviour resulting in the smaller final settlement. The settlement is around 30% of the initial length of the column, thus transformation between deformed and undeformed configurations becomes of utmost importance.
(b) The response of the ALE problem under small deformation (close to infinitesimal strain (around 2%)) is very close to the linear poroelastic case, which is a source of verification. Although the final settlements are almost equal the ALE reaches the final settlement slightly earlier. (a) The pore pressure in the ALE case vanishes faster than the linear one indicating a faster transition from the transient state to the steady-state, which could also be concluded from Figure  13a where the ALE steady-state maximum settlement is reached quicker.  gradients at different times are shown, respectively, in Figures 17, 18, 19, and 20b. Comparing Figures 17, 18b, and 19 indicates that the pore pressure results in an increase in hydraulic conductivity in all directions (since it is a hydrostatic pressure) while the macroscopic displacement gradient ∇ X u (0) 22 (which is negative) increases K 22 but decreases K 11 .
The macroscale deformation and pore pressure result in the microscale deformation (represented by ∇ Y u (1) ) whose spatial profiles at different and Ω h ∇ X · w L dV , where w L is the relative fluid velocity of the linear case.  times are provided in Figure 20. The macroscale and microscale displacement gradients constitute the leading order deformation gradient which can be studied via Figure 21. Furthermore, the overall hydraulic response of the medium, namely, the relative fluid velocity of the ALE case is compared with the linear one in Figure 22a showing faster drainage and transient to steady-state transformation. An equivalent hydraulic response in reference configuration can be envisaged using the transformed relative fluid velocity shown in Figure 22b.
Finally, for completeness of the ALE-linear comparison, we plot the micro-macro tangents, namely, the tensors M and Q which are also used in the calculation of the linear poroelastic model parameters. The initial values of these tensors (M 2222 = −0.23, M 2211 = −0.076, and Q ii = −0.113) are used in the linear case. The significant variations in the latter components show the (a) Although shortly after the application of the Neumann BC the pore pressure spatial profile of ALE and linear cases are relatively close to each other they diverge at longer times, mainly due to the faster drainage in the ALE case.
(b) The spatial pore pressure profiles of ALE and linear cases are closer when the applied Neumann BC is small. Fig. 17: The variations of pore pressure along the ordinate direction (Y-axis) at three different times, one shortly after the application of the load, one close to the steady-state, and one in between showing the transition from application of the BC until the final condition. We highlight that at the steady-state of this problem the pore pressure approaches zero everywhere in the space.
importance of considering different sources of nonlinearity including the nonlinear material model, the large deformation, and how their negligence could be misleading. ) (see also [22] and references therein) the dimensional variables can be obtained resembling the brain tissue via

Conclusion
Starting from the ALE formulation of the FSI problem at the pore level, we develop the governing PDEs for homogenised porohyperelastic problems under finite strain. This is accompanied by a standard fluid RVE problem (to determine the hydraulic conductivity) and a hyperelastic solid problem with pore pressure and deformation-dependent Neumann B.C. together with a body force, driven by the macroscale displacement gradient. The latter RVE problem provides, firstly, the field of microscale response to study the local phenomena and, secondly, the average microscale displacement gradient tensor. The latter is required for the homogenised governing PDEs in an online form to determine the macroscopic mechanical and hydraulic responses of the medium. Since solving the solid RVE problem via DNS for every numerical quadrature point in every time increment is very time-consuming, we employ an ANN as a real-time surrogate model for the RVE solid problem. A reliable adaptive sampling algorithm is introduced to provide (a) Agreeing with Figure 13a, the settlement profile of ALE diverges from the linear case. It reaches steady-state considerably sooner (at t = 38) than the linear case. Furthermore, the final settlement of the ALE case is around 20% smaller.
(b) Although the final settlements of both cases under small deformation are, more or less, the same, the settlement of the ALE case, during the transient state, is slightly more which is due to faster fluid drainage.  : Since the compression test is confined, the main element of the macroscale displacement gradient is ∇ X u (0) 22 . We highlight that the other elements are negligible compared to this one.
an optimal training dataset for ANN parameter tuning using the Adam optimisation algorithm.
A sensitivity analysis is carried out by varying the pore pressure and macroscale displacement gradient as the inputs of the RVE problem and the microscopic response as the outputs. A considerable deformation and strain energy density concentration is observed at the interface of the RVE pores with a magnitude much greater than the average values. We show that, considering the fully nonlinear RVE problem, several microscale deformations that are neglected in Eulerian/linear poroelastic formulation can take significant values at large deformations. Furthermore, we show that the transformed hydraulic conductivity varies nonlinearly due to the pore pressure and macroscopic deformation.
Following the incremental weak formulation and numerical implementation of the macroscale problem (which completes the solution cycle/approach), we perform a confined compression/consolidation problem as a proof-of-concept multiscale test for the proposed method. Although this problem results in a uniaxial macroscale displacement gradient tensor we show that the microscale deformation and, consequently, the leading order deformation gradient tensor are three-dimensional. As an implementation verification and comparison means, the same problem is also solved via linear poroelasticity using the initial effective properties derived from the same initial model parameters. As expected from the strain stiffening of the neo-Hookean material model the settlement of the nonlinear ALE case is smaller than the conventional linear poroelastic case. The importance of employing the present method is more evident by studying the hydraulic response of the medium which shows a much faster 33 is similar to ∇ Y u 11 ) showing the combined effects of macroscale displacement gradient and pore pressure on microscale displacement gradient tensor.
transition from transient to steady-state (faster fluid drainage). One reason for the latter is the sizeable increase in the transformed hydraulic conductivity due to the macroscale and microscale deformations. Furthermore, the spatial variations of different parameters with their interactions are studied in detail. Finally, the model response is dimensionalised using the typical/general characteristic values in the brain tissue and soil mechanics applications showing one of the advantages of solving the problems non-dimensionally.
The present methodology is applicable in a wide range of scenarios from biological applications to soil and rock mechanics which provide numerous future directions in higher dimensions. From a numerical point of view, the solution strategy can be improved by translating the ANN into Unified Form Language (UFL) to avoid linearising it within the time increments. Furthermore, the methodology can be extended to include local fracture/damage and path-dependent RVE solid material (e.g. viscoelastic) since the response is determined using a strain energy density function.

A Expansion of general fields
In this section, we apply multi-scale expansion to some basic fields, namely, F, J = det(F), and G plus the expanded form of Nanson's formula that are frequently used in ALE-FSI formulation. Assuming n = 1 which results in two-scale expansion, see Equation (20), the deformation gradient may be expanded as F = F (0) + F (1) = (∇ X + 1 ∇ Y )(u (0) + u (1) ) + I = ∇ X u (0) + ∇ X u (1) + 1 ∇ Y u (0) (a) Although, as shown in Figure 19, the maximum of ∇ X u       . The initial values of these tangent tensors determine the poroelastic properties of the theory of linear poroelasticity such as effective elasticity tensor, Biot coefficient, and Biot modulus. The considerable deviation from the initial values shows the deep effects of nonlinearity in finite strain poroelasticity.