Draft PR: one-sided fluid-shell interaction

[WeiyiVirtonomy]

Description

For one-sided fluid-shell interaction problems (e.g. blood flowing in elastic shell vessel), the shell particles cannot prevent penetration of fluid particles due to incompleteness of kernel support.

TODO

• Modifying dW_ijV_j and e_ij of shell neighbors by counting in the contribution of 'ghost particles' in the cut off radius of fluid particles

[WeiyiVirtonomy]

@BenceVirtonomy @FabienPean-Virtonomy @JohnVirtonomy @YuVirtonomy I'm working on this branch

[WeiyiVirtonomy]

Test on a 2D channel flow: Penetration can be prevented if we directly add dummy particles to the neighborhood:

But if we use the equivalent dW_ijV_j and e_ij of real shell particle and dummy particles, extremely high velocity in y direction appears near the wall, and some particles fly out of the computational domain. It seems that the force from the wall is incorrect. Magnitude of velocity:

vel_y:

The calculation of equivalent dW_ijV_j and e_ij might be improper. Also, in BaseFluidShellIntegration1stHalf<Integration1stHalfType>::interaction, we have Real p_in_shell = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_shell_external_acceleration);. Maybe the change in r_ij should also be taken into consideration?

[Xiangyu-Hu]

Test on a 2D channel flow: Penetration can be prevented if we directly add dummy particles to the neighborhood: But if we use the equivalent dW_ijV_j and e_ij of real shell particle and dummy particles, extremely high velocity in y direction appears near the wall, and some particles fly out of the computational domain. It seems that the force from the wall is incorrect. Magnitude of velocity: vel_y:

The calculation of equivalent dW_ijV_j and e_ij might be improper. Also, in BaseFluidShellIntegration1stHalf<Integration1stHalfType>::interaction, we have Real p_in_shell = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_shell_external_acceleration);. Maybe the change in r_ij should also be taken into consideration?

There must a big bug there, because the dummy particle is simply the degenerated case of the new model. The baseline model ( if the new summation not activated) should reproduce results using the dummy particle.

[DrChiZhang]

Test on a 2D channel flow: Penetration can be prevented if we directly add dummy particles to the neighborhood: But if we use the equivalent dW_ijV_j and e_ij of real shell particle and dummy particles, extremely high velocity in y direction appears near the wall, and some particles fly out of the computational domain. It seems that the force from the wall is incorrect. Magnitude of velocity: vel_y: The calculation of equivalent dW_ijV_j and e_ij might be improper. Also, in BaseFluidShellIntegration1stHalf<Integration1stHalfType>::interaction, we have Real p_in_shell = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_shell_external_acceleration);. Maybe the change in r_ij should also be taken into consideration?

There must a big bug there, because the dummy particle is simply the degenerated case of the new model. The baseline model ( if the new summation not activated) should reproduce results using the dummy particle.

This issue may due to the viscous force calculation. Did you use the code from this draft PR. https://github.com/Xiangyu-Hu/SPHinXsys/pull/184

[WeiyiVirtonomy]

Test on a 2D channel flow: Penetration can be prevented if we directly add dummy particles to the neighborhood: But if we use the equivalent dW_ijV_j and e_ij of real shell particle and dummy particles, extremely high velocity in y direction appears near the wall, and some particles fly out of the computational domain. It seems that the force from the wall is incorrect. Magnitude of velocity: vel_y: The calculation of equivalent dW_ijV_j and e_ij might be improper. Also, in BaseFluidShellIntegration1stHalf<Integration1stHalfType>::interaction, we have Real p_in_shell = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_shell_external_acceleration);. Maybe the change in r_ij should also be taken into consideration?

There must a big bug there, because the dummy particle is simply the degenerated case of the new model. The baseline model ( if the new summation not activated) should reproduce results using the dummy particle.

This issue may due to the viscous force calculation. Did you use the code from this draft PR. #184

Right, I pulled codes from this PR. Is there any update in viscous force calculation?

[WeiyiVirtonomy]

