Multiscale homogenization for dual porosity time-dependent Darcy–Brinkman/Darcy coupling and its application to the lymph node

We study the coupling between time-dependent Darcy–Brinkman and the Darcy equations at the microscale subjected to inhomogeneous body forces and initial conditions to describe a double porosity problem. We derive the homogenized governing equations for this problem using the asymptotic homogenization technique, and as macroscopic results, we obtain a coupling between two Darcy equations, one of which with memory effects, with mass exchange between phases. The memory effects are a consequence of considering the time dependence in the Darcy–Brinkman equation, and they allow us to study in more detail the role of time in the problem under consideration. After the formulation of the model, we solve it in a simplified setting and we use it to describe the movement of fluid within a vascularized lymph node.


Introduction
The flow of fluids through porous media is a fundamental process with wide-ranging applications in various fields, such as hydrogeology and biology.In particular, porous media with dual porosity, where fluid flow occurs in different compartments with different pore structures, have garnered growing attention because of their widespread occurrence in both natural and engineered environments [1][2][3][4][5][6][7].
In this work, we present a multiscale model for dual-porosity porous media using the asymptotic homogenization technique [8,9] that couples a time-dependent Darcy-Brinkman equation [10] with a Darcy equation [8,9,11,12] which describes the fluid flow of the blood vessels inside the node, as our starting point.Here, we assume that the Darcy equation depends on time only parametrically.The theoretical justification of the Darcy-Brinkman equation is less warranted from a multiscale perspective compared to the Darcy equation [8,[13][14][15][16][17][18][19]; however, since its structure is halfway between the Darcy and the Stokes equations, it can serve as a valid starting point for a multiscale formulation [1], and grants the ability to specify boundary conditions in more detail [20,21].Moreover, the Darcy-Brinkman equation is also well defined in a time-dependent setting [14,22,23].
We consider a multiscale volume load that drives the coarse scale fluid flow acting on both the Darcy and Darcy-Brinkman problems [24].These forces can occur from the use of electromagnetic fields, in magnetorheological fluids and in electrolytes that permeate non-uniform tissues [24].Moreover, we consider a multiscale initial condition of the time-dependent problem that we take into account.Such initial conditions can arise for the flow of fluids in porous media with drug injection [25] or when as an initial condition we have another multiscale motion.
The derived model comprises porous media flow with memory effects (e.g.[8,[26][27][28][29]), which are encoded in the effective hydraulic conductivity; thus, the nature of the obtained differential problem is intrinsically different with respect to [1].The latter is represented by a convolution in time between the time-dependent hydraulic conductivity (obtained by solving the differential problem at the cell level) and the macroscopic pressure gradient.The macroscopic Darcy equation with memory effect that we obtain is new and differs from previous works on this subject such as [8,[26][27][28][29] because we consider the effect of inhomogeneous body forces, a multiscale initial condition and the fluid exchange with another phase (described here by the Darcy equation).Using asymptotic homogenization, we pass from an equation at the microscale where the present time is necessary for deducing the future [30] to a macroscopic equation where the present and the past are necessary to deduce the future [31].The memory effects appear to be a natural framework when dealing with the homogenization of this kind of problems [30,31].
Among the countless applications of dual-porosity porous media, the main motivation that guided us to propose and study this model is the application of the latter to the fluid flow of a lymph node, an essential part of the lymphatic system.The lymphatic system is composed of a network of vessels, capillaries and organs [32].The interstitial fluid is drained by the lymphatic capillaries and it is called lymph once inside the lymphatic system.Inside the lymphatic vessels, there is a series of one-way valves which prevent retrograde flows.The part of the vessel between two valves is called lymphangion.Lymphangion walls are innervated with sympathetic and parasympathetic nerves [33][34][35] and can perform rhythmic contractions, making lymph transport a time-dependent flow [33,[36][37][38].The lymphatic system is an integral part of the immune system thanks to the lymph nodes.The lymph node is vital for immune defence, housing B and T cells that circulate in the body to protect against infection.B cells generate antibodies to fight antigens, while stimulated B cells transform into plasma or memory cells.Antigenpresenting cells capture antigens and present them to T cells in the lymph nodes, activating the adaptive immune response.Lymph transport has an important function in transporting immune cells, proteins, cancer metastasis, drugs and so on [39][40][41].Furthermore, the movement of fluid through lymph nodes promotes the expression of chemokines, establishing a chemokine gradient that directs the movement of immune cells into the node.Changes in lymph transport play an important role in several pathologies.Elevated fluid flow augments the growth rate and sensitivity to drugs in specific forms of lymphomas [42].Inadequate lymph transport can lead to a condition called lymphoedema, where excess interstitial fluid accumulates in the tissues, and it is often caused due to an impairment of lymph nodes [36,43].From a mechanical perspective, the key characteristics of the lymph node include the lymphoid compartment (LC), which is a porous bulk region and forms the parenchyma of the lymph node, and the subcapsular sinus (SCS), a thin channel near the wall that surround the LC and that allows free-fluid flow [33].The fluid can enter the LC from the SCS through a conduit system network created by fibroblastic reticular cells (FRCs) [44][45][46] which forms the main porous region of the lymph node.The lymph node is a highly vascularized organ, and inside the LC there are plenty of blood vessels that allow fluid and substance exchange [47][48][49], making the lymph node an important connection region between the lymphatic system and the blood system.
To the best of our knowledge, only a few models in the literature try to describe the fluid dynamical aspects of lymph node behaviour [50][51][52].In [53,54], it is explored how the fluid flow in the lymph node is influenced by its internal structure, using an image-based model to establish a relationship between the greyscale values of the images and the permeability of the lymph node tissue, using this as the royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 permeability for the steady Darcy equation.In [47], the authors conduct a parameter sensitivity analysis by using a computational flow model based on a mouse popliteal lymph node, coupling a steady Darcy-Brinkman equation in the LC with a steady Navier-Stokes equation in the SCS.Grebennikov and others [45,46] study the lymph flow through the conduit system network using an object-oriented computational algorithm to generate the three-dimensional geometry of the FRC graph network.In [55], the fluid flow within the lymph node is simulated by using a microfluidic device that mimics the microenvironment of a lymph node.Birmingham et al. [39] focused on the lymph node's SCS fluid dynamics thanks to a microfluidic platform, evaluating how physiological flow patterns impact the adhesion of metastatic cancer cells (thanks to the wall shear stress values).In [56], a numerical method is developed to simulate the fluid flow in the lymph node using boundary integral equations, and then in [57] the authors provide an artificial neural network model to describe the lymph node drainage function.In [1,58,59], we have the first explicit models that describe the fluid flow in the lymph node in simplified geometries.In particular, in [58,59], they found a divergence-free explicit and numerical solution in a time-dependent setting, with a very idealized geometry for [58] and a spherical geometry in [59].Girelli et al. [1] describe the blood vessel drainage function in the lymph node considering the multiscale nature of the latter in a steady setting, obtaining a rigorous mathematical model using the asymptotic homogenization technique; thanks to this approach, they were able to describe the fluid flow inside the conduit system network (formed by FRC) and inside the blood vessels networks.
This work addresses a crucial limitation characterizing the work in [1], as it takes into account the time-dependent character of the pulsatile flow which takes place in the lymph node due to the pulsation of the lymphangions [36,60,61] and its upscaling onto the macroscale.In fact, we present a multiscale model using the asymptotic homogenization technique [8,9] that couples a time-dependent Darcy-Brinkman equation [10], which describes the fluid flow inside the conduit system network formed by FRC, with a Darcy equation [8,9,11,12], which describes the fluid flow of the blood vessels inside the node, as our starting point.The Darcy equation depends on time only parametrically because the time-dependency we take into account is given by the lymphangion pulsation and affects mainly the lymph flow inside the FRC network.This model focuses on the porous region of the lymph node (LC) and the fluid exchange occurring solely between the node and the blood vessels within this specific area [47][48][49].It is crucial to emphasize that, as highlighted in [1], our analysis begins with the Darcy-Brinkman and Darcy equations.This implies that variations in pore-scale geometry have been effectively averaged out, eliminating the necessity for precise details about the intricate and challenging-to-describe microstructure geometry of the lymph node.
In §2, we show the starting equations and boundary conditions of our problem, based on the balance equations of continuum mechanics.In §3, we employ the asymptotic homogenization technique and in §4, we find the macroscopic averaged equations.In §5, we find the explicit solution to our problem in spherical geometry under the hypothesis of axisymmetry with respect to the azimuthal angle and isotropy of the porous medium using the Fourier transform and using the explicit solution that we already found in [1].In §6, we describe the numerical simulation that we use to solve the problem in the cell domain.Finally, in §7, we solve the explicit solution we found in §5 using lymph node physiological data obtained from the literature, allowing us to study in more detail the behaviour of the lymph inside the lymph node in a time-dependent setting.

Formulation of the starting problem
We consider the domain V ¼ V v < V m , where Ω m and Ω v denote regions of the matrix and the vessels, respectively.A sketch of a portion of the three-dimensional domain at hand comprising the two phases is provided in figure 1.We describe the fluid flow in Ω v via Darcy's Law supplemented by a body force as in the work [24], namely In order to model the unsteady fluid flow within the matrix phase Ω m , we exploit the following Darcy-Brinkman equation [1,10,24]: royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 For γ = v, m, v g is the velocity of the fluid in the region V g , p g is the pressure in the region V g , b g is the inhomogeneous external force density in the region V g , Kg ðxÞ is the hydraulic conductivity tensor in the region V g , μ e is the effective viscosity, and ρ 0 is the fluid density.For γ = v, m, we suppose that the hydraulic conductivity tensor is symmetric and positive definite: , 8 a = 0 : a Á Kg ðxÞa .0: ð2:3Þ Our starting point consists of both Darcy and Darcy-Brinkman equations, assuming that the pore structure is homogenized in both compartments.In this scenario, the hydraulic conductivity tensor Kg ðxÞ functions as a means to represent the essential microscale geometric information.
The matrix and the vessels are coupled via the following interface conditions: where G ¼ @V m > @V v , n is the outer normal to Ω m , τ is any tangential vector to Γ, pðtÞ is a function that depends only on time and α is a constant that depends on the physico-chemical properties of the interface.The first interface condition of (2.4) is related to the normal component of the velocity; if we take p ¼ sðp m À p v Þ, we obtain the well-known Starling equation [62,63], used to describe the fluid exchange between two regions separated by a porous membrane, where σ is the Staverman reflection coefficient, π v the osmotic pressure of Ω v and π m the osmotic pressure of Ω m .For simplicity, in this work, we assume that the osmotic pressures π v and π m depend only on time, but they can depend also on space [3].The parameter L p is typically experimentally measured and depends on the geometry and porosity/leakage of the vessels' walls Γ.We consider this particular interface condition because we aim to apply this model to the fluid flow of a lymph node, but our formulation remains valid for other applications too.The second equation of (2.4) is the Beavers-Joseph-Saffman boundary condition [64], an interface condition on the tangent component of the velocity introduced in [65,66].This interface condition introduces a slip velocity, which is proportional to the normal component of the velocity gradient of the fluid near the boundary.This interface condition is often used to describe the connection between free-fluid and porous regions, as well as an interface condition for dual-porosity media [2,67].Moreover, linking the Darcy-Brinkman equation with the Darcy equation necessitates conditions reliant on the actual velocity.Therefore, the Beavers-Joseph conditions offer a more comprehensive grasp of the underlying physics compared to merely imposing a zero tangent velocity.The initial condition is where v m,0 (x) need to satisfy r Á v m,0 ðxÞ ¼ 0 and must be compatible with the interface conditions (2.4).The non-dimensional form of the Darcy equation, the Darcy-Brinkman equation, and the interface conditions for the domain V g , with γ = m, v, are as follows.We denote with a prime symbol the following non-dimensional quantities: royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 denoting the separation between two vascularized regions.Instead of delving into the intricate particulars of individual vessels, the vascular network region is conceptualized as a geometric domain denoted as Ω v , consisting of interconnected cylinders with radius r c (illustrated in figure 1), where the Darcy equation is applicable.Accordingly, d is precisely characterized as the distance between two neighbouring cylinders within this model representation.Substituting (2.6) into (2.2),we obtain, by neglecting the primes where As d represents the fine-scale length in our problem, the most natural choice as a representative fine-scale conductivity, which we denote as K ref , resides, given its physical dimensions, in choosing d 2 =m, i.e.
Substituting these relations into (2.1) and (2.8), we obtain the non-dimensional equations We also need to non-dimensionalize the interface conditions (2.4).In particular, we adopt the same distinguished limit embraced in [3,68] which ensures that the blood flux stays finite when the lengthscale separation that exists in the system becomes more and more pronounced.This is equivalent to assuming where L p ¼ L p mL 2 =d 3 [3].Moreover, we close our problem using periodic boundary conditions at the boundary @V n G.

The multiscale formulation
The goal of this section is to derive a macroscale model for the continuum system of equations (2.10), (2.5), (2.11) and (2.12), using the asymptotic homogenization technique [8,9,69].We assume that there is a clear separation between the spatial fine scale d and the coarse scale L, which means that the quantity e defined in (2.7) is small: To achieve spatial scale decoupling, we introduce a fresh local variable where x represents the coarse-scale spatial coordinates and y represents the fine-scale spatial coordinates: they have to be considered independent in a formal way.We assume that p g , v g , K g and b g (where γ = m, v) depend on both x and y.We assume two main hypotheses concerning the geometry of the multiscale problem: local periodicity and macroscopic uniformity.Local periodicity means that p g , v g , K g and b g are y-periodic.By making this assumption, we are able to focus our study on a restricted portion of the fine-scale domain.Instead, macroscopic uniformity means ignoring geometric variations within the cell structure and inclusions concerning the coarse-scale variable x, i.e. the microstructure is unique.Then we may consider the royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 utilization of a single periodic cell, denoted as Ω γ , for each macroscale point x, and ð3:2Þ The differential operator becomes ð3:3Þ we now employ a power series expansion for e for the variables of our problem as follows (where γ = m, v): and he @v e m @t ðx, y, tÞ ¼ Àer x p e m ðx, y, tÞ À r y p e m ðx, y, tÞ and if we substitute them into the interface conditions (2.12) and the initial condition (2.5), we obtain and v e m ðx, y, 0Þ ¼ v e m,0 ðx, yÞ in V m : ð3:11Þ

Coefficients of order e 1
If we collect the terms of order e 1 from equations (3.8) to (3.11), we have The ansatz (3.28) solves the problem (3.26)-(3.27)(it is a solution up to a y-constant function), provided that the auxiliary vector and scalar fields h v and hv solve the following cell problems: and v ðx, y, tÞ] Á n on G Â ½0, T: To ensure the solution uniqueness, we need to impose that hh v ðx, y, tÞi Vv ¼ 0 and h hv ðx, y, tÞi Vv ¼ 0, where the average operator hð †Þi Vg is defined as We apply the Fourier transform defined by fðtÞe À2pitv dt ð3:34Þ to system (3.33), and we obtain (note that the initial condition can be written as dðtÞv ð0Þ m,0 ðx, yÞ, where δ(t) is the Dirac delta) We have that (3.39) and (3.40) are the unique solutions of the auxiliary Darcy-Brinkman problem (3.33) provided that the auxiliary second rank tensor Qg , the auxiliary vectors ĥg , qg and the auxiliary scalar function ĥg solve the following cell problems: royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 If we apply the inverse Fourier transform F À1 to (3.41) and (3.42), we obtain where e 2 is a small number >0, so that h m ðx, y, tÞ ¼ @ h m @t ðx, y, tÞ, Q m ðx, y, tÞ ¼ @ Q m @t ðx, y, tÞ; ð3:46Þ substituting these equations into (3.43) and integrating over time, we have

Derivation of the macroscopic model
We now apply the average operator (3.31)-(3.48)We apply the average operator and the divergence theorem to equation (3.20) and, following the same computations as before, we obtain and using the divergence theorem and the interface conditions (3.23) m ðx, tÞ À p ð0Þ v ðx, tÞ À pðtÞ], ð4:6Þ where we used the fact that n v = −n; hence Equation (4.11) is the well-known diffusion problem related to the Darcy equation, where we can find additional terms that are related to the fluid exchange between phases and the multiscale forces [1,24].When the multiscale force b v is zero, the unique solution hv ðx, y, tÞ of the system (3.30) is zero, and in this case, equation (4.11) becomes the Darcy equation with fluid exchange between phases as derived and solved in [2,3,68].On the other hand, equation (4.10) is the Darcy equation diffusion problem with memory effect, with additional terms related to the multiscale forces [24], the fluid exchange between phases and the multiscale initial condition.We note that if the multiscale force b m and the initial condition v m,0 are both zeros, the unique solution qm of the system (3.44) is zero.Taking into account the time dependence of the problem, we obtain a very different model from the previous one [1-3,24,68].However, even when ignoring the contributions related to the external volume loads and royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 the initial condition, the final model that we have obtained differs from the one in [8,26] due to the coupling and the fluid exchange between the Darcy equation with memory effect and the classical Darcy equation and due to the Darcy-Brinkman type cell problem which needs to be solved to compute the hydraulic conductivity h@ Q=@ti Vm for the matrix compartment Ω m .Now we dimensionalize equations (4.1), (4.3), (4.4) and (4.7).We write the equations in dimensional form because, with regard to their application to the lymph node, this makes it easier to interpret the results and compare them with other studies on the same topic.We have where |Ω| is the lymph node total volume, jV tot m j is the phase m total volume, jV tot v j is the phase v total volume and S tot is the total vessel surface.Hence equations (4.1), (4.4), (4.5) and (4.8) are (in dimensional form)

Explicit solution
In this section, we explicitly solve the problems given in §4.
For simplicity, we suppose v m,0 = 0, b m = 0, b v = 0; hence the unique solution of problems (3.30) and (3.44) is zero.Moreover, we assume the isotropy of both the porous media, which means Vv and d 2 =mhQ m ðy, tÞi Vm , respectively.We have that K v is constant in space due to the geometry and the hypotheses used, and it is found solving the cell problem (3.29) using COMSOL Multiphysics (see §6 for more details about the numerical simulations of the cell problem).We have that the dimensional value for the vessel hydraulic conductivity K v d 2 /μ is computed using the Kozeny-Carman formula, i.e. ð1=c 0 ÞðjV tot v j=S tot Þ 2 ; see table 1  We find Qm ðy, vÞ by solving the cell problem (3.41) using COMSOL Multiphysics, with α = 1 (see §6 for more information about the numerical simulations of the cell problem and figure 2 for the results).
As we are working in the frequency domain, we apply the Fourier transform (3.34) to (4.9)-(4.9).We consider a spherical domain Ω with radius R, denoting by r the radial coordinate, θ the polar coordinate and ϕ the azimuthal angle.Moreover, we assume axisymmetry with respect to the azimuthal angle ϕ.Hence, we obtain the following problem:  Figure 2. The Fourier transform of the dimensional average hydraulic conductivity in mm 3 s mg À1 calculated numerically with physiological values of the lymph node (table 1).

Cell problem numerical simulations
In this section, we discuss the numerical simulations used to find the solutions of the cell problems (3.29), (3.30), (3.41) and (3.42).We can see the geometry of the cell problems using the lymph node physiological data found in [1, appendix B] and summarized in table 1 in figure 1.
As we mentioned in §5, we suppose u m,0 = 0, f m = 0, f v = 0, which means that the solutions of the problems (3.30) and (3.42) are zeros.Moreover, we suppose that both porous media are isotropic, which means For the solution of the cell problem in the phase Ω v , we refer to [1, appendix C].For the phase Ω m , we solve the problem (3.41) with α = 1 using the steady Brinkman equation module in COMSOL Multiphysics.We use the PARDISO solver and the P 3  2 À P 1 finite-element discretization for the velocity and the pressure, respectively.Moreover, we perform a parametric sweep analysis varying ω to obtain the solution for different frequencies.After that, we apply the average operator to obtain h Qm ðy, vÞi Vm .We obtain the dimensional hydraulic conductivity K m ðvÞ shown in figure 2. If we take the values of the permeability h Qm ðy, vÞi Vm and we perform an inverse Fourier transform (with the MATLAB built-in command ifft), we find the dimensional average hydraulic conductivity shown in figure 3.As expected, we have a permeability that decreases in time [8].The decreasing imaginary part of h Qm ðy, vÞi Vm in figure 2 is connected to the decrease in time of the permeability (due to the memory effect of the macroscopic Darcy's Law).If we compare the ifft result with the solution of the cell problem (3.45), we have a maximum relative error of 0.188%.Moreover, we have that the real part of the permeability is similar to the one found in the steady case in [1].We believe these values to be realistic because, at time zero, we have the same permeability values we found in the steady setting [1].Additionally, as time increases, the permeability tends to zero very fast, reflecting its royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 diminishing influence on the flow dynamics at the time we are studying our flow, a characteristic commonly observed in memory-type equations [8].
To study in more detail the mesh of the previous solution, we perform an adaptive mesh refinement study.After this process, we compare the two K m ðvÞ that we found, and we obtain, for all the values of ω, a maximum relative error of 0.06%.
To test the Brinkman equation module in Comsol with complex numbers, we run some simulations in a cylindrical geometry and we compare the solution to the following explicit solution (found in [1]): Taking the average, we obtain where l c is the cylinder length, Ω cyl is the cylinder volume, K Ã is the permeability, rc is the cylinder radius, μ Ã is the viscosity, J 0 and J 1 are the Bessel functions of the first and second kind, respectively.We run the numerical simulations and compare the solution to the explicit one (6.2) with the following (nondimensional) data: l c = 1, rc = 7.7 × 10 −3 , K Ã = 6.66 × 10 −6 + i, 6.66 × 10 −6 + 10i, 6.66 × 10 −6 + 100i, μ Ã = 1.We found that the maximum relative error between the numerical and the explicit solution is 0.67% for the real part and 1.9% for the imaginary part.

Study of the fluid flow in a lymph node
In this section, we study the fluid flow inside the lymph node using the model results obtained in the previous sections of this work.In particular, we use the explicit solution of §5 applied to the porous region of the node (the LC) because all the lymph node vascularization is situated in this region [47][48][49].Moreover, we assume that the lymph node has a spherical shape, and we use physiological lymph node data inspired by a mouse popliteal lymph node [1,47,59].The explicit solution of §5 is implemented in MATLAB.To compute the inverse Fourier transform, we use the MATLAB built-in command ifft, that is, the inverse fast Fourier transform.
We have that Ω v and Ω m are the portions of the domain that represent the blood vessels and the FRC network, respectively.Choosing pðtÞ ¼ sðp m ðtÞ À p v ðtÞÞ, we describe the fluid exchange between the phases mentioned above using the well-known Starling equation [62,63].Moreover, all the physiological data used are summarized in table 1, and the meaning of these data is explained in [1, appendix B].  royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 As we can see from the previous here the hydraulic conductivity of the phase Ω m depends on time due to the memory term in the Darcy equation (4.1).This is very different from the previous works on the lymph node [1,47,53,54,[57][58][59], as it allows for a comprehensive analysis of the lymph's time behaviour through a rigorous homogenization method.
Using physiological data obtained from the literature on the lymph node, we find that the parameter that governs the time dependency of our problem is η ≈ 0.1, in line with the Womersley number found in the lymphatic system [36].
At the boundary, we have the following boundary conditions: p m ðR, z, tÞ ¼ p m ðz, tÞ and p v ðR, z, tÞ ¼ p v , where p m ðz, tÞ is a general function in ζ and t, and p v is the mean of the blood vessel pressure data (in general constant).
As we mentioned in [1], the variation with respect to z ¼ cos u of the boundary condition p m is essential to mimic the pressure distribution in the SCS.The precise pressure distribution of the SCS is not well known so, to circumvent this lack of data in the literature and inspired by [45], we take a linear relation between the pressure and z ¼ cos u.This linear relation connects the maximum value of the pressure p m,max ¼ 3:9 mmHg % 5:2 Â 10 5 mPa with the minimum value of p m,min ¼ 3 mmHg % 4 Â 10 5 mPa, taken from the numerical results of [47].Moreover, to describe the pulsatile inflow, we take a pressure distribution that varies over time in the following way: We can see the behaviour of the boundary condition (7.1) with respect to ζ at different times in figure 4. We use the following physiological parameters: σ = 0.88, π v − π m = 1.02 × 10 6 mPa, L p ¼ 5:475 Â 10 À10 mm s À1 mPa À1 and p v ¼ 6:66 Â 10 5 .We can see the behaviour of the interstitial pressure p m over time in figure 5.As we can see, the minimum of the interstitial pressure p m increases and moves from the centre of the node to the lower part of the node.This is due to the fact that, as time passes, the boundary pressure distribution (7.1) remains linear but the maximum of the pressure increases, and this effect combines with the effect of the fluid exchange between phases.This phenomenon highlights the importance of the time dependency of the flow inside the node.We note that the Darcy equation linearly relates the fluid discharge to the pressure gradient, so the lymph moves accordingly to the pressure, which means that a sink term is represented by a lower pressure region in the node.
As expected and in accordance with [1], increasing L p results in a decrease of p m and an increase of p v at the centre of the node.We can see this behaviour of the interstitial pressure p m in the contour plot of figure 6 at time t = 1 s.As we can see, the results are in line with the steady one [1], which means that the increase of L p decreases and moves the minimum of p m towards the centre of the node.royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 in [1]; from this value to higher values of p v , the maximum of p m starts to move towards the centre of the node.This behaviour is accordance with the findings in [1].Moreover, we can have flow inversion for only certain values of the time; for instance, we can see figure 9, where we plot the blood vessel pressure p v at different times with p v ¼ 1:35 Â 10 6 mPa.For the initial times, we have a flow inversion, which means that the fluid moves from phase Ω v to phase Ω m , resulting in a region of minimum pressure at the centre of the node for p v .As time passes, the boundary condition (7.1) increases the pressure in the upper region of the lymph node (near θ = 0), and this results in a region of higher pressure zone near θ = 0 and the region of lower pressure zone moves near θ = π.
Increasing Δπ means an increase in the concentration difference between the blood vessels and the FRC phase, and this results in a decrease of the minimum of p m and moves it to the centre of the node, and an increase of the maximum of p v .We can see a contour plot of p m at t = 1 s varying Δπ in figure 10.
From the fact that the boundary condition p m ðz, tÞ can be chosen as we want, we solve our explicit solution using the solution we found in our previous paper [59] using the stream function approach.We refer to [59] for more details about the computations and the data of this pressure distribution; this solution is calculated in a div-free setting, where we fix p m ðR 2 , À 1Þ ¼ 6:18 Â 10 5 mPa and with a time pulsation of the form ð1 À cosðptÞÞ=2.In figure 11, we can see the boundary pressure distribution.The fast increment near ζ = 1 represents the inlet condition, and the fast decrement near ζ = −1 represents the outlet condition.In figures 12 and 13, we can see the pressure and the velocity distribution over time using the parameters p v ¼ 1:06 Â 10 6 mPa, π v − π m = 1.02 × 10 6 mPa, L p = 5.475 × 10 −11 mm s −1 mPa −1 , and the above-mentioned boundary pressure.As we can see, for initial times, we have a similar pressure distribution to that we found in the previous case but, as time passes, it becomes more evident the higher pressure near the inlet and the lower pressure near the outlet.We can see this behaviour for the velocity magnitude too, where at t = 1 s, we can see a higher difference between the inlet-outlet velocity and the velocity at the centre of the node.
Moreover, we have that the velocities that we found inside the LC are in agreement with the literature [47,55,[79][80][81].We found a higher velocity (at time t = 1 s) with respect to the steady case [1]; moreover, we have that the pressure change (at time t = 1 s) between the upper region (near the inlet condition) and the lower region (near the outlet condition) is higher here too.
Finally, we have performed a parameter sensitivity analysis and found that on varying the average hydraulic conductivity, the resulting pressure (interstitial and blood pressure) exhibits only a relatively royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 small variation.In particular, a variation of 1, 2, and 5% of the average hydraulic conductivity corresponds to the relative variations of the interstitial pressure of 0.012, 0.0238 and 0.0578%.

Conclusion
In this paper, we the asymptotic homogenization technique in a time-dependent setting, starting from the equations and interface conditions (2.1), (2.2), (2.4), assuming both local periodicity and macroscopic uniformity.This model is a non-trivial extension of our previous model [1]; here, we have considered a timedependent Darcy-Brinkman equation, which results in a Darcy equation with memory at the macroscale for the phase Ω m .The time dependency is considered for the pulsation behaviour of the lymph [36,37,61,82] and allows us to study the fluid behaviour inside the lymph node in more detail [39,40,59].With lymph node physiological data, the characteristic time is η ≈ 0.1: this is in agreement with the Womersley number found for the lymphatic system [36].
This model has been designed to be applied to the flow of lymph within the lymph node, but the derivation of the model has been intentionally kept as general as possible to be applicable to other problems as well.
After the derivation of the macroscopic equations that described the time-dependent fluid flow ( §4), in §5, we have found the explicit solution of the proposed model in a spherical domain, using the computations of the solution that we found in [1] and the properties of the Fourier transform.Subsequently, we have studied the fluid and pressure distribution within the lymph node using the royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983 above-mentioned explicit solution and physiological data inspired by an idealized spherical mouse popliteal lymph node.Incorporating the temporal component into the proposed model has allowed us to study how lymph pulsations influence these quantities within a lymph node.Regarding this multiscale formulation, our primary focus has been directed towards the porous region of the lymph node (the LC).Moreover, we have placed particular emphasis on studying the fluid exchange that occurs between the interstitial space of the lymph node and the blood vessels, which are exclusively present in this specific part of the node [47][48][49].We have analysed how various parameters influence fluid absorption and pressure (i.e.velocity) over time, and the obtained results align with findings documented in the literature [1,47,[73][74][75]83].
Let us now explore considerations that could enhance the model in the future.In this work, we solely focused on the LC, but it would be interesting to couple this model with the fluid flow in the SCS to describe the fluid flow in the entire lymph node.
The Beavers-Joseph-Saffman interface condition (2.4) used in this work is initially established in a two-dimensional setting, and expanding it to three dimensions presents a significant and complex challenge [2,[84][85][86][87][88].Moreover, it would be interesting to study how the physico-chemical properties of the interface affect the solution [67].
When we applied our model to the lymph node, we assumed negligible forces to simplify our model and due to a lack of data regarding them in the literature.However, it is important to note that in practical situations, such forces can play a significant role, especially when using electromagnetic fields, as seen in [89,90] within the context of cancer hyperthermia.Hence, when we have access to physiological data, it becomes crucial to account for the impact of inhomogeneous volume loads, as discussed in [24].We can say the same for the velocity initial condition, which we have assumed to be zero.In general, we can have more general fluid behaviour as an initial condition, in particular in drug delivery applications [25].
In our model, we have taken into account the time variation of the concentration of protein inside the node, but we have supposed it to be constant due to a lack of precise data.Incorporating the temporal and spatial behaviour of protein concentration, as discussed in [3], and together with physiological data, would be a compelling and valuable direction to explore.
In this work, we have assumed a rigid porous matrix to keep the model as simple as possible and due to the lack of biological data regarding this problem; a possible extension of this model is to take into account a deformable matrix that interacts with the lymph flow inside the node, for example by considering the modelling framework developed in [91,92].
Finally, due to our explicit analysis, we assumed a spherical lymph node.In general, lymph nodes have a more complex geometry, typically an ellipsoid shape [47,53,54].If we can access more realistic data regarding their shape through medical imaging [53,54], our model could be used for more realistic numerical simulations to make physiologically meaningful predictions in the future.

Figure 1 .
Figure 1.A sketch of a periodic portion of the domain comprising the vessels Ω v (a) and matrix Ω m (b).

Figure 3 .
Figure3.The inverse Fourier transform of the dimensional average hydraulic conductivity in mm 3 s mg À1 calculated numerically using the command ifft to the solution plot in figure2with physiological values of the lymph node (table1).

Figure 7 .Figure 8 .
Figure 7.The variation of the interstitial pressure p m (in mPa) varying the mean blood vessel pressure p v (in mPa) at time t = 1 s, with the parameters in table 1, π v − π m = 1.02 × 10 6 mPa, L p = 5.475 × 10 −10 mm s −1 mPa −1 and the boundary conditions

Table 1 .
Physiological and estimated parameters; see[1, appendix]for more details.¼ cos u, M(ω) = M v + M m (ω), I n is the modified Bessel function of the first kind of order n, P n is the Legendre polynomial of the first kind of order n, and (the following relations are obtained from [1, appendix A]) royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231983