ANN-aided incremental multiscale-remodelling-based finite strain poroelasticity

Mechanical modelling of poroelastic media under finite strain is usually carried out via phenomenological models neglecting complex micro-macro scales interdependency. One reason is that the mathematical two-scale analysis is only straightforward assuming infinitesimal strain theory. Exploiting the potential of ANNs for fast and reliable upscaling and localisation procedures, we propose an incremental numerical approach that considers rearrangement of the cell properties based on its current deformation, which leads to the remodelling of the macroscopic model after each time increment. This computational framework is valid for finite strain and large deformation problems while it ensures infinitesimal strain increments within time steps. The full effects of the interdependency between the properties and response of macro and micro scales are considered for the first time providing more accurate predictive analysis of fluid-saturated porous media which is studied via a numerical consolidation example. Furthermore, the (nonlinear) deviation from Darcy’s law is captured in fluid filtration numerical analyses. Finally, the brain tissue mechanical response under uniaxial cyclic test is simulated and studied.


Introduction
Poroelasticity deals with the mechanical and hydraulic responses of the materials consisting of an elastic porous solid matrix interacting with viscous fluid percolating its pores. This field covers a wide range of applications from soil and rock mechanics to soft biological tissues such as the brain and cancerous ones [1,2]. There are different theories developed for poroelastic problems, including the phenomenological Biot's theory [1] and microstructurederived equations in [3]. The former considers only the homogenised form to be fitted into experiments while the latter is based on asymptotic homogenisation providing analytical relationships, in the form of PDEs to be solved in the cell domain, between the microscale and effective proper-ties. The cell problems, then, can be solved using numerical tools such as Finite Element (FE) method as in [4][5][6]. Representative solutions of cell problems computed this way can provide suitable information to create data sets which can be used in the context of Artificial Intelligence, as done in [7] to address the effective response of heterogeneous poroelastic materials.
One considerable limitation of using asymptotic homogenisation is that the upscaling procedure is only straightforward for infinitesimal deformations which restricts us to the choice of linear poroelasticity (in the sense of linear constitutive equations). On the other hand, the phenomenological approaches such as Biot's poroelastic theory [1] and biphasic mixture theory [8] have been developed to predict the nonlinear poroelastic behaviour of soft tissues [9][10][11][12][13][14]. These theories, however, do not provide sufficient links between micro and macro scales, which considering the complexity and multiscale and multiphysics nature of poroelastic problems, is crucial for a realistic analysis of the scenarios of interest. The fluid flow in poroelastic media is usually described via linear Darcy's law [15], which again is valid in the cases with infinitesimal fluid velocity and solid deformation. Nonlinear fluid flow in poroelastic media has been considered as another challenge which, so far, is addressed via cubic and quadratic corrections (in terms of averaged velocity or Reynolds number) to linear Darcy's law in a rigid porous media [16][17][18][19]. However, not only another parameter (apart from hydraulic conductivity) is to be assumed/determined for the correction term but also the solid deformation does not play any role in the fluid flow nonlinearity using such models.
In this study, we aim at overcoming the mentioned hurdles to achieve a more accurate model describing mechanical and hydraulic behaviour of poroelastic media under large deformation by means of an ANN-aided incremental multiscale-remodelling-based method for numerical analysis of such complex problems. Incremental analysis of structures has long been used in the numerical analysis of "pathdependent" problems such as soil mechanics, plasticity, large (geometric) deformations, etc. [20][21][22][23][24][25]. In such frameworks, three configurations are considered in the deformation history, namely, reference configuration which is the original state at the beginning of the analysis, the current deformed configuration at time t + Δt, and the intermediate configuration at time t just before the deformed configuration. In other words, the intermediate configuration is at the beginning of the time increment, which ends at the deformed configuration. The path from the former to the latter is usually linear, while the transformation from the reference configuration to the deformed configuration can be nonlinear. Although, in linear poroelastic problems, the constitutive equations are linear the nonlinearity of mechanical and hydraulic response arises from solid-fluid interaction. In fact, as a poroelastic problem is generally history/path-dependent, it is usually performed in an incremental numerical framework regardless of being linear [26] or nonlinear [12], however, the full potential of such frameworks have not been exploited. In particular, this framework allows the rearrangement of material properties in each time increment (remodelling) without considerable loss of efficiency if the updated properties are available in real-time. A similar approach is exploited in poroplastic models in [14] in which the structural reorganisation is assumed to obey a phenomenological flow rule driven by stress. In the field of heterogeneous media, the microstructural evolution is captured based on asymptotic homogenisation in [27] which is limited to solving the cell problems at each time increment and spatial point to update the macroscopic properties imposing a high computational cost. A similar framework is also applied in the field of poroelasticity (from a theoretical viewpoint), partially in [28] and fully in [29], where the bulk growth is also addressed. However, due to the high computational complexity imposed by the upscaling and localisation, it might be challenging for the real-world scenarios involving non-linear constitutive response of the constituents.
The key to overcoming this problem could be found in Artificial Intelligence (AI). Data-driven modelling has been successfully applied in a wide range of problems from marketing [30] to computational analysis of mechanics problems. The well known Artificial Neural Networks (ANNs) [31] has been successfully applied in the field of computational mechanics by providing a powerful interpolation mean. This method is inspired by human brain architecture, acting as a transfer function by providing accurate outputs from given inputs. An ANN consists of some hidden layers each adopting a number of neurones (see, Fig. 2 for ANN architecture) which includes a weight and a bias to be tuned during an optimisation procedure called ANN training. The number of hidden layers and the neurones inside each are chosen based on the user's experience. Having a large number of neurones and layers could lead to inefficiency of training and output calculation (so-called feed-forward procedure) while having a small number of neurones can lead to a loss of accuracy. The training procedure, which determines the network's accuracy, is based on an optimisation procedure minimising a distance function from, in the simplest case, a priori provided training dataset consisting of a number of exact outputs with their corresponding inputs. In computational mechanics, ANNs are employed to complement standard approaches such as FE [7] or, in specific cases, to carry out the whole computations [32,33]. In poroelasticity, ANNs are employed for effective model parameter identification in order to efficiently solve complex problems with spatially dependent porosity and solid matrix properties by the present authors in [7]. The training dataset in the latter study was acquired by solving a certain number of cell problems. We employ ANNs for fast computation of both effective model parameters and coefficients required for localisation procedure.
Here, we integrate the mentioned techniques and develop an ANN-informed incremental computational methodology for numerical analysis of poroelastic media under finite deformation considering rearrangement of microstructural properties (porosity and material properties of the solid matrix) and their effects on the effective coefficients of homogenised system of PDEs. The following steps are taken for this purpose: • First, a certain number of cell problems are solved in order to provide a training dataset with inputs being the microscopic properties and the outputs being the tensors required for effective properties calculation as well as localisation procedure. • A suitable ANN is designed and trained via the provided training dataset replacing the time-consuming FE analysis of the cell problems for the calculation of the effective properties of the medium.
• Then, the time domain is discretised, ensuring infinitesimal deformation in each time increment followed by the space discretisation of macroscale domain for FE analysis. • The homogenised system of PDEs obtained via asymptotic homogenisation multiscale analysis for poroelastic media [3] (linear poroelasticity) is employed for numerical analysis of macroscopic response in one increment. • With the help of ANN, this is followed by localisation analysis which calculates the average microscale solid matrix deformation due to the macroscale mechanical and hydraulic response. • The porosity and solid matrix material properties, using a specific strain energy function (here, neo-Hookean), are updated based on the current microscopic deformation (remodelling). • The ANN, again, is employed to carry out upscaling for rearranged effective coefficients identification based on the current (updated) microscopic properties. • The new properties of macroscale system of PDEs are used for the numerical analysis of the next time increment starting from the fourth step.
This methodology builds on the successful integration of FE (for homogenised response approximation) and Artificial Neural Network (ANN) (for upscaling) in [7]. The latter study concerns complex spatially dependent distribution of the underlying microscale properties such as porosity and solid matrix Poisson ratio under infinitesimal strain, while the present method, in addition, addresses the problem of finite strain poroelasticity with ANNs for upscaling, as well as localisation. Here, through the kinematics of the solid matrix, the current configuration is obtained based on which the deformation dependent porosity, as well as linearisation of neo-Hookean strain energy function (remodelling) is carried out and the nonlinear poroelastic (or poro-hyperelastic) model is achieved. The governing equations of the mentioned framework are presented in details in Sect. 2. In Sect. 3, the described methodology is employed, first, for the analysis of a simple poroelastic problem in order to study the micro-macro interdependency of deformation and properties highlighting the feasibility and importance of the present methodology. Then, Darcy's experiment is simulated showing that employing this framework, there is no need for any correction to Darcy's law as the nonlinear fluid flow is captured automatically by updating the hydraulic conductivity. Last but not least, a cyclic uniaxial test, similar to the experiments on the brain tissue in [2], is carried out showing the viability of the simulation of complex real-world problems by the present computational framework. The promising results of the numerical examples open several directions for the future works, which together with discussion and conclusions are provided in Sect. 4. Fur-thermore, a part of the previously developed homogenisation procedure is placed in "Appendix".