Test on a 2D channel flow: Penetration can be prevented if we directly add dummy particles to the neighborhood: But if we use the equivalent dW_ijV_j and e_ij of real shell particle and dummy particles, extremely high velocity in y direction appears near the wall, and some particles fly out of the computational domain. It seems that the force from the wall is incorrect. Magnitude of velocity: vel_y: The calculation of equivalent dW_ijV_j and e_ij might be improper. Also, in BaseFluidShellIntegration1stHalf<Integration1stHalfType>::interaction, we have Real p_in_shell = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_shell_external_acceleration);. Maybe the change in r_ij should also be taken into consideration?

There must a big bug there, because the dummy particle is simply the degenerated case of the new model. The baseline model ( if the new summation not activated) should reproduce results using the dummy particle.

There was a bug in the loop. It should work now The result after fixing the bug:

[DrChiZhang]

Test on a 2D channel flow: Penetration can be prevented if we directly add dummy particles to the neighborhood: But if we use the equivalent dW_ijV_j and e_ij of real shell particle and dummy particles, extremely high velocity in y direction appears near the wall, and some particles fly out of the computational domain. It seems that the force from the wall is incorrect. Magnitude of velocity: vel_y: The calculation of equivalent dW_ijV_j and e_ij might be improper. Also, in BaseFluidShellIntegration1stHalf<Integration1stHalfType>::interaction, we have Real p_in_shell = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_shell_external_acceleration);. Maybe the change in r_ij should also be taken into consideration?

There must a big bug there, because the dummy particle is simply the degenerated case of the new model. The baseline model ( if the new summation not activated) should reproduce results using the dummy particle.

This issue may due to the viscous force calculation. Did you use the code from this draft PR. #184

Right, I pulled codes from this PR. Is there any update in viscous force calculation?

I remember some new algorithm is introduced in that part, which is not fully investigated. if it works now, it should correct.

[WeiyiVirtonomy]

@DongWuTUM We can update the formulation of shell curvature in this file: tutorials/markdown/one-side-fluid-shell-interaction.md Thanks : )

[WeiyiVirtonomy]

I checked $\sumj \nabla W{ij} V_j$ (dW_ijV_j e_ij, j includes all fluid and shell particles in the neighborhood) of the fluid particles, which should be 0. Near the shell boundary, dW_ijV_je_ij is non-zero, which indicates that the equivalent value is too large. Cylinder: The magnitude of dW_ijV_j*e_ij near the wall is about 520 near the wall, and the non-zero value caused by insufficient support at inlet and outlet is about 1500.

FDA nozzel:

[Xiangyu-Hu]

I checked ∑j∇WijVj (dW_ijV_j * e_ij, j includes all fluid and shell particles in the neighborhood) of the fluid particles, which should be 0. Near the shell boundary, dW_ijV_j_e_ij is non-zero, which indicates that the equivalent value is too large. Cylinder: The magnitude of dW_ijV_j_e_ij near the wall is about 520 near the wall, and the non-zero value caused by insufficient support at inlet and outlet is about 1500. FDA nozzel:

I am fully understand what you mean. Could you make it for specifically?

[WeiyiVirtonomy]

There were some mistakes in the parameters of FDA nozzle. It works now after correcting those parameters. The velocity distribution is very close to that of the case with volumetric solid boundary

Shell boundary:

Solid boundary:

Axial velocity:

I also changed neighborhood.W_ij_ to the summation of the real particle and ghost particles, which is used in density summation. It seems that this makes the simulation more stable

[WeiyiVirtonomy]

I checked ∑j∇WijVj (dW_ijV_j * e_ij, j includes all fluid and shell particles in the neighborhood) of the fluid particles, which should be 0. Near the shell boundary, dW_ijV_j_e_ij is non-zero, which indicates that the equivalent value is too large. Cylinder: The magnitude of dW_ijV_j_e_ij near the wall is about 520 near the wall, and the non-zero value caused by insufficient support at inlet and outlet is about 1500. FDA nozzel:

