Simulating the Hydrodynamic Nature of Porosity

My name is August young, I'm a PhD candidate at Duke University, and I am simulating the dynamic nature of porosity using ND Solve. So we are going to be discussing the hydrodynamic nature of porosity. And this is really quite foundational. Normally porosity is just assumed to be a constant and any porous media, and there is a term called effective porosity, and that just means what part of that empty void space in your porous media is conducive to flow? Well, we're going to show that that effective porosity is flow dependent. Meaning as you increase your flow speed, you're changing your effective porosity. So why do we care about this redefinition of porosity? Why should we focus on it? Well, if we think about groundwater flows and contaminated aquifers, those remediation projects are extremely costly and timely. They can run on the order of magnitude of a century and they can cost anywhere from tens of millions to a few multiples of that dollars. So we're looking at trying to be more effective with our use of remediation techniques. So if we now scale down and look at our pore space, since we're concerned with porosity, we can divide our pore space into two different zones. They're called the mobile and immobile zones. So if we look at our magnification A, we see some nice through flow. We have well connected channels, no large obstructions. And if we look at magnification B, we see that there's a vortex sitting in a little cavity, and that little cavity does not participate in the through flow. So that volume of that pore space could just be blocked off and it would be doing the same thing it's doing here in terms of flow. So when we think about the effective porosity of a medium, we're thinking about all of the mobile zones, what is conducive to flow. And like I just discussed, that recirculation that you saw in that little cavity is the dead end pore. In that dead end pore, contaminants hang out when you just run a steady flow field within your groundwater, your medium, like an aquifer, you are just relying on diffusion to remove those contaminants from that cavity. And that's a very slow mechanism. Because of that, oftentimes most of these remediation schedules stop before all your contaminants are removed. And you have something called contaminant rebound as diffusion continues to move those contaminants back into your through channels. So those dead end pores are really what we need to target if we want to be thorough about removing contaminants and tech. Typically, the dead end pour is idealized by this very idealistic geometry where we have a square through channel or a rectangular screw through channel and a square dead end pore hanging off of that through channel. And they're the same dimensionally, so we've non dimensionalized the whole problem not working in units anymore. So we can communicate across different scales. And originally when geologists and fluid dynamicists looked at these pores, they said we're going to make a really crude approximation and we're going to say that there is absolutely no flow within that dead end pore, meaning it's stagnant. But we know that that is not true because there's a boundary between our dead end pore and our through flow. And that would violate the boundary conditions that you have flow along the top side and nothing moving along the bottom side. So that oversimplification has to be changed. And the most common way that this is resolved is that there's a geometric boundary condition between our through channel and our cavity. And that boundary is referred to as a separatrix. So the and within that dead end pore space we just see a cavity taking up that entire volume. But we know that this approximation is false. So the improvement on this is that that separation between the through channel and the dead end pore is defined by a fluidic boundary condition. So there's actually impingement of that through channel flow into the cavity. And that boundary between the through channel and and the dead end pore is known as the separatrix. And this is the model that we will proceed with. So how do we model the dead end poor? Well, we use the super function, as I like to call it ndsolve. And before we can do that, we do need to do a few things. But we'll get to the beauty of why we use NDSOLVE here. So we define our solver mesh, we just discretize our domain and using our two element mesh function, and we look at our boundary nodes so that we can satisfy our boundary conditions as well as our internal mesh. And then we also use a refinement region within the area that we expect that fluidic separation to be that separatrix so that like we can really see what's going on there. Since our definition of that location is going to be dependent on how well we refine our flow field in that area. So the governing equations that we'll use to define this flow field are the Navier Stokes equation and the equations and the continuity equation. We make a few assumptions. We say that our flow is incompressible, steady state, two dimensional. And the result of this we have these three equations right here that we are able to code into Mathematica with almost no, no lines of, of coding. If you were to try to put Python or Matlab, it would be hundreds of lines. And the boundary conditions that we enforce is just that. Within that idealized pore space, we don't have any flow into the pore or into the walls or along the walls. We input a parabolic flow profile and we just set a arbitrary pressure condition at the outlet. So the flow solver that we were in lives in an external package. And we're going to demonstrate what that flow solver outputs, which is a stream plot. We tell it the Reynolds number that we would like. So Reynolds number is just a non dimensional way to describe velocity. They're directly proportional and we just modify our velocity with some fluid parameters. So when I mentioned Reynolds number, think velocity. So what we get when we tell it we want a Reynolds Number of 1 within a lower and upper bound. And this last parameter just describes how many grid points we would like in our domain so that the stream plot will include those grid points. We get this nice stream plot and it very clearly shows where there is recirculatory flow in our domain and where there's through flow. So how do we get this result in four lines of code, which could of course be one if I just go back into notebook form and condense it, which I have expanded here, so it will all be shown on one screen. So what we do is we just say our outputs are going to be equivalent to what ndsolve gives us. We input our governing equations, our boundary conditions, our three variables, and the region that we'd like it to solve over and we get that beautiful stream plot. So now we need to describe where that separatrix is. And in order to do that, we have another module that will take our stream plot and then it then overlays an interpolating function and it allows the user to move these points and determine where this interpolating function should live. Now obviously this is not correct, so we'd say, okay, it looks like we need to be about in this space. Now, of course, the way that this program works, as I previously mentioned, the location of the separatrix is dependent on how you generate the stream plot, what your working precision is, et cetera, but more or less what we're Trying to show is the movement of this. And it's not going to be the exact separatrix that we're characterizing, but it's close enough. So then we save the points that we've selected and we output it to our data matrix and then we input those points into an interpolating function and we get a nice mathematical description of what our separatrix looks like and where is located. Now that we have that description of the separatrix, we can actually show where it is in that dead end pore space and we can quantify things like the maximum penetration of that separatrix into the dead end pore. And so that's something that was important to do because previously this, this work had been done for a very small range of Reynolds numbers. And so we needed to verify our numerical results against previous results that were generated using a commercial flow solver. And with our interpolating function we can also do things like now integrate to find our mobile and immobile areas, we can define our effective porosity within the pore space. So here are our previous results in blue. I apologize, that's mislabeled. Our previous results are in yellow and our new results are in blue. And you can see that there's decent agreements between our new results and our old. We have about a 2% difference at most between our data points. And we find that that separatrix that I was talking about really doesn't move too much below a Reynolds number of one and but above a Reynolds number of one, which is what you will see when you start to induce flow field in a vertical circulation. When you're cleaning up in those remediation projects, we see a substantial decrease in that sub matrix depth actually starts to move up within our dead end pore space. So now that we're starting to think about ERs in that previous model, we see that in this regime of Reynolds numbers, our separatrix location of it actually allows for flow and about 15 to 20% of that dead end pore area that that geometric model says isn't there. So while this might seem like not a lot, we do have to think about the fact that we are looking at a single pore and then scale it up to the tens and hundreds of meters that we would be looking at within these remediation systems. So if we for instance, instead of just looking at that single Reynolds number that we looked at, we looked at Reynolds numbers ranging that whole range of laminar flow, which is when we have that nice smooth flow without any turbulence, we See substantial movement of that separatrix that you can illustrate with a stream plot. So here we have a Reynolds number of one. We increase our order of magnitude by one, we increase, we multiply that by five, and you see that the separatrix is almost now at that geometric boundary of the pore space. And by the time we get to a Reynolds Number of 100, which is about the cap of where we would be for laminar flow and porous media, we are basically at that geometric definition. So if we then just look at where that separatrix is for those stream plots, you can see here, there's a nice straight blue line that, that defines what that geometric boundary is. We see that by the time we get to the end of our laminar flow regime, those two models are in agreeance. But when you're closer to those low Reynolds numbers which you see in natural groundwater flows, we're not in agreeance with that previous model at all. And if we zoom out, this is what it looks like in our entire pore space. So if we wanted to then describe the porosity of this geometry and look at what is conducive to through flow, AKA define our effective porosity, we would then look at how we could fit a function, a relationship between our effective porosity and our Reynolds numbers. And when we go across our entire laminar flow regime, so that's a number of around 0 to 100, we find that there's a really great cubicle fit that describes the relationship between effective porosity and velocity. And if we wanted to look back at the fluidic and geometric boundary conditions, that for instance, the Reynolds Number of 1, we see a 20% error in the effective porosity of the dead end pore geometry. And it's not until we get to a Reynolds number of around 50 that the error in this definition is less than 1. So that that Reynolds number is really uncommon even in your induced flows until you get right to the well where you're inducing your flow. So what we see here is a substantial error that propagates out. I'm 5% might seem low, but again, think about the spatial scale that we're working on. So what we've learned here is that the flow models that describe porosity traditionally underestimate the porosity of the dead end pore space when we use that geometric boundary condition, especially when we're at those lower Reynolds numbers, where we have impingement, deep impingement of that separatrix into the dead end pore space. And we also found that there's a nice cubic relationship that describes porosity as a function of velocity. And porosity, again I would like to reemphasize, is typically assumed to just be constant. And we're able to model this in a fraction of the code that would have been required in a low level programming language like matlab, Python, C languages that are far more common in research for fluid dynamics. So our conclusion is that if we want to better estimate the cost and timeframe of a remediation project, the porosity of the medium should be defined as flow dependent, since that is required if we're going to accurately model our velocity fields. And we can quite easily do that in Map Mathematica. That about wraps it up. So if anybody would like to ask any questions, I can field them now. So we have seeded this geometry with points and put this in the time domain. But this is an easy way to solve the problem without being in the time domain. Yes, yeah, yes. So that is also something. So this is kind of the foundation and that is an application that we have actually used with this. So there's two mechanisms that you will induce. If you start, if you put an oscillatory flow into that geometry, you will actually get that vortex that recirculates to completely eject. And then you will also get, when you. So if you ramp up your flow at the top of your sine or cosine, whatever it is you want to put in, you'll get a deep sweep where you have 100% of your through flow scrapes along that boundary. So you have maybe some really, really small recirculating vortices in the corner, but they're, they're very negligible compared to the overall volume of the dead end force base. So that will remove somewhere around like 70% of the contaminants and of, you know, of that period that we're looking at. And then the vortex ejection will remove not all of the contaminants that are in your pore, but the substantial amount, the remaining. So those two mechanisms will change your effective porosity of that dot and pore space from 0 to 1, which is quite more than the 20% that we were looking at here. Yeah, yes. Yeah. So you would be looking at your pulse. How can you optimize your pulse for certain matrices and shapes and et cetera. Yeah, yeah, thank you. Right, and so that movement of those contaminants is dependent on where that diffusion occurs. Like across that separatrix is what's moving, diffusion is what's Moving those contaminants. And so here we're characterizing exactly where that separatrix is. So some of them are. So they all behave differently and they do different things. So some of them will be adsorbed into whatever material it is that we're trying to remediate. And some of them won't be. And they'll, they will just be hanging out in those cavities and they might be aggregating or they might be equally diffused before we start to induce a velocity within there. But let me pull up an image of a vertical circulation well, and we can see how we would move contaminants within the groundwater matrix. So here we have a vertical circulation well, what we're doing is we are pumping. Think about dipole flow. So you have a source and a sink, and they're ones pumping out and one's pulling in. And up here we have our source and then we have our sink. So as we're pushing, our water comes out of this source and gets pulled into our sink. So our contaminants are effectively moving along this dipole flow pattern. And we're capturing them right here at the well bore. No, no, no. We're starting with contaminated water. Yeah. And we're trying to remove our contaminants. Yeah, yeah, exactly. So that's something that would have to be done application by application. So looking at certain mediums, and of course this would also be an approximation, but if you got a representative volume of your medium using image analysis, try to determine about what your pores look like. Do you, if you're working with rocks, you have fractures mostly. But if you're looking at like granular media, are your media rounded, Are they jagged? So washed or unwashed, those are big differences. We can start to say, well, it looks like we have a lot of small crevices that are small cavities that are on the same order of magnitude as our through channels. That's an important distinction you can start to make. And then you can say, is it washed? Are they rounded? Is it unwashed? Are they jagged? Do they look like triangles? So you'd have to do some sort of analysis on your media to really make a more accurate approximation of what's happening in the pore space. But this has actually been so as deep sweeps and vortex ejections that we were talking about that oscillatory flow actually has been applied to unwashed media. So the jagged cavities. And it did clean up the same amount of as the case of steady flow. And it was about 12% quicker for that. But the issue with that model was that there are preferential flow paths along the boundary. Wasn't actually like a real vertical circulation while out in practice, I'm not sure if that has been applied at all. I've been asked a lot instead of about extraction, about actually seeding porous media with certain bacteria or food for the bacteria to clean up these cavities. But I have not heard much about extraction of oil. So unless somebody has written your code for you for your specific application and that other language, you're going to need to do a lot of modification. So those external packages that I used, I wrote, but I was using functions that were very well documented, so we can pull that up. This might look like a lot of code. It's really not, I promise, is like one Navier Stokes equation. Like you are working with a substantial chunk of code. So this is our flow solver that I talked about right here. All I'm doing is setting up the region and I'm applying a mesh to that region and I'm defining that refinement region that I wanted to use to be able to really zoom in on that separatrix. Down here I'm just setting up some equations that are going to describe my boundary conditions. And here is the governing equations which I showed you before here. Mathematica is great. It already has my type of boundary condition that I want to apply. So all I needed to do is follow the format I needed to apply those. And here's, like I said, the one line of code that I expanded into four, just presentation format that actually solves my whole system. Here I generate my grid points that I want my stream plot to include. And this is me generating the stream plot. So that's the entire flow solver. Here's the code that I use to describe the separatrix. So that's that interpolating function. And then finally we have the module where I let the user describe the location of the separatrix by placing those points. Yeah, it's super condensed and I've done this for different geometry dead end pores. And all you're doing is just instead of saying I want to test this one Reynolds number, I just put a table in and run through a range of Reynolds numbers. So when I go through, and I've actually this, this code has been simplified just for that single geometry that we looked at. When I go through and I apply this to a range of geometries, this does not get any longer. Like the pore spaces. Yeah, yes, exactly. So that first are approximately this geometry everywhere. And then you could go in and you could sample different. If you expect that the formation that you're looking at has, for instance, different layers. So maybe we have rock and then soil and etc. Because that's much more realistic when we're looking at actual groundwater flows. Yes. Well said. Okay, thank you for your time, everybody.