Governing equations
Let us consider an arbitrary poroelastic domain Ω ∈ R 3 which consists of the solid matrix Ω s and interstitial fluid Ω f such that Ω = Ω s ∪ Ω f . We assume that the poroelastic medium has the average pore size d which is infinitesimal compared to the medium size L (e.g. as shown in Fig. 1). As the exploited analytical homogenisation and localisation are valid for infinitesimal strain, the time discretisation is the first matter to consider, clarifying the applicability of the multiscale methodology.
The time-dependent state U t (at time t) of the problem is advanced to the time instant t + Δt as follows: In one increment, U t+Δt is calculated from the intermediate configuration (configuration at time t), the time increment Δt, and solution increment ΔU. As the intermediate configuration and Δt are fixed, the problem decreases to finding ΔU and consequently U t+Δt via At this stage, for the sake of simplicity, we drop Δ notation, however, we emphasise that the following equations and mechanical relationships are considered to be within one increment (i.e. from time t to t + Δt), thus, applying the infinitesimal deformation assumptions.
Considering the mentioned condition, in one increment, the solid compartment is linear elastically interacting with incompressible Newtonian fluid percolating its pores with no-slip boundary condition on the solid-fluid interface. The corresponding equations are for solid, fluid, and interface conditions, respectively, where u indicates the displacement of solid phase, v represents the fluid velocity, C is the fourth rank solid matrix elasticity tensor, and μ is the interstitial fluid dynamic viscosity. τ and σ are the solid Cauchy and fluid viscous stress tensors, respectively. We highlight that the mentioned variables are at the physical scale and multiscale analysis is not applied, so far. Furthermore, the symmetric gradient operator ξ is defined as

Non-dimensionalisation
Poroelasticity has a wide range of applications from soil and rock mechanics to different biological tissues (such as brain tissue). It also has a complex interdependency between the units of measurement due to its multiscale and multiphysics nature which is a significant source of ambiguity leading to a misunderstanding of the model parameters and response [19]. In this study, before we proceed further, we take advantage of Non-dimensionalisation procedure which is similar to the one in [34]. For each example, we adopt four formally independent characteristic values for macroscale length L, average microscale cell dimension d, fluid viscosity μ c , and force f c such that where x and y are macroscale and microscale spatial variables, F is (any) force, and μ is the interstitial fluid dynamic viscosity. All other parameters in this study are nondimensionalised with respect to the mentioned independent characteristic values as where r indicates the non-dimensional parameter. K, M, and C are the hydraulic conductivity, Biot modulus, and effective (drained) elasticity tensor which are to be defined in Sect. 2.2.

