Numerical solution to stationary viscous flow with a Darcy term

[Irid]

Hi,
I'm trying to solve the flow profile inside an inhomogeneous porous material between two parallel moving plates (essentially Couette flow with a deviation), and I model my system by the following equations:
$\nabla^2 \mathbf{u} = p(x,y,z) \mathbf{u}\\ \nabla \cdot \mathbf{u} = 0$

where p(x,y,z) is a given function describing the resistance of the porous material. The boundary conditions on u are of Dirichlet type. I am looking for a numerical solution to this problem. I have tried to solve the viscous equations by a relaxation technique (Jacobi method, etc.) and then impose the incompressibility by some other means such as:
1) Pressure projection method. The problem I encounter is that I have a hard time solving the pressure Poisson equation with Neumann boundary conditions (grad of pressure is the velocity at the boundary).
2) Penalty function. To my first equation I add a penalty term proportional to the local divergence. This effectivelly couples the flow in various directions, but the total of the divergence keeps increasing and my solution never settles to some stationary state.

My main problem is to enforce the incompressibility (zero divergence). Any ideas? I think this is not a difficult problem, I just can't seem to find any examples in the literature...

 

[Chestermiller]

Hi,
I'm trying to solve the flow profile inside an inhomogeneous porous material between two parallel moving plates (essentially Couette flow with a deviation), and I model my system by the following equations:
$\nabla^2 \mathbf{u} = p(x,y,z) \mathbf{u}\\ \nabla \cdot \mathbf{u} = 0$

where p(x,y,z) is a given function describing the resistance of the porous material. The boundary conditions on u are of Dirichlet type. I am looking for a numerical solution to this problem. I have tried to solve the viscous equations by a relaxation technique (Jacobi method, etc.) and then impose the incompressibility by some other means such as:
1) Pressure projection method. The problem I encounter is that I have a hard time solving the pressure Poisson equation with Neumann boundary conditions (grad of pressure is the velocity at the boundary).
2) Penalty function. To my first equation I add a penalty term proportional to the local divergence. This effectivelly couples the flow in various directions, but the total of the divergence keeps increasing and my solution never settles to some stationary state.

My main problem is to enforce the incompressibility (zero divergence). Any ideas? I think this is not a difficult problem, I just can't seem to find any examples in the literature...

I've had considerable experience modeling problems for flow in porous media, but I'm not able to understand your description. Please describe the physical problem, including boundary conditions. What are the equations and boundary conditions in component form? I don't see any time term, so I suppose this is steady flow.

 

[Irid]

I have two parallel infinite plates, one of which is moving at constant speed. This gives rise to Couette flow. Now at the bottom (stationary) plate I adsorb a porous particle (we may consider high porosity case only, i.e. the particle is soaked with the liquid). The porosity is a given function p(x,y,z) and goes to one rapidly away from the particle. So essentially we have a Couette flow with a perturbation near the origin. My initial attemps have flow profiles such as shown in this graph:

https://docs.google.com/file/d/0B6xidb7BFzkxMzhBRU1hTDNNRm8/edit?usp=sharing

The problem is that the incompressibility constraint is not satisfied, and the divergence keeps growing with every iteration. Maybe I have to use staggered grids to enforce incompressibility?..

 

[Irid]

I've had considerable experience modeling problems for flow in porous media, but I'm not able to understand your description. Please describe the physical problem, including boundary conditions. What are the equations and boundary conditions in component form? I don't see any time term, so I suppose this is steady flow.

Yes, I am in a steady-state condition, so the flow is laminar, but with a distortion due to inhomogeneous porosity, localized near the origin. The boundary conditions are exactly the same as for Couette flow, so
u(x,y,z=0) = 0 (stationary plate)
u(x,y,z=H) = u0 (moving plate)

v,w = 0 at the plates.
on the sides, x,y --> +/- infty, I should recover Couette flow:
u(z) = u0 z/H
while v,w are zero again.

 

[Chestermiller]

I have two parallel infinite plates, one of which is moving at constant speed. This gives rise to Couette flow. Now at the bottom (stationary) plate I adsorb a porous particle (we may consider high porosity case only, i.e. the particle is soaked with the liquid). The porosity is a given function p(x,y,z) and goes to one rapidly away from the particle. So essentially we have a Couette flow with a perturbation near the origin. My initial attemps have flow profiles such as shown in this graph:

https://docs.google.com/file/d/0B6xidb7BFzkxMzhBRU1hTDNNRm8/edit?usp=sharing

The problem is that the incompressibility constraint is not satisfied, and the divergence keeps growing with every iteration. Maybe I have to use staggered grids to enforce incompressibility?..

Is this a 2D flow, or is it 3D. What is the shape of the particle? From the figure, it looks like half cylinder into the paper. I assume inertial is not important. Have you solved the problem for a non-porous particle? What is the boundary condition you use at the surface of the porous particle?

 

[Irid]

Is this a 2D flow, or is it 3D. What is the shape of the particle? From the figure, it looks like half cylinder into the paper. I assume inertial is not important. Have you solved the problem for a non-porous particle? What is the boundary condition you use at the surface of the porous particle?

To start off, I would be happy to solve a 2D system, but eventually I hope to be able to tackle the full 3D problem. The particle is irregular-shaped, here is how the density profile looks like:

https://docs.google.com/file/d/0B6xidb7BFzkxb0p1UVM2azNQNmc/edit?usp=sharing

For now I just assume that the porosity is linearly related to the above density. There is no sharp boundary between the particle and the free liquid. You may assume that the whole space is inhomogenous porous medium, but keep in mind that the porosity goes to 1 (no particle) far away from the origin.

 

