Modeling Fluid Dynamics

Hello, everybody. My name is Oliver and I work in the numerics group. And today I'll be talking about fluid dynamics. So let's get started. So this talk is in a series of other PDE modeling talks. So this initial slide has some links to the documentation of various fields of physics that we've written monographs for that provide some functionality to model in those domains with partial differential equations. So we've done acoustics, mass transport, heat transfer. In version 13, we've added solid mechanics. We've continued to extend that until this version, version 14. There was just a talk about that by Alessandro, who explained the new additions to hyper elasticity and solid mechanics. Earlier this morning, Sebastian gave a talk about electromagnetics. I'm going to talk about fluid dynamics. And a little later today, Ricardo is talking about the Schrodinger equation. So maybe with the exception of the last one, this is the list of fields of physics that we currently have PDE models for. That doesn't mean that NDSolve can't solve other equations. It's just that for these fields of physics, we have a framework to generate the equations and monographs to assist you in making these models. Okay. The monograph for the fluid flow PD component is called Laminar flow. This is a link. You can later download the presentation and follow that link to the documentation. This is going to be available in version 14. Right, so. So the new function that we've introduced is fluid flow PD component, and we're going to see examples of that in just a second. What can this function do? Well, it solves a compressible Navier Stokes equation, it can solve the incompressible Navier Stokes equation, you can solve non Newtonian flow, and you can couple this PD operator to other fields of physics to create multiphysics examples. Right, So a super quick recap of the Navier Stokes equation. We start with a compressible version of this. The dependent variables of the Navier Stokes equation are this bold V, which is a vector for the velocity U, V and W in 3D, and P is the pressure that are the U. So the velocity and the pressure are the primary dependent variables, the yield variables and. And we have a tau, which is a viscous stress tensor. And the second equation, the first equation is a momentum equation, and the second equation is a continuity equation. The viscous stress tensor is mainly dependent on the dynamic viscosity and the strain rate. And the strain rate is the change of velocity, the gradient of the velocity field. So that's for the compressible version. Now, the Incompressible version looks very similar. The thing is, now that the density is constant, rho is constant. And that means that this term here simplifies a bit. So this term disappears. And also here the constant is eliminated. Now. Now what's going to become important later in the talk is the following. So if you look at this, P is only present in this divergence. So only the derivative of P is here. And in this equation or in this, in the incompressible version, P is not present at all. So we're going to encounter some problems. I'm going to show you how we solve these type of things. Okay, let's start by looking at a first example. So this is a famous benchmark type of problem. It's called the lift driven cavity example. And you have a rectangular region. On top of this region we have a fluid flow field in the positive X direction. The V velocity is zero. All other walls, U and V, the velocities are set to zero. And we have a small pressure point, the reference pressure point, which we set to zero down here. Now, the way to use this function, the fluid flow PD component, is as follows. You set up variables U, V and P, the primary dependent variables, and they are dependent of X and y. It's a 2D problem. And the parameter for this case is just a Reynolds number. A Reynolds number is a ratio. It's an inertial force versus a viscous force. This is a very simple setup. You just have to specify this Reynolds number. We'll see other ways to set this up in a second, but for now we just use this Reynolds number. Fluid flow PD component then gets the vars and pars and will set the right hand side to zero. And it's going to generate an equation for us that is the Navier Stokes equation in this case. Okay, we want some boundary conditions. I've mentioned them earlier. So we have a Dirichlet condition on top. This is why y equals 1. We have a top condition where the u velocity is 1. So a fluid flow with 1 meter per second in that direction and 0 in the y direction and all other areas are set to zero. So the wall is set to zero. And we have a pressure point at this point here. Now, the first thing that's maybe a bit odd is that you see Dirichlet conditions here. And the reason is that we don't yet have fluid flow specific boundary conditions or named fluid flow conditions. And the reason for that is that we plan to extend the fluid flow capacities that we have a bit further and before committing to some specific boundary conditions, we first like to understand better better how these extensions should work so that the boundary conditions that we come up with then fit nicely in all scenarios that we'd like to tackle. So for now we're using dirty condition. But that's not a disadvantage Then solving the problem is you plug it into ndsolve the pde, the boundary conditions, dependent variables, you specify the region and you have to specify a method option. We'll get. Plot this then. So here is the solution of the U and V velocity. We see some reflow areas, recirculation areas here and here. So stream plot is one way to visualize the solutions. Then now I'd like to get back to this method option. So here I've evaluated the same example I've given, I've commented this out and we get a load of errors. And the reason is that as I mentioned earlier, P is only present in its derivative form. So only the derivative of the pressure is here and here. We don't have a P at all in this equation. The way to stabilize this in the finite element method is to say that the fluid velocities are interpolated by second order interpolations, U and V and P. The pressure is only interpolated in the first order. That's the way to deal with this. And in finite element speak. This is a Taylor Hood type of element. The first version you see here is for the stationary case. And this is how you'd specify this for the transient case. We may in the future. This is not in 14.0 yet. We may in the future maybe have something like this. The fluid flow options that will automatically generate you the correct options for this setup. Or we may go ahead and like wrap everything in a PDE model that you'll pass to ND solve. But again, this is still future. For now you'll have to give these method options to stabilize the solution of the Navier Stokes equation. Okay, parameter specification. We've seen the Reynolds number. You could equally well say I have a specific material like glycerol here. This is an entity or. Or you can specify the dynamic viscosity and the mass density manually. All of these ways work for boundary conditions. Typically we have inflow conditions, so it's a Dirichlet condition where you specify the U and V at the inflow. An outflow condition is often done with specifying the pressure or setting it to zero. Wall boundary conditions, like a no slip condition is when you have U and V at the wall set to zero. And the slip boundary condition is that you have one, two components set to zero and the other one is free to move. There is another type of boundary condition that's called a traction boundary condition. And basically what that means is this component that's in the divergence is the Neumann value that's been put in here and here I've extended it completely. This is a somewhat complicated boundary condition. It's useful nevertheless, and I have an extensive example of how to set these traction boundary conditions up. This is just to let you know that we can also handle these type of conditions. Right. Let's look at a time dependent problem. Now the difference now is that U, V and P get this time variable t and here's another T. And then you'll set up the fluid flow component essentially in the same way. Now in this example we have a long pipe, maybe that has an obstacle. This hole is slightly off center. We'll have an inflow condition on this left hand side. Flow profile is going to come in. We have a wall here and up here. And at the end we have an outflow condition. And we're going to see that in a second again. So setup is U V and P time dependent. I specify a mass density of 1 dynamic viscosity of 10 minus 3. And the operator is then just calling fluid flow PD component, it generates the operator. And next we set up the boundary conditions. So here the inflow condition specifies the U velocity. So the velocity in the x direction and the V velocity is zero. Now this expression here is interesting. So this part here is a profile, an inflow profile. It's a parabolic profile. But we have a function here that we're going to talk about in a second that's going to ramp up this inflow profile. Then we have for x smaller than or larger than 0 and smaller than length, we have the wall condition. And at the outflow we have P equals to zero and the initial conditions are all set to zero. Now because the initial conditions are all zero, this entire system is at rest. And what we'd like to do is we slowly like to ramp up this inflow condition. So I've written a small ramp function that for up until what, 2, 2.5 is 0 and then slowly goes up to 1 and at about 6 or 5 it has reached the value of 1 and that is going to scale the inflow profile. So when we time integrate this, we have the operator set to zero. The boundary conditions, initial conditions, UVMP, the other dependent variables, we simulate from T0 to 10 and this highlighted thing is just the method option that we talked about earlier. So 120 seconds is what we need to solve this. And what you see now is the velocity in the X direction. So what we see is a nice Von Karma on the vortex street. So you'll just see how this behaves in this fluid flow problem. Okay, Right. So the next thing I'd like to show is non Newtonian fluid flow. So in non Newtonian fluid flow we have a nonlinear stress strain relation. And what that means is that this we have an apparent viscosity. So this dynamic viscosity now is an apparent viscosity and that is a function of a shear rate. The shear rate is computed from the strain rate. So this is basically the implementation that's in the code. So take the strain rate, tens are contracted, multiply by two, take the square root literally what's written here. And that's the function that computes the shear rate. So let's look at an example of that. Right before we do that, there are several non Newtonian fluid flow models that we've implemented. Power law, Coral Bingham and more. All of these models are available or visible to you in the documentation. You can see how we've implemented them. And the reason for that is if you want to make your own custom model, you can do so, and these can serve as guides to do that, to extend them. So implementation code is available and you can also add your own non Newtonian fluid flow model if you like to do so. Right. So as an example, I have a power law flow. So here the apparent viscosity is some formula. The details are not really relevant now, but the way to specify this is then a fluid dynamics material model, the power law. You specify the exponent, minimal shear rate, all these different parameters that these material models have. And that's the way how you specify them. In this example, I have a pipe that's opening is becoming a bit wider. On the left hand side we'll have an inflow and on the right hand side an outflow. And the rest is all a wall condition. What I'm going to do now is we're going to solve two, two equations. So one is a normal Newtonian flow and the other one is a non Newtonian flow. You see here the parameters of the power law model are the same as the normal fluid flow parameters. I've just added the fluid dynamics material model with some specifics of these for these parameters. And then we create two fluid flow PDEs to solve. And this is now we plug them into ND solve and this is now the solution. So here we have the flow profile at a certain position in the Newtonian case. And this is the same flow boundaries. And now with a non Newtonian flow, you see the difference between these two types of flows. So that's something we can also do. Another important thing is energy transport. So if you have a fluid flow, you'd like to know how much energy it transports. So meaning coupling fluid flow with heat transfer, you can couple these equations with mass transport and other things. It's not a problem. I'm just focusing on heat transfer here. The way to do this is by the Buzines approximation. And we're going to have a small look at that, what that is. So the Buzines approximation assumes that the density is constant and except in this buoyancy term. So this is the term where the density is not constant. I've marked this here in red. And all other densities are reference densities which are constant. This term is modeled by a reference density. Again, 1 minus alpha. Alpha is a thermal coefficient T minus 2 ref. That's basically how this is done. We're going to see this in more detail in a second. So what we do now is this is a flow. What we're going to do is a flow driven by heat difference. So we have a rectangular region where we have heat on the bottom, cold at the top, and we'll see how that temperature difference is going to drive a fluid flow in this region. So this is the mesh for this. Now we have two variables, the fluid variables U, V and P and the heat variables like the temperature time dependent. Again, the next step is to set up the parameters for this model. So I've used symbolic parameters for now. I'll show you later why I've done that. So the reference density, the dynamic viscosity, specific heat capacity, thermal conductivity. This is the heat equation here that we're solving. And here is the key where the magic happens. So the heat convection velocity is this term here. This velocity term is exactly the same velocity from the Navier Stokes equation. So it's coupled back into the heat transfer equation. It's a fully coupled model. And this is now how you set up the equations you call fluid flow. PD component, heat transfer component. And this adds this Businesk approximation that I've defined down here again here, T smaller one means that for time smaller than one, we use this as a scaling factor again to slowly ramp this up like in the previous example. Boundary conditions, everything for the fluid flow, everything is a Wall. So U and V are set to zero and everything true means everything on the boundary. And we set a small pressure point, a reference pressure point here in this lower corner. The temperature are a heat temperature condition, hot at the bottom and cold at the top. And we collect these all in a variable. Okay, then. Now I've specified a list of parameters. What cold, hot, the reference temperature, the density and stuff. What they all are the initial conditions. And then I create the model. So t end is 500 seconds, the model. And then I replace all symbolic parameters at this point with these actual material parameters. We'll plug this into ndesolve. Now we have U, V, P and T. So the temperature is now an additional dependent variable. This takes a little more than two minutes on my laptop. And here is the solution of this. So now you see how the heat is slowly stirring up the fluid. And we'll see nice rotational circles here of how the fluid is trying to convect inside this container, if you will. Okay. Right. So one last thing I'd like to mention is compressible flow. So the difference now is that the density is no longer constant. What that means is in this case we have a constant. So this symbol is constant. And the equation, it's easy seen in the mass transport, in the continuity equation here, rho is constant. But if we now specify the mass density to be a function of space and space or could also be time, then the equation down here is going to change. So the fluid flow PD component automatically understands what whether you need a compressible or an incompressible Navier Stokes equation is going to generate you the correct thing to do that you can use for your equations. Okay, that's what I wanted to show. Basically we provide functionality to do 2D, 3D compressible, incompressible, Newtonian, non Newtonian fluid flow. All of this has its standard reference page, of course, and there's an entire monograph focusing on a laminar flow. If you have questions, you can always reach me@wolfram rubenco.wolfram.com My email is not a secret, so I enjoy contacts with customers. One question I have for you is where would you like to see this go in the future? So I mean, 14.0 is going to be new for you. I'm already thinking about 14.1 and maybe you have some ideas of what you'd like to see in 14.1. Maybe turbulence modeling, like Reynolds averaged turbulence model might be interesting. Let me know what you think. So I'm going to end the presentation here. Now and I'm going to see if I have questions. Stop sharing screen. Yeah. Okay, so let me see. Great development. Thank you. Thank you, I appreciate that. Okay, so I'm going to stick around a little longer here. If someone has questions, I mean, if you want to chat with me about PDE modeling, numerics, anything like that, ping me. It's always a pleasure talking with customers and we'll see if we can. If I can help you if you have a specific problem. But it's fine if you also want to just chat with me. No problem. Okay, so thank you. Thank you for the nice comment. Okay, first question. Would you recommend to use the free fluid flow component for investing Corvette Taylor Stream with heat differences? Are they designed to handle parameter parameter verification to look at the emergence of coherent structures? I don't know. I've never tried this, so I'm a little bit hesitant to say yes because I simply don't know what could potentially come up there. Maybe give it a try and if you run into issues, contact me and I'll try to figure out how to solve these things. So I can't quite say whether it's possible or not as I haven't tried this myself. Okay, thank you for the presentation. Thank you. Can you outline how these Anton asks, can you outline how these functions can be applied to simulating the traveling of oil slicks in ocean sea? Oh, Anton, Probably I can't because I'm lacking knowledge about how oil slick in sea ocean is modeled. I don't know too much about that, so it's probably something I'll have to dig into and see what typically typical modeling approaches are. Okay. Another question from Eric Brown is parameter fitting against experimental data. Well, fitting is principally of course possible in Mathematica you'd use find fit or something like that to use to create a fit. Either a symbolic fit and that you could use as a parameter in. I don't know, the viscosity. I assume that should definitely be possible. Or you could use an interpolating function for the viscosity, but a symbolic fit for a viscosity would be beneficial because then the internal compiler can generate efficient code. It's also going to work with an interpolating function, but it probably would take a little bit longer to solve then. But I think this is doable. If I understand this question correctly. I think it's doable. Okay, so do you think we can use neural networks instead with some of the basis functions? Okay, yeah. This is a famous question and Stephen asked similar questions. Okay, so neural Nets, I think they can be useful in engineering, but not. At least not yet, maybe as a solver for partial differential equations. Neural nets, I think, are useful engineering. When you don't have an exact physical understanding of how something works, you can use them to train them and provide a heuristic model. In some scenarios, I think that's where they're useful, and I think that's where they're used, for example, in turbulence modeling, because we don't understand how turbulence works, really. I mean, to the deepest level there, neural nets are used to generate heuristic models because we simply don't understand the physics. My stance on this is I personally try or like to understand how things work by understanding how the physics behind things work. And neural nets is something I see, at least in the engineering domain, as a last resort type of thing. If we don't understand how something works, we can use a neural net to help us model this. But the goal, in my opinion, should always be to understand the physics better than try to use a neural net. At least that's my stance on this. But I'm pretty sure this is subject to debate, and many people would disagree with me. And that's fine, of course, but that's just my stance on neural nets in PDE modeling. Okay, Anton, I don't quite understand this last question you have. Okay, but you agree with me that the physics understanding is good. That. So that that's a good thing then? I guess. Okay, so if there are no other questions. Okay, so thanks a lot for attending. Enjoy the rest of your conference. And if you want to talk about PDEs, numerics, anything like that, ping me. I'm always keen on interacting with people. Have a good day. Bye now.