search

Xiangyu-Hu / SPHinXsys

SPHinXsys provides C++ APIs for engineering simulation and optimization. It aims at complex systems driven by fluid, structure, multi-body dynamics and beyond. The multi-physics library is based on a unique and unified computational framework by which strong coupling has been achieved for all involved physics.
https://www.sphinxsys.org/
Apache License 2.0
259 stars 199 forks source link

Fluid shell interaction with contribution of ghost particles #465

Closed WeiyiVirtonomy closed 2 months ago

WeiyiVirtonomy commented 6 months ago

Same as PR #453 https://github.com/Xiangyu-Hu/SPHinXsys/pull/453

WeiyiVirtonomy commented 6 months ago

Observation: weird shell particle position

I made a test on Turek's fsi benchmark case. Under some conditions, the shell can have weird deformation:

image image

The positions of shell particles are quite jaggy, and fluid particles fill in the spaces between two shell particles. I used several sets of parameters to see what caused this problem.

It seems that the problem might be related to the ratio between fluid spacing and shell thickness. It happens when thickness>fluid spacing. @DongWuTUM Do you have any idea what might cause shell particles to have this weird distribution?

Xiangyu-Hu commented 6 months ago

It looks like Hourglass mode.

DongWuTUM commented 6 months ago

The results seem to have the instability of hourglass modes. However, it occurs only when the shell thickness is different from the fluid resolution. You can check if the contact force is correct or not because the hourglass modes may occur if the contact force is enormous or the distribution of force is extremely uneven. You can also try the hourglass control algorithm by "Dynamics1Level stress_relaxation_first_half(plate_body_inner, 3, true)“ referring to 2d plate case.

WeiyiVirtonomy commented 6 months ago

The results seem to have the instability of hourglass modes. However, it occurs only when the shell thickness is different from the fluid resolution. You can check if the contact force is correct or not because the hourglass modes may occur if the contact force is enormous or the distribution of force is extremely uneven. You can also try the hourglass control algorithm by "Dynamics1Level stress_relaxation_first_half(plate_body_inner, 3, true)“ referring to 2d plate case.

Thanks a lot for the help. Yes, there were quite large and uneven forces from fluid before the hourglass modes occurred. Actually, I also met this when fluid resolution is equal to shell resolution, but it will soon disappear.

image

I have already set hourglass control to true, but this problem still exists. In the 2d aortic valve, the hourglass modes become so serious that the simulation soon crashes. Is there any other way to solve this problem?

DongWuTUM commented 6 months ago

Try to increase hourglass_controlfactor? But I don't think it is a good idea. Seems like quite large and uneven forces are the problems. Could we find a method to make the force smooth?

WeiyiVirtonomy commented 6 months ago

Try to increase hourglass_controlfactor? But I don't think it is a good idea. Seems like quite large and uneven forces are the problems. Could we find a method to make the force smooth?

I increased the coefficient to 0.008. It seems that the hourglass modes are suppressed to a lower level, but as you said, the best way is to make sure fluid forces on the shell have a smoother distribution. Maybe we can try to average the fluid force in the inner neighbor of a shell particle, but I'm not sure if this is a good idea.

image
Xiangyu-Hu commented 6 months ago

Please make sure the reason, or where is the exact problem, before new modification. Because very often some issue is merely due to a bug.

WeiyiVirtonomy commented 6 months ago

I'll not modify hourglass coefficient at the moment. The problem in 2D aortic valve mainly comes from U_max. It looks good after changing U_max to 2*U_f. I'll try to tune parameters of fsi2 as well

image
WeiyiVirtonomy commented 6 months ago

3D aortic valve:

The simulation can run with larger U_max and physical coefficient, but the particle distribution of shell near the boundary is not so good

image image
WeiyiVirtonomy commented 6 months ago

Kernel completeness check

Fluid particles are relaxed using RelaxionStepComplex this time. The initial $\sumj W{ij}V_j$ at this corner is too large (>1).

image

Very large $\sumj dW{ij}Vj\mathbf{e}{ij}$:

image

Even for the cylinder, $\sumj dW{ij}Vj\mathbf{e}{ij}$ is larger than that of the solid wall. The repelling force given by the shell boundary might be too strong. And I doubt the use of curvature, though volume is increased or decreased accordingly, $\sumj dW{ij}Vj\mathbf{e}{ij}$ is almost the same as the one using constant volume.

Hourglass of Turek's benchmark when thickness of shell>fluid_resolution

I tried to use damping or increase U_max, but it doesn't help. Sometimes using damping can make the condition even worse

WeiyiVirtonomy commented 6 months ago