Homogenisation
In this section, we outline the mathematical two-scale homogenisation based on non-dimensionalised parameters which is developed and employed in several studies in the literature [3,6,34]. For the reader's convenience and as we will refer to the notations later on, we provide more details of this procedure in Appendix A.
Removing the prime symbol (indicating nondimensionalised variables), for the sake of simplicity, all the Eqs. (2)-(8) will be written the same except for Eq. (5) which adopts the form in Ω f (11) where = d/L 1 is the length scale separation parameter. The coefficient 2 appears due to the scaling of Stokes flow in porous media (see e.g. [3,34]).
An infinitesimal length scale separation parameter allows to make use of asymptotic homogenisation technique. We assume that x and y ( y = x/ ) indicate two formally independent macroscale and microscale spatial variables, respectively.
Then, the differential operators transform into two independent ones for macroscale and microscale as (12) and, similarly, which separates the problem into two-scales. Applying the transformation (12) and (13) into Eqs. Equating the coefficients of 0 and 1 , by replacing every field (stress, displacements etc.) by its power of series representation, ψ (x, y) = ∞ l=0 ψ (l) (x, y) l , yields zero-th and first order equations (i.e. by choosing l = 0, 1) (Eqs. (66)-(79) in "Appendix") which are the principal ones of the upscaling process based on which the leading order solid displacement u (0) and hydrostatic pressure p (0) are microscale independent (locally constant) and could be referred to as macroscale variables. The zero-th (leading) and first order variables are shown by r(0) and r(1) , respectively.
As the relative fluid velocity and the solutions of the cell problems (will be defined later) are not constant at the microscale, we apply the integral average to exploit the acquired equations in the homogenisation process and macroscale. Assuming y-periodicity, no-growth limit [34], and the following ansatzes The system of PDEs governing the macroscale problems (in homogenised domain Ω h ) in a time increment is obtained as where, indicate effective stress, macroscopic strain tensor, and relative fluid velocity. This system of PDEs is characterised by some coefficients, namely, effective elasticity tensor, Biot coefficient and modulus, and hydraulic conductivity, that are defined, respectively, as where, C and φ indicate the elasticity tensor of solid matrix and porosity (i.e. volume fraction of fluid phase), respectively. The fourth rank tensor M and the second rank tensors Q and W are the solutions of three systems of PDEs resulting from the multiscale procedure (for more details regarding the auxiliary cell problems please see "Appendix" and [5]).

Remark 1
Although the solution of the auxiliary cell problems are obtained by solving a linear elastic type problem via FEM, considering Eq. (88), the second rank tensor ξ y (a) has the dimension L 2 f c meaning that Biot modulus M has the same dimension as C which is correct according to Eq. (20). We highlight that, consequently, the vector a has the characteristic dimension d L 2 f c which agrees with the solid ansatz (15).

Localisation
In this section, based on the asymptotic homogenisation, we develop an analytical localisation and properties rearrangement procedure in order to obtain the microscopic deformation due to the macroscale mechanical and hydraulic response. The ansatz (15)- (17), in fact, are the bridge between microscale and macroscale. As we mentioned earlier (and from Eq. (67) in "Appendix") u (0) is locally-constant meaning that it can only cause rigid body motions which does not affect the microscale deformation gradient tensor. The same condition holds for p (0) (see Eqs. (68) and (69) in "Appendix").
Taking the microscopic spatial gradient of Eq. (15) reads which, taking into account that ε (0) , and p (0) are microscale spatially independent (locally constant), reduces to Thus, taking the integral average from both sides yields where, ε (1) = ξ y u (1) is exploited to update solid matrix Young's modulus and Poisson's ratio, as well as porosity.
The coefficients of macroscale system of PDEs are computed based on the averaged values of strain and velocity components over their corresponding phase in the cell problem. Taking the integral average of Eq. (17) we can write which, considering the uniqueness condition (83), results in p (1) f = 0 and highlights that p (1) does not play a role in material properties update as they are considered in an average sense. Consequently, ε (1) is the only parameter that plays a role in this procedure so we call it microscopic strain tensor.

Microscopic properties update
At this stage, we flash back to the increment notation Δ (so ε (1) becomes Δε (1) ) for clarification in updating the stiffness/tangent tensor of the solid phase based on the current strain tensor ε (1) t+Δt . The latter can be calculated via ε (1) (1) .
We make use of the Lagrangian-based formulation to update the solid matrix mechanical properties. For the sake of clarification, we highlight that this formulation is only used for determining the updated stiffness tensor and can be replaced by other methods.
Assuming that the pore section remains circular (so one can consider that the average microscopic displacement gradient tensor is symmetric), the deformation gradient tensor F and Green strain tensor E can be calculated via where I is the second rank identity tensor. We highlight that ε (1) is the integral average over the entire cell which could be transformed to the integral average over the solid matrix by using the coefficient |Ω| |Ω s | which is equal to 1 1−φ in the case of the assumed unit cube cell. As the fluid part is assumed incompressible, then φ, V i s , and φ i are related by and where, J := det(F) is the total solid matrix volume change at the point. V i s and φ i are, respectively, the solid and fluid phase volume fractions in the reference configuration.
In this study, we consider isotropic compressible neo-Hookean strain energy function, which, can be written [35] where C 10 and D 1 are material parameters. We determine the latter parameters based on initial Young's modulus E i and Poisson's ratio ν i via The components of the deformation dependent tangent matrix for this material are given by The fourth rank tensor C t i jkl (note that the superscript t, here, does not indicate the time) is the deformation dependent tangent matrix calculated based on ε (1) t+Δt . We take this tensor as the current stiffness tensor for the next time increment. In other words, for each increment, we use the rearranged stiffness tensor based on the intermediate configuration. Having sufficiently small strain increments one can for the next increment. In this study, for the sake of simplicity, we assume that the solid matrix remains isotropic under the deformation so that its Young's modulus and Poisson's ratio can be updated by

