 Hello dear listeners, I'm presenting a work dedicated to, dedicated to, I'm presenting a work dedicated to massively parallel SPH meso modeling of shock-loaded spherical particles. I'll start with experimental data and mathematical modeling problems. Here is the experimental setup we are going to describe. Octogen plane wave generator on the bottom, then above with copper and dental plates. And dental plate, dental is tungsten alloy. And dental plate, there are some notches, some sheen in which golden spheres, very small golden spheres, are placed or the notch is empty or it's filled with groove plate. When octogen is detonated, it's expected that the shock wave goes through copper, dental, come to a layer of spheres and the spheres go up where their velocity is measured by PDV probes. So what do the experiment, experimentators obtain? They obtain the following plots of velocity from time dependence from PDV probes. And what do we see? We see that when a shock breaks out, there are some fragments of material with different velocities and different decay appear. That's the problem, so what is really measured in the experiments if it's not the initial spheres? So we need to solve the following problems. What physical processes take place in the layer of golden spheres during a shock wave propagation? Do the particles preserve their shape or if they do not, then what's the composition of ejected material? And does cumulation or fermentation take place? And to answer all these questions, we need to conduct direct numerical modeling on metascale as is. So we directly model golden spheres, dental plate, shock wave propagation in the layer, and etc. So how to do it? So modeling method choice and computational domain layout. There are some modeling complications because we have non-trivial geometry, many, many spheres. Many samples, contacts and shares flows are expected. If we choose traditional mesh methods, we obtain high accuracy, they will research. If we choose a layer view, we will have to use specific methods of contact and free surface tracking. If we use Lagrange point of view, then we will have to use mesh reconstruction. But we have an option to use Lagrange and meshless methods. So we do not need mesh, it means that we do not need to store connectivity between elements, it may change freely. In these methods, contact surfaces attract automatically. They are computationally hard, but are simply program implemented. Of all the variety of these meshless methods, we choose smooth particle hydrodynamics and its variant, which we call contact-smooth particle hydrodynamics. Contact because it contains even problem solution on contact between particles. This is a particle method with limited interaction length between particles, and that is why can be paralyzed effectively. On the slide you can see the figures that show what's the difference between mesh and particles. The mesh everything is connected, the particles are separated and they can freely mix. So that was the computational domain layout in real life on the left. And on the right you can see how we do model this. We took a lot of spheres, placed them on a thin layer of the plate, and hit all the system in the solid wall to obtain the needed shock wave. The data on velocity is got from the experimental work. Each particle, each SPH particle carries question of state idealistic plasticity model and liquid solid phase transition. There may be misunderstanding between golden particles and SPH particles. They are like this placed in a one single golden sphere. This is how they are traditionally figured like balls. So one golden particle, which flows out from the experimental step, consists of many SPH particles. There are several computational domain layout styles we used. So the first is, I call it plain pack. It contains three million of SPH particles, looks like this. Then SPH looks like this. Of course, simulation is in 3G. And also the new results we obtained recently is a modeling of irregular pack, as it is, of course, in the experiment. How do we do it? We randomly place the spheres. There are critical boundary conditions here, here and here. We place them randomly, and then we just pack them using our code, like direct modeling of ballistical packing. This simulation uses 8 and 2 million particles. There are some sides view of the sample pack. And it's clearly seen that there is so many SPH particles that all the modeling needs a parallel approach. So what do we do here? We use parallelization technique for distributed memory system. We call it Werner dynamic domain decomposition or VD3. So what's the idea? First, definition of the Werner diagram. If we have a set of points on some sample, 1.2, 0.3 and etc. Then each point makes its own diagram cell. Points in which, so every point in this area, are closer to this point than to any other point. So here points are all closer to this point, here all these points are close to this point. And all these cells comprise the Werner diagram. What's the idea of the composition? If this sample is not just a sample, it's our SPH sample filled with particles, and we put the Werner diagram on it, we say that each of these points correspond to one process. So one process, two processes, and one process. And we let these points move with SPH particles flow and load in sub-domains. How it looks like. So for example, process one has very high load. We just move the point generator out, and the boundaries move also, and some particles from area one move to adjacent areas. And that is why the load is reduced, and so the simulation is quite balanced. We show here the strong scalability test result for sample static problem of plate at rest, containing 52 million particles from 100 to 1000 processes. So it's clearly seen that scalability is quite good, and it is enough to effectively model the problem. Some remark, one process contains well up to one to hundred of thousands. Now the results of the mesoscale modeling. First, regular packing results, then spec. When a shockwave comes to layer of spheres, the accumulation effects appear here. So we see these eddy structures, typical for with my mescophan stability, and they continue to progress in time. When a shockwave comes to the surface of the pack, jets appear that freely expand into vacuum. We do not model air, only gold, tungsten, and vacuum. The overall figure looks like this, from here to here. So the structure, the initial structure is completely lost. And some jets continue to propagate into vacuum. Oh, here are the results for plane pack. And here you can see that the jets that do appear at the contact between plate and spheres, they appear and they have free space to flow into. So the jet head doesn't hit into the particles that are situated before it. It just continues to freely expand. And we see that in such a situation, huge jets appear, which heads are here, and this is the jet body. The next result is real pack, the whole simulation. I suppose that nothing can be understandable from this video, so let's have a closer view. When a shockwave comes to the surface between spheres and plate, also eddy structures are formed, but they are not symmetric, not so right. But the picture is the same as for regular dense pack. On the right-hand side you can see the pressure distribution picture. Behind the shockwave front, quite all the material is melted, except for some spots of solid centers of golden particles. When a shockwave comes to the surface of the layer, also jets appear from the upper layer of the pack. And as in the dense pack situation, they continue to freely propagate. So that's the side view. Quite all the material is melted except for some spots. And that's the more closer view of the initial film. So here is what is left from the transmitter or plate in our simulation. That's the lower layer of spheres deformed. This is the main mass. The main mass is situated here, and these are remains of the jets. So some plots. Mass per area and velocity dependence here is mass per area of a small slice as the jets propagate this X. So the main mass has a velocity from, well, 2 kilometers per second and more. So in the experiment we observed that after a shockwave breaks out of the layer, then the distribution of velocity is, well, is quite in a good agreement with what we can see here. Well, this is what we obtained. We've explained what happens in the layer. We developed a highly efficient SPH code for parallel mass scale modeling of materials in extremes, and this code allows us to reveal the underlining mechanisms determining complex flow evolution. We were the first who showed that, in the problem that we demonstrate here, contains ejecta. We have found that ejecta characteristics and with Myrmeshkov instabilities effect are developed, and they depend on initial layer packing configuration that we use in our modeling. Well, the last slide showed that the main mass velocity distribution of ejecta is in a good agreement with experimental data provided. And the last word. So the code that we use, it provides the information that is impossible to obtain from the experiment. That's all. Present the talk. So there are five minutes for discussion. Did you use an excluded volume potential in your simulation so that the spheres don't overlap? We used Vendland kernel as the approximating function. It is standard for SPH simulations. Yes or no? You asked about potential, yeah? Two spheres occupy the same place. Yes, of course, yeah. Does your model allow two spheres to be at the same place at the same time? No, the SPH approximation does not allow this. If the pressure becomes high, they come close. Could you tell me what a golden sphere is? You call it a sphere, golden. What's golden about the sphere? Well, yeah, actually it's a ball. And you call it, so all balls are going at the same time. Yeah. So do you ever do the simulations with pro-late shapes? What is pro-late shape? Well, it's like a football. So there's an aspect ratio of a different one. No, we've done it only for right shape. So it's not the same for speaker, yeah? Hi, I'm Rabbi Bonifu, but I'm not Dr. Vieng. Not Dr. Vieng. Just a correction. Which one is which one? Okay, hello everyone. My name is Bobak Rabbi Bonifu, as you already heard. I'm coming from the University of Lille, and I'm working on melting the event by thermal convection. Together with Dr. Calzaverini and Dr. Silvia Hirata. Melting and freezing convection is happening in larger scale in environment, in phenomena like lava lakes, magma chambers, and the one that we are interested in, in melting of ice in North Pole Arctic. Why we are interested in? Because melting of ice results in formation of ponds, and these ponds are darker than surrounding environment. They absorb more heat, and the rate of melting is increasing. The fact is that during last 20 years, we have several models, large scale models, for predicting the environment in Arctic, but they all failed. We consider that we assume that this problem is due to this formation of ponds, because the scale of these ponds are much smaller than the large scales for simulation of the Arctic. They cannot be resolved, so there should be some kind of tuning to have the effect of these ponds in the simulations. What is happening, in fact, is that the cycle of melting of ice starts during the springtime when the snow layer at the top starts to melt and forms a layer of water on the top of the ice. During the summertime, the rate of melting grows, and sometimes we have some channels that connects the water inside the ponds to the ocean, the open water of the ocean. During winter season, we have again the formation of ice. Sometimes these ponds completely solidify, sometimes we have only a cap of ice on top of it. The governing fact that these ponds are governed by cold water at the bottom due to being in touch with ice, warmer water at the surface because of warm air, warmer environment, and solar radiation, and the fact that density profile of water is inverse between 0 to 4 degrees. So these ponds are prone to have convection. We created a model to analyze this system. This model consists of ice. At the beginning we start to hit the system with temperature T0 greater than the temperature Tm, which is the melting temperature of solid. We fix the height of the system as Hmax. When we start the simulation, a layer of liquid appears at the bottom, and at some height the convection starts. The governing equation for the system is Navier-Soch's equation, of course together with temperature equation and the boundary conditions. For the liquid part of the melting, we have the part that is in touch with solid layer, and it's not regular. At the bottom it's flat and the lateral boundary condition, which respectively, at the top, we have a Stefan boundary condition, which is known together with 0 velocity. At the bottom constant temperature, and for the lateral wall, we have periodic boundary condition. It is known that if we take the Stefan boundary condition into temperature equation, we can bring it into single-domain formulation. Considering the length of the system, the total length of the system, as the length scale, we can dimensionalize our system, which will be dependent to three global parameters of a Stefan number, which is the specific heat to Latin heat times temperature difference, and it is defined the rate of the melting. Prandt number, which is already known, I guess everybody knows here, and Rayleigh, which is known also as free convection or natural convection, and it defines if it is below some critical value, the main temperature is transferred by conduction, and when it is higher than the critical value, the temperature is transferred mainly due to convection. For the materials that we know already, the Stefan problem for rocks in the form of lava or magma is in order 1 to 10, and Prandt is in order 10 to the power of 4 to 10 to the power of 8. For water, which we are interested in, Stefan number is in order of 10 to the power of minus 2, and Prandt in order of 10. Okay, there is another parameter, which is important in our computation, and it is the average height of the melting, which is defined where the interface, the solid-liquid interface is located. Respectively, we have another parameter, which is Rayleigh effective, and is defined based on this melting front height. In this way, Rayleigh effective at time 0 is 0, of course, and at time when the total solid is melted, Rayleigh will be Rayleigh max. In a system of melting, it is governed by heat. Either it's coming inside the system or going out of the system. The heat, if it is normalized, it is known as Nusselt number, Nusselt effective coming in, Nusselt effective in, and Nusselt effective out. Writing it in a form of equation, it is a difference of temperature at the bottom, I mean Nusselt in is the gradient of temperature at the bottom, and Nusselt out is the rate of changes of melting front divided by a Stefan number. Conservation of energy requires that what is coming in with the energy coming inside the system should be equal to the energy going out of the system, plus the temperature increasing in the liquid layer. We know that this temperature increase in the liquid layer is always greater than 0. It means that Nusselt out is always smaller than Nusselt in. For the conductive case, the solution is known. It is known as a Stefan problem. The height of interface is proportional to the square root of time, and the profile of temperature is in form of error function. Likewise, Nusselt effective in and out can be computed numerically, and as you can see, it is independent of time. The problem of melting has been addressed quite a lot. I mentioned some of them, which is related to our project. For instance, Peter and Bijan, they started performing research on melting, lateral melting. They are heating the system from the side. Also, they addressed the problem of heating numerically. They use the same setup. They are heating the system from the bottom, but in different contexts. This is for the simulation of magma and lava lakes. Experimentally, we have Davies who performed an experiment on the melting of ice, and he focused mainly on the shape of the interface and what is happening to the interface. We have Hill who performed on Gilles Seroll, which is quite different to ice for high railing, and we have Yen. Again, he performed experimental on melting of ice, but in a different condition. They started to melt the ice from the bottom. It's a different configuration to pounds. They focused on the interface and the shape of the interface, and they saw something they called as shark skin pattern in the interface. The question is that can we connect heat flux to the global parameters of the system? Is there any characteristic for the roughness or the length of the domes that are happening or roles that are appearing in the interface? Do they have some kind of feedback? Which one has effects on the other one? The interface on the flow or the flow on the interface? The solution for the conductive regime, the analytical solution for the height of the system, the height of the interface is unknown, but the fact that we have a relation between Nusselt and the average height of the system can be used to derive some scaling relations. Similar to Navier-Ressault's Rayleigh-Bernard system, one can assume that Nusselt is a power scaling of Rayleigh, Prandtl, and Stefan, likewise for the interface, average of interface height, depending on time, Prandtl and Stefan, and if we work the relation in this way, we find this equation. There are some observations in order from this equation. First of all, we saw that Nusselt is not dependent on time. It's constant during conductive case, so alpha is zero in that case. We see that the height of the interface is proportional to the square root of time, which we already know from the Stefan problem. In Malkus' scaling, similar to Rayleigh-Bernard, if the alpha is one-third, we see that the velocity of the interface is constant, and if we can see the alpha one-half in the ultimate regime similar to Rayleigh-Bernard system, the acceleration of the interface will be constant. For the numerical part, we created a code. This code is based on lattice Boltzmann. In 2D and 3D, it can be used for simulating Rayleigh-Bernard system or melting in different times, of course. The governing part for the melting is, at the beginning, the melting fraction is zero. We compute velocity and temperature based on the lattice Boltzmann method. From this temperature and velocity, we compute enthalpy. We derive linear interpolation of enthalpy for computing the liquid fraction, and this liquid fraction for next time can be computed by finite difference. We can compute these loops several times until we get enough accuracy, but it's known that even one iteration is enough to have a solution in a conductive case similar to accurate enough to analytical solution. However, in the solid part, we have to cancel the velocity, unwanted velocity. By the method of penalization force, we can remove unwanted velocity appearing due to the computation. Before going further, let me show you a movie to see what is happening qualitatively in the system. First, this system starts at 1.10. We heated up the system from the beginning at the beginning so that the interface is moving parallel. It's flat, and then we see the convection appearing in the system. This convection goes further and further. The cells are getting bigger and bigger. Finally, at some point, it starts to merge. For instance, at this point, you see that we have only four roles previously we had three, and here, no, we have merging of two cells. This way, we cannot appreciate what is happening in the system. So let's look at some snapshot of this configuration. At the beginning, as I told you, we have the regular steady convection appearing in the system. This is the first destabilization in the system. At some point, these also destabilize again, and the roles start to merge together, the roles merging completely, and now we have less roles here. Again, it goes further and further until we have only one role, and it's confined only by lateral boundary conditions. We try not to get to that point because the effect of lateral boundary condition also on the system has been observed. So we try to only have one parameter at a time for our simulations. Similar to Rayleigh-Benard system, we can look at the new cells coming in the system versus Rayleigh for convective melting. If we put Rayleigh-Benard system, which we computed with the same code, the parameters for convective melting and Rayleigh-Benard are the same. The only thing that is changing for Rayleigh-Benard system is the height of the system, similar to convective melting. You are seeing that they are in good agreement with each other, and at high Rayleigh regime, they are converging to the same point. Already I told you about Olvarava and the performance simulation, but they use different parameters. I put their computation on our computation only for the sake of comparison, and what we see here is that first, the onset of convection for convective melting is delayed compared to Rayleigh-Benard system. The heat flux is slightly larger than convective melting than Rayleigh-Benard system. At the ultimate regime, they are getting closer to each other, and even it is inconsistent with other computation, although they use different parameters. And only for observation, I put the slope of this line for Rayleigh-Benard system, which you see is less than one-third here. Likewise, we can look at Reynolds number, which is defined at the average velocity times the height of the interface divided by viscosity. You see the similar behavior for Rayleigh-Benard system and convective, and why we have this similarity? This is a question. If you consider the velocity of the interface as this, if it is much smaller than the average velocity of the liquid part, the system has enough flexibility to stabilize this and behaves like a Rayleigh-Benard system. From this, we get this inequality, and if we bring Nusselt out to the right-hand side and plot it, we see that this parameter is always greater than one, and at the beginning of convection is in order of 10. Also, we have 3D system. It is similar to... It was a movie, but it didn't play. Similar to 2D scenario, we have the system of convection, but the freedom here is higher. We put a 2D convective melting Rayleigh-Benard and 3D on top of each other to compare. As you can see, the heat flux in 3D is higher than 2D and even higher than Rayleigh-Benard 3D. For the sake of comparison with experimental, we chose Grossman-Lose theory, and we plotted on top of our plots, and you can see that the Rayleigh-Benard system is really good, it's comparable in a good condition, and the message is that Nusselt melting in 3D is higher than 2D and is higher than in 3D. Another thing that we can look in is the morphology of the interface. You see there is some kind of polygonal hexagonal pattern in the interface, especially at the onset of convection. It has been set for experimental that is kind of hexagonal. It's difficult to conclude this kind of hexagonal pattern. There are different kinds of polygonal that I can see here at least. For the co-elixification, we can use two different parameters. One is the height of the roles in the interface, which can be identified as a deviation, and the length of the roles as a correlation function. We put 2D and 3D here. I try to wrap it up as soon as possible because there is some more information here. At the beginning, it is zero because the interface is flat. It follows with a jump, and it says more or less constant. Remember that it's been scaled by the height of the interface, and then we have the second jump. It is where the roles start to merge together, and it more or less says in the same regime. It means that the aspect ratio of the cells stays in the same value. However, for a deviation of interface, for 3D, you see that the depth of the role is much higher than 2D. The message here is that in 3D, the path to single-dome shape cell is much smoother than 3D, and the cells are much deeper in 3D. Another parameter that can be analyzed is the effect of a cell phone. We did it only for 2D at the moment. You see that the shape, the behavior of the system based on different cell phone numbers is more or less the same, so the effect of a cell phone is negligible. It's not that huge. We computed that this is of order of 10 to the power of minus 2, to the power of minus 0.1. The other fact that we can see is the relation between nu-celled in and nu-celled out. It more or less stays in a constant value, and it's really close to the solution of the conductive case. This is the good news because it says that there is a relation between heat coming in and out of this system, and it can be used for different application in tuning the global system, global simulations. And for the summary, as I told you already, the Rayleigh-Bernard, for the large Rayleigh at the ultimate regime is comparable to melting system, and the fact that this kind of power law is possible to be used with the parameter alpha to the power of 0.25. Why they are comparable? We talked about it because the velocity of the average inside the fluid layer is much higher than the velocity of the interface. We saw this kind of behavior, the polygonal behavior, polygonal shape of the solid liquid interface. The domes in 3D are much deeper, and also it's going smoother to create a single dome shape. And there is a weak dependency between Nusselt and SFN parameter, SFN number. For the future work, we plan to analyze the system with cavity and see the effect of cavity, any form of cavity. At the moment, we put circular cavity and box in our simulation. The effect of bulk heating, which is also important is due to radiation and it is what is happening, in fact, in reality. And for the last part, we want to put shear velocity on the interface. It can be realized at having wind or some kind of current in the system and having some kind of force in convection, melting like open channel that is happening in Arctic, in melting of ice. Thank you very much. We already done it. So the solution for analytical part is available, not from the convective regime. In the conductive case, we can create a liquid layer at the beginning and start a simulation from there. And the main reason we did it, in fact, is that at the beginning of simulation, the velocity of interface is infinite, theoretically. And there is some kind of bursting in this simulation. By putting this kind of initial liquid, we cancel the effect of infinite velocity. Yes? Does it work now? We have to communicate with them and we put our result in them and they are willing to do that? No, at the moment I cannot say it is working perfectly or not because nobody has done it yet. It's melting faster than our prediction. It's melting from the top because the density profile of water is inverse. We can have the same set of melting from the bottom and have the same effect. So a lot of the heat transfer, the surface area is increased enormously, which would account for a faster rate. So if you model the water underneath as a multi-phase flow, a particulate with ice and water, then you might capture enough of the heat transfer mechanism to account for the melting. In fact, it is happening. Pardon? It is happening, in fact. First of all, we needed to understand what is happening in convection in the normal way. What you are talking about is what is happening in multiple...