It seems that the old formulation (W_ij = W_ijthickness, dW_ijV_j=dW_ijV_jthickness) performs better when shell thickness > resolution_fluid. There are still fluid particles trapped inside shell at the beginning and the hourglass modes persist throughout the simulation, but most of the time fluid particles are kept away from the shell. What's more, when I use ghost particles, the displacement of the beam is always very small, while this one is closer to the volumetric solid results.

image image

What I think now is that when shell_thickness > fluid_resolution, this simple integration gives a larger repelling force than the ghost particle method. The first one keeps fluid particles nearly 0.5thickness away from the shell, while the second one can only keep them about 0.5 resolution_fluid away. Let's say thickness = 2 dp_fluid=2 dp_shell. The force of shell given by first formulation is F0 dp thickness=2 F0 dp dp, while the one given by ghost particle integration is $F0 dp dp+F1 dp * dp$, since the ghost particle is already 1dp more away from fluid particle, F1 will be much smaller than F0.

image

I still have no idea where the large ununiform force from fluid comes from.

image

For volumetric solid at the same time, the force is balanced from two sides. It looks like some shell particles receive force like the upper surface, some like the lower surface.

image
WeiyiVirtonomy commented 6 months ago

@Xiangyu-Hu I have updated the documentation. https://github.com/Xiangyu-Hu/SPHinXsys/tree/feature/one_sided_fluid_shell_interaction/tutorials/fluid_shell_interaction I would really appreciate it if you could check the first section when you have time.

DongWuTUM commented 6 months ago

@Xiangyu-Hu I have updated the documentation. https://github.com/Xiangyu-Hu/SPHinXsys/tree/feature/one_sided_fluid_shell_interaction/tutorials/fluid_shell_interaction I would really appreciate it if you could check the first section when you have time.

@WeiyiKong321 Weiyi, when previewing the 'one-side-fluid-shell-interaction.md' file on GitHub, some formulas are not displaying. Is there a good way to properly compile the Markdown file and ensure correct preview on GitHub?

WeiyiVirtonomy commented 6 months ago

@Xiangyu-Hu I have updated the documentation. https://github.com/Xiangyu-Hu/SPHinXsys/tree/feature/one_sided_fluid_shell_interaction/tutorials/fluid_shell_interaction I would really appreciate it if you could check the first section when you have time.

@WeiyiKong321 Weiyi, when previewing the 'one-side-fluid-shell-interaction.md' file on GitHub, some formulas are not displaying. Is there a good way to properly compile the Markdown file and ensure correct preview on GitHub?

There was no problem in the preview of vscode, not sure why they can't display properly on GitHub :( I've modified those formulas, could you please take a look again?

Xiangyu-Hu commented 6 months ago

I have read the description and thinking that we should handle one-side and two side cases differently. For one side, the integral should include all the ghost particle within cut-off radius. However, for both side case, we should only take into account the ghost particles within the thickness, because there are still fluid particles on the other side of the shell.

WeiyiVirtonomy commented 6 months ago

I have read the description and thinking that we should handle one-side and two side cases differently. For one side, the integral should include all the ghost particle within cut-off radius. However, for both side case, we should only take into account the ghost particles within the thickness, because there are still fluid particles on the other side of the shell.

I'm not sure if it will help when shell thickness =2 * fluid resolution. When the cut-off radius <= thickness, it will be the same the one-sided formulation, but we still have the hourglass modes

WeiyiVirtonomy commented 6 months ago

I have read the description and thinking that we should handle one-side and two side cases differently. For one side, the integral should include all the ghost particle within cut-off radius. However, for both side case, we should only take into account the ghost particles within the thickness, because there are still fluid particles on the other side of the shell.

If it's integrated within the thickness, then maybe it's better to use Gaussian integration. By the way, do you think if there would be any method to improve the calculation of fluid forces on the shell?

WeiyiVirtonomy commented 5 months ago

Shell self contact test

image
Xiangyu-Hu commented 5 months ago

Shell self contact test image

Great!

DongWuTUM commented 5 months ago

Shell self contact test image

Great. Thanks.

WeiyiVirtonomy commented 5 months ago

Some questions about coding style

I'm not sure which one would be better, writing new classes for shell self-contact density summation and force, or just use SelfContactDensitySummation and ContactDensitySummation.

These two classes take SurfaceContactRelation as the input variable, if I want to reuse them, then ShellSelfContactRelation must be a child class of SurfaceContactRelation, which seems unnecessary. Also, as these classes use Vol_ and mass_, I need to modify dW_ijV_j to dW_ijV_j*particle_distance_*particle_distance_ in NeighborBuilderShellSelfContact. If you change those variables to ParticleVolume and ParticleMass someday, then NeighborBuilderShellSelfContact will also need modification, which might be inconvenient.

Xiangyu-Hu commented 5 months ago

Some questions about coding style

