The details of the numerical methods and validation can be found in our previous study. In short, NS and PNP were solved in a block iterative manner, and iterated until self-consistency within each timestep was achieved (flowchart, Figure ). For the NS equation, the non-linear convection term was linearized. A 2 nd order accurate backward differential formula (BDF 2) was used for discretizing the time dependent terms. Following standard strategy, a variational multiscale method (VMS) was used for the NS equation. For the PNP equation, a streamline upwind Petrov-Galerkin (SUPG) style stabilizer was adopted for the convection and electric migration terms. In addition, we utilize an octree based adaptive mesh generator that enables massive parallelization and complex geometry. No-slip boundary conditions for fluid flow were applied for all solid surfaces in the domain.