ANN application for localisation and homogenisation
In order to obtain the coefficients of the system of PDEs governing the macroscale mechanical and hydraulic response (e.g. Biot modulus and coefficient), the tensors M, Q, and K are required which could be determined by solving the auxiliary cell problems (defined via Eqs. (83)-(84) in "Appendix"). On the other hand, M and Q are necessary to carry out the localisation to update the solid matrix elasticity tensor and the cell porosity. We notice that due to the inhomogeneous/une-qual strain and pore pressure distribution in the space and time the mentioned cell properties are spatially and temporally dependent imposing the need to solve the cell problems at each spatial point and time increment which is not feasible (or too expensive). In order to address the above mentioned issue, we employ Artificial Neural Networks (ANNs) [31] in a similar way that is used in a former study of the current authors [7]. However, here, instead of computing the coefficients of macroscale system of PDEs, we make use of the capacity of parallel ANNs to calculate M, Q, and K so that we carry out the localisation procedure using the same outputs which means increasing the efficiency of the framework. This network can be written as a micro-macro scale transfer function i.e.
which, assuming solid matrix isotropy and geometric rotational invariance with respect to three orthogonal axes, reduces to [6,26] Considering that throughout the homogenisation procedure the parameters E and μ are non-dimensionalised and the problems are linear elastic solid and newtonian fluid types one can reduce the input size of the ANNs to At this point we highlight that, according to Eq. (20) and Remark 1, if the nondimensional solid matrix Young's modulus E is different from its initial value (which is the case when including material remodelling), as it is excluded from ANN inputs, one should consider its effects on Q 11 multiplying it by E i E .
In this study, we make use of five parallel networks each responsible for only one of the outputs in order to decrease the complexity while increase the accuracy of the training procedure. Figure 2 schematically shows the configuration of this type of ANNs while, mathematically, the feed-forward process which calculates the output based on the inputs through hidden layers can be expressed via where a (i) k is the value of the k-th neurone in the (i)-th layer, (i) i ∈ (0, 1, .., L) is the number of hidden layers, k and j indicate the number of the neurones in the (i)-th and (i − 1)th layers, respectively. The activation function ReLU (x) = max(0, x) determines whether a neurone is active (adopts nonzero value).
At this stage, in order to achieve accurate outputs, we need to train the networks which is viable, in the simplest case, by means of a training dataset consisting of sufficiently rich "exact" outputs each corresponding to distinct inputs. The later information, in this case, is provided by solving the auxiliary cell problems via FE method. We solve the cell problems for 50 different porosities φ covering the interval from 0.082 to 0.783 with fine steps and for each porosity 50 different Poisson's ratio ν are considered so that the role of both variables are well reflected.
In the training procedure, the ANNs are fed with the same inputs (φ and ν) as the one in the FE problems (the target ones). Then, based on the later target, the L 2 error is computed and an optimiser (in this case Adam optimiser is chosen [36]) to minimise the difference between the ANNs output and the target values (more details are available widely in the literature, e.g. [7,37,38]).

Incremental nonlinear algorithm
In strongly coupled solid-fluid problems such as the poroelastic models, fine discretisation of both space and time is required even assuming infinitesimal deformations. The latter is achieved by means of the incremental solution via Finite Element Method in which the total time of the problem is divided into sufficiently small increments Δt so that accurate numerical temporal integration/differentiation can be achieved. Here, the simple linear time integration scheme is used. The time discretisation described via Eq. (1) can be rewritten for the macroscopic poroelastic case as We note that, in total, the network size is even smaller than the one in [7] but, although the relationship is more complex here, the results are more accurate In other words, the force vector and fluid relative velocity (which are history dependent) are calculated based on the intermediate configuration (configuration at time t), the time increment Δt, and displacement and pore pressure increment Δu (0) and Δp (0) , respectively. As the intermediate configuration and Δt are fixed, the problem decreases to finding Δu (0) and Δp (0) for which the infinitesimal deformation assumption is valid, thus, the employed mathematical homogenisation and localisation are applicable. Consider a vector of DOF increments that must yield balance of linear momentum and the mass conservation provided, respectively, by Eqs. (18) and (20) which can be approximated using the weak formulation. The latter can be obtained by integrating mentioned equations with respect to the corresponding volume and multiplying them with arbitrary test functions. In this case, two test functions δu (0) and δ p (0) which represent the variation of the macroscale displacement and pore pressure are required. Using Eqs. (49) and (50) and applying the divergence theorem the weak form of the governing equations to be solved in each increment is which rt or r t show the situation of r at the intermediate configuration (at time t). We highlight that the first terms inside each integral are fixed and known which form the residuals of the intermediate configuration and could be factorised. However, in this study, in order to avoid the accumulation of the residuals (numerical errors), we consider them in the formulation and implementation [39].
Then, by substituting the constitutive equations (19) and (21) into (51) The coefficientsC t , Q α t , M t , and K t are the updated ones based on the updated microscopic material and geometrical parameters on the intermediate configuration.
As in [7], free drainage can be introduced on a surface by enforcing as a part of the fourth term on the right hand side of Eq. (52).
Here, Δp and Δx are the differences in pressure and distance between the surface of free drainage (S f d ) and the environment. We highlight that the latter refers to the surrounding space where the experiment is being carried out, in particular, at/near to the free drainage surface. If the free drainage condition is not introduced, an impermeable boundary condition is imposed.