I'm not sure which one would be better, writing new classes for shell self-contact density summation and force, or just use SelfContactDensitySummation and ContactDensitySummation.

These two classes take SurfaceContactRelation as the input variable, if I want to reuse them, then ShellSelfContactRelation must be a child class of SurfaceContactRelation, which seems unnecessary. Also, as these classes use Vol_ and mass_, I need to modify dW_ijV_j to dW_ijV_j*particle_distance_*particle_distance_ in NeighborBuilderShellSelfContact. If you change those variables to ParticleVolume and ParticleMass someday, then NeighborBuilderShellSelfContact will also need modification, which might be inconvenient.

@WeiyiKong321 I would prefer not change the local dynamics but introduce new body relations. The reason is that we have so many different local dynamics and only handful body relations. Also the idea of the onside shell model is more on geometric other than a dynamic method.

WeiyiVirtonomy commented 5 months ago

Some questions about coding style

I'm not sure which one would be better, writing new classes for shell self-contact density summation and force, or just use SelfContactDensitySummation and ContactDensitySummation. These two classes take SurfaceContactRelation as the input variable, if I want to reuse them, then ShellSelfContactRelation must be a child class of SurfaceContactRelation, which seems unnecessary. Also, as these classes use Vol_ and mass_, I need to modify dW_ijV_j to dW_ijV_j*particle_distance_*particle_distance_ in NeighborBuilderShellSelfContact. If you change those variables to ParticleVolume and ParticleMass someday, then NeighborBuilderShellSelfContact will also need modification, which might be inconvenient.

@WeiyiKong321 I would prefer not change the local dynamics but introduce new body relations. The reason is that we have so many different local dynamics and only handful body relations. Also the idea of the onside shell model is more on geometric other than a dynamic method.

Thanks, I got it. May I change the dyanamic classes from SelfContactDensitySummation(SelfSurfaceContactRelation &self_contact_relation) to SelfContactDensitySummation(BaseInnerRelation &self_contact_relation)? I think it looks not so good to make ShellSelfContactRelation a child class of SelfSurfaceContactRelation, since there's no clear inheritance relation between them.

Xiangyu-Hu commented 5 months ago

Some questions about coding style

I'm not sure which one would be better, writing new classes for shell self-contact density summation and force, or just use SelfContactDensitySummation and ContactDensitySummation. These two classes take SurfaceContactRelation as the input variable, if I want to reuse them, then ShellSelfContactRelation must be a child class of SurfaceContactRelation, which seems unnecessary. Also, as these classes use Vol_ and mass_, I need to modify dW_ijV_j to dW_ijV_j*particle_distance_*particle_distance_ in NeighborBuilderShellSelfContact. If you change those variables to ParticleVolume and ParticleMass someday, then NeighborBuilderShellSelfContact will also need modification, which might be inconvenient.

@WeiyiKong321 I would prefer not change the local dynamics but introduce new body relations. The reason is that we have so many different local dynamics and only handful body relations. Also the idea of the onside shell model is more on geometric other than a dynamic method.

Thanks, I got it. May I change the dyanamic classes from SelfContactDensitySummation(SelfSurfaceContactRelation &self_contact_relation) to SelfContactDensitySummation(BaseInnerRelation &self_contact_relation)? I think it looks not so good to make ShellSelfContactRelation a child class of SelfSurfaceContactRelation, since there's no clear inheritance relation between them.

Ok.

WeiyiVirtonomy commented 5 months ago

Hydrostatic fsi test

I build the wall boundary also as shell, so that I can use the fluid damping to run this hydrostatic test.

shell_resolution=fluid_resolution=thickness

When the thickness of elastic gate is 0.05, the maximum displacement is 1.0e-4, while the volumetric solid result 1.2e-4. When the thickness is 0.025, the maximum displacement is 6.7e-4, solid result is 7.8e-4. The force on shell from fluid is about 916.

image

shell_resolution=fluid_resolution=0.5 * thickness

However, when the shell resolution = 0.5 * thickness = 0.025, the maximum displacement becomes 1.7e-4, which is too large. The maximum force on shell from fluid is about 914.

shell_resolution=0.5 fluid_resolution=0.5 thickness

This time the deformation becomes 1.1e-4. The force on shell from fluid is about 585.

It seems that the resolution and distance do have some influence on the result.

Xiangyu-Hu commented 5 months ago

Hydrostatic fsi test

I build the wall boundary also as shell, so that I can use the fluid damping to run this hydrostatic test.

shell_resolution=fluid_resolution=thickness

When the thickness of elastic gate is 0.05, the maximum displacement is 1.0e-4, while the volumetric solid result 1.2e-4. When the thickness is 0.025, the maximum displacement is 6.7e-4, solid result is 7.8e-4. The force on shell from fluid is about 916. image