[Chestermiller]

To start off, I would be happy to solve a 2D system, but eventually I hope to be able to tackle the full 3D problem. The particle is irregular-shaped, here is how the density profile looks like:

https://docs.google.com/file/d/0B6xidb7BFzkxb0p1UVM2azNQNmc/edit?usp=sharing

For now I just assume that the porosity is linearly related to the above density. There is no sharp boundary between the particle and the free liquid. You may assume that the whole space is inhomogenous porous medium, but keep in mind that the porosity goes to 1 (no particle) far away from the origin.

For this porous particle, how does the permeability vary. Also, I'm sure you know that the differential equations for flow inside a porous medium are different mathematically from the differential equation for viscous flow in a non-porous medium. How do you make the transition?

 

[Irid]

For this porous particle, how does the permeability vary. Also, I'm sure you know that the differential equations for flow inside a porous medium are different mathematically from the differential equation for viscous flow in a non-porous medium. How do you make the transition?

I assume that the permeability is inversely proportional to the density as shown in the figure. Eventually I might consider some other constitutive equation, but this is a minor issue I believe.
In my situation, the porosity is very high, so for the most part I have liquid with a small perturbation due to a little bit of porous medium, so we are computing this small perturbation to the usual viscous steady flow.

The biggest problem I have is that I can't find a numerical scheme which ensures zero divergence of my flow. I solve the viscous flow equation by relaxation and then at each iteration I add a correction proportional to the grad of pressure from the previous iteration. The speed converges, while the divergence of speed keeps increasing with every step, and so the solution never settles to some stable situation

 

[Chestermiller]

Let me make sure I understand correctly. The differential equations you are solving are:

$$v_x=-\frac{k(x,z)}{μ}\frac{\partial p}{\partial x}$$
$$v_z=-\frac{k(x,z)}{μ}\frac{\partial p}{\partial z}$$

Is this correct?

Chet

 

[Irid]

The equations are

$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial z^2} = k(x,z) u$

$\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial z^2} = k(x,z) v$


$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial z} = 0$

 

[Chestermiller]

The equations are

$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial z^2} = k(x,z) u$

$\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial z^2} = k(x,z) v$


$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial z} = 0$

I don't quite understand your model structure, but I do understand the mathematics.

Have you tried solving this using stream function and vorticity? That would enable you to automatically satisfy the continuity equation, and would be pretty easy to set up and solve. Of course, that approach wouldn't work in 3D, but at this point you seem desperate to get an answer. Would you be willing to try that?

Chet

 

[Irid]

How easy would that be? I can't find examples of something similar having been done. I get two coupled 3rd order differential equations if I introduce vorticity.

 

[Chestermiller]

How easy would that be? I can't find examples of something similar having been done. I get two coupled 3rd order differential equations if I introduce vorticity.

$$u=\frac{\partial \psi}{\partial z}$$
$$v=-\frac{\partial \psi}{\partial x}$$
$$\frac{\partial \omega}{\partial z}=k\frac{\partial \psi}{\partial z}$$
$$\frac{\partial \omega}{\partial x}=k\frac{\partial \psi}{\partial x}$$
where
$$ω=\nabla^2 \psi$$
Take the partial of the first equation with respect to z and add the partial of the second equations with respect to x:

$$\nabla^2 \omega=kω+\frac{\partial k}{\partial x}\frac{\partial \psi}{\partial x}+\frac{\partial k}{\partial z}\frac{\partial \psi}{\partial z}$$

The last two equations would be solved simultaneously for the two unknowns (they shouldn't be combined). You would have to work out how to express the boundary conditions numerically. These equations are very ellyptic looking. I think the best way to solve them would be to use a linear equation solver (undoubetedly banded). You would have to work out how to order the unknowns in the solution vector so that the band width is as narrow as possible.

Chet

 

[Irid]

$$u=\frac{\partial \psi}{\partial z}$$
$$v=-\frac{\partial \psi}{\partial x}$$
$$\frac{\partial \omega}{\partial z}=k\frac{\partial \psi}{\partial z}$$
$$\frac{\partial \omega}{\partial x}=k\frac{\partial \psi}{\partial x}$$
where
$$ω=\nabla^2 \psi$$
Take the partial of the first equation with respect to z and add the partial of the second equations with respect to x:

$$\nabla^2 \omega=kω+\frac{\partial k}{\partial x}\frac{\partial \psi}{\partial x}+\frac{\partial k}{\partial z}\frac{\partial \psi}{\partial z}$$

The last two equations would be solved simultaneously for the two unknowns (they shouldn't be combined). You would have to work out how to express the boundary conditions numerically. These equations are very ellyptic looking. I think the best way to solve them would be to use a linear equation solver (undoubetedly banded). You would have to work out how to order the unknowns in the solution vector so that the band width is as narrow as possible.

Chet

Hmm.. I'll give it a try. Right now I'm looking forward to learn some FreeFEM, because it's going to be too time consuming to write all the scripts myself..

 

[Chestermiller]

Hmm.. I'll give it a try. Right now I'm looking forward to learn some FreeFEM, because it's going to be too time consuming to write all the scripts myself..

I think I could code this and have it running in about a day. It might even be possible to solve the equations using relaxation, which would shorten the time even more.

I'm very confident that this formulation will solve your problem (since I've solved many, many, equation sets almost exactly like this one). If you would like some advice in how the set up the boundary conditions for the numerical scheme, I'd be glad to help. Of course, away from the boundary, you would be using central differences.

Chet