I am fully understand what you mean. Could you make it for specifically?

I was thinking that $\nabla 1=\sumj \nabla W{ij} V_j \approx \mathbf{0}$ should be fulfilled for each fluid particle with the contribution of ghost shell particles, but it looks like it's not a zero vector. Also, could you please explain a bit more about the integration you mentioned in last week's meeting? I still don't quite understand that. Thanks a lot.

[Xiangyu-Hu]

There were some mistakes in the parameters of FDA nozzle. It works now after correcting those parameters. The velocity distribution is very close to that of the case with volumetric solid boundary

Shell boundary:

Solid boundary:

Axial velocity:

I also changed neighborhood.W_ij_ to the summation of the real particle and ghost particles, which is used in density summation. It seems that this makes the simulation more stable

Good. and Yes. the density summation you have mentioned should gives better results.

[WeiyiVirtonomy]

Some problems:

1. It seems that shell resolution cannot be larger than fluid resolution. This is the Hagen flow in a rigid shell wall boundary when resolution_wall > resolution_fluid:

   

2. Test of the flow in an elastic tube fixed at the inlet and outlet. Fluid particles penetrate the wall just behind the constraint part (might be fixed if the shell particles are relaxed):

   

   When setting resolution_shell = 0.5*resolution_fluid, the penetration is delayed, but shell and fluid will properties diverge to NaN later.

I hope that we can discuss these problems tomorrow.

[Xiangyu-Hu]

The shell simulation has hourglass artifacts. Did you update the shell normal direction at each time step?

[WeiyiVirtonomy]

Check $\sumj W{ij}V_j \approx 1$ and $\sumj \nabla W{ij}V_j \approx \mathbf{0}$:

Volumetric solid wall, with relaxed fluid particles:

• $\sumj W{ij}V_j$

The inner fluid particles without contact to the wall, the total W_ijV_j is about 0.79, while the value of particles near the wall is about 0.94. I'm not sure about the reason, but the total number of inner and wall contact particles is larger near the boundary.

• $\sumj \nabla W{ij}V_j$

The gradient summation near the wall is much higher than the inner region. The magnitude is about 476.

Shell, with relaxed fluid particles:

• $\sumj W{ij}V_j$

The summation for inner particles is 0.79, for particles near the wall is 0.89.

• $\sumj \nabla W{ij}V_j$

  The gradient summation for particles near the wall is about 285.

Volumetric solid, fluid particles not relaxed

• $\sumj W{ij}V_j$

The summation for inner particles is 0.82, for particles near the wall varies from 0.81~0.98.

• $\sumj \nabla W{ij}V_j$

The gradient summation for inner particles is 0, for particles near the wall is about 595.

Shell, fluid particles not relaxed

• $\sumj W{ij}V_j$

The summation for inner particles is 0.82, for particles near the wall varies from 0.76~0.94.

• $\sumj \nabla W{ij}V_j$

The gradient summation for inner particles is 0, for particles near the wall is about 515.

To conclude, for cylinder, the summation of W_ijV_j of volumetric solid is closer to 1 than the shell, and the difference is larger when the fluid particles are relaxed. The magnitude of dW_ijV_je_ij of volumetric solid is also higher. I think this means that the contribution of "ghost" shell particles is still smaller than that of volumetric solid.

[Xiangyu-Hu]

When you relax the fluid particles, did you use the wall particles as boundary particles and compute the force from them?

[WeiyiVirtonomy]

The problem of unequal shell thickness and shell resolution:

Take shell thickness=2 shell_resolution of a cylinder wall boundary as an example. The closest distance between a shell particle and a fluid particle is not 0.5 (fluid_resolution+shell_resolution), but 0.5 (fluid_resolution+shell thickness), so that the radius of fluid domain is expanded by 0.5shell thickness behind the inlet buffer. I think that's why the accuracy of fluid velocity becomes lower in this case.

[WeiyiVirtonomy]

When you relax the fluid particles, did you use the wall particles as boundary particles and compute the force from them?

