Mesoscopic Approach Research Articles

A Mesoscopic Approach for the Numerical Simulation of a Mass Concrete Structure Construction Using Post-Cooling Systems

This study introduces an innovative numerical approach to simulate the construction of large concrete structures incorporating post-cooling systems employing the finite element method (FEM). The proposed methodology integrates critical construction parameters, including temperature control mechanisms, while accounting for concrete hydration and environmental conditions. Compared to traditional discrete models, this approach provides similar accuracy with substantially reduced computational costs, enhancing predictive capabilities in the thermal analysis of mass concrete. The method was applied to simulate the construction of a water intake pillar at the Tocoma hydroelectric plant in Venezuela, where the simulated results closely matched in situ temperature measurements. The findings highlight the method’s efficiency and accuracy in simulating post-cooling systems, offering a practical solution for improving the safety and cost-effectiveness in large-scale concrete construction projects.

Thermally induced oscillatory rarefied gas flow inside a rectangular cavity

Thermally induced oscillatory rarefied gas flow inside a two-dimensional rectangular cavity is investigated based on the hybrid macro-/mesoscopic scheme. The effects of the Knudsen (Kn) numbers and the oscillation frequency of lid temperature on the flow parameters are analyzed. The Shakhov model equation is solved numerically based on the mesoscopic approach in the near-wall region, and the macroscopic approach is adopted in the bulk flow region to reduce the computational cost. To close the numerical iteration procedure, the velocity distribution functions serving as the pseudo boundary between macroscopic and mesoscopic methods are reconstructed using the high-order Hermite polynomials. Numerical simulations demonstrate that the temperature profile at the central vertical of the cavity predicted by the hybrid method is in good agreement with results from the mesoscopic method, with a maximum error of 0.23%. In addition, the computational memory cost can be saved up to about 69.91%. The hybrid approach is able to capture the nonlinear phenomenon in the thermally induced oscillatory rarefied gas flow under high Kn numbers, where the horizontal velocity no longer obeys the law of periodic oscillating cosine function, and the rise time of the horizontal velocity is much longer than the fall time. The thickness of the viscous penetration layer and the disturbed region increases as the Kn number increases and decreases as the Strouhal number increases.

Critical surface adsorption of confined binary liquids with locally conserved mass and composition

Close to a solid surface, the properties of a fluid deviate significantly from their bulk values. In this context, we study the surface adsorption profiles of a symmetric binary liquid confined to a slit pore by means of molecular dynamics simulations; the latter naturally entails that mass and concentration are locally conserved. Near a bulk consolute point, where the liquid exhibits a demixing transition with the local concentration as the order parameter, we determine the order parameter profiles and characterise the relevant critical scaling behaviour, in the regime of strong surface attraction, for a range of pore widths and temperatures. The obtained order parameter profiles decay monotonically near the surfaces, also in the presence of a pronounced layering in the number density. Overall, our results agree qualitatively with recent theoretical predictions from a mesoscopic field-theoretical approach for the canonical ensemble.

The diffuse interface description of fluid lipid membranes captures key features of the hemifusion pathway and lateral stress profile.

Topological transitions of lipid membranes are ubiquitous in key biological processes for cell life, like neurotransmission, fertilization, morphogenesis, and viral infections. Despite this, they are not well understood due to their multiscale nature, which limits the use of molecular models and calls for a mesoscopic approach such as the celebrated Canham-Helfrich one. Unfortunately, such a model cannot handle topological transitions, hiding the crucial involved forces and the appearance of the experimentally observed hemifused intermediates. In this work, we describe the membrane as a diffuse interface preserving the Canham-Helfrich elasticity. We show that pivotal features of the hemifusion pathway are captured by this mesoscopic approach, e.g. a (meta)stable hemifusion state and the fusogenic behavior of negative monolayer spontaneous curvatures. The membrane lateral stress profile is calculated as a function of the elastic rigidities, yielding a coarse-grained version of molecular models findings. Insights into the fusogenic mechanism are reported and discussed.

A mesoscopic approach for nanoscale evaporation heat transfer characteristics