Implementation verification
Due to the complexity of the methodology, the verification procedure is not a straightforward task. In this case, we divide the framework into simpler cases and verify them step by step.
-The upscaling at the beginning of each iteration is carried out via ANN which is verified by producing random test dataset and observing high accuracy shown in Fig. 3 and as is done in [7]. -At each iteration a linear consolidation problem is solved.
This part is verified agreeing Terzaghi's analytical solution for a column of poroelastic material as in [6,26]. -In order to verify the general incremental nonlinear analysis based on the specified strain energy function we perform a uniaxial test a linear elastic model, a typical nonlinear model, and the remodelling based incremental nonlinear analysis. The results are provided in Fig. 4. -Finally, in order to verify the integrity of the presented incremental nonlinear poroelastic framework we compare the maximum settlement of a column of poroelastic material under compression computed by the presented incremental nonlinear model with the linear one shown in Fig. 5. Furthermore, the results provided in the next section are also helpful for verification purposes. The associated model geometry, material properties, BCs etc. is the same as the first numerical example.

Numerical examples
Employing the present AI-assisted incremental computational method based on remodelling and asymptotic homogenisation/localisation we are able to analyse nonlinear poroelastic problems. In [7] the initial spatially dependent porosity and solid matrix Poisson ratio was considered. However, they were assumed as constants and the method was not valid at large deformations. Here, we consider a nonlinear elastic material for the solid matrix, which is valid at finite strain. Thus, the porosity depends on the solid matrix deformation. One major advantage of this method over the other methods in the field of nonlinear poroelasticity is that this framework, for the first time, considers the full macro and micro scale response/properties interdependencies. Hence it responds to several questions, in this study, including how are the microscale, as well as effective, properties affected by the macroscale mechanical and hydraulic response? Is an appropriate deviation from the overall linear Darcy's law for fluid flow in porous media (under finite deformation) enforced or fractional/exponential Darcy's law is needed? Why the preconditioning effects (residual strains), hysteresis response, and Mullins effects [40] appear in poroelastic media such as brain tissue under cyclic loading? And, finally, how important is (or when is it important) to employ such a nonlinear method? We highlight that the Mullins effect refers to the situation at which the constitutive behaviour of the material depends on the loading/deformation history, particularly, in the case of cyclic tests, which can be seen as a pseudoelastic behaviour. Furthermore, hysteresis is the dependence of the system's response on its history (could be inelastic), which exists in rate-dependent and path-dependent problems such as plasticity, viscoelasticity and poroelasticity. In the latter case, this behaviour is due to the Fluid-Solid Interaction and the time/path-dependent nature of the problem due to fluid drainage or absorption [41]. The numerical examples performed in this section respond to the above-mentioned questions.

Uni-axial consolidation test
Let us start from a simple example which has been exploited many times through the history of fully or partially fluidsaturated porous media. In this test, a column of poroelastic material is compressed by a constant mechanical pressure with free drainage BC on top (please see [7] for a representative sketch of the problem). The model is given sufficient time to approach the steady-state where the time-dependent Although this test is very simple, it is helpful to gain a better insight and to observe basic facts during the consolidation procedure. In order to focus on the relevant novelties (although, without considerable additional cost, spatially dependent initial properties could be enforced), we embraced initially (macroscopic) homogeneous problems. Let us assume initial isotropic solid matrix Young's modulus and Poisson's ratio, respectively, equal to E i = 15e6 [Pa] and ν i = 0.3 with initial porosity of φ i = 0.3. The mechan- We, first, visualise the macroscale response together with microscale and effective properties in order to understand their interdependency. Figure 6 shows the mac-roscale response of the model along the axial axis (length) at different times with Fig. 6a, b being the settlement and interstitial pore pressure, respectively. The study of the spatial profile of settlement (which is the major displacement component in this study) and pore pressure are important and the starting point to understand the profile of micro and macro scale properties during the analysis. The spatial distribution of the properties during the analysis are monitored in the uniaxial consolidation test. As it is shown in Figs. 7 and 8, due to the response-driven remodelling, both microscopic and effective properties of the medium become spatially dependent instantly after load application. Considering initially homogeneous medium allows us to show how the spatial distribution of the characterising properties follow the macroscopic hydraulic and mechanical responses. Considering Eq. (31), we notice that shortly after the load application, in the major part of the length, pore pressure profile is the variable enforcing the variation of the micro and macro scale properties, while the spatial gradient of displacement affects only the part near the top of the model. This is due to the sharp decrease in pore pressure which is a consequence of the free drainage BC at top. In fact, the profiles of the properties at short times, mostly, should follow the variations of a high pore pressure while at longer times they are dominated by the solid strain.
At small times, the porosity distribution shown in Fig.  7a (which is similar to the studies focused on the porosity variation such as [42]) reveals sharp spatial variations starting from the values below the initial porosity at the top of the model due to the solid strain (which shows a volume increase in the solid matrix) reaching to the values greater than the initial one under the influence of pore pressure. The latter highlights the fact that under the positive pore pressure, due to external compressive load, the porosity increases, causing a contraction in the solid matrix. However, as time passes, the pore pressure decreases and the porosity profile follows the macroscale solid strain showing an obvious effect of consolidation, which means higher volume fraction of the solid matter. We highlight that the direction of a positive pore fluid hydrostatic pressure acting on the solid-fluid interface is from the centre of the pore towards the solid domain.
The solid matrix material properties (Young's modulus and Poisson's ratio) are updated by the linearisation of the strain energy function given by Eqs. (38) and (41) based on the microscale kinematics. Figure 7b shows that considering the porosity variation, the solid matrix Young's modulus is smaller than its initial value when there is an expansion in On the other hand, this profile is not the same in the case of Poisson's ratio shown in Fig. 7c. In fact, the latter decreases in both contraction and expansion as it is lower than its initial value in all the model. As the profile of Poisson's ratio and Young's modulus are not varying in the same direc-tion, in order to monitor the overall material compressibility we compute and provide the variation of the bulk modulus which includes both Young's modulus and Poisson's ratio and considered directly as a measure of the compressibility (and consequently strain-softening or hardening). Figure 7d reveals that the solid domain properties variation obeys the given neo-Hookean material model, which could be served as a verification of the FE implementation. Fig. 8 The profiles of the effective poroelastic properties calculated based on the macroscopic mechanical and hydraulic response. We highlight that, due to the nonlinear and asymptotic nature of upscaling, the reader may not generalise the shown results and conclusions for the problems with different ranges of initial values At this stage, exploiting the instant upscaling computation via ANN, the fourth rank tensor M, the second rank tensor Q, and the hydraulic conductivity K are obtained. The poroelastic propertiesC,α, K, and M used in the homogenised model (see Eqs. (18)-(21)), are calculated via Eqs. (25)- (28). Furthermore, the effective Young's modulus, Poisson's ratio, and shear modulus are calculated based on cubic symmetricC via Eqs. (42)- (43) replacing C i j byC i j . Figure 8 shows the coefficients of macroscale system of equations that characterise the macroscale mechanical and hydraulic response of the homogenised model. The profile of effective Young's modulus demonstrated in Fig. 8a is derived from displacement and pore pressure profiles and the resultant interaction of varying porosity and solid matrix stiffening/softening. In other words, as demonstrated in the literature [5,6], increasing porosity at constant solid matrix properties results in a decrease in effective Young's modulus, while, considering remodelling of the solid material properties, it causes higher solid matrix Young's modulus. The latter, in turn, at constant porosity, increases the effective Young's modulus. A decrease in porosity has inverse effects. The result of this complex behaviour can only be determined using the provided AI-assisted multiscale framework. It worth to note that, despite the complex behaviour, effective solid matrix Young's modulus stands at a higher value than its initial value throughout the time and space in the problem.
Effective shear modulus in Fig. 8b has a similar complexity as effective Young's modulus although its profile is, somehow, smoother. In fact, similar to effective Young's modulus, it adopts values higher than the initial one resulting in strain stiffening. Effective Poisson's ratio provided in Fig.  8c has a less complex profile as it depends only on the solid matrix Poisson's ratio and porosity. As an overall description, effective Poisson's ratio decreases, increasing the pore pressure and increases, due to the consolidation at the steadystate.
The fluid drainage and consolidation procedure is considerably influenced by Biot modulus, Hydraulic conductivity, and Biot coefficient. Biot modulus determines the rate of change in pore pressure due to the macroscopic solid and fluid volume change. Figure 8d shows that this parameter exhibits an increase due to the positive pore pressure while it decreases due to the negative macroscale solid strain. The considerable decrease in hydraulic conductivity and Biot coefficient in Figs. 8e, f is due to the decreasing porosity at the top of the model as well as, in the case of Biot coefficient, the resultant change in the solid matrix material properties. We highlight that Biot coefficient acts as the coupling term between solid and fluid and hydraulic conductivity relates the spatial gradient of pore pressure to the fluid velocity relative to the solid displacement rate. The overall effects and impor-tance of the variation in the effective poroelastic properties are reflected in the macroscopic response of the medium.
As a comparison mean and to understand the influences of the interplay between macroscale response and macroscale/microscale properties, the settlement and pore pressure profiles of both linear Biot consolidation and the present method are provided in Fig. 9. The latter shows a higher polynomial degree of nonlinearity during the transient state while the difference between them is closer to a coefficient change approaching the steady-state.
Applying divergence theorem and having small time incrementation the fluid drainage in one time increment from unit section surface area could be computed via We can also compute the total fluid drainage until m th increment via Δε n v f (55) Figure 9d shows the normalised fluid drainage (w.r.t the maximum value of the linear case) Noteworthy is that all the mentioned properties in both scales are spatially dependent during the transient state while they become homogeneous at the steady-state. This highlights the fact that considering the described material response-properties interdependency is particularly more important in transient problems, e.g. cyclic loading. The effects of spatially dependent solid matrix material properties and porosity is thoroughly discussed in [7].

