Extending square droplet example at high Re to droplet impact on a solid surface

[nilot-pal]

Hi all,
I am experiencing a problem in SPHinXsys which I will attempt to describe in a single issue. However, if you want me to create additional issues, let me know and I will do likewise.
I was trying to extend the square droplet deformation example to circular droplet impact at high Re and We nos. For this, I was first trying to get the square droplet example to work at high Re. Although the equations in SPHinXsys are not non-dimensionalized, my naive understanding of the calculation of Re is as follows:
$Re = \rho_w U_{max} l/\mu_w$
where $\rho_w$ = density of heavier fluid (here water)
$U_{max}$ = maximum anticipated flow speed
$l$ = characteristic length = smoothing length = 1.3 * particle ref. spacing
$\mu_w$ = dynamic viscosity of heavier fluid
For the square droplet example, this is how I calculated $U_{max}$:
$U_{max} = max (v1,v2,v3)$
where v1, v2 and v3 are found using the balance of inertial, viscous and surface tension forces as follows:
v1 = $\sqrt(\sigma/(\rho_w*l))$; v2 = $\sigma/\mu_w$; v3 = $\mu_w/(\rho_w*l)$
Before setting up the parameters to get a high Re, I changed each from it's default value one at a time to look at the effect on the droplet behavior. Note that I've been using the acoustic Riemann solver with $\eta = 3$ for these tests. My results were as follows:

$DL$	$l=1.3*DL/40$	$\rho_w$	$\mu_w$	$\sigma$	$U_{max}$	$Re$	$Result$
2.0		1.0	0.2	1.0	$${\color{red}5.0}$$	1.625	12.12.2023_12.38.50_REC (screenrec.com)
$${\color{red}5.4e-3}$$		1.0	0.2	1.0	$${\color{red}1139.6}$$	1.0	17.12.2023_13.51.28_REC (screenrec.com)
2.0		$${\color{red}1.0e3}$$	0.2	1.0	$${\color{red}5.0}$$	1625	12.12.2023_12.44.42_REC (screenrec.com)
2.0		1.0	$${\color{red}1.0e-3}$$	1.0	$${\color{red}1.0e3}$$	65e3	12.12.2023_13.16.56_REC (screenrec.com)
2.0		1.0	0.2	$${\color{red}7.3e-2}$$	$${\color{red}3.07 }$$	1.0	12.12.2023_12.50.42_REC (screenrec.com)

The red colored values indicate they've been changed from their default values. In the result column, the red particles are water, blue particles are air and white is solid. One thing I noted in all these simulations is that the advective and acoustic time steps are equal. Also, for high Re cases, the particles just disappear. For low Re cases, the droplet behavior is quite stable (for the 2nd case, the droplet changes shape and then doesn't move at all. High dissipation, maybe?)
My questions are:

1. Is the acoustic Riemann solver not dissipative enough for high Re flows? A high $U_{max}$ means a high artificial sound speed and hence high dissipation. In the 2nd and 4th cases in my table above, $U_{max}$ is high but the droplet behaves normally only when the Re is low. As an aside, I had also tried the dissipative Riemann solver for the square droplet example (Re = 7155) and got the following result:
   30.11.2023_21.09.35_REC (screenrec.com)
2. The constraint I put on $U_{max}$ for the experiments shown in the table is:
   $U_{max} = max (v1,v2,v3)$
   For a droplet impact experiment, I thought my constraints would be the Re and We nos. I have done one such experiment with the following values of parameters:
   $\rho_w = 1.0$, $U_{max} = 715.5$, $l = 1.0$, $\mu_w = 0.2$, $\sigma = 3938$, $Re = 7155$, $We = 260$. The reasons I took these values are as follows:
   1. From my table, I thought that a high value of $\rho_w$ or a low value of $\mu_w$ don't work well.
   2. A small value of $l$ means a very small time step and hence longer time for simulation.
   3. The only parameters left were $U_{max}$ and $\sigma$ and hence their values were dictated by the Re and We nos.
   My result was as follows:
   https://screenrec.com/share/79XJkLjaVx
   Here, the droplet first spreads on the bottom wall and then the film moves upwards along the side walls. Till the formation of this film, the droplet seems to behave as expected (spreading at high Re and We is recorded in literature). After that it seems like there is not enough dissipation to stop this film from moving along the side walls.
   What might be causing this issue?