I'm not sure if I understood you correctly. Did you mean computing the summation of W_ij V_j and dW_ijV_j e_ij of boundary particles? I used the summation of all particles in inner configuration and contact configuration.

[Xiangyu-Hu]

When you relax the fluid particles, did you use the wall particles as boundary particles and compute the force from them?

I'm not sure if I understood you correctly. Did you mean computing the summation of W_ij V_j and dW_ijV_j e_ij of boundary particles? I used the summation of all particles in inner configuration and contact configuration.

I mean when you relax the fluid particles, the relaxation force from wall particle should be considered too.

[WeiyiVirtonomy]

When you relax the fluid particles, did you use the wall particles as boundary particles and compute the force from them?

I'm not sure if I understood you correctly. Did you mean computing the summation of W_ij V_j and dW_ijV_j e_ij of boundary particles? I used the summation of all particles in inner configuration and contact configuration.

I mean when you relax the fluid particles, the relaxation force from wall particle should be considered too.

I was not using RelaxationStepComplex, but the particle distribution looks still very weird after using this class. Fluid particles cluster near the boundary, and there's a gap between wall particles and fluid particles. Do you have any examples of 3D complex relaxation?

[Xiangyu-Hu]

When you relax the fluid particles, did you use the wall particles as boundary particles and compute the force from them?

I'm not sure if I understood you correctly. Did you mean computing the summation of W_ij V_j and dW_ijV_j e_ij of boundary particles? I used the summation of all particles in inner configuration and contact configuration.

I mean when you relax the fluid particles, the relaxation force from wall particle should be considered too.

I was not using RelaxationStepComplex, but the particle distribution looks still very weird after using this class. Fluid particles cluster near the boundary, and there's a gap between wall particles and fluid particles. Do you have any examples of 3D complex relaxation?

You need to know how to use it correctly.

[WeiyiVirtonomy]

When you relax the fluid particles, did you use the wall particles as boundary particles and compute the force from them?

I'm not sure if I understood you correctly. Did you mean computing the summation of W_ij V_j and dW_ijV_j e_ij of boundary particles? I used the summation of all particles in inner configuration and contact configuration.

I mean when you relax the fluid particles, the relaxation force from wall particle should be considered too.

I was not using RelaxationStepComplex, but the particle distribution looks still very weird after using this class. Fluid particles cluster near the boundary, and there's a gap between wall particles and fluid particles. Do you have any examples of 3D complex relaxation?

You need to know how to use it correctly.

It seems that I can't use the exact shape of fluid block in defineComponentLevelSetShape, otherwise the level set correction will wrongly constrain particles near the wall boundary. If level_set_correction is set to false, then the particles near inlet and outlet won't be correctly constraint. Would you mind telling me who is in charge of relaxation, so that I might ask them about how to set the level set shape?

[Xiangyu-Hu]

I think that you need some to check how the code works.

[WeiyiVirtonomy]

Test on 2D channel flow with an elastic gate fails

https://github.com/Xiangyu-Hu/SPHinXsys/tree/documentation/tests/user_examples/test_2d_simplified_valve_shell

I'm testing a 2D simplified valve model, but the simulation crashes soon. Some singularities with extremely high velocity appear near the beam. The beam then folds into a strange shape. Flow field becomes very messy, and some fluid particles fly away.

This happens no matter if I take the ghost particle contribution into consideration or not. @ChiZhangatTUM I remember that SPHinXsys had a similar test case of fluid-shell interaction before https://www.youtube.com/watch?v=VpzV8dhn2Vs, which worked well. It would be helpful if someone could tell me what were the parameters you used.

[WeiyiVirtonomy]

I have created a new branch from master without codes from the old fluid-shell interaction branch (shell_fluid_interaction.h and shell_fluid_interaction.cpp are kept though, since I need these classes to calculate force from fluid on shell). I'll create a new PR and abandon this one later. https://github.com/Xiangyu-Hu/SPHinXsys/tree/feature/one_sided_fluid_shell_interaction