Deviation from Darcy's law
The applicability of Darcy's law [15] into poroelastic problems in asymptotic range (at small Reynolds number) has been under debate since its publication. In [19] it is mentioned that the linear fit to the experimental data yields a negative pressure gradient when the flow rate is zero (e.g. −∇ p = −0.745773 + 0.640324q where q is the flow rate) which means that there is a need to nonlinear corrections to Darcy's law. In this study, we model the Darcy's filtration experiment for a wide range of pressure differences (ΔP) reaching to a maximum ΔP max = E i /4 (where E i is the initial solid matrix Young's modulus). We calculate the relative fluid velocity when the simulation reaches the steady-state We calculate the dimensionless deviation from Darcy's law via Fig. 9 A comparison between the final response of the poroelastic medium calculated via the presented incremental nonlinear method (labeled as "NL") and the outcome of the linear poroelasticity (labeled as "L") where K i and r max are the initial hydraulic conductivity and the maximum value of r throughout the experiment, respectively. v r f is the relative fluid velocity computed via the presented incremental nonlinear method in which the solid matrix properties and porosity (consequently hydraulic conductivity) are updated at each increment. Figure 10a shows that v r f deviates considerably from the relative fluid velocity computed via Darcy's law (with the initial hydraulic conductivity) at pressure differences higher than 0.2ΔP max which agrees with the information provided in the literature (see, e.g. [43]). The non-dimensional deviation value y in Eq. (56) at different pressure drops is provided in Fig. 10b. The latter provides a mean to compare the results provided via the presented method with the linear and quadratic deviations (i.e. y = (ΔP/ΔP max ) and y = (ΔP/ΔP max ) 2 , respectively) which shows a great match with the quadratic one. The latter was also concluded in several studies such as in [16,19].  Figure 10c reveals the fact that the average cell dimension plays an important role in the deviation from the linear Darcy's law highlighting that the more the average cell dimension, the more the deviation. One could also conclude that the more the cell average dimension, the less is ΔP limit at which the nonlinear model starts to, considerably, deviate from the linear. This can be seen in the definition of dimensional hydraulic conductivity which is directly proportional to d 2 . The results provided in Fig. 10c is, qualitatively, in agreement with the effect of the soil grain size on deviation from Darcy's law provided by means of experimental results of [44].
From the microscale cell problems which are solved to determine the hydraulic conductivity and the dimensionalisation procedure we notice that, apart from average cell dimension and viscosity, porosity is another important parameter in this problem which is considerably affected by solid matrix Poisson's ratio (ν) at a specific pore pressure distribution. The reason is that ν plays a major role in the solid matrix volume change leading to a change in porosity. Figure 10d shows the influence of solid matrix Poisson's ratio on deviation from linear poroelasticity.
From our observations, we conclude that considering the medium remodelling after each time increment, there is no need for the formerly proposed corrections to Darcy's law.

