Characterization of Cellular Damage Induced by the Bubble Bursting Phenomenon

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Characterization of Cellular Damage Induced by the Bubble Bursting Phenomenon Autor: Pablo Matute Conejero Tutor: D. Alfonso Miguel Gañán Calvo Cotutor: D. Jose María López-Herrera Sánchez Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Characterization of Cellular Damage Induced by the Bubble Bursting Phenomenon Autor: Pablo Matute Conejero Tutor: D. Alfonso Miguel Gañán Calvo Catedrático de Mecánica de Fluidos Cotutor: D. Jose María López-Herrera Sánchez Catedrático de Mecánica de Fluidos Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado: Characterization of Cellular Damage Induced by the Bubble Bursting Phenomenon Autor: Pablo Matute Conejero Tutor: D. Alfonso Miguel Gañán Calvo Cotutor: D. Jose María López-Herrera Sánchez El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos A mi madre y abuela, referentes fundamentales en mi vida, por su constante apoyo, generosidad y ejemplo. A mi padre, que me ve desde el cielo. A mi tía Ana, por su cercanía y apoyo constante. A quienes me han acompañado durante estos cuatro inolvidables años, compartiendo retos, risas, desesperación y muchos buenos momentos. A todos ellos, mi más sincero agradecimiento. Pablo Matute Conejero Grado en Ingeniería Aeroespacial Sevilla, 2025 I Resumen E l fenómeno conocido como ’Bubble Bursting’ ha demostrado jugar un papel importante en diversos escenarios. Entre todos ellos, uno particularmente interesante es su capacidad para provocar la muerte celular en biorreactores. El objetivo de este estudio es caracterizar correctamente la letalidad de una burbuja que explota en la superficie del agua, teniendo en cuenta una serie de mecanismos que jamás se han considerado, además de los que ya se contabilizan. Para ello, se han llevado a cabo una serie de simulaciones, variando el número de Ohnesorge, mediante Basilisk y Python. III 1 Introduction. The Bubble Bursting Problem. O ne of the most developed fluid mechanics fields in the past 50 years has been microfluidics. Driven by the latest discoveries of nanotechnology, this discipline has proven to be crucial for numerous technical and scientific applications. Countless examples of this can be provided, covering disciplines from bioengineering or medicine to manufacturing industry or the aerospace sector. In the present, microfluidics keeps expanding itself and will continue to provide significant advancements for humanity, opening the door to new opportunities in the research and development of technologies with a direct impact on everyday life. Among all phenomena related to it, the bubble bursting problem deserves special attention. Bubble bursting consists in the collapse of an air bubble, in this case next to the water surface, and its evolution in time. At first sight this may not appear to be a very significant problem; nevertheless, it has a direct impact on life on Earth, as it is responsible for the water cycle. [5] The cavity left by the bubble leads to an explosion that generates aerosols that transport into the atmosphere salts, organic matter, microorganisms and other compounds. These act as cloud condensation nuclei and promote their formation, being then essential [13]. Besides, the rupture of the bubble has a significant impact on life surrounding the cavity. Pressure gradients and other forces involved in the process can affect microorganisms. In marine ecosystems, it has proven to have consequences in the ecosystem’s dynamics and biogeochemical processes of the marine life. For all these reasons, the study of the bubble bursting phenomena does not only contribute to the understanding of microfluidics, but also it allows us to understand the atmospheric and biological processes involved. 1.1 The Bubble Bursting Problem In the following lines, we will provide a brief description of the problem as well as the tools and procedures used for its resolution. Whenever a bubble rises through a liquid and gets to the surface, a gas-liquid interface, the interplay between different physical forces leads to an explosion. This explosion is nothing more than a jet that reaches a very high velocity and emits droplets. In this study we will consider that the inner air bubble has identical properties as the exterior air and that the bubble attaches to the interface without modifying it, as shown in Figure 1.2. Although this is not very precise, as what really happens is that there is a deformation of the water surface, it offers a good approximation. This will be further analyzed in the subsequent chapters. 1 2Chapter 1. Introduction. The Bubble Bursting Problem. Figure 1.1 Bubble bursting process on a bubble attached to the underwater surface, evolution on the inside and outside is shown. Extracted from [15]. Figure 1.2 Bubble considered at this study. As mentioned, when the bubble adheres to the gas-liquid interface, a thin liquid layer is formed. Due to thermal fluctuations, the thin water film that separates the inner bubble air from the exterior breaks. The new interface will develop a rapid evolution of its geometry, searching for a new state of equilibrium. Figure 1.3 Evolution of the water-air interface. 1.1.1 Global Parameters Bond number This phenomenon is guided by vastly different forces that, however, are fully defined by two dimensionless numbers: Bond and Ohnesorge. Bond number quantifies the relation between gravitational forces and surface tension. Bo =ρgL2 σ(1.1) In our case we will focus on a Bo«1, meaning that capillary effects will practically determine the evolution of the system over gravity. Since gravity then is not very significant, we can assume that the initial geometry of the bubble is nearly a sphere – a circle in 2D. As a simplification, it will be considered completely spherical. 1.1 The Bubble Bursting Problem 3 ρgR2 σ≪1R2≪0.072 1000 ×9.81 R≪2.7mm If the bubble radius (it’s characteristic length) is much less than 2.7 mm, this should be a good approximation. Ohnesorge number Ohnesorge number measures the importance of viscous forces in front of a combination between inertial forces and surface tension. Oh =µ √ρσL(1.2) Viscous forces oversee damping the effects of inertia and surface tension in the dynamics of the system, thereby modifying the evolution of the phenomena. Furthermore, Oh will determine the bubble radius if every other variable is fixed. R=µ Oh ·√ρσ 2 Therefore, depending on Bond and Ohnesorge numbers, the system will produce a different outcome. 1.1.2 The importance of Oh & the jet phenomena As said before, the collapse of the bubble is followed by a jet, which is responsible for the droplets (Figure 1.4). Figure 1.4 Jet and first droplet emitted for Oh =0.00833 ; (R=200µm). . Numerical data show that around Oh = 0.03, the jet that follows the collapse reaches maximum velocity [ 6 ]. This is due to a balance between inertia, capillarity and viscous forces, in such a way that the energy available for jet formation is maximized. Another interesting Ohnesorge number is 0.052, after which no liquid spout will be emitted into the air [9]. 4Chapter 1. Introduction. The Bubble Bursting Problem. 1.2 Basilisk To run our different simulations, we will use Basilisk. Basilisk is an open-source program, designed for Linux, based on the C programming language. Its main purpose is to serve as an efficient advanced numerical method solver. To achieve this, Basilisk uses adaptive Cartesian grids, meaning it refines the grid where a higher accuracy is needed - in this study, fluid interfaces. We will further explain different methods and functions that Basilisk handles. 1.2.1 Basilisk’s Quadtree/Octree Adaptive Mesh Refinement (AMR) AMR is the core of Basilisk’s computational efficiency; it is the procedure that it uses to determine where to increase the grid refinement. Often, grids acquire a greater precision by raising the number of cells in the whole domain; however, AMR adjusts the computational grid resolution only in the regions where a finer detail is required, such as areas with steep gradients, turbulence or the interface between two fluids - this is specifically our case. Depending on whether a 2D or 3D scenario is run, Basilisk will use a quadtree or octree structure, respectively. 1. 2D Quadtree The domain starts as a single cell. If a higher refinement is needed, the cell is divided into four other cells. 2. 3D Octree The domain stats as a cubic cell. If a higher refinement is needed, the cube will be divided into another eight cubic cells. By defining an error threshold, the program will decide where to increase the grid level. This is constantly checked at every time step and if the error is below a second error threshold, meaning that such a level of refinement is no longer needed, the cells are merged back into a larger cell. To grant that each cell can compute its derivatives, halo cells (ghost cells) are introduced. These are nonphysical (artificial) cells that serve as communication between different level refinement grids, for interpolations or for boundary conditions. Explained very briefly, if a cell is at the border of the domain it lacks at least a neighbour, thus to secure boundary conditions ghost cells are used. For example, if we want to impose in a 2D domain that the velocity of the flow when it enters the domain is u=1.0 a halo cell is created and stores this value, allowing the first real cell to calculate its derivatives. In this sense, if a cell neighbors a greater or smaller level grid, the fine cell estimates the value by using interpolations with adjacent cells. Figure 1.5 Tree-Grid structure. Picture extracted from [19]. More information is available there. 1.2 Basilisk 5 1.2.2 VOF Method The goal of the VOF Method is to determine the interface between two fluids using a volume fraction function, “f”, which varies between 0 and 1. The whole domain of study is divided into many cells, in which the following advection equation is measured. ∂f ∂t+u·∇f=0(1.3) If “f” is equal to 0 or 1, it means that that cell is filled by one or another fluid; if it takes values in between, that means that the cell contains the interface. Two different techniques can be used then to reconstruct the interface: Standard VOF or PLIC-VOF, this last one used by Basilisk. PLIC-VOF (Piecewise Linear Interface Calculation) is an evolution of the Standard VOF, solving one of its main problems, the diffusion of the interface when advecting “f”. This means that the whole cell is considered to be the interface, which clearly is not. PLIC-VOF calculates the shape of the interface inside the cell following this procedure. 1. The vector n, normal to the interface is calculated. n=∇f |∇f|(1.4) 2. A scalar "d" is calculated such that the volume fraction of the reconstructed interface matches the given VOF fraction f inside the cell. Af luid =f Acell =f(∆x)2Af luid(d) = f(∆x)2(1.5) Being ∆xthe length of the cell. 3. Finally the vector which contains the shape of the interface, xis calculated. n·x=d(1.6) In practice, the bisection method is applied to equation Rn·x−f(∆x)2=0 and n·x=d is searched in an iterative process. Note that a local reference system has been used, setting as the origin the cells lower left corner. 1.2.3 Navier-Stokes centered Navier-Stokes solver navier-stokes/centered.h is a Basilisk function used for incompressible flows where density varies (multiphase fluids). Its main objective is to solve the equation 1.7 through the Helmholtz-Hodge projection method. ∂u ∂t+u·∇u=−∇p+ν∇2u,∇·u=0(1.7) This involves, first, predicting the velocity without taking into account the pressure. 6Chapter 1. Introduction. The Bubble Bursting Problem. ρu∗−un ∆t=−ρun·∇un+µ∇2un+f Following this, incompressibility is imposed ∇·un+1=0 , leading to a Poisson equation for pressure: ∇·1 ρ∇p=1 ∆t∇·u∗ Lastly, the velocity field is corrected to achieve the incompressibility condition. un+1=u∗−∆t ρ∇p This has proven to be computationally efficient for multi-phase flows, where VOF is needed. 2 Bubble modelling & simulations T his chapter will focus on the explanation of the simulations carried out, providing an analytical demonstration of the hypothesis made. 2.1 Bubble 2.1.1 Fluid properties First of all, two fluids are present in this set of simulations: water and air. The properties of both fluids correspond to their tabulated value at 20ºC. Table 2.1 Density and dynamic viscosity of water at different temperatures. Temperature (°C) Density (kg/m3) Dynamic viscosity (kg/(m·s)) 0 999.82 0.001792 5 1000.00 0.001520 10 999.77 0.001308 15 999.19 0.001139 20 998.29 0.001003 25 997.13 0.000891 30 995.71 0.000798 35 994.08 0.000720 Table 2.2 Density and dynamic viscosity of air at different temperatures. Temperature (°C) Density (kg/m3) Dynamic viscosity (kg/(m·s)) 0 1.292 1.71 ×10−5 5 1.269 1.73 ×10−5 10 1.247 1.76 ×10−5 15 1.225 1.80 ×10−5 20 1.204 1.82 ×10−5 25 1.184 1.85 ×10−5 30 1.164 1.88 ×10−5 35 1.146 1.91 ×10−5 7 8Chapter 2. Bubble modelling & simulations Another relevant parameter of the fluid is the surface tension, which depends not only on temperature, but also on the purity of the liquid, the atmospheric pressure and the nature of the adjacent liquid or gas. Table 2.3 shows surface tension for pure water, in contact with air and at atmospheric pressure depending on its temperature. Table 2.3 Surface tension of water at different temperatures . Temperature (°C) Surface tension (N/m) 0 0.07564 5 0.07494 10 0.07422 15 0.07349 20 0.07275 25 0.07200 30 0.07124 35 0.07047 2.1.2 Bubble geometrical model In this section, we will justify the initial geometry of the bubble. In the first instance, we shall ignore the deformation that occurs when the bubble reaches the surface. This can be proven by comparing floatability forces and surface tension forces. If we consider a generic bubble with radius R, we can estimate the floatability force as: Ffloat = (ρ1−ρ0)gV (2.1) Being the energy then: Efloat = (ρ1−ρ0)gV R(1−sinα), which will have to be the same as the resistance that the surface tension opposes: fσ=σcosα(2.2) Eσ=Z Zspherical cap fσ(α)dS =Z2π 0Zπ 2 α0 fσ(α)R2cosαdαdφ=Z2π 0Zπ 2 α0 σR2cos2αdαdφ(2.3) Guided by the scheme in Figure 2.2, we can estimate that: Ef loat =Eσ Eσ=σπR2π 2−α0−1 2sin(2α0) (ρ1−ρ0)g4 3πR4(1−sinα0) = σπR2π 2−α0−1 2sin(2α0) 2.1 Bubble 9 Figure 2.1 Deduction of the portion of the bubble that rises above the water. Deformation due to pressure has not been taken into account. 4(ρ1−ρ0)g 3σR2=1 1−sinα0 (π 2−α0−1 2sin(2α0)) sinα0=h R;sin(2α0) = 2sin(α0)q1−sin2(α0);sin(2α0) = 2h Rs1−h R2 4(ρ1−ρ0)g 3σR2=1 1−h R  π 2−arcsinh R−h Rs1−h R2  If we scan over a range of values of R the next graph is obtained (2.2). It is clear from the results that, under the assumption of completely no deformation, bubbles with a radius not greater than 1mm adjust well enough to this approximation. This bubble will rise 1% of the radius ( 0.5% of the diameter) above the water. While acknowledging that this constitutes a significant oversimplification, increasingly higher with the bubbles’ radii, other studies that precisely determine the full bubble shape, such as [ 18 ], have ultimately reached the same conclusion (Figure 2.3). Hence, taking into account the deformation induced by gravity (Bond number) and what was just exposed, a bubble with a radius not greater than 1mm can be considered a perfect circle - a sphere in 3D, like the one shown in Figure 2.4 - without introducing significant errors. 2.1.3 Initial conditions considered As it is shown in Figure 2.4 the bubble already starts opened at the top. As this happens in a very short period of time, compared to the explosion of the cavity, it should not be a source of significant errors. Then, driven by the internal pressure of the bubble, the gas is expelled very quickly and a series of vortices appears. This internal pressure has been modeled following the Young-Laplace equation. 16 Chapter 2. Bubble modelling & simulations can reach sizes up to 15µm and studies such as [ 2 ] measure the Young’s modulus of the whole cell at 0.7kPa. Therefore, our goal is to estimate lethal thresholds for the work of the pressure gradients and the viscous forces. •Lethal criterion for pressure gradients. 1. Mechanical Energy Stored in the Membrane The elastic energy per unit volume stored in a CHO cell under a tangential strain ε is given by the following formula: w=1 2Esε2[J/m2] (2.7) where Esis the Young’s modulus of the cell N/m2 2. Total Mechanical Power If this energy is deposited within a short exposure time ∆t , the total instantaneous power is: ˙ W=w ∆t=1 2Esε2·1 ∆t(2.8) 3. Numerical Estimation –Es=700 N/m2 –ε=0.6(60% strain) –∆t=10−6s Substituting into the expression: ˙ W=1 2·700 10−6·(0.6)2(2.9) =1.26 ×108(2.10) ≈108W/m3(2.11) •Lethal criterion for ∇·(τ′·u) The remaining term of the energy equation 2.6 can be divided into two others if expanded, one being the traditional EDR and the other that accounts for the work done by viscous forces. ∇·(τ′·u) = τ′:∇u+u·∇·τ′ A very interesting way of setting a lethal threshold for this case is relating the shear stress developed in the fluid with these two terms and measuring directly that shear stress, which is the physical variable that is directly involved with lysis. Nevertheless, we will use the traditional EDR ( =τ′:∇u )criterion for CHO cells (threshold =108 ) for the whole expression, due to the fact that it has proven to characterize correctly a lethal volume when no other mortal mechanisms are involved. 2.2 Measuring impact on surrounding life 17 Despite it may seem we are ignoring the contribution of the term u·∇·τ′ , if we make a simple deduction of a mortal threshold its result is again 108 , since the mortal mechanism is similar to the pressure gradient’s one, thus making it a reasonable criterion. In practice, u·∇·τ′ is usually overshadowed by τ′:∇u and −∇·(Pu) , hence playing a minor role (Figure 3.28. Code 2.3 Event that saves and calculate variables related to the impact on surrounding life. event save_simulation (t += 0.001; t <= TTF) { scalar EDR_Walls[], EDR_propio[]; vector grad_p[]; foreach() { grad_p.x[] = (p[1,0] - p[-1,0]) / (2. * Delta); grad_p.y[] = (p[0,1] - p[0,-1]) / (2. * Delta); double mu = f[] * (mu1) + (1.0 - f[]) * (mu2); double dvdz = (u.x[1,0] - u.x[-1,0]) / (2.0 * Delta); double dvdr = (u.x[0,1] - u.x[0,-1]) / (2.0 * Delta); double dudz = (u.y[1,0] - u.y[-1,0]) / (2.0 * Delta); double dudr = (u.y[0,1] - u.y[0,-1]) / (2.0 * Delta); double d2udz2 = (u.y[1,0] - 2.*u.y[] + u.y[-1,0]) / (sq(Delta)); double d2dvdr2 = (u.x[0,1] - 2.*u.x[] + u.x[0,-1]) / (sq(Delta)); double d2dvdz2 = (u.x[1,0] - 2.*u.x[] + u.x[-1,0]) / (sq(Delta)); double d2udrdz = (u.y[1,1] - u.y[-1,1] - u.y[1,-1] + u.y[-1,-1]) / (4.0 * sq(Delta)); double d2dvdrdz = (u.x[1,1] - u.x[-1,1] - u.x[1,-1] + u.x[-1,-1]) / (4.0 * sq(Delta)); double tau_zz = 2.0 * mu * dudz; double tau_rr = 2.0 * mu * dvdr; double tau_rz = mu * (dudz + dvdr); double tau_thetatheta = 2.0 * mu * u.y[] / y; double Pot_dis_vis = tau_rr * dudr + tau_zz * dvdz + tau_thetatheta * u.y[] / y + tau_rz * (dudz + dvdr); double Pot_grad_p = - (u.y[] * grad_p.y[] + u.x[] * grad_p.x[]); double Pot_F_vis = mu * ( u.y[] * (2.0*(d2dvdr2 + (1./y)*dvdr - u.y []/sq(y)) + d2dvdz2 + d2udrdz) + u.x[] * (d2dvdrdz + (1./y)*dudz + (1./y)*dvdr + d2udz2)); EDR_Walls[] = Pot_dis_vis; EDR_propio[] = Pot_dis_vis + (Pot_grad_p > 0. ? Pot_grad_p : 0.) + ( Pot_F_vis > 0. ? Pot_F_vis : 0.); } char fname[256]; sprintf(fname, "snapshot-%g.csv", t); 18 Chapter 2. Bubble modelling & simulations FILE *fp = fopen(fname, "w"); fprintf(fp, "x,y,u.x,u.y,f,EDR_propio,EDR_Walls\n"); foreach() fprintf(fp, "%g,%g,%g,%g,%g,%g,%g\n", x, y, u.x[], u.y[], f[], EDR_propio[], EDR_Walls[]); fclose(fp); } Code 2.3 shows the event in Basilisk in which the variables are calculated and saved. It is important to note, in light of the impact on surrounding life, that the work of pressure gradients and viscous forces has only been accounted for when they are positive, as otherwise they could produce values that would overshadow other lethal mechanisms that may occur within the fluid. That is, pressure gradients will only be accounted for when they are contrary to the movement of the fluid, and the work of viscous forces will only be accounted for when they provide energy to the cell. 3 Data Analysis I n this chapter we will go over the results of the simulations that took place. From now on the historically used formula for EDR =τ:∇U will be named EDR_Walls and our own developed formula will be called EDR_propio =−∇·(Pu)+∇·(τ′·u). 3.1 Introduction to data analysis During the process of investigation on how to compute the lethal zone, several findings have been made; usually through trial and error. In total, three different methods were undertaken to fulfill this task, each one with its own pros and cons. •Method 1: Grid Interpolation and Backtracking with mask. •Method 2: Precision Augmented Interpolation and Backtracking with no mask. •Method 3: Basilisk based Adaptative Grid and Forward Tracking with mask. The justification of developing three methods lies in the pursuit of higher accuracy and rigor, along with accordance to what experimental data shows. 3.1.1 Determination of an dimensionless lethal threshold It is important to know that, in every simulation conducted, the results have no dimensions, since the code has been non-dimensionalized based on the Ohnesorge number. Therefore, as it will have to be used in every method, it is crucial to determine the dimensionless threshold for each Ohnesorge number. Since every mortal boundary has been estimated at 108it is a very simple calculation: [EDRSI] = kg m·s3 EDRadim =EDRSI ·T3 ρ·R2;kg/(m·s3)·s3 kg/m3·m2=kg/m kg/m=1 T=rρR3 σ 19 20 Chapter 3. Data Analysis 3.2 Method 1: Grid Interpolation and Backtracking with mask The procedure followed to determine the lethal zone consists of two different phases. First of all, we have run a set of simulations for different Ohnesorge numbers, saving variables of interest. Then, the lethal zone will be determined by making a backward particle tracking based on the velocity field. Full code can be found in Appendix A. The methodology followed has been to evaluate from the final time step back to the initial one every coordinate saved in the simulation and search for those that had recorded a mortal EDR value. Whenever a ’particle’ records lethality, it is then tracked until t=0 is reached in the simulation. This is computationally very simple; however, what was previously an advantage became a hindrance: Basilisk’s AMR (the adaptive grid). As a consequence, we were forced to construct a regular grid and to interpolate the velocity field based on the data that was saved with the adaptive grid, losing some precision and adding a considerable computational cost. 3.2.1 Particle Tracking Table 3.1 Saved variables in each time step in Method 1. x y u.x u.y f EDR_propio EDR_Walls ... ... ... ... ... ... ... −4.84375 0.15625 1.6062 ×10−5−5.88018 ×10−71−1.369 ×10−13 −1.369 ×10−13 ... ... ... ... ... ... ... Note: in the table above, f is the variable from the VOF Method. f=1 means that the cell is fully submerged in water. For every time step computed, the particles that fulfill f(x,y,tn) = 1and EDRi(x,y,tn)>ε∗ are stored L(ti) = {(xj,yj)}N j=1. Then, every mortal point is backtracked based on the velocity field, in what effectively is a Lagrangian tracking. x(n−1) j=x(n) j−ux(x(n) j,y(n) j,tn)·∆t(3.1) y(n−1) j=y(n) j−uy(x(n) j,y(n) j,tn)·∆t(3.2) As said before, the determination of the velocity field has been a complex task, in which we had to find a compromise solution between accuracy and computational cost. The approach taken was to process the data with Python instead of Basilisk C and to use some of its libraries. Specifically, to reconstruct the velocity field, we have used: scipy.interpolate.griddata((x_data, y_data), u_data, (x_j, y_j)) Given a discrete field, this function implements a Delaunay triangulation based on the nearest points available, the ones that were saved in the simulation, thereby creating a regular grid. Due to the high precision needed in the surroundings of the bubble, a double stage grid has been developed, offering a standard resolution of 2048 ×2048 in the whole domain that doubles for x∈[−5,1],y∈[0,4] . This has proven to be crucial as the vortex generated in the gas phase during the collapse of the bubble can add significant errors in the Lagrangian backtracking of the positions. This becomes particularly relevant for particles near the interface, which are, as recorded in this study, very important, ultimately leading to some inaccuracies. 3.2 Method 1: Grid Interpolation and Backtracking with mask 21 Code 3.1 Function that generates the two-stage grid. def generar_malla_no_uniforme(res): nx_fino = int(res * 0.75) nx_grueso = res - nx_fino x1 = np.linspace(XMIN, XFOCO, nx_fino, endpoint=False) x2 = np.linspace(XFOCO, XMAX, nx_grueso) x_grid = np.concatenate([x1, x2]) ny_fino = int(res * 0.7) ny_grueso = res - ny_fino y1 = np.linspace(YMIN, YFOCO, ny_fino, endpoint=False) y2 = np.linspace(YFOCO, YMAX, ny_grueso) y_grid = np.concatenate([y1, y2]) return x_grid, y_grid Figure 3.1 The two stages of the grid. 22 Chapter 3. Data Analysis Figure 3.2 Grid. During the backtracking, at each instant, the velocity field assigned to the particle corresponds to the one that the cell in which the particle is located has. These cells are nothing more than the adjacent space to each point of the two-stage grid. This is surely the reason for the inaccuracies found in the analysis of the results. Code 3.2 Functions that determine the particles’ velocity. def buscar_indices(x, y): ix = np.searchsorted(XGRID, x) - 1 iy = np.searchsorted(YGRID, y) - 1 ix = np.clip(ix, 0, len(XGRID) - 1) iy = np.clip(iy, 0, len(YGRID) - 1) return ix, iy def integrar_trayectoria(puntos, dumps): for archivo in reversed(dumps): df = pd.read_csv(archivo) u, v = cargar_velocidad(df) x = puntos[:, 0].numpy() y = puntos[:, 1].numpy() ix, iy = buscar_indices(x, y) puntos[:, 0] += u[iy, ix] * DT puntos[:, 1] += v[iy, ix] * DT dentro = (puntos[:, 0] >= XMIN) & (puntos[:, 0] <= XMAX) & \ (puntos[:, 1] >= YMIN) & (puntos[:, 1] <= YMAX) puntos = puntos[dentro] if puntos.shape[0] == 0: break return puntos 3.3 Method 2: Precision Augmented Interpolation and Backtracking with no mask 23 The function np.searchsorted is the one in charge of this. 3.2.2 Final Volume Once known the position (xj,yj) of every lethal point at t=0 , since we are working in an axisymmetric domain with Xas the axis of revolution, the following formula is used: V=∑ (i,j) lethal mask 2πyi,j·∆x·∆y However, in order to avoid counting the same region more than once, a procedure that we will refer to as ’mask’ has been used. This mask consists of calculating in which cell the particle finally is. Then, the coordinates and characteristic of the cell are used in the previous formula. If two or more particles end in the same cell it will only be computed once thus not counting the same region twice or more. 3.2.3 Pros and cons Although it achieves a surprisingly good reconstruction of the initial shape of the bubble, particles near the top water interface depict a great numerical imprecision (Figure 3.5). Furthermore, when backtracking, the same particle can be counted lethal more than once if it stays in the lethal zone for more than one saved instant, as new lethal particles are generated in each instant without considering those that have been generated in previous instants. Yet, this problem is solved with the final mask. That said, the values recorded are much less than those recorded by Walls et al in [20]. 3.3 Method 2: Precision Augmented Interpolation and Backtracking with no mask In this case, the procedure followed is exactly the same as in the previous method; however, a higher precision in the results, as well as a greater accord to experimental data, has been pursued. This has been achieved through an additional interpolation procedure and with some changes in the reconstruction of the mortal volume. 3.3.1 Precision Augmented Particle Tracking When it comes to the determination of the velocity of a specific particle, method 1, showed great inaccuracies. Therefore, a bilinear interpolation has been implemented in method 2, based on the regular grid, allowing us to achieve a higher level of accuracy at first instance. Code 3.3 Function that implements bilinear interpolation. def interpolar_bilineal(x, y, campo, xgrid, ygrid): ix, iy = buscar_indices(x, y) x1, x2 = xgrid[ix], xgrid[ix + 1] y1, y2 = ygrid[iy], ygrid[iy + 1] dx = (x - x1) / (x2 - x1) dy = (y - y1) / (y2 - y1) v00 = campo[iy, ix] v10 = campo[iy, ix + 1] v01 = campo[iy + 1, ix] 24 Chapter 3. Data Analysis v11 = campo[iy + 1, ix + 1] return (1 - dx) * (1 - dy) * v00 + dx * (1 - dy) * v10 + (1 - dx) * dy * v01 + dx * dy * v11 3.3.2 Final Volume Again, once known the position (xj,yj) of every lethal point at t=0 , the following formula is used 1: V=∑ (i,j) lethal 2πyi,j·∆x·∆y In this case, no mask has been used for the reason that previous data showed a very low lethal volume in comparison to [ 20 ] and experimental data. Thus, in the pursuit of achieving concordance with these results, we chose to eliminate it. This means that in this method, the volume is calculated considering every lethal particle. 3.3.3 Pros and cons of this method The ’precision augmented backtracking’ shows a high precision, understanding this as that the particles are one next to each other as it is expected; however, it does not achieve a greater accuracy, when it comes to the initial shape of the bubble, than Method 1. Nevertheless, the results of EDR_Walls are of the same order of magnitude as those recorded in [20]. For this method, counting the particles more than once in the backtracking is possible due to the elimination of the mask. That is for sure the reason for the higher levels of mortal volume. 3.4 Method 3: Basilisk based Adaptative Grid and Forward Tracking with mask Lastly, a third method was implemented to mitigate the issues encountered in the analysis of the results of Method 2. In this one, the same procedure that Walls et al uses to reconstruct the lethal zone is implemented; however, in this case, an adaptive positioning of particles through the fluid has been developed, following a similar approach to Basilisk’s AMR. The particles have been tracked forward, in contrast with previous methods, saving their dimensionless EDR values through the simulation. This guarantees that lethal zones (lethal particles that record mortal EDR values) are counted just once. Then, based on the spatial distribution of the particles, again using a mask, the mortal volume is computed. 3.4.1 Tracking Table 3.2 Saved variables in each time step in Method 1. x y EDR_propio EDR_Walls u.x u.y Lev f ... ... ... ... ... ... ... ... −4.843 0.156 −1.36 ×10−13 −1.36 ×10−13 1.60 ×10−5−5.88 ×10−77 1 ... ... ... ... ... ... ... ... 1∆xand ∆y are calculated through the dimensions of the grid with the higher resolution, as data in the previous method showed that all the lethal points will surely end up there. 3.4 Method 3: Basilisk based Adaptative Grid and Forward Tracking with mask 25 Note: A new set of variables was needed to implement this method, thus running new simulations. In particular, Lev, which indicates the level of the Basilisk grid in that place, was crucial to develop an adaptive positioning of the particles. Firstly, an initial grid was developed, assigning each cell of the Basilisk’s AMR at t=0 a particle. Code 3.4 Method 3 intial mesh. def generar_malla_inicial(path_csv, L0=10.0): df = pd.read_csv(path_csv) df = df[(df[’f’] == 1) & (df[’y’] < 4) & (df[’Lev’].between(5, 12))] x_list, y_list, h_list, tipo_list = [], [], [], [] for _, row in df.iterrows(): x = row[’x’] y = row[’y’] lev = int(row[’Lev’]) h = L0 / (2 ** lev) r = h / 4 # radio de dispersión local if lev >= 7: for angle in [0, 2*np.pi/3, 4*np.pi/3]: dx = r * np.cos(angle) dy = r * np.sin(angle) x_list.append(x + dx) y_list.append(y + dy) h_list.append(h) tipo_list.append(3) elif lev >= 5: for angle in [0, np.pi/3, 2*np.pi/3, np.pi, 4*np.pi/3, 5*np. pi/3]: dx = r * np.cos(angle) dy = r * np.sin(angle) x_list.append(x + dx) y_list.append(y + dy) h_list.append(h) tipo_list.append(6) x_part = np.array(x_list, dtype=np.float32) y_part = np.array(y_list, dtype=np.float32) h_part = np.array(h_list, dtype=np.float32) tipo = np.array(tipo_list, dtype=np.uint8) print(f" Malla refinada generada: {len(x_part)} partículas.") return x_part, y_part, h_part, tipo Note: This part of the code belongs to the latest version, meaning that, in this case, there is not just one particle in each cell, but three or six, depending on the level of the cell. A higher level 32 Chapter 3. Data Analysis Figure 3.11 Non-dimensional lethal volume in each case as a function of the Laplace number. Displayed in logaritmic scale. 3.5.2 Total volume Now, another set of graphs is presented, depicting this time the lethality in absolute terms. As a function of the Ohnesorge number Figure 3.12 Dimensional lethal volume as a function of the Ohnesorge number. Graphs obtained as a result of method 1, 2 and 3 are displayed from left to right. 3.5 Results 33 Figure 3.13 Dimensional lethal volume as a function of the Ohnesorge number. Graphs obtained as a result of method 1, 2 and 3 are displayed from left to right in logaritmic scale. The results in Figure 3.12 are remarkable. The balance between pressure gradients, viscous dissipation, work done by the viscous forces and the bubbles’ volume leads to a maximum in lethality, in absolute terms, at Oh ≈0.008 or R≈210µm for Method 1, Oh ≈0.01 or R≈150µm for Method 2 and a maximum at Oh ≈0.006 or R≈380µmfor Method 3. This has never been studied before, nor recorded. In fact, Donald E. Spiel states in [ 16 ] that “Microorganisms are known to be killed by bursting bubbles. The smaller the bubble, the more lethal. ... The centripetal acceleration of the opening cap of, say, a 1mm-diameter bubble is about 58,500g’s! The caps of smaller bubbles... are subjected to even larger accelerations.” As a function of the bubbles radii Figure 3.14 Dimensional lethal volume as a function of the bubbles radii. Graphs obtained as a result of method 1, 2 and 3 are displayed from left to right. Figure 3.15 Dimensional lethal volume as a function of the bubbles radii. Graphs obtained as a result of method 1, 2 and 3 are displayed from left to right in logaritmic scale. 34 Chapter 3. Data Analysis 3.5.3 Lethal zone morphology Method 1. Lethal zone at t=0 The reconstruction of the lethal zone has as an outcome the following images. These images do not appear to correspond with the results obtained in the graphs. This is mainly due to the way the lethal zone was plotted: an attempt was made to generate a closed area instead of plotting each point individually. For the following methods, only the points will be plotted. Figure 3.16 Lethal zone at t=0. Method 1 . 3.5 Results 35 Figure 3.17 Lethal zone at t=0. Method 1 . It can be noted in these images that no protrusion like recorded in [ 20 ] appears. This will be addressed in the next subsection. Method 2. Lethal zone at t=0 The Lagrangian reconstruction at t=0 in Method 2 has yielded the following images. Figure 3.18 Lethal zone at t=0. Method 2 . 36 Chapter 3. Data Analysis Figure 3.19 Lethal zone at t=0. Method 2 . Figure 3.20 Lethal zone at t=0. Method 2 . As it can be noticed, the results differ from those obtained in [ 20 ]. While Walls et all predicts a protrusion of the lethal zone beneath the bubble, that should correspond to those particles that will become part of the jet, our study records in Method 2 that the particles come from the lateral part of the bubble. As a demonstration of this last statement, we present in Figure 3.21 the Lagrangian tracking of 3.5 Results 37 different particles, adjacent to the bubble, forwardly in time. In addition, it is presented in Figure 3.22 a comparison between previous studies and the actual one. Figure 3.21 Forward-in-time Lagrangian tracking of several particles around the bubble reveals evidence of an upward extending tail that spreads along the bubble’s lateral region is observed . 38 Chapter 3. Data Analysis Figure 3.22 Comparison of the protrusion recorded in [ 20 ] (left) and the tail recorded in the actual one (right). This implies that the most exposed cells are those situated at the side of the bubble and not below. This makes certain cell cultures more likely to suffer greater damage than others, such as interfacial cell cultures. Method 3. Lethal zone at t=0 Conversely, method 3 provides a lethal zone in t=0 closer to what Walls et al describes. Figure 3.23 Lethal zone at t=0 based on EDR_propio. Method 3. 3.5 Results 39 Figure 3.24 Lethal zone comparison at t=0 based for EDR_propio and EDR_Walls. Method 3. Figure 3.25 Lethal zone at t=0 based on EDR_propio. Method 3. In this case, the tail recorded in method 2 does not appear, having instead a subtle protrusion. Temporal evolution of the lethal zone Furthermore, we have been able to define accurately how the lethal zone develops over time. This has been computed directly with Basilisk C, not needing, therefore, any interpolation. As it is visible in Figure 3.26 and 3.27, when our own EDR formula is used, the volumes where the fluid is mortal become bigger, which is what is expected. However, what is really interesting is that in Figure 3.27 the lethal area spreads through the jet, while in Figure 3.26 it does not. This is in accordance with what experimental data depict in [ 3 ], where it says that, after a set of experiments, most of the cells that died went through the jet. 40 Chapter 3. Data Analysis Figure 3.26 Temporal evolution of the lethal zone for Oh =0.00833 , calculated with the historical EDR formula "EDR Walls" in this case. Figure 3.27 Temporal evolution of the lethal zone for Oh =0.00833, calculated with our formula. In addition, it is shown in Figure 3.28 and 3.29 the spatial distribution of lethal mechanisms within the fluid when the jet is formed, which is the most critical instant for cell damage (Figure 3.4). 3.5 Results 41 Figure 3.28 Spatial distribution of the lethal mechanisms for Oh =0.00833. Figure 3.29 Spatial distribution of EDR-propio for Oh =0.00833. 48 Chapter A. Used Codes double *tau_max, int *plano) { double traza = tau_rr + tau_zz; double det = tau_rr * tau_zz - tau_rz * tau_rz; double discriminante = sq(0.5 * traza) - det; double sqrt_term = (discriminante > 0.0) ? sqrt(discriminante) : 0.0; double lambda1 = 0.5 * traza + sqrt_term; double lambda2 = 0.5 * traza - sqrt_term; double lambda3 = tau_tt; double d12 = fabs(lambda1 - lambda2); double d13 = fabs(lambda1 - lambda3); double d23 = fabs(lambda2 - lambda3); double maxdiff = d12; *plano = 0; if (d13 > maxdiff) { maxdiff = d13; *plano = 1; } if (d23 > maxdiff) { maxdiff = d23; *plano = 2; } *tau_max = 0.5 * maxdiff; } int main(int argc, char *argv[]) { if (getenv("OH")) { Oh = atof(getenv("OH")); } size (L0); DT = HUGE; LEVEL = 12; origin (-L0/2., 0.); init_grid (1 << 5); rho2 = 12e-4; mu2 = Oh/55; rho1 = 1.; mu1 = Oh; f.sigma = 1.; fprintf (stderr, " Oh = %g Level = %d rho2 = %g \n", mu1, LEVEL, rho2) ; run(); A.1 Basilisk 49 } double geometry(double x, double y) { double C1 = sq(x + R1 + 2*L2) + sq(y) - sq(R1); double C2 = sq(x + Rc) + sq(y - L3) - sq(Rc); double D1 = - x - 1e-8; double D2 = y - L3; double D3 = - x - (2*Rc); double D1D2 = min(D1, D2); double D1D2D3 = max(D1D2, D3); double D1D2D3C1 = min(D1D2D3, C1); double D1D2D3C1C2 = min(-D1D2D3C1, C2); return -D1D2D3C1C2; } event init (t = 0) { if (!restore (file = "dump")) { double eps = 0.05; refine ( sq(x + R1 + 2*L2) + sq(y) > sq(R1-eps) && sq(x + R1 + 2*L2) + sq(y) < sq(R1+eps) && level < LEVEL); refine (x < eps && x > -(2*Rc+eps) && level < LEVEL); fraction(f, geometry(x,y)); } } event init_pressure (t = 0) { double delta = 2.0 * Rc; foreach() { if (f[] == 0) { if (x < -delta) { p[] = 0.0; } else if (x >= -delta && x < 0) { double tanh_val = tanh(20*(x+delta/2)); p[]=-1-tanh_val; } else if (x > 0) { p[] = -2.0; } } else { p[] = -2.0; } } } event adapt (i++) { fprintf(stderr,"instante = %g dt = %g LEVEL = %d\n", t, dt, LEVEL); scalar f1[]; foreach() f1[] = f[]; 50 Chapter A. Used Codes adapt_wavelet ({f1}, (double[]){1.e-3}, minlevel = 7, maxlevel = LEVEL ); } event save_simulation (t += 0.001; t <= TTF) { scalar EDR_Walls[], EDR_propio[]; vector grad_p[]; foreach() { grad_p.x[] = (p[1,0] - p[-1,0]) / (2. * Delta); grad_p.y[] = (p[0,1] - p[0,-1]) / (2. * Delta); double mu = f[] * (mu1) + (1.0 - f[]) * (mu2); double dvdz = (u.x[1,0] - u.x[-1,0]) / (2.0 * Delta); double dvdr = (u.x[0,1] - u.x[0,-1]) / (2.0 * Delta); double dudz = (u.y[1,0] - u.y[-1,0]) / (2.0 * Delta); double dudr = (u.y[0,1] - u.y[0,-1]) / (2.0 * Delta); double d2udz2 = (u.y[1,0] - 2.*u.y[] + u.y[-1,0]) / (sq(Delta)); double d2dvdr2 = (u.x[0,1] - 2.*u.x[] + u.x[0,-1]) / (sq(Delta)); double d2dvdz2 = (u.x[1,0] - 2.*u.x[] + u.x[-1,0]) / (sq(Delta)); double d2udrdz = (u.y[1,1] - u.y[-1,1] - u.y[1,-1] + u.y[-1,-1]) / (4.0 * sq(Delta)); double d2dvdrdz = (u.x[1,1] - u.x[-1,1] - u.x[1,-1] + u.x[-1,-1]) / (4.0 * sq(Delta)); double tau_zz = 2.0 * mu * dudz; double tau_rr = 2.0 * mu * dvdr; double tau_rz = mu * (dudz + dvdr); double tau_thetatheta = 2.0 * mu * u.y[] / y; double Pot_dis_vis = tau_rr * dudr + tau_zz * dvdz + tau_thetatheta * u.y[] / y + tau_rz * (dudz + dvdr); double Pot_grad_p = - (u.y[] * grad_p.y[] + u.x[] * grad_p.x[]); double Pot_F_vis = mu * ( u.y[] * (2.0*(d2dvdr2 + (1./y)*dvdr - u.y []/sq(y)) + d2dvdz2 + d2udrdz) + u.x[] * (d2dvdrdz + (1./y)*dudz + (1./y)*dvdr + d2udz2)); EDR_Walls[] = Pot_dis_vis; EDR_propio[] = Pot_dis_vis + (Pot_grad_p > 0. ? Pot_grad_p : 0.) + ( Pot_F_vis > 0. ? Pot_F_vis : 0.); } char fname[256]; sprintf(fname, "snapshot-%g.csv", t); FILE *fp = fopen(fname, "w"); fprintf(fp, "x,y,u.x,u.y,f,EDR_propio,EDR_Walls\n"); foreach() A.2 Post-processing 51 fprintf(fp, "%g,%g,%g,%g,%g,%g,%g\n", x, y, u.x[], u.y[], f[], EDR_propio[], EDR_Walls[]); fclose(fp); } This code is an adaptation of the code provided by D. José María López-Herrera Sánchez. The initial pressure field and everything related to the lethality has been developed for this final degree project. A.2 Post-processing A.2.1 Lethal zone reconstruction at t=0. Method 1. Code A.2 Lethal Zone Reconstruction. Backtracking. Method 1. import pandas as pd import numpy as np import torch from glob import glob import matplotlib.pyplot as plt import os # ======================== # CONFIGURACIÓN # ======================== DEVICE = "cpu" DT = -0.001 RES = 2048 XMIN, XMAX = -5.0, 5.0 YMIN, YMAX = 0.0, 10.0 XFOCO, YFOCO = 1.0, 4.0 IMG_DIR = os.path.join("..", "imagenes") os.makedirs(IMG_DIR, exist_ok=True) # ======================== # UMBRAL ADAPTATIVO # ======================== def calcular_umbral_adimensional(oh): mu = 1e-3 # Pas rho = 998.29 # kg/ m sigma = 0.07275 # N/m EDR_fisica = 1e8 # W/m # CORRECTA fórmula: R = mu^2 / (Oh^2 * rho * sigma) R = (mu**2) / (oh**2 * rho * sigma) factor = np.sqrt(sigma**3 / (rho**3 * R**5)) return EDR_fisica / factor, R 52 Chapter A. Used Codes # ======================== # REJILLA NO UNIFORME # ======================== def generar_malla_no_uniforme(res): nx_fino = int(res * 0.75) nx_grueso = res - nx_fino x1 = np.linspace(XMIN, XFOCO, nx_fino, endpoint=False) x2 = np.linspace(XFOCO, XMAX, nx_grueso) x_grid = np.concatenate([x1, x2]) ny_fino = int(res * 0.7) ny_grueso = res - ny_fino y1 = np.linspace(YMIN, YFOCO, ny_fino, endpoint=False) y2 = np.linspace(YFOCO, YMAX, ny_grueso) y_grid = np.concatenate([y1, y2]) return x_grid, y_grid XGRID, YGRID = generar_malla_no_uniforme(RES) # ======================== # FUNCIONES # ======================== def extraer_puntos_letales(dumps, umbral): puntos = {"EDR_propio": [], "EDR_Walls": []} for archivo in dumps: df = pd.read_csv(archivo) if not {"x", "y", "f"}.issubset(df.columns): continue for var in puntos: if var in df.columns: mask = (df[var] > umbral) & (df["f"] == 1.0) pts = df.loc[mask, ["x", "y"]].values.tolist() puntos[var].extend(pts) for clave in puntos: puntos[clave] = torch.tensor(puntos[clave], dtype=torch.float32) return puntos def cargar_velocidad(df): ugrid = griddata_np(df, "u.x", XGRID, YGRID) vgrid = griddata_np(df, "u.y", XGRID, YGRID) return torch.tensor(ugrid), torch.tensor(vgrid) def griddata_np(df, columna, xgrid, ygrid): from scipy.interpolate import griddata puntos = df[["x", "y"]].values A.2 Post-processing 53 valores = df[columna].values X, Y = np.meshgrid(xgrid, ygrid, indexing="xy") interp = griddata(puntos, valores, (X, Y), method="linear", fill_value=0) return interp def buscar_indices(x, y): ix = np.searchsorted(XGRID, x) - 1 iy = np.searchsorted(YGRID, y) - 1 ix = np.clip(ix, 0, len(XGRID) - 1) iy = np.clip(iy, 0, len(YGRID) - 1) return ix, iy def integrar_trayectoria(puntos, dumps): for archivo in reversed(dumps): df = pd.read_csv(archivo) u, v = cargar_velocidad(df) x = puntos[:, 0].numpy() y = puntos[:, 1].numpy() ix, iy = buscar_indices(x, y) puntos[:, 0] += u[iy, ix] * DT puntos[:, 1] += v[iy, ix] * DT dentro = (puntos[:, 0] >= XMIN) & (puntos[:, 0] <= XMAX) & \ (puntos[:, 1] >= YMIN) & (puntos[:, 1] <= YMAX) puntos = puntos[dentro] if puntos.shape[0] == 0: break return puntos def reconstruir_mascara(puntos): mask = np.zeros((len(YGRID), len(XGRID)), dtype=bool) x = puntos[:, 0].numpy() y = puntos[:, 1].numpy() ix, iy = buscar_indices(x, y) mask[iy, ix] = True return mask def graficar(mask, nombre, etiqueta): X, Y = np.meshgrid(XGRID, YGRID) plt.figure(figsize=(6, 6)) plt.contourf(X, Y, mask, levels=[0.5, 1], colors=["#fc8d62"], alpha =0.8) plt.contour(X, Y, mask, levels=[0.5], colors="black", linewidths =0.8) plt.xlabel("x") plt.ylabel("y (radial)") plt.title(f"Zona t=0 ({etiqueta})") 54 Chapter A. Used Codes plt.axis("equal") plt.grid(True, linestyle="--", linewidth=0.5, alpha=0.5) plt.tight_layout() archivo = os.path.join(IMG_DIR, f"zona_t0_{nombre}_{etiqueta}.png") plt.savefig(archivo, dpi=300) plt.close() print(f" Imagen guardada: {archivo}") def volumen_revolucion(mask): dx = np.diff(XGRID).mean() dy = np.diff(YGRID).mean() Y, _ = np.meshgrid(YGRID, XGRID, indexing="ij") return np.sum(2 * np.pi * Y[mask] * dx * dy) # ======================== # MAIN # ======================== if __name__ == "__main__": folder = os.path.basename(os.getcwd()) try: oh_val = float(folder.replace("Oh_", "")) UMBRAL, R = calcular_umbral_adimensional(oh_val) print(f" Oh = {oh_val:.5f} R = {R:.6e} m") print(f" Umbral EDR adimensional usado: {UMBRAL:.2f}") except Exception as e: print(f" Error al interpretar Oh desde el nombre de la carpeta: {e}") UMBRAL = 92.59 dumps = sorted(glob("snapshot-*.csv")) puntos_dict = extraer_puntos_letales(dumps, UMBRAL) for clave, etiqueta in [("EDR_propio", "propio"), ("EDR_Walls", " walls")]: puntos = puntos_dict[clave] if puntos.shape[0] == 0: print(f" No hay puntos letales para {clave}") continue puntos_t0 = integrar_trayectoria(puntos, dumps) mask = reconstruir_mascara(puntos_t0) nombre = folder.lower() graficar(mask, nombre, etiqueta) volumen = volumen_revolucion(mask) print(f"[VOLUMEN-{etiqueta.upper()}] Volumen del sólido de revolución: {volumen:.6f} u") A.2.2 Lethal zone reconstruction at t=0. Method 2. A.2 Post-processing 55 Code A.3 Lethal Zone Reconstruction. Augmented Precision. Method 2. import time import pandas as pd import numpy as np import torch from glob import glob import matplotlib.pyplot as plt import os from scipy.interpolate import griddata import argparse # ======================== # ARGUMENTOS # ======================== parser = argparse.ArgumentParser() parser.add_argument("--input_dir", required=True, help="Carpeta de entrada con snapshots CSV") parser.add_argument("--output_images_dir", required=True, help="Carpeta de salida para imágenes") parser.add_argument("--output_volumes_dir", required=True, help="Carpeta de salida para volúmenes") parser.add_argument("--sufijo", required=True, help="Sufijo para identificar los archivos de salida") args = parser.parse_args() INPUT_DIR = args.input_dir IMG_DIR = args.output_images_dir VOL_DIR = args.output_volumes_dir SUFIJO = args.sufijo os.makedirs(IMG_DIR, exist_ok=True) os.makedirs(VOL_DIR, exist_ok=True) # ======================== # CONFIGURACIÓN # ======================== DEVICE = "cpu" DT = -0.001 RES = 2048 XMIN, XMAX = -5.0, 5.0 YMIN, YMAX = 0.0, 10.0 XFOCO, YFOCO = 1.0, 4.0 # ======================== # FUNCIONES AUXILIARES # ======================== def calcular_umbral_adimensional(oh): mu = 1e-3 # Pas 56 Chapter A. Used Codes rho = 998.29 # kg/ m sigma = 0.07275 # N/m EDR_fisica = 1e8 # W/m R = (mu**2) / (oh**2 * rho * sigma) factor = np.sqrt(sigma**3 / (rho**3 * R**5)) return EDR_fisica / factor, R def generar_malla_no_uniforme(res): nx_fino = int(res * 0.75) nx_grueso = res - nx_fino x1 = np.linspace(XMIN, XFOCO, nx_fino, endpoint=False) x2 = np.linspace(XFOCO, XMAX, nx_grueso) x_grid = np.concatenate([x1, x2]) ny_fino = int(res * 0.7) ny_grueso = res - ny_fino y1 = np.linspace(YMIN, YFOCO, ny_fino, endpoint=False) y2 = np.linspace(YFOCO, YMAX, ny_grueso) y_grid = np.concatenate([y1, y2]) return x_grid, y_grid XGRID, YGRID = generar_malla_no_uniforme(RES) def buscar_indices(x, y): ix = np.searchsorted(XGRID, x) - 1 iy = np.searchsorted(YGRID, y) - 1 ix = np.clip(ix, 0, len(XGRID) - 2) iy = np.clip(iy, 0, len(YGRID) - 2) return ix, iy def griddata_np(df, columna, xgrid, ygrid): puntos = df[["x", "y"]].values valores = df[columna].values X, Y = np.meshgrid(xgrid, ygrid, indexing="xy") interp = griddata(puntos, valores, (X, Y), method="linear", fill_value=0) return interp def extraer_puntos_letales(dumps, umbral): puntos = {"propio": [], "walls": []} total_archivos = len(dumps) print(f" Procesando {total_archivos} snapshots para puntos letales ...") for idx, archivo in enumerate(dumps): df = pd.read_csv(archivo) if not {"x", "y", "f"}.issubset(df.columns): A.2 Post-processing 57 continue if "EDR_propio" in df.columns: mask = (df["EDR_propio"] > umbral) & (df["f"] == 1.0) pts = df.loc[mask, ["x", "y"]].values.tolist() puntos["propio"].extend(pts) if "EDR_Walls" in df.columns: mask = (df["EDR_Walls"] > umbral) & (df["f"] == 1.0) pts = df.loc[mask, ["x", "y"]].values.tolist() puntos["walls"].extend(pts) if idx % (total_archivos // 10 + 1) == 0: porcentaje = (idx + 1) / total_archivos * 100 print(f" Progreso extracción: {porcentaje:.1f}%") for clave in puntos: puntos[clave] = torch.tensor(puntos[clave], dtype=torch.float32) return puntos def integrar_trayectorias(puntos_dict, dumps): resultados = {} total_dumps = len(dumps) print(f" Integrando trayectorias hacia t=0 para todas las zonas...") for zona, puntos in puntos_dict.items(): if puntos.shape[0] == 0: print(f" No se encontraron puntos para {zona}. Saltando integración.") resultados[zona] = torch.empty((0,2)) continue puntos = puntos.clone() for idx, archivo in enumerate(reversed(dumps)): df = pd.read_csv(archivo) u = torch.tensor(griddata_np(df, "u.x", XGRID, YGRID)) v = torch.tensor(griddata_np(df, "u.y", XGRID, YGRID)) x = puntos[:, 0].numpy() y = puntos[:, 1].numpy() ix, iy = buscar_indices(x, y) puntos[:, 0] += u[iy, ix] * DT puntos[:, 1] += v[iy, ix] * DT dentro = (puntos[:, 0] >= XMIN) & (puntos[:, 0] <= XMAX) & \ (puntos[:, 1] >= YMIN) & (puntos[:, 1] <= YMAX) puntos = puntos[dentro] 64 Chapter A. Used Codes letal_walls = max_walls > EDR_adim celda_ids = [] for x, y, hi in zip(x0, y0, h): x_r = round(x, 12) y_r = round(y, 12) h_r = round(hi, 12) celda_ids.append((x_r, y_r, h_r)) celda_ids = np.array(celda_ids) celdas = defaultdict(list) for i, key in enumerate(celda_ids): celdas[tuple(key)].append(i) Vp, Vw = 0.0, 0.0 centros_x_propio, centros_y_propio, h_propio = [], [], [] centros_x_walls, centros_y_walls, h_walls = [], [], [] for key, indices in celdas.items(): x_c, y_c, h_val = key A_celda = h_val**2 if np.any(letal_propio[indices]): Vp += 2 * np.pi * y_c * A_celda centros_x_propio.append(x_c) centros_y_propio.append(y_c) h_propio.append(h_val) if np.any(letal_walls[indices]): Vw += 2 * np.pi * y_c * A_celda centros_x_walls.append(x_c) centros_y_walls.append(y_c) h_walls.append(h_val) with open(directorio_salida / "volumen_letal.txt", "w") as f: f.write(f"Volumen letal por EDR_propio: {Vp:.6e} m^3\n") f.write(f"Volumen letal por EDR_Walls: {Vw:.6e} m^3\n") f.write(f"Total celdas evaluadas: {len(celdas)}\n") f.write(f"Celdas letales (propio): {len(centros_x_propio)}\n") f.write(f"Celdas letales (walls): {len(centros_x_walls)}\n") def graficar_celdas(x, y, h, nombre): sizes = (np.array(h) / max(h))**2 * 50 plt.figure(figsize=(10, 10)) plt.scatter(x, y, s=sizes, c=’red’, alpha=0.4) plt.xlabel("x") plt.ylabel("y") plt.title(nombre) plt.xlim(-5, 0) plt.ylim(0, 4) A.3 Other codes 65 plt.gca().set_aspect(’equal’, adjustable=’box’) plt.tight_layout() plt.savefig(directorio_salida / f"{nombre}.png", dpi=600) plt.close() if centros_x_propio: graficar_celdas(centros_x_propio, centros_y_propio, h_propio, " zona_letal_EDR_propio") else: print(f" No hay celdas letales por EDR_propio para { directorio_entrada.name}") if centros_x_walls: graficar_celdas(centros_x_walls, centros_y_walls, h_walls, " zona_letal_EDR_Walls") else: print(f" No hay celdas letales por EDR_Walls para { directorio_entrada.name}") print(f" Procesado {directorio_entrada.name} {directorio_salida. name}") # ======================= # Bucle principal # ======================= if __name__ == "__main__": raiz_entrada = Path("datos procesamiento") raiz_salida = Path("resultados") raiz_salida.mkdir(exist_ok=True) for carpeta in sorted(raiz_entrada.iterdir()): if carpeta.is_dir() and carpeta.name.startswith("Oh_"): match = re.match(r"Oh_(\d+\.\d+)", carpeta.name) if not match: print(f" Nombre de carpeta inválido: {carpeta.name}") continue oh = float(match.group(1)) carpeta_salida = raiz_salida / carpeta.name procesar_zona_letal(carpeta, carpeta_salida, oh) The codes above were used for the reconstruction of the lethal zone at t=0. Assistance from large language models (LLMs) was employed during the development of these specific codes, contributing significantly to reducing execution times. All outputs were closely reviewed and validated by the author. A.3 Other codes 66 Chapter A. Used Codes Code A.6 LAUNCHER: Lethal Zone Reconstruction. Grid and Particle Tracking. Method 3.. import subprocess from pathlib import Path csv_root = Path("csv") output_root = Path("datos procesamiento") output_root.mkdir(exist_ok=True) # Buscar carpetas tipo Oh_* carpetas_oh = sorted([d for d in csv_root.iterdir() if d.is_dir() and d. name.startswith("Oh_")]) for carpeta in carpetas_oh: nombre = carpeta.name carpeta_salida = output_root / nombre carpeta_salida.mkdir(parents=True, exist_ok=True) print(f"\ n Procesando {nombre}") resultado = subprocess.run([ "python", "MallaHistEDR3_modular.py", "--input", str(carpeta), "--output", str(carpeta_salida), "--dt", "0.001", "--nproc", "8" ]) if resultado.returncode != 0: print(f" Error procesando {nombre}") else: print(f" Completado {nombre}") Note: In case anyone is interested in running these codes, the only requirement is to create a folder and place them all there. Inside this main folder, there must be another named ’datos procesamiento’. This one should contain additional folders, each named after the Ohnesorge number corresponding to a simulation Oh_0.* . Each of these must have several .csv files with the relevant variables for each time step saved in the simulation. Appendix B Basilisk Installation Guide B.1 Introduction The following guide aims to facilitate the installation and use of Basilisk, focusing on the interface and basic commands more than on its programming. Basilisk is a free software program based on solving partial differential equations using adaptive Cartesian meshes, that is, solving in an environment, cell by cell, the equations that are given to it. This makes it very useful for simulations in fields such as fluid mechanics. Everything that is going to be explained is collected in the different entries of the page of Basilisk. B.2 System requirements •Operating System: Linux •RAM Memory: 4GB •Free disk space: 10GB (independent from that needed for simulations) B.3 Installation B.3.1 Download Basilisk has been designed for the Linux operating system and it is only possible to operate the software from it. Therefore, it is strictly necessary to have it installed on the PC; however, both Mac and Windows offer a simple way to run Linux without having to create DualBoots or use a different computer. This guide covers only the Windows installation process. B.3.2 Installation process for Windows Enabling the Linux subsystem for Windows First of all, we should enable WSL (Windows Subsystem for Linux). To do so, we must open the Control Panel and go into Uninstall a program, below Programs. 67 68 Chapter B. Basilisk Installation Guide Figure B.1 Control Panel. Then, inside Turn Windows features on or off, we must select the option where it says Linux Subsystem for Windows. Figure B.2 Turn Windows features on or off. B.3 Installation 69 Figure B.3 Linux Subsystem for Windows. Ubuntu installation Next thing to do is downloading Ubuntu from Microsoft Store. There are many options available for it, the best one is the one that says just Ubuntu, without anything else. Terminal initialization Once installed, we should proceed to execute the app 1 . As soon as it is opened, you will be asked for a user and a password. After the initial setup, what is shown in Figure B.4 should appear, possibly followed by additional messages. Figure B.4 Ubuntu Terminal. Basilisk installation For the installation of Basilisk, the instructions on the Installation section of the program website are going to be followed; therefore, it is highly recommended to read them carefully. This section will also address potential issues that may arise during the process and how to solve them 2 Below a set of codes will be shown that must be written in the Ubuntu terminal, one after another. First, it must be written: sudo apt install darcs make gawk darcs clone http://basilisk.fr/basilisk 1It is not advisable to run Ubuntu as an administrator if you are not familiar with it, as it may damage the computer 2 It is very likely that unexpected errors will occur during the installation. Artificial Intelligence such as ChatGPT or Grok can be very useful in resolving them -simply attaching a screenshot of the error, accompanied by a small description should work 70 Chapter B. Basilisk Installation Guide This way, the program will be installed, and to update it, it should be enough to type: cd basilisk darcs pull After entering the first sentence, it is possible that we may encounter the error displayed in Figure B.5. Figure B.5 Darcs error during the installation. To solve this we will install darcs using cabal, a program that facilitates the building, packaging and distribution of libraries. The following sentences must be entered: sudo add-apt-repository universe sudo apt update sudo apt install build-essential ghc cabal-install cabal update Then, we will add a line of code to one of the generated files. This will make the use of darcs easier, as it will not be necessary to specify, whenever it is called, that it is stored in the cabal folder. It should be written: nano ~/.bashrc Add the next line at the end of the file that will open: export PATH="$HOME/.cabal/bin:$PATH" Once written, press: Ctrl + O, then Enter and Ctrl + X. source ~/.bashrc Lastly, it must be written. B.3 Installation 71 sudo apt install zlib1g-dev sudo apt install g++ cabal install darcs Thus, darcs will be installed and the only command left to complete the installation is the following one. The update process is the same as previously shown, as if no error had occurred. darcs clone http://basilisk.fr/basilisk Useful libraries for Basilisk For the correct use and display of Basilisk, a number of programs are needed, which will vary depending on the needs of what is programmed in Basilisk. Here are the basic and most common: emacs,gnuplot,bview and ImageMagick. •Emacs. Emacs is the software from which Basilisk codes will be programmed. Its use will be explained later. To install it: sudo apt-get install emacs emacs & •Gnuplot. Gnuplot allows users to create graphs from data. sudo apt install gnuplot •bView. bView is another more sophisticated visualization program that allows multiple functionalities that will be discussed in detail later. •ImageMagick. Finally, ImageMagick is a .ppm files viewer very useful for post-processing. 3 sudo apt install imagemagick 3 In the Basilisk Installation Guide, from their own web, it is recommended to execute the next sentece, although with what is already installed is enough: sudo apt install gnuplot imagemagick ffmpeg graphviz valgrind gifsicle pstoedit 72 Chapter B. Basilisk Installation Guide B.4 Use of Basilisk As explained in the introduction, the next part of the guide will focus on basic commands, the Basilisk interface, and post-processing, but will also cover aspects of programming. To understand how Basilisk works and its different functions, it is very helpful to visit the program page. B.4.1 Linux use To navigate efficiently, it is essential to know how Linux basic commands work. Knowing how to navigate between directories, create them and delete them is necessary. The LinuxCommand website provides a good explanation of everything one needs to know. B.4.2 Emacs use. Programming and compiling. Emacs is the core of Basilisk, where we will program. To start it, it should be written: emacs & A window, shown in Figure B.6, should appear. Figure B.6 Emacs initial interface. First of all, by clicking in File in the left top corner and then in Visit New File a space, like one shown in Figure B.7, will pop up. There we can write our first program. The superior bar, where it says Name, is where the name of our program should be written, followed by a .c.Test.c, for example. Next, another tab will open (Figure B.8), where we will finally code. Every other aspect is very intuitive; however, note that before compiling a program, it must be saved by clicking in Save. To edit a file, instead of pressing Visit New File, click on Open File. B.4 Use of Basilisk 73 Figure B.7 Emacs interface 2. Figure B.8 Emacs interface 3. Once written a program, we must compile, as it is common in every language, in the Ubuntu terminal. If any mistake is made in the code, Ubuntu will show an error in the terminal. It is important to know which libraries are used in the program for a successful compilation. Generally speaking, a Basilisk code should be compiled the following way. qcc -O2 -g -Wall programa.c -o programa -lm B.4.3 BView2D Basilisk allows the user to use an easy way of visualizing variables saved during a simulation through bview2D. If Basilisk has been installed just as explained here, these steps are to be followed to open the interface. Bview2D is opened with the next sentence: List of Codes 2.1 Implemented code for setting initial pressures 11 2.2 Code for boundary conditions 12 2.3 Event that saves and calculate variables related to the impact on surrounding life 17 3.1 Function that generates the two-stage grid 21 3.2 Functions that determine the particles’ velocity 22 3.3 Function that implements bilinear interpolation 23 3.4 Method 3 intial mesh 25 3.5 Code fragment that shows how the velocity field is interpolated in Method 3 26 A.1 Bubble Bursting Code 47 A.2 Lethal Zone Reconstruction. Backtracking. Method 1 51 A.3 Lethal Zone Reconstruction. Augmented Precision. Method 2 55 A.4 Lethal Zone Reconstruction. Grid and Particle Tracking. Method 3. 60 A.5 Lethal Zone Reconstruction. Lethal Volumes, EDR and Images. Method 3. 63 A.6 LAUNCHER: Lethal Zone Reconstruction. Grid and Particle Tracking. Method 3. 66 81 Bibliography [1] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and Peter Walter, Molecular biology of the cell, 6 ed., Garland Science, 2014. [2] Emrah Celik, Midhat H. Abdulreda, Dony Maiguel, Jie Li, and Vincent T. Moy, Rearrangement of microtubule network under biochemical and mechanical stimulations, Methods 60 (2013), no. 2, 195–201. [3] J. J. Chalmers, Y. Liu, M. Aucoin, J. Schultz, M. Glazman, B. Fiechtner, and A. J. Sinskey, Quantification of damage to suspended insect cells as a result of bubble rupture, Biotechnology Progress 10 (1994), no. 1, 40–46. [4] Dennis E. Discher, Paul Janmey, and Yu li Wang, Tissue cells feel and respond to the stiffness of their substrate, Science 310 (2005), no. 5751, 1139–1143. [5] Alfonso M. Gañán-Calvo, The ocean fine spray, Proceedings of the National Academy of Sciences (PNAS) (2022), Manuscript compiled on February 3, 2022. [6] Alfonso M. Gañán-Calvo and José M. López-Herrera, On the physics of transient ejection from bubble bursting, Journal of Fluid Mechanics 929 (2021), A12. [7] Lionel Guillou, Avin Babataheri, Michael Saitakis, Armelle Bohineust, Stéphanie Dogniaux, Claire Hivroz, Abdul I. Barakat, and Julien Husson, T-lymphocyte passive deformation is controlled by unfolding of membrane surface reservoirs, Molecular Biology of the Cell 27 (2016), no. 22, 3574–3582, Epub 2016 Sep 7. [8] Wei-Shou Hu, Claudia Berdugo, and John J. Chalmers, “the potential of hydrodynamic damage to animal cells of industrial relevance: current understanding”, Biotechnology Advances 29 (2011), no. 4, 448–460. [9] Ji San Lee, Byung Mook Weon, Su Ji Park, Jung Ho Je, Kamel Fezzaa, and Wah-Keat Lee, Size limits the formation of liquid jets during bubble bursting, Nature Communications 2 (2011), 367. [10] Oliver McRae, Peter L. L. Walls, Venkatesh Natarajan, Chris Antoniou, and James C. Bird, Elucidating the effects of microbubble pinch-off dynamics on mammalian cell viability, Biotechnology and Bioengineering 121 (2024), no. 2, 406–415. [11] B. Neunstoecklin, M. Stettler, T. Solacroup, H. Broly, M. Morbidelli, and M. Soos, Determination of the maximum operating range of hydrodynamic stress in mammalian cell culture, Journal of Biotechnology 194 (2015), 100–109, Epub 2014 Dec 18. 83 84 Bibliography [12] Thomas S. Penfield, Jing Xiong, and Linhong Ma, Mechanical heterogeneity in the plasma membrane of living cells, Biophysical Journal 114 (2018), no. 8, 1909–1920. [13] Eugenio Rastelli, Cinzia Corinaldesi, Antonio Dell’Anno, Marco Lo Martire, Silvestro Greco, Maria Cristina Facchini, Matteo Rinaldi, Colin O’Dowd, Darius Ceburnis, and Roberto Danovaro, Transfer of labile organic matter and microbes from the ocean surface to the marine aerosol: an experimental approach, Scientific Reports 7(2017), 11475. [14] W. Rawicz, K. C. Olbrich, T. McIntosh, D. Needham, and E. Evans, Effect of chain length and unsaturation on elasticity of lipid bilayers, Biophysical Journal 79 (2000), no. 1, 328–339. [15] Daniel Shaw, Qi Li, Janine Nunes, and Luc Deike, Ocean emission of microplastic, PNAS Nexus 2(2023). [16] Donald E. Spiel, The droplets produced by individual bubbles bursting on a sea water surface, pp. 291–325, Springer Netherlands, Dordrecht, 1999. [17] Shinya Tanaka, Kei Takemura, Takuya Takamatsu, and Masaru Mochizuki, Effects of stretching speed on mechanical rupture of single cells, Biophysical Reviews and Letters 16 (2021), no. 4, 1–9. [18] Yoshiaki Toba, Drop production by bursting of air bubbles on the sea surface (ii): Theoretical study on the shape of floating bubbles, Journal of the Oceanographical Society of Japan 15 (1959), no. 3, 121–130, Received September 4, 1959. [19] Antoon van Hooft, The tree-grid structure in basilisk, n.d. [20] Peter L. L. Walls, Oliver McRae, Venkatesh Natarajan, Chris Johnson, Chris Antoniou, and James C. Bird, Quantifying the potential for bursting bubbles to damage suspended cells, Scientific Reports 7(2017), no. 15102. [21] Fang Yuan, Chen Yang, and Pei Zhong, Cell membrane deformation and bioeffects produced by tandem bubble-induced jetting flow, Proceedings of the National Academy of Sciences 112 (2015), no. 51, E7039–E7047, Epub 2015 Dec 9. [22] Martin Černý, Martin Gregor, and Vladimír Šubr, Mapping mechanical properties of cellular membranes via high-resolution techniques, Scientific Reports 12 (2022), no. 1, 1453.