[DrChiZhang]

Hi Nilot, thx for sharing the results.
Would do you please provide:
1, the published reference of the case you tested;
2, the code you used for test.
Then, we can reproduce and analyse you issue.
thx

[nilot-pal]

@ChiZhangatTUM, sure.

1. For the droplet impact problem, I was using the following reference:
   Yang, X., & Kong, S. C. (2018). 3D simulation of drop impact on dry surface using SPH method. International Journal of Computational Methods, 15(03), 1850011.
   Note that they have used a non dimensional time step of 1e-7$\times V/D$, I however used the time step as follows:
   $\Delta t$ = min($\Delta t_1$, $\Delta t_2$, $\Delta t_3$, $\Delta t_4$, $\Delta t_{ad}$, $\Delta t_{ac}$)
   where $\Delta t_1$ = 0.25$\min_{b} \left( \frac{h}{|a_{b}|} \right)^{\frac{1}{2}}$
   $\Delta t_2$ = 0.25$\frac{h}{c_s}$
   $\Delta t_3$ = 0.25$\left( \frac{h}{|g|} \right)^{\frac{1}{2}}$
   $\Delta t_4$ = 0.25 $ \min \left( \frac{\min(\rho_k, \rho_l) \times h^3}{2\pi a} \right)^{0.5}$
   There is also this confusion whether they actually had air particles in the simulation. However the experimental paper they used for validation is:
   Rioboo, R., Marengo, M., & Tropea, C. (2002). Time evolution of liquid drop impact onto solid, dry surfaces. Experiments in fluids, 33(1), 112-124.
   Above paper describes the following experimental setup and it looks like there was no vacuum.
   

2. Code for the impact problem (currently in 2D before moving to 3D):
   This is from an older version of SPHinXsys which I modified to incorporate the additional time steps above.
   https://github.com/nilot-pal/SPHinXsys-nilo
   Folder: tests/user_examples/test_2d_droplet_impact_on_solid
   A small detail: For calculating Re, I used the diameter of the droplet which is 2.0. This is different from what I was assuming for the square droplet deformation. Prior to this, I was using the smoothing length in the expression for Re, which gave me the parameter values as:
   $\rho_w = 1.0$, $U_{max} = 22015.38$, $l=\frac{1.3}{20}$, $\mu_w = 0.2$, $\sigma = 121169.239$, $Re = 7155$, $We = 260$.
   however the droplet started disintegrating even before spreading on the bottom wall.
   Hence I also have this question: How is Re calculated in SPHinXsys?

[nilot-pal]

Hi all,
Since I am not an expert in Riemann solvers, I made an attempt to simulate the impact problem by incorporating into SPHinXsys the WCSPH with artificial viscosity and surface tension models from the paper. My new results are as follows:
2D: 01.01.2024_13.48.58_REC (screenrec.com)
3D: 02.01.2024_11.04.47_REC (screenrec.com)
While you analyze the results, I've some questions and concerns about the correct implementation of the models from the paper in SPHinXsys, which I can share with you on this thread or in a face to face meeting if you are available to meet anytime soon. Let me know whichever way you want.

[DrChiZhang]

Hi all, Since I am not an expert in Riemann solvers, I made an attempt to simulate the impact problem by incorporating into SPHinXsys the WCSPH with artificial viscosity and surface tension models from the paper. My new results are as follows: 2D: 01.01.2024_13.48.58_REC (screenrec.com) 3D: 02.01.2024_11.04.47_REC (screenrec.com) While you analyze the results, I've some questions and concerns about the correct implementation of the models from the paper in SPHinXsys, which I can share with you on this thread or in a face to face meeting if you are available to meet anytime soon. Let me know whichever way you want.

Hi Nilot-pal, thx for sharing your results.
Is there any comparison between the results obtained with Riemann-based and Artificialviscosity-based WCSPH on this case?

[DrChiZhang]

1. 3D simulation of drop impact on dry surface using SPH method.

In the Ref, the time step seems VERY small, compared with that of SPHinXsys. You can try to change the time step size of SPHinxsys, to see if the simulation can be stable.