Uniaxial cyclic test on brain tissue
Detailed understanding of complex and important mechanical properties and response of brain tissues provides us with vital accurate predictions for the design of treatment protocols and assessment of risk. For this purpose, numerical analysis with appropriate models and accurate parameters can provide crucial information on the effects of external loads. The application of these external loads can be slow or fast and their effects can be short or long term. It is seen that the brain tissue mechanical response differs considerably in each case. The complexity of this tissue behaviour is well studied in [45] introducing several challenges in this field, however, an appropriate poroelastic modelling of brain tissue was missing. The mechanical response of brain tissue under uniaxial experimental tests at free drainage from, respectively, sides and top of the model in [2] provides the first direct evidence of poroelastic behaviour of brain parenchyma which is proven by comparing the results of uniaxial consolidation experiment with Terzaghi's analytical solution. However, the numerical/analytical models for modelling the uniaxial cyclic experiments are not considered via poroelastic models. According to the latter study, under uniaxial cyclic loading, "brain tissue exhibits a peculiar nonlinear mechanical behaviour, exhibiting hysteresis, Mullins effect and residual strain, qualitatively similar to that observed in filled elastomers".
In this subsection, we consider the application of the presented framework into the mechanical modelling of biological soft tissues which are typically considered as nonlinear problems. In order to focus on the role of the novel incremental nonlinear poroelastic methodology for mechanical response analysis of the medium, we, again, choose the simple neo-Hookean hyperelastic model for the solid matrix properties rearrangement. We show, for the first time, that the poroelastic nature of brain tissue considering microscopic properties rearrangement due to the macroscopic mechanical and hydraulic response plays a major role (if not the only reason) in the appearance of the mentioned peculiar behaviour.
First, one should pay attention to the characteristic values of the parameters of this application. For example, due to small dimension of the tissue, it is more convenient if the macroscale average length scale is millimetres (L = 10 −3 [m]) so that one can also avoid numerical errors due to very small dimensions of the elements. One can perform the calculations based on dimensionless parameters calculated via (10) in order to avoid ambiguity. The dynamic viscosity of cerebrospinal fluid CSF [46] is assumed as the characteristic one μ c = 10 −3 [Pa s]. Average microscopic periodic cell dimension d = 20 −6 [m] is chosen from the literature [47][48][49] and unit force F c = 10 −3 [N] is adopted to reach an appropriate unit pressure. For example, characteristic time and pressure are t c = 1[s] and P c = 10 3 [Pa], respectively, and, assuming solid matrix initial Poisson's ratio ν i = 0.3 [−] and Young's modulus E i = 13.5 × 10 3 [Pa], the dimensionless initial Young's modulus is E i = 13.5 [−] which is the same as the one assumed in ANN training.
We first make use of the conventional linear poroelastic model, providing a comparison means to understand to what extent the provided incremental nonlinear framework is important in mechanical modelling of brain tissue. Figure 11b shows that, provided appropriate parameters, even the linear poroelastic model shows the hysteresis behaviour (indicating the energy dissipation) of the tissue without introducing viscoelasticity, damage, etc. The preconditioning effect, however, is considerably smaller than the empirically captured values in [2]. We highlight that the cycle time T cycle = 108 [s] is chosen to introduce, more or less, the same loading conditions as in the latter study. Employing the presented incremental method for this brain tissue deformation example, a considerable difference is observed. The mechanical response of the model shown in Fig. 11a exhibits the effect of preconditioning to a larger extent which agrees with the experimental results provided in [2] for brain tissue (and in [50,51] for fibrous connective tissues such as ligament) highlighting the importance of considering remodelling and micro-macro interaction of nonlinear poroelastic media under finite deformation. The preconditioning effects causes the cycles to move to the right hand side of the plot which means that the deformed length of the model increases after each cycle which was again observed in the latter study [2]. Furthermore, the hysteresis behaviour is slightly different from the linear poroelastic one and is closer to the experimentally observed ones.
After the first cycle, there is a considerable cycle residual strain (the gap between the beginning of the loading and the end of the unloading path in one cycle) which is related to the fluid absorption due to the model deformation. As the fluid flow inside the medium requires time, we expect that increasing the loading rate the cycle residual strain and consequently, the preconditioning effects decreases and the medium responds similarly to an elastic one. Figure 12 highlights that not only the preconditioning effects but also the hysteresis behaviour are due to fluid flow and its interaction with the solid phase. Furthermore, it is shown that by increasing the loading rate, the minimum stretch decreases, which means a decrease in the medium compressibility. In this case, the incompressible interstitial fluid does not percolate through the porous medium due to the insufficient time The model response shows that the latter is more accurate, particularly, in case of soft poroelastic tissues. We highlight that, in the nonlinear case, the material model of the solid phase of the cell is assumed neo-Hookean with initial elastic properties adopted from the linear one in which case static fluid phase filling the pores could be assumed (undrained case). Considering Eqs. (19) and (20) the model deformation follows the undrained elasticity tensor (for more details see [3,5,34]) which expected to be far less compressible than the drained medium [5].
On the other hand, decreasing the strain rate, the given time for a small deformation increases which, in turn, decreases the resultant pore pressure driven by the solid deformation (see Eq. (20)) approaching the overall material response to the undrained case. The latter is expected to be considerably more compressible and softer as without Eq. (20) the pores could be considered as void spaces. Figure 13 shows that increasing the cycle time, the dissipated energy due to the hysteresis response as well as the preconditioning effects and residual strains increases up to a maximum point (see Fig. 13a, b) from which they decrease approaching zero at very long times shown in Figs. 13c, d. The latter shows the undrained response of the porous media under cyclic loading which is, compared with Fig. 12d, considerably softer and more compressible (note the difference in the applied load).
In the experiments carried out in [2] a behaviour similar to what is called Mullins effect [40] together with residual strain is observed. In [2] this behaviour is modelled by means of the pseudo-elastic model provided in [52] neglecting the solid-fluid interaction, multiscale nature of the problem, and the observed residual strain. Figure 14 shows that the present modelling methodology results in a macroscopic mechanical response similar to what is observed in the experiments, however, based on simple compressible neo-Hookean material employed in the present multiscale and multiphysics methodology. In this figure, the model deformation under both cyclic and monotonic loading is provided so that the dependence of the model response on the previous cycles or the deformation history is clear. As there is no term in the material model accounting for the Mullins effect, we conclude that this type of response is due to the nonlinear poroelastic nature of the medium. Furthermore, due to the chosen material model which imposes strain-softening in tensile deformation and strain-hardening is compressive one, as expected, a greater residual strain is captured under tensile loading (shown in Fig. 14a) compared with the compression test (shown in 14b).
The authors highlight that, although several key features and aspects of the brain tissue response to the external loads are addressed, the accurate modelling and simulation of this complex tissue requires further study and effort focused, for example, on the solid phase material model as well as the cell geometry properties.

