# Sink Term

 Register Blogs Members List Search Today's Posts Mark Forums Read

 December 5, 2023, 03:42 Sink Term #1 New Member   Join Date: Sep 2020 Posts: 25 Rep Power: 5 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

 December 11, 2023, 13:37 #2 Super Moderator   Hans Bihs Join Date: Jun 2009 Location: Trondheim, Norway Posts: 377 Rep Power: 17 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 __________________ Hans Bihs Team REEF3D www.reef3d.com

 December 12, 2023, 00:17 #3 New Member   Join Date: Sep 2020 Posts: 25 Rep Power: 5 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<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

 December 12, 2023, 02:40 #4 Member   Felix S. Join Date: Feb 2021 Location: Germany, Braunschweig Posts: 85 Rep Power: 6 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

 December 12, 2023, 08:02 #5 Super Moderator   Hans Bihs Join Date: Jun 2009 Location: Trondheim, Norway Posts: 377 Rep Power: 17 Darcy porosity is implemented in B240 onwards. __________________ Hans Bihs Team REEF3D www.reef3d.com

December 24, 2023, 05:55
#6
New Member

Join Date: Sep 2020
Posts: 25
Rep Power: 5
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.

Quote:
 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