shell_resolution=fluid_resolution=0.5 * thickness

However, when the shell resolution = 0.5 * thickness = 0.025, the maximum displacement becomes 1.7e-4, which is too large. The maximum force on shell from fluid is about 914.

shell_resolution=0.5 fluid_resolution=0.5 thickness

This time the deformation becomes 1.1e-4. The force on shell from fluid is about 585.

It seems that the resolution and distance do have some influence on the result.

How far from the reference solution?

WeiyiVirtonomy commented 5 months ago

Question about acc_prior caused by force from fluid

The pressure force on solid from fluid is:

Real face_wall_external_acceleration = (acc_prior_k[index_j] - acc_ave_[index_i]).dot(e_ij);
Real p_in_wall = p_k[index_j] + rho_n_k[index_j] * r_ij * SMAX(Real(0), face_wall_external_acceleration);
Real u_jump = 2.0 * (vel_k[index_j] - vel_ave_[index_i]).dot(n_[index_i]);
force += (riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] - (p_in_wall + p_k[index_j]) * e_ij) * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n];

For hydrostatic problem, riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] and rho_n_k[index_j] * r_ij * SMAX(Real(0), face_wall_external_acceleration) are much smaller than p_k[index_j]. Neglecting these parts, the force is only related to p_k[index_j]* e_ij * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n].

The pressure of fluid particles in contact with solid particles is all about rhof g H, so the force is determined by `Vol[index_i] contact_neighborhood.dW_ijVj[n] e_ij `. The acceleration on the shell particle is then:

$$\frac{1}{\rho_s}\sumjdW{ij}Vj e{ij}$$

Assume that shell particle distance dp_s is the same as fluid particle distance dp_f, then the relative distance between shell and fluid particles q=r_ij/h_f are almost the same.

image

This term will be proportional to inv_h of fluid, i.e. proportional to particle distance of fluid:

$$\sumj dW{ij}Vje{ij} \propto \sum_j \frac{1}{dp_f^3}dW1(q{ij})dpf^2e{ij}=\sum_j \frac{1}{dp_f}dW1(q{ij})e_{ij}$$

Then the acceleration prior of the shell particle will be dependent on dp_f: $$\frac{1}{\rho_s}\sumjdW{ij}Vj e{ij} \propto \frac{1}{dp_f} \sum_j dW1(q{ij})e_{ij}$$

That's why when dp_f=dp_s=0.5*t=0.025, the displacement of shell becomes nearly twice that when dp_f=dp_s=0.05. Under the same conditions, the solutions change due to the changes in fluid resolution, which is surely wrong.

One way I can think of to fix the problem right now is to create a new kernel with a smoothing length of $dp_{avg} = 0.5(dp_f+dp_s)$. The volume of a shell particle will be modified as $A_j dp_{avg}$, instead of $A_j t$, but the mass is still the actual mass $\rho_s A_j * t$. Then the acceleration becomes:

$$\frac{V_i}{m_i}\sumjdW{ij}Vj e{ij} \propto\frac{dp_{avg}}{rhos*t}\frac{1}{dp{avg}} \sum_j dW1(q{ij})e_{ij} \propto \frac{1}{t}$$

@DongWuTUM Dong, what do you think about this?

DongWuTUM commented 5 months ago

Question about acc_prior caused by force from fluid

The pressure force on solid from fluid is:

Real face_wall_external_acceleration = (acc_prior_k[index_j] - acc_ave_[index_i]).dot(e_ij);
Real p_in_wall = p_k[index_j] + rho_n_k[index_j] * r_ij * SMAX(Real(0), face_wall_external_acceleration);
Real u_jump = 2.0 * (vel_k[index_j] - vel_ave_[index_i]).dot(n_[index_i]);
force += (riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] - (p_in_wall + p_k[index_j]) * e_ij) * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n];

For hydrostatic problem, riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] and rho_n_k[index_j] * r_ij * SMAX(Real(0), face_wall_external_acceleration) are much smaller than p_k[index_j]. Neglecting these parts, the force is only related to p_k[index_j]* e_ij * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n].

The pressure of fluid particles in contact with solid particles is all about rhof g H, so the force is determined by `Vol[index_i] contact_neighborhood.dW_ijVj[n] e_ij `. The acceleration on the shell particle is then:

1ρs∑jdWijVjeij

Assume that shell particle distance dp_s is the same as fluid particle distance dp_f, then the relative distance between shell and fluid particles q=r_ij/h_f are almost the same. image

This term will be proportional to inv_h of fluid, i.e. proportional to particle distance of fluid:

∑jdWijVjeij∝∑j1dpf3dW1(qij)dpf2eij=∑j1dpfdW1(qij)eij