We propose a mesoscopic approach for investigating steady/transient nanoscale evaporation heat transfer characteristics, whose kinetic nature makes it an ideal tool to model the evaporation kinetics at the liquid-vapor interface and the microscopic wall-fluid interactions. Simulation of the evaporation of a flat nanoscale thin liquid film demonstrates its capability of resolving kinetically limited evaporation and liquid film adsorption. For the simulation of an evaporating meniscus in a nanochannel, our proposed mesoscopic approach considers the conservation of mass, momentum and energy of both liquid and vapor phases, thereby naturally incorporating both the multi-dimensional heat transfer effect and vapor transport resistance in nano-confined space. We demonstrate that solid-fluid interaction plays dominant roles on the interfacial transport during nanoscale evaporation, and vapor transport resistance has nonnegligible influences on the evaporation heat/mass transport in nano-confined spaces. By handling statistical behavior of molecule ensembles, this mesoscopic approach not only preserves microscopic physics but also has high computational efficiencies, enabling the simulation of space and time domains on the practical application level. This work paves the way for quantitative investigations of steady/transient nanoscale evaporation heat transfer characteristics using mesoscopic approach.

Experimental and mesoscopic simulation study of size effect on splitting tensile performance of spontaneous combustion coal gangue concrete

The spontaneous-combustion coal gangue concrete (SCGAC) demonstrates different mechanical behavior compared to ordinary concrete owing to its distinct porous structure and the presence of fibrous materials in its coarse aggregates. Presently, research on the macroscopic mechanical properties of SCGAC is reasonably comprehensive; however, it falls short of fully understanding its size effect. This study examines the splitting tensile strength of cubic SCGAC specimens with three strength grades, i.e., C30, C35, and C40, and three sizes, i.e., 100 mm, 150 mm, and 200 mm; proposes a mesoscopic modeling approach based on the computer vision; analyzes the intrinsic mechanisms of splitting failure and the size effect of different strength grades on the SCGAC by utilizing mesoscopic simulation methods. Subsequently, an analysis of the size effect on the splitting tensile strength of SCGAC is conducted based on classical size effect laws. The results reveal the following: The method proposed in this study for mesoscopic modeling is recognized for its low cost, high efficiency, and accurate simulations, fulfilling the required standards for mesoscopic mechanical analysis. SCGAC, as a heterogeneous material, demonstrates significant size effect phenomena that become more pronounced with enhancements in strength grades. The Carpinter fractal size effect law is ideally suited to characterize the size effect patterns in SCGAC. This paper enhances the comprehension of the size effect of SCGAC through a combination of macroscopic experiments and mesoscopic simulations, thereby facilitating its engineering applications.

Mesoscopic model of cycling trip energy expenditure based on operating modes

Cyclist physical exertion is largely ignored in quantitative travel analysis, partly due to a lack of appropriate tools. Microscopic models of second-by-second energy expenditure based on equations of motion are data intensive and cannot be applied to hypothetical routes (such as needed for route choice modelling). Macroscopic models of aggregate energy expenditure based on a fixed assumed energy intensity are insensitive to traveller, trip, and contextual factors that are relevant for behavioural research and policy analysis (such as bicycle type or trip purpose). Building on concepts from motor vehicle emissions analysis, this paper proposes a mesoscopic approach to model cycling trip energy expenditure based on the distribution of travel time across discrete states of motion (“operating modes”) for different classes of traveller and trip (“model segments”). We aim to answer two key questions for model implementation: 1) which variables most effectively classify trips into model segments and 2) what operating mode definition most consistently characterizes cycling energy expenditure within model segments? We also evaluate the precision of the mesoscopic model relative to cycling energy estimates from microscopic and macroscopic models. Applied to a dataset of naturalistic cycling trips in Vancouver, Canada, the proposed mesoscopic model with six model segments based on 3 segmenting variables (rider gender, electric-assist bicycle, and high or low speed tier) explains up to 28 % of the variance in trip-level energy estimates from the microscopic model (within around 35 W of the microscopic estimates, on average). Further research to develop cycling trip energy models for general application is discussed.

Colloid and nanoparticle-driven phase behavior in weakly perturbed nematic liquid crystals

We study theoretically weakly colloidal particle or nano-particle perturbed nematic liquid crystal (LC) order. We consider cases where the particles are essentially homogeneously dispersed in LC matrix. A mesoscopic Landau-de Gennes approach is used where we focus on the temperature variations in the amplitude order parameter field for different particle volume concentrations ϕ. In the dilute regime, we describe the nematic order in terms of interface and volume amplitude order fields. This modeling suggests that the temperature-driven order-disorder crossover exhibits a two-step transformation. In the case that particles locally support LC order, on decreasing temperature firstly the interface order is established, followed by the build-up of volume order. The latter appears via the 1st order phase transition. On the other hand, the interface order growth is formed either gradually or via the 1st order interfacial transition. In the regime of relatively high values of ϕ the system can be effectively described in terms of a single amplitude order parameter field. In this case, the temperature-driven phase behaviour is qualitatively similar to the interface order temperature response in the diluted regime. We analytically express phase transition temperatures and determine critical points in the temperature-interface interaction phase space. Key novelties of the paper are the determination of conditions for particle-LC interface-driven discontinuous change of nematic ordering and nanoparticle-enabled critical point condition in pressure-driven order-disorder phase transformations. Experimental evidences supporting the predicted phase behavior are presented.