Conclusions and future work
In this study, ANNs, a standard multiscale approach (asymptotic or two-scale homogenisation), and an incremental FE analysis framework are integrated to provide a novel and robust computational method for analysis of multiphysics and multiscale poroelastic problems under global deformations that are large enough to invalidate the assumptions of infinitesimal strain theory. This framework considers the full interdependency between the homogenised hydraulic and mechanical response and microscopic (and, consequently, macroscopic) properties rearrangement (multiscale remodelling) which encompasses different sources of nonlinearity, namely, solid-fluid interaction as well as the nonlinear governing equations (due to updating effective coefficients). The whole procedure is based on non-dimesionalised variables and properties so that it is applicable in a wide range of problems in different scales such as soil and biomechanics.
Applied into the soil mechanics (in the sense of the characteristic values), Terzaghi's consolidation and Darcy's fluid Fig. 13 The effects of increasing T cycle . We highlight that all other properties are the same as the reference case shown in Fig. 11a filtration tests are reconstructed numerically under finite deformation providing detailed information on the combined effects of the solid deformation and pore pressure on the microscopic and effective properties of the medium, as well as its overall response. Agreeing with the experimental results in the literature, for example, the fluid filtration example shows that, employing the presented method, there is no need to the corrections to Darcy's law for fluid flow in poroelastic media. Furthermore, the average cell dimension (grain size) and initial Poisson's ratio nonlinearly increase and decrease the deviation from Darcy's law, respectively. In Terzaghi's test, the distribution (in space and time) of the mechanical and hydraulic response, as well as microscopic and effective properties are monitored until the steady-state is reached. In general, we observe spatially and time-dependent distribution of the properties and variables highlighting the complexity of the problem. Among many other results, we show that, in this problem, the porosity of the medium under mechanical pressure increases instantly after loading (due to the resultant pore pressure) while it decreases during the transient state, as a consequence of the consolidation and fluid drainage, reaching a constant value smaller than the initial one at the steady-state. Moreover, the captured dependency of the effective properties on the changes in microscopic properties (porosity and solid matrix material properties) due to the solid matrix contraction/expansion and its effects on the overall response of the media highlights the importance of employing the presented method.
Furthermore, we simulate a uniaxial cyclic test on the brain tissue as a poroelastic medium (which is known for its complex and peculiar mechanical response). The test is characterised to mimic the conditions of the experimental tests so that we can evaluate the results. The hysteresis and preconditioning effects together with a behaviour similar to the observed response in Mullins tests are captured for the first time only based on the poroelastic nature of the tissue. The effects of the loading/deformation rate is thoroughly studied as an important source of ambiguity in brain tissue properties identification. Although the numerical results agree with the experimental ones in several aspects, more effort is required for accurate modelling of brain tissue which can be the subject of a future study.
Apart from the development of the mentioned method, several aspects of fundamental issues in the context of poroelasticity have been addressed for the first time in this study. However, there are more challenges to be overcome in the future studies which can be divided into three directions. Firstly, an application-specific and robust nonlinear material law for the solid matrix, able to accurately describe the complex poroelastic media such as brain tissue, is to be cho-sen/discovered. This, in turn, requires a robust and specific framework for microscopic model parameter identification using AI (similar to poroelastography but considering solid matrix properties). Secondly, the future works could be oriented to remove some simplifying assumptions such as pore section and solid matrix material, respectively, remaining circular and isotropic under the deformation during the analysis. Addressing the nonlinear cell problems, in the context of poroelasticity and composites, introduced in [27][28][29] will be considered as a potential future direction. Last but not least, based on the observations using the presented computational framework, phenomenological homogenised material models can be adopted/introduced considering the mentioned interdependencies and interactions between the micro and macro properties and responses. In other words, as soft poroelastic tissues may not be easily available, one can partially replace the experiments with the numerical simulations using the present method to gain a deep insight of the required elements in the material models for poroelastic media. 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://creativecomm ons.org/licenses/by/4.0/.
where Ω s , Ω f , Γ , A, and a represent the solid and fluid domains, their interface, a third rank tensor and a vector, respectively. n is the inward unit vector normal to the solidfluid interface Γ , and M := ξ y (A), Q := ξ y (a), P is an auxiliary vector that encodes microscale information by relating the first order hydrostatic pressure to ∇ x p (0) ( p (1) = −P ·∇ x p (0) ), and I is the second rank identity tensor.