Then the acceleration prior of the shell particle will be dependent on dp_f: 1ρs∑jdWijVjeij∝1dpf∑jdW1(qij)eij

That's why when dp_f=dp_s=0.5*t=0.025, the displacement of shell becomes nearly twice that when dp_f=dp_s=0.05. Under the same conditions, the solutions change due to the changes in fluid resolution, which is surely wrong.

One way I can think of to fix the problem right now is to create a new kernel with a smoothing length of dpavg=0.5∗(dpf+dps). The volume of a shell particle will be modified as Aj∗dpavg, instead of Aj∗t, but the mass is still the actual mass ρs∗Aj∗t. Then the acceleration becomes:

Vimi∑jdWijVjeij∝dpavgrhos∗t1dpavg∑jdW1(qij)eij∝1t

@DongWuTUM Dong, what do you think about this?

Have you tried the last equation? I have carefully read what you wrote. But I am still confused about "The force is determined by e_ij * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n]." The force is also related to 'p_k[index_j]', right? When particle distance is changed to smaller, the 'p_k[index_j]' is also smaller. I may have misunderstood. If it is also convenient for you, we can have a discussion through Zoom.

WeiyiVirtonomy commented 5 months ago

Question about acc_prior caused by force from fluid

The pressure force on solid from fluid is:

Real face_wall_external_acceleration = (acc_prior_k[index_j] - acc_ave_[index_i]).dot(e_ij);
Real p_in_wall = p_k[index_j] + rho_n_k[index_j] * r_ij * SMAX(Real(0), face_wall_external_acceleration);
Real u_jump = 2.0 * (vel_k[index_j] - vel_ave_[index_i]).dot(n_[index_i]);
force += (riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] - (p_in_wall + p_k[index_j]) * e_ij) * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n];

For hydrostatic problem, riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] and rho_n_k[index_j] * r_ij * SMAX(Real(0), face_wall_external_acceleration) are much smaller than p_k[index_j]. Neglecting these parts, the force is only related to p_k[index_j]* e_ij * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n]. The pressure of fluid particles in contact with solid particles is all about rhof g H, so the force is determined by `Vol[index_i] contact_neighborhood.dW_ijVj[n] e_ij `. The acceleration on the shell particle is then: 1ρs∑jdWijVjeij Assume that shell particle distance dp_s is the same as fluid particle distance dp_f, then the relative distance between shell and fluid particles q=r_ij/h_f are almost the same. image This term will be proportional to inv_h of fluid, i.e. proportional to particle distance of fluid: ∑jdWijVjeij∝∑j1dpf3dW1(qij)dpf2eij=∑j1dpfdW1(qij)eij Then the acceleration prior of the shell particle will be dependent on dp_f: 1ρs∑jdWijVjeij∝1dpf∑jdW1(qij)eij That's why when dp_f=dp_s=0.5*t=0.025, the displacement of shell becomes nearly twice that when dp_f=dp_s=0.05. Under the same conditions, the solutions change due to the changes in fluid resolution, which is surely wrong. One way I can think of to fix the problem right now is to create a new kernel with a smoothing length of dpavg=0.5∗(dpf+dps). The volume of a shell particle will be modified as Aj∗dpavg, instead of Aj∗t, but the mass is still the actual mass ρs∗Aj∗t. Then the acceleration becomes: Vimi∑jdWijVjeij∝dpavgrhos∗t1dpavg∑jdW1(qij)eij∝1t @DongWuTUM Dong, what do you think about this?

Have you tried the last equation? I have carefully read what you wrote. But I am still confused about "The force is determined by e_ij * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n]." The force is also related to 'p_k[index_j]', right? When particle distance is changed to smaller, the 'p_k[index_j]' is also smaller. I may have misunderstood. If it is also convenient for you, we can have a discussion through Zoom.

The pressure of fluid particles will be sth like rho_f g (DH - 0.5 * resolution_fluid), but since DH is much larger than fluid resolution, they will be nearly constant. It would be great if we could have a zoom meeting, but we can also discuss about it on Friday's meeting if that suits you better.

DongWuTUM commented 5 months ago

We can have a meeting now if it is Okay for you. The Zoom link is the same as before, and when you join, I will receive an email.

WeiyiVirtonomy commented 5 months ago

Normal smoothing tested on T-shaped pipe

I made some changes to the one-sided NeighborBuilder. Now it will only calculate ghost particles when the normal direction of shell is in the opposite direction of e_ij, so particles in domain 1 cannot see the ghost particle of the vertical shell boundary, while particles in domain 2 cannot see ghost of horizontal shell boundary. Therefore, there's no more overlapping of ghost particles.

image

I tried to smoothen shell normal direction by the codes in ShellNormalPrediction: smoothed_[index_i] /= temp_[index_i].norm() + TinyReal;

