[paulathikalam]

Hi,

I'm attempting to incorporate a sink term into the momentum equation.

For this, I modified the porosity term in the iowave_velsource.cpp file.

Is this procedure correct? For example, will the following code work if I want to include a derivative and specific constants?

porousterm= 0.5 * a->u(i, j, k) * fabs(a->u(i, j, k)) + 0.5 * ((a->u(i+1, j+1, k+1) - a->u(i-1, j-1, k-1)) / p->dt*2);

for others directions i have kept it as

porousterm= 0;

I'm not familiar with modifying such scripts, so this approach could be incorrect. Can anyone weigh in on this?

Paul

[valgrinda]

Just fyi, porous media (Darcy or VRANS) is already included in the CFD model. Please check the B 240 options onwards.

Have a look at this paper: https://reef3d.files.wordpress.com/2...cej_porous.pdf

[paulathikalam]

Dear Prof. Hans,

Thanks for the reply.

I had actually gone through*the paper before, and the reason I did not want to use the VRANS was due to the length scale constraints l<<r<<L. This condition is not true in the case I'm attempting to work with.

I'd like to have a sink term like the one in the VRANS equation. (vrans_veg_velsource.cpp)

Fi = p->B309*0.25*PI*pow(D_val,2.0)*N_val*((a->u(i,j,k) - un(i,j,k))/p->dt);
Fd = 0.5*Cd_val*N_val*D_val*a->u(i,j,k)*fabs(a->u(i,j,k));

But I need it in iowave_velsource.cpp.

the porous term modification in iowave_velsource.cpp works fine without the derivative but the moment when i plug in the derivative,
((a->u(i+1, j+1, k+1) - a->u(i-1, j-1, k-1)) / p->dt*2)
the simulation crashes with velocity exceeding N61 error. I'd be grateful if you could help me with this.

Thanks

Paul

[Felix_Sp]

The Darcy options are also in this study https://www.sciencedirect.com/science/article/pii/S2214166915000028 or here: https://reef3d.files.wordpress.com/2...3d_oe-2015.pdf.

Why do you want to calculate the time derivative of the node at (i+1,j+1,j+1), so the node in the top, right, back, and the node (i-1,j-1,j-1), so the node in the bottom, left, front? You know that, i, j, k are the indices of the grid node right? That is why in the vrans_veg_velsource.cpp there is an u(i, j, k) and an un(i, j, k) (the velocity of the prior time step; see vrans_veg.cpp).

Why don't you just set B295 to 0, B 308 to 0 and Cd/n/D of B 309 to a value that results in 0.5*Cd_val*N_val*D_val being 0.5, so 1, 1, 1. Also set B 309 to a value that leads to p->B309*0.25*PI*pow(D_val,2.0)*N_val being your Cm value (which is one at the moment). Should result in the thing you want.

Or just change Fi in the vrans_veg_velsource.cpp to your liking and set B295 to 0, B 308 to 0. So sth like Fi = 0.5*a->u(i, j, k)*fabs(a->u(i, j, k))+0.5*((a->u(i, j, k) - un(i, j, k))/p->dt*2). And Fd = 0.

Hope this helps somehow
Cheers

[valgrinda]

Darcy porosity is implemented in B240 onwards.

[paulathikalam]

Thank you Felix and hans.

Why do you want to calculate the time derivative of the node at (i+1,j+1,j+1), so the node in the top, right, back, and the node (i-1,j-1,j-1), so the node in the bottom, left, front? You know that, i, j, k are the indices of the grid node right? That is why in the vrans_veg_velsource.cpp there is an u(i, j, k) and an un(i, j, k) (the velocity of the prior time step; see vrans_veg.cpp).

As said before, I am still learning how things work.

Why don't you just set B295 to 0, B 308 to 0 and Cd/n/D of B 309 to a value that results in 0.5*Cd_val*N_val*D_val being 0.5, so 1, 1, 1. Also set B 309 to a value that leads to p->B309*0.25*PI*pow(D_val,2.0)*N_val being your Cm value (which is one at the moment). Should result in the thing you want.

This suggestion would work.

Originally Posted by Felix_Sp

The Darcy options are also in this study https://www.sciencedirect.com/science/article/pii/S2214166915000028 or here: https://reef3d.files.wordpress.com/2...3d_oe-2015.pdf.

Why do you want to calculate the time derivative of the node at (i+1,j+1,j+1), so the node in the top, right, back, and the node (i-1,j-1,j-1), so the node in the bottom, left, front? You know that, i, j, k are the indices of the grid node right? That is why in the vrans_veg_velsource.cpp there is an u(i, j, k) and an un(i, j, k) (the velocity of the prior time step; see vrans_veg.cpp).

Why don't you just set B295 to 0, B 308 to 0 and Cd/n/D of B 309 to a value that results in 0.5*Cd_val*N_val*D_val being 0.5, so 1, 1, 1. Also set B 309 to a value that leads to p->B309*0.25*PI*pow(D_val,2.0)*N_val being your Cm value (which is one at the moment). Should result in the thing you want.

Or just change Fi in the vrans_veg_velsource.cpp to your liking and set B295 to 0, B 308 to 0. So sth like Fi = 0.5*a->u(i, j, k)*fabs(a->u(i, j, k))+0.5*((a->u(i, j, k) - un(i, j, k))/p->dt*2). And Fd = 0.

Hope this helps somehow
Cheers