Mesoscopic Kinetic Approach of Nonequilibrium Effects for Shock Waves.

A shock wave is a flow phenomenon that needs to be considered in the development of high-speed aircraft and engines. The traditional computational fluid dynamics (CFD) method describes it from the perspective of macroscopic variables, such as the Mach number, pressure, density, and temperature. The thickness of the shock wave is close to the level of the molecular free path, and molecular motion has a strong influence on the shock wave. According to the analysis of the Chapman-Enskog approach, the nonequilibrium effect is the source term that causes the fluid system to deviate from the equilibrium state. The nonequilibrium effect can be used to obtain a description of the physical characteristics of shock waves that are different from the macroscopic variables. The basic idea of the nonequilibrium effect approach is to obtain the nonequilibrium moment of the molecular velocity distribution function by solving the Boltzmann-Bhatnagar-Gross-Krook (Boltzmann BGK) equations or multiple relaxation times Boltzmann (MRT-Boltzmann) equations and to explore the nonequilibrium effect near the shock wave from the molecular motion level. This article introduces the theory and understanding of the nonequilibrium effect approach and reviews the research progress of nonequilibrium behavior in shock-related flow phenomena. The role of nonequilibrium moments played on the macroscopic governing equations of fluids is discussed, the physical meaning of nonequilibrium moments is given from the perspective of molecular motion, and the relationship between nonequilibrium moments and equilibrium moments is analyzed. Studies on the nonequilibrium effects of shock problems, such as the Riemann problem, shock reflection, shock wave/boundary layer interaction, and detonation wave, are introduced. It reveals the nonequilibrium behavior of the shock wave from the mesoscopic level, which is different from the traditional macro perspective and shows the application potential of the mesoscopic kinetic approach of the nonequilibrium effect in the shock problem.

Development and application of hybrid AIMD/cDFT simulations for atomic-to-mesoscale chemistry.

Many important chemical processes involve reactivity and dynamics in complex solutions. Gaining a fundamental understanding of these reaction mechanisms is a challenging goal that requires advanced computational and experimental approaches. However, important techniques such as molecular simulation have limitations in terms of scales of time, length, and system complexity. Furthermore, among the currently available solvation models, there are very few designed to describe the interaction between the molecular scale and the mesoscale. To help address this challenge, here, we establish a novel hybrid approach that couples first-principles plane-wave density functional theory with classical density functional theory (cDFT). In this approach, a region of interest described by abinitio molecular dynamics (AIMD) interacts with the surrounding medium described using cDFT to arrive at a self-consistent ground state. cDFT is a robust but efficient mesoscopic approach to accurate thermodynamics of bulk electrolyte solutions over a wide concentration range (up to 2M concentrations). Benchmarking against commonly used continuum models of solvation, such as SMD, as well as experiments, demonstrates that our hybrid AIMD-cDFT method is able to produce reasonable solvation energies for a variety of molecules and ions. With this model, we also examined the solvent effects on a prototype SN2 reaction of the nucleophilic attack of a chloride ion on methyl chloride in the solution. The resulting reaction pathway profile and the solution phase barrier agree well with experiment, showing that our AIMD/cDFT hybrid approach can provide insight into the specific role of the solvent on the reaction coordinate.

Biaxial Structures of Localized Deformations and Line-like Distortions in Effectively 2D Nematic Films

We numerically studied localized elastic distortions in curved, effectively two-dimensional nematic shells. We used a mesoscopic Landau-de Gennes-type approach, in which the orientational order is theoretically considered by introducing the appropriate tensor nematic order parameter, while the three-dimensional shell shape is described by the curvature tensor. We limited our theoretical consideration to axially symmetric shapes of nematic shells. It was shown that in the surface regions of stomatocyte-class nematic shell shapes with large enough magnitudes of extrinsic (deviatoric) curvature, the direction of the in-plane orientational ordering can be mutually perpendicular above and below the narrow neck region. We demonstrate that such line-like nematic distortion configurations may run along the parallels (i.e., along the circular lines of constant latitude) located in the narrow neck regions of stomatocyte-like nematic shells. It was shown that nematic distortions are enabled by the order reconstruction mechanism. We propose that the regions of nematic shells that are strongly elastically deformed, i.e., topological defects and line-like distortions, may attract appropriately surface-decorated nanoparticles (NPs), which could potentially be useful for the controlled assembly of NPs.