The results with and without normal smoothing doesn't show much difference.

No normal smoothing

image

With normal smoothing

image
WeiyiVirtonomy commented 5 months ago

Hydrostatic fsi test

I build the wall boundary also as shell, so that I can use the fluid damping to run this hydrostatic test.

shell_resolution=fluid_resolution=thickness

When the thickness of elastic gate is 0.05, the maximum displacement is 1.0e-4, while the volumetric solid result 1.2e-4. When the thickness is 0.025, the maximum displacement is 6.7e-4, solid result is 7.8e-4. The force on shell from fluid is about 916. image

shell_resolution=fluid_resolution=0.5 * thickness

However, when the shell resolution = 0.5 * thickness = 0.025, the maximum displacement becomes 1.7e-4, which is too large. The maximum force on shell from fluid is about 914.

shell_resolution=0.5 fluid_resolution=0.5 thickness

This time the deformation becomes 1.1e-4. The force on shell from fluid is about 585. It seems that the resolution and distance do have some influence on the result.

How far from the reference solution?

I changed possion ratio to 0.49, so that we can compare to the analytical solution.

Analytical solution: $$\delta_{max} = \frac{pL^4}{384EI}=7.3e-5$$

Shell with resolution_shell = resolution_fluid = shell_thickness = 0.05: 9e-5, about 23% error

Volumetric solid with resolution_shell = resolution_fluid = thickness/4 = 0.125: 6.9e-5, about 5% error

DongWuTUM commented 5 months ago

Hydrostatic fsi test

I build the wall boundary also as shell, so that I can use the fluid damping to run this hydrostatic test.

shell_resolution=fluid_resolution=thickness

When the thickness of elastic gate is 0.05, the maximum displacement is 1.0e-4, while the volumetric solid result 1.2e-4. When the thickness is 0.025, the maximum displacement is 6.7e-4, solid result is 7.8e-4. The force on shell from fluid is about 916. image

shell_resolution=fluid_resolution=0.5 * thickness

However, when the shell resolution = 0.5 * thickness = 0.025, the maximum displacement becomes 1.7e-4, which is too large. The maximum force on shell from fluid is about 914.

shell_resolution=0.5 fluid_resolution=0.5 thickness

This time the deformation becomes 1.1e-4. The force on shell from fluid is about 585. It seems that the resolution and distance do have some influence on the result.

How far from the reference solution?

I changed possion ratio to 0.49, so that we can compare to the analytical solution.

Analytical solution: δmax=pL4384EI=7.3e−5

Shell with resolution_shell = resolution_fluid = shell_thickness = 0.05: 9e-5, about 23% error

Volumetric solid with resolution_shell = resolution_fluid = thickness/4 = 0.125: 6.9e-5, about 5% error

The resolution of Fluid-shell interaction is smaller than Fluid-solid.

WeiyiVirtonomy commented 5 months ago

Curvature change leads to leakage

Fluid particles tends to escape from the corner, possibly because the volume of ghost particles is reduced too much due to negative curvature. When adding ghost particle only for those satisfying displacement.dot(n_j) <= 0:

image

It's better when adding ghost particle for every shell particle

image

No curvature calculation, add ghost particle of each shell particle:

image

No curvature calculation, add ghost particle only for those satisfying displacement.dot(n_j) <= 0: tend to penetrate

image

I think we need to carefully choose a threshold close to one for the criterion displacement.dot(n_j) <= threshold

WeiyiVirtonomy commented 5 months ago

Use ghost particles for shell-fluid contact

I think we should mimic the volumetric solid, and calculate the force from fluid on ghost shell particles as well. The force is then averaged over the thickness. The algorithm is similar to the fluid-shell contact:

Assuming that the normal direction points from fluid to shell, the position of the ghost shell particle is: $$r_{dummy} = r_i + n_i * dp_s$$ Then we calculate the summation of W_ij and eij: $$\sum{dummy} dW_{ij} Vj e{ij}$$

The process is repeated until the distance between dummy shell particle and fluid particle j has exceeded the cut-off radius of fluid kernel.

The force on shell particle: force += (riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] - (p_in_wall + p_k[index_j]) * e_ij) * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n];

Here, Vol_[index_i] the equal to shell particle area * dp_s. I think it would be better to use curvature information to calculate an equivalent volume of the shell, but I haven't tested it yet.

The prior acceleration: acc_prior_[index_i] = force / particles_->ParticleMass(index_i);

Here, the mass is the actual mass of shell, i.e. area * thickness. The total force from fluid is thus averaged over the shell thickness.

Test results

The method is tested on the hydrostatic fluid-shell interaction.

dp_s < t (t=0.05, Dam_H = 2)

dp_s > t (t=0.0125, Dam_H = 1):

Discussion