In SPHinXsys, the Re number is determined by the viscosity with predefined characteristic length and velocity.

[DrChiZhang]

1. 3D simulation of drop impact on dry surface using SPH method.

In the Ref, the time step seems VERY small, compared with that of SPHinXsys. You can try to change the time step size of SPHinxsys, to see if the simulation can be stable.

In SPHinXsys, the Re number is determined by the viscosity with predefined characteristic length and velocity.

FYI, in SPHinXsys, the U_max or U_ref should be 10 times of the characteristic velocity, following the weakly-compressible assumption.

[nilot-pal]

@ChiZhangatTUM , well thank you for the replies. So I did a quick check with the Riemann solver and am getting similar results for the 2D case (I am sure the same can be inferred for 3D).
result in 2D with Riemann solver: 09.01.2024_23.38.28_REC (screenrec.com)
I've my own naming convention when it comes to artificial viscosity and Riemann solver based codes, which I am sharing here:

Whether using artificial viscosity or Riemann solver based code, a comparison of results for 3D droplet impact among experiments, the Iowa State paper and SPHinXsys for a particular set of values of Re and We nos. is attached below:

My goal is to reproduce the plot from experiments and/or the Iowa State paper, but currently you can see the results do not match beyond a certain value of the non-dimensional time. There are more details which I am finding difficult to explain without a meeting.

[Xiangyu-Hu]

@ChiZhangatTUM , well thank you for the replies. So I did a quick check with the Riemann solver and am getting similar results for the 2D case (I am sure the same can be inferred for 3D). result in 2D with Riemann solver: 09.01.2024_23.38.28_REC (screenrec.com) I've my own naming convention when it comes to artificial viscosity and Riemann solver based codes, which I am sharing here: Whether using artificial viscosity or Riemann solver based code, a comparison of results for 3D droplet impact among experiments, the Iowa State paper and SPHinXsys for a particular set of values of Re and We nos. is attached below:

My goal is to reproduce the plot from experiments and/or the Iowa State paper, but currently you can see the results do not match beyond a certain value of the non-dimensional time. There are more details which I am finding difficult to explain without a meeting.

@nilot-pal Have you ever test the oscillation frequency of the droplet to make the surface tension parameter is correctly implemented?

[nilot-pal]

@Xiangyu-Hu, I understand your concern and hence want to clarify it. I think there are two different surface tension models used in the Iowa State paper (Yang, X., & Kong, S. C. (2018). 3D simulation of drop impact on dry surface using SPH method. International Journal of Computational Methods, 15(03), 1850011) . One of them (let's call it S1) is for the droplet-wall interaction and the other (let's call it S2) is for treatment of the air-water interphase. Whether the authors actually used S2 is not clear from their paper, as their problem setup doesn't show air or the properties table doesn't contain physical properties of air. So, the results that I shared are from the simulation of a water droplet (no surrounding air) impacting a solid wall (with S1 implemented). Now when I implemented S1 for the droplet oscillation problem (refer #378), I get the following plots for the X and Y coordinates of the center of mass of upper left quarter:

This isn't the expected behavior which renders the simple S1 model insufficient to capture the air-water interface.
Coming to model S2, it is based on the following paper:
Zhang, M. and Deng, X.-L. [2015] “A sharp interface method for SPH,” J. Comput. Phys.
302, 469–484.
Since you already have a robust surface tension model (color function gradient based) in your code, I will now attempt to simulate the droplet impact problem with:

1. S1 model for water-wall interaction
2. SPHinXsys surface tension model for air-water interaction
   Let me know if you have any concerns or suggestions.

[Xiangyu-Hu]

Reference in ne

@nilot-pal I am quite confused. Are you using the new surface tension model in SPHinXsys or you are still using the old one?

[Xiangyu-Hu]

@nilot-pal Also, in the reference paper you referred, they are not using an explicit surface tension at al. The surface tension is implicit given from the EOS, therefore, they obtain the surface tension after the simulation by measure the pressure drop.

[Xiangyu-Hu]

@nilot-pal Sorry that I missed that the new reference paper. It is using a explicit surface tension model.
It looks like it is a single phase simulation