Therapy-induced modulation of tumor vasculature and oxygenation in a murine glioblastoma model quantified by deep learning-based feature extraction

Glioblastoma presents characteristically with an exuberant, poorly functional vasculature that causes malperfusion, hypoxia and necrosis. Despite limited clinical efficacy, anti-angiogenesis resulting in vascular normalization remains a promising therapeutic approach. Yet, fundamental questions concerning anti-angiogenic therapy remain unanswered, partly due to the scale and resolution gap between microscopy and clinical imaging and a lack of quantitative data readouts. To what extend does treatment lead to vessel regression or vessel normalization and does it ameliorate or aggravate hypoxia? Clearly, a better understanding of the underlying mechanisms would greatly benefit the development of desperately needed improved treatment regimens. Here, using orthotopic transplantation of Gli36 cells, a widely used murine glioma model, we present a mesoscopic approach based on light sheet fluorescence microscopic imaging of wholemount stained tumors. Deep learning-based segmentation followed by automated feature extraction allowed quantitative analyses of the entire tumor vasculature and oxygenation statuses. Unexpectedly in this model, the response to both cytotoxic and anti-angiogenic therapy was dominated by vessel normalization with little evidence for vessel regression. Equally surprising, only cytotoxic therapy resulted in a significant alleviation of hypoxia. Taken together, we provide and evaluate a quantitative workflow that addresses some of the most urgent mechanistic questions in anti-angiogenic therapy.

Comprehensive experimental and numerical validation of Lattice Boltzmann fluid flow and particle simulations in a child respiratory tract

The numerical simulation of inhaled aerosols in medical research starts to play a crucial role in understanding local deposition within the respiratory tract, a feat often unattainable experimentally. Research on children is particularly challenging due to the limited availability of in vivo data and the inherent morphological intricacies. CFD solvers based on Finite Volume Methods (FVM) have been widely employed to solve the flow field in such studies. Recently, Lattice Boltzmann Methods (LBM), a mesoscopic approach, have gained prominence, especially for their scalability on High-Performance Computers.This study endeavours to compare the effectiveness of LBM and FVM in simulating particulate flows within a child’s respiratory tract, supporting research related to particle deposition and medication delivery using LBM. Considering a 5-year-old child’s airway model at a steady inspiratory flow, the results are compared with in vitro experiments. Notably, both LBM and FVM exhibit favourable agreement with experimental data for the mean velocity field and the turbulence intensity.For particle deposition, both numerical methods yield comparable results, aligning well with in vitro experiments across a particle size range of 0.1–20 µm. Discrepancies are identified in the upper airways and trachea, indicating a lower deposition fraction than in the experiment. Nonetheless, both LBM and FVM offer invaluable insights into particle behaviour for different sizes, which are not easily achievable experimentally.In terms of practical implications, the findings of this study hold significance for respiratory medicine and drug delivery systems — potential health impacts, targeted drug delivery strategies or optimisation of respiratory therapies.

A large-scale hybrid micro- and mesoscopic simulation approach for railway operation

In this study, the existing agent-based simulation framework MATSim is extended towards a module for the simulation of trains and their interactions with the infrastructure. The proposed railway simulation module railsim offers several advantages over existing railway planning tools: railsim is open-source, easily adaptable, allows for a microscopic, mesoscopic or hybrid micro-mesoscopic simulation setup and is fully integrated in an existing simulation framework. Trains have a spatial dimension along several infrastructure segments. Trains accelerate and decelerate based on vehicle- and infrastructure-specific attributes as well as the available capacity ahead of the vehicle. The interaction of vehicles as well as the modeling of the infrastructure is adapted to the railway system. Trains are enabled to reserve infrastructure elements (fixed block segments, moving blocks), and reservation requests interact with each other. The illustrative Switzerland case study reveals that an overall unique capacity of 3 trains per section and 10 trains per station ensures most trains reach their destinations successfully, while lower capacities lead to insufficient infrastructure and train disruptions.