With ghost particles, the results are in the same range no matter what dp_s and dp_f we use, but there are still some problems.

If we look at case 1~4, the maximum displacement seems to be getting smaller as the resolution is refined, instead of converging to a constant value. One possible explanation is that the volume of shell particle is not corrected by curvature, so the total force becomes smaller when the layer of ghost particle increases, but in this case the beam is nearly flat, curvature shouldn't make much difference. I'll check if using an equivalent r_ij could help.

Another problem is that when the fluid particle gets too close to the shell particle in the initial step, it seems that fluid receives a large repelling force. In most cases the initial fluid particle relaxation should not be a problem, and this seems more likely to be caused by the instability of inviscid fluid, instead of fluid-shell contact.

image image

I'm not sure if this method makes sense. I hope we can have some further discussion tmr.

WeiyiVirtonomy commented 5 months ago

Suggestion

I think it would be better to change the ParticleMass() to mass_ in BasePressureForceAccelerationFromFluid:

Real face_wall_external_acceleration = (force_prior_k[index_j] / mass_k[index_j] - force_ave_[index_i] / particles_->ParticleMass(index_i)).dot(e_ij);

force_ave_ is defined as the area force for shell, so the accave of shell should be force_ave_[index_i] / particles_->mass_[index_i].

If you change it to particles_->mass_[index_i], I can reuse classes in fluid_structure_interaction.h for force from fluid on shell, instead of writing new classes.

Xiangyu-Hu commented 5 months ago

Suggestion

I think it would be better to change the ParticleMass() to mass_ in BasePressureForceAccelerationFromFluid:

Real face_wall_external_acceleration = (force_prior_k[index_j] / mass_k[index_j] - force_ave_[index_i] / particles_->ParticleMass(index_i)).dot(e_ij);

force_ave_ is defined as the area force for shell, so the accave of shell should be force_ave_[index_i] / particles_->mass_[index_i].

If you change it to particles_->mass_[index_i], I can reuse classes in fluid_structure_interaction.h for force from fluid on shell, instead of writing new classes.

Actually, I prefer to use the particle mass function because it is misleading. Anyway, we can discuss it tommorow.

WeiyiVirtonomy commented 5 months ago

About the calculation of shell curvature using the cut-off radius of fluid

I'm not sure what would be the most efficient way to do this. Should I create a new CellLinkedList of shell with the cut-off radius of fluid, or try to use the same CellLinkedList but changing the SearchDepth?

DongWuTUM commented 5 months ago

About the calculation of shell curvature using the cut-off radius of fluid

I'm not sure what would be the most efficient way to do this. Should I create a new CellLinkedList of shell with the cut-off radius of fluid, or try to use the same CellLinkedList but changing the SearchDepth?

If the cut-off radius is changed, the CellLinkedList should be changed accordingly.

WeiyiVirtonomy commented 5 months ago

The curvature calculated from fluid cut-off radius looks good at the 90 corner.

image

Without using B_matrix, the curvature of particles near the corner is quite close, and the corner curvature is slightly smaller than the 2 particles beside it. This doesn't cause problem for fluid flow, but not sure if the inaccuracy will cause problem at sharper corner.

image

If we want to use B matrix correction for curvature calculation of an elastic shell, then a new kind of B must be registered. I'll continue to test after back from vacation.

Test on 60 degree angle:

image
DongWuTUM commented 5 months ago

The curvature calculated from fluid cut-off radius looks good at the 90 corner. image

Without using B_matrix, the curvature of particles near the corner is quite close, and the corner curvature is slightly smaller than the 2 particles beside it. This doesn't cause problem for fluid flow, but not sure if the inaccuracy will cause problem at sharper corner. image

If we want to use B matrix correction for curvature calculation of an elastic shell, then a new kind of B must be registered. I'll continue to test after back from vacation.

Thanks. Have a nice holiday!

WeiyiVirtonomy commented 4 months ago

Test on flow around a flat balloon

Fluid resolution dp_f = 0.15, shell resolution dp_s = 0.075. The disthe upper and lower layers is 1 * dp_f.

image

We can't use correction matrix B for this case, otherwise the curvature of the flat part will be calculated to 0:

image

When B matrix is not used, the curvature of the upper and lower layer is 9.1 and 10.1 respectively.

Theoretically, there should be 1 ghost particle in total:

image

The length of the 1st ghost particle generated by the upper shell particle is (1-9.1 0.075) dp_s = 0.32 dp_s, while the length of the 2nd is already below 0. That generated by the lower shell particle is (1-10.1 0.075) dp_s = 0.24 dp_s. Thus, the length of the ghost particle generated in total is 0.56 dp_s, which is smaller than the theoretical value of 1 dp_s.

Although some fluid particles seem to be too close to the shell, the fluid field looks still ok.