Optimal management of full train load services in the shunting yard: A comprehensive study on Shunt-In Shunt-Out policies

A key aim of the European Union is to double freight rail traffic by 2050 in order to reduce pollution emissions and alleviate congestion by shifting traffic from roads to rail networks. To accomplish this objective, it is crucial to minimize emissions and high costs associated with shunting yard operations while maintaining an acceptable level of service. This research paper introduces a new MINLP model that optimizes the shunt-in and shunt-out (SISO) operations of wagons in a shunting yard that handles full train load services. Additionally, an efficient Multi-Objective Dijkstra Algorithm (MDA) is proposed to handle simultaneous shunt-out operations in a multi-train scenario. The MINLP model takes a mesoscopic approach, aiming to minimize the number of SISO operations while satisfying strategic and tactical objectives such as wagon fleet size, operational costs, and shunting locomotive emissions. Several versions of the mathematical model are described, each employing different Shunt-In (SI) policies with varying criteria for wagon selection and strong goal orientation. The Multi-Objective Dijkstra Algorithm determines the Multi-Objective Shortest Path between mandatory shunts, considering both the time required for shunting and clustering costs. It provides information on the clusters of wagons that need to be shunted out. To assess the effectiveness of the MINLP model, real train timetables for freight trains in the Bettembourg Eurohub Sud Terminal (Luxembourg) are used, and various KPIs related to tactical and strategic objectives are evaluated. Furthermore, the performance of the Multi-Objective Dijkstra Algorithm is compared to the Shunt-Out sub-model, considering average computation time and solution quality in relation to the MINLP model. The computational results demonstrate that the criteria for wagon selection have a significant impact on the analyzed KPIs.

An efficient greedy heuristic for the real-time train platforming problem

This paper focuses on a special class of real-time railway traffic management problem: efficiently reallocating trains to platforms at stations in case of slight schedule changes or major disruptions. The need for a real-time response with a high level of quality makes this problem particularly challenging.To address this issue, the authors propose a mesoscopic approach that involves preprocessing the data to determine feasible routes and other disruption parameters. They develop a greedy interchange heuristic to solve the mesoscopic real-time train platforming problem, providing high-quality routing and timing decisions within the computational time constraints of real-time management problems.The performance of the proposed heuristic is analyzed through case studies using both synthetic and realistic scenarios from the Spanish railway traffic system. For large instances of the Atocha-Cercanías station case study, the solutions are generated from 5 to 10 times faster by the heuristic algorithm than by the exact method. The authors conclude that the proposed heuristic is a promising solution for real-time train platforming problems.

Mesoscopic discrete modeling of compression and fracture behavior of concrete: Effects of aggregate size distribution and interface transition zone

This paper presents a mesoscopic discrete modeling approach for studying compression and fracture behavior of concrete. A modified piecewise linear sieve curve was proposed to generate irregularly shaped aggregates following specific particle size distributions, and then a random particle model and the parallel bond model were implemented within a 3D framework of discrete element model (DEM). To characterize the interface transition zone (ITZ) which is generally the weak region in concrete, a reduction factor was employed to adjust the strength and consider the effect of ITZ on the mechanical response of concrete. The numerical modeling was calibrated and validated by comparing simulation results with experimental data from compression, splitting tensile, and dog-bone tests. The results demonstrate that the random aggregate model provides accurate representations of the concrete mesostructure, and the discrete model can effectively capture the effects of aggregate size distribution and ITZ on mechanical properties and crack propagation of concrete.

Consistent lifting relations for the initialization of total-energy double-distribution-function kinetic models.

A lifting relation connecting the distribution function explicitly with the hydrodynamic variables is necessary for the Boltzmann equation-based mesoscopic approaches in order to correctly initialize a nonuniform hydrodynamic flow. We derive two lifting relations for Guo etal.'s total-energy double-distribution-function (DDF) kinetic model [Z. L. Guo etal., Phys. Rev. E 75, 036704 (2007)1539-375510.1103/PhysRevE.75.036704], one from the Hermite expansion of the conserved and nonconserved moments, and the second from the O(τ) Chapman-Enskog (CE) approximation of the Maxwellian exponential equilibrium. While both forms are consistent to the compressible Navier-Stokes-Fourier system theoretically, we stress that the latter may introduce numerical oscillations under the recently optimized discrete velocity models [Y. M. Qi etal., Phys. Fluids 34, 116101 (2022)10.1063/5.0120490], namely a 27 discrete velocity model of the seventh-order Gauss-Hermite quadrature (GHQ) accuracy (D3V27A7) for the velocity field combined with a 13 discrete velocity model of the fifth-order GHQ accuracy (D3V13A5) for the total energy. It is shown that the Hermite-expansion-based lifting relation can be alternatively derived from the latter approach using the truncated Hermite-polynomial equilibrium. Additionally, a relationship between the order of CE expansions and the truncated order of Hermite equilibria is developed to determine the minimal order of a Hermite equilibria required to recover any multiple-timescale macroscopic system. Next, three-dimensional compressible Taylor-Green vortex flows with different initial conditions and Ma numbers are simulated to demonstrate the effectiveness and potential issues of these lifting relations. The Hermite-expansion-based lifting relation works well in all cases, while the Chapman-Enskog-expansion-based lifting relation may produce numerical oscillations and a theoretical model is developed to predict such oscillations. Furthermore, the corresponding lifting relations for Qi etal.'s total energy DDF model [Y. M. Qi etal., Phys. Fluids 34, 116101 (2022)10.1063/5.0120490] are derived, and additional simulations are performed to illustrate the generality of our approach.

Comparative Evaluation of Road Pricing Schemes: A Simulation Approach (Australian Perspective)

Road network pricing and congestion charging continue to be debated as efficient instruments to address traffic congestion and emissions. For cities where the schemes have not been implemented yet, the impacts of these schemes are typically evaluated using transport simulation models to understand the impacts and design effective solutions before the schemes are deployed. This paper considers a simulation approach for the city of Melbourne in Australia to investigate the potential impacts of road network pricing on reducing private vehicle travel, road congestion, and vehicle emissions. The study uses a dynamic traffic simulation model developed for Melbourne using the AIMSUN modeling tool, which was extended for modeling road user pricing and congestion charging, including considerations and formulations of distance-based, delay-based, joint-distance-and-delay-based, and cordon-based schemes under low-cost, medium-cost, and high-cost regimes. The study’s contributions also include an extension of the modeling framework to include public transport options to allow for providing travelers with the option of choosing an alternative mode of transport if they do not wish to pay. A mesoscopic stochastic route choice modeling approach was adopted to examine the impact of road pricing inside a nominated charging zone within the network. The results showed it would be possible to achieve a reduction of 11% in vehicle count, a 20% reduction in travel time, a 13% reduction in emissions, and a 3% increase in travel speed within the proposed pricing zone under a high-cost pricing scenario. The results also showed a significant reduction in emissions resulting from shifting drivers who are not willing to pay the congestion charge to public transport. When 20% of car drivers shifted to public transport, carbon emissions were reduced by up to 30% and network performance improved by 45%, compared to the baseline scenario without pricing. The findings of this research provide important directions for policymakers in deciding on the type and scope of charging schemes to use and how these could reshape transportation taxation systems by moving away from taxes on vehicles through registration fees and towards user-pay taxations where travelers pay for the amount of travel they do or the pollution and emissions they are responsible for.

A new numerical mesoscopic scale one-domain approach solver for free fluid/porous medium interaction

A new numerical continuum one-domain approach (ODA) solver is presented for the simulation of the transfer processes between a free fluid and a porous medium. The solver is developed in the mesoscopic scale framework, where a continuous variation of the physical parameters of the porous medium (e.g., porosity and permeability) is assumed. The Navier–Stokes–Brinkman equations are solved along with the continuity equation, under the hypothesis of incompressible fluid. The porous medium is assumed to be fully saturated and can potentially be anisotropic. The domain is discretized with unstructured meshes allowing local refinements. A fractional time step procedure is applied, where one predictor and two corrector steps are solved within each time iteration. The predictor step is solved in the framework of a marching in space and time procedure, with some important numerical advantages. The two corrector steps require the solution of large linear systems, whose matrices are sparse, symmetric and positive definite, with M-matrix property over Delaunay-meshes. A fast and efficient solution is obtained using a preconditioned conjugate gradient method. The discretization adopted for the two corrector steps can be regarded as a Two-Point-Flux-Approximation (TPFA) scheme, which, unlike the standard TPFA schemes, does not require the grid mesh to be K-orthogonal, (with K the anisotropy tensor). As demonstrated with the provided test cases, the proposed scheme correctly retains the anisotropy effects within the porous medium. Furthermore, it overcomes the restrictions of existing mesoscopic scale one-domain approaches proposed in the literature.