image image

I haven't tested on smaller distances since there will be some problems with normal direction, but we can still check the curvature:

image

The distance is now 1.02 dp_s. The curvature of the flat part is 10.8 and 9.7, which doesn't differ much from the 2 dp_s case, while there shouldn't be any ghost particle in this case.

It seems that the curvature is not accurate enough when B correction is not used, but when B is used, it becomes the real curvature.

WeiyiVirtonomy commented 4 months ago

Test on flow around an elastic balloon

An elastic balloon is placed in the flow. The inlet velocity gradually increases to the maximum value from 0 to 1 sec, and maintains constant thereafter. The balloon is stiff enough that no self-contact occurs.

image

At 0.5s:

image

At 1s, the curvature correction is working properly:

image

At 1.5s:

image

At 2.5s:

image

The balloon nearly recovers to the initial shape when the flow becomes stable:

image

Update:

Test on fluid-shell interaction + shell self-contact

image
WeiyiVirtonomy commented 4 months ago

Discussion: use different mean curvature for shell self-contact and fluid-shell interaction

For shell self-contact problem, the curvature must be the actual curvature with B matrix correction. I would suggest to use two separate variables "TotalMeanCurvature" and "AverageTotalMeanCurvature" for self-contact and fluid-shell interaction, and two different classes to calculate them. What do you think about this?

Maybe we even don't need to adjust the volume of ghost particles for self-contact problem?

DongWuTUM commented 4 months ago

Thank you! Then, you could try not to use the B for curvature calculation. The results seem Okay. For the figure below, how do you get the normal direction, and could you try to get the correct one first? Or you could generate the particle by directly giving the particle position and normal direction. image

The results of the test on flow around an elastic balloon are very good. I suggest no adjustment on the volume of ghost particles for self-contact problem. But could you try both and check which one is better?

WeiyiVirtonomy commented 4 months ago

Fail on 3D pipe with sharp angle

The pipe has a sharp angle of 10 degree. Fluid resolution = radius / 20 = 0.03, shell resolution = 0.5 * fluid resolution = 0.015.

image

Very high velocity singularities near the corner, though the total kernel is not higher than the other particles. There's also no obvious particle repelling.

image

I'll try to add a fillet at the corner and see if the result will be better.

WeiyiVirtonomy commented 4 months ago

Thank you! Then, you could try not to use the B for curvature calculation. The results seem Okay. For the figure below, how do you get the normal direction, and could you try to get the correct one first? Or you could generate the particle by directly giving the particle position and normal direction. image

The results of the test on flow around an elastic balloon are very good. I suggest no adjustment on the volume of ghost particles for self-contact problem. But could you try both and check which one is better?

I was using shell particle relaxation from SPHinXsys. It's a bit difficult to assign normal for this case, I'll try on a 3D case later and use virtonomy's surface remeshing to get the normal.

I haven't seen much difference in self-contact with or without curvature change yet, but not sure what it will be like when the deformed shape has very large curvature.

DongWuTUM commented 4 months ago

Fail on 3D pipe with sharp angle

The pipe has a sharp angle of 10 degree. Fluid resolution = radius / 20 = 0.03, shell resolution = 0.5 * fluid resolution = 0.015. image

Very high velocity singularities near the corner, though the total kernel is not higher than the other particles. There's also no obvious particle repelling. image

I'll try to add a fillet at the corner and see if the result will be better.

At the corner, there should be high-velocity particles. You are right; the sharp corner may cause problems (singularity).

WeiyiVirtonomy commented 4 months ago

Fail on 3D pipe with sharp angle

The pipe has a sharp angle of 10 degree. Fluid resolution = radius / 20 = 0.03, shell resolution = 0.5 * fluid resolution = 0.015. image Very high velocity singularities near the corner, though the total kernel is not higher than the other particles. There's also no obvious particle repelling. image I'll try to add a fillet at the corner and see if the result will be better.

At the corner, there should be high-velocity particles. You are right; the sharp corner may cause problems (singularity).

dp_f = 0.4, dp_s = 0.2:

Fillet of radius 0.2

Compared with solid wall:

image

At the same time, unlike solid wall, shell wall still has velocity near the corner:

Shell:

image

Solid:

image

At 10s, the singularities seem to disappear:

Shell:

image

Solid:

image

Fillet of radius 0.3

With a fillet of radius 0.3, the result is closer to that of solid:

image image

Comparison of wall particle distribution at the corner:

image

I still don't know if the singularities are caused by geometry or shell ghost particles. Since W_ijV_j summation near corner seems to be normal, I think it's more likely to be caused by geometry.

DongWuTUM commented 4 months ago

It seems Okay when the fillet is applied on the sharp corner. I think it is the geometry problem.