MATHEMATICA MONTISNIGRI
Vol LVIII (2023)

CONTINUAL MODELING OF PROCESSES OF HOMOGETEROGENEOUS MELTING AND FRAGMENTATION OF METALBY ULTRASHORT LASER PULSE
V.I. MAZHUKIN ⃰, M.M. DEMIN, A.V. SHAPRANOV,
E.N. BYKOVSKAYA
Keldysh Institute of Applied Mathematics of RAS, Moscow, Russia*Corresponding author. E-mail: vim@modhef.ru

DOI: 10.20948/mathmontis-2023-58-6
Summary. The mechanism of homo-heterogeneous melting of aluminum with subsequentfragmentation of the melt under the influence of an ultrashort laser pulse is implemented
within the framework of a nonequilibrium continuum model with dynamically adapting
computational grids. The interaction of slow heterogeneous and fast homogeneous melting isinvestigated. Melt parameters and features of mechanical fragmentation of the liquid phase of
the metal due to negative pressure developing in the hydrodynamic unloading wave aredetermined.
1

INTRODUCTION

The study of dynamic fragmentation in shock-loaded metals and the assessment of thegeometric and kinematic properties of the obtained fragments is an urgent problem for both
fundamental [1] and applied science [2]. Among the dynamic fragmentation processes, thespallation fracture of solid materials has become one of the widely studied phenomena for along time [3], [4]. In this case, the spallation was defined as a rupture in the body due to themagnitude of the stress exceeding the strength of the substance. The main mechanism ofdestruction is a compression pulse propagating in a solid sample. This wave reflects from thefree surface and after interaction with an incident unloading wave creates tensile stresses thatcan lead to damage inside the material, ranging from small voids and cracks to complete
destruction and ejection of the spalled material.
Recently, there has been an increasing interest in such phenomena in the liquid phase,
which develop after partial or complete melting of the sample, for example, as a result of theaction of super-powerful ultrashort laser pulses on metals, which are additionally subjected tocompression or stretching. Experimental and theoretical [8]-[10] studies have made it possible
to determine a relatively narrow range of fluence values in which spallative ablation from theirradiated surface occurs due to unloading phenomena in the melt. Dynamic disintegration
leads to the formation of a cloud of liquid melt droplets released into the environment at ahigh speed. The release of melt droplets indicates the important role of melting mechanisms in
the formation of a liquid phase in a metal under the influence of ultrashort laser pulses.
Within the framework of classical thermodynamics, melting and solidification of solids isdefined as a phase transformation of the first kind occurring at a certain temperature Tm,
which corresponds to the equality of the free energies of the solid and liquid phases at thesame pressure. The equality of free energies indicates that the solid and liquid phases are in2020 Mathematics Subject Classification: 70-10, 76-10, 82C21.
Key words and Phrases: Laser ablation, Overheating, Homo-heterogeneous melting, Spallation, Fragmentation,
Nonequilibrium continuum model, Dynamic adaptation.

80

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskayaequilibrium. Rapid phase transformations of the first kind, including melting/solidification,
are associated with the emergence of metastable superheated/undercooled states. Recently, inmany theoretical and experimental works, increased attention has been paid to fundamental
problems in which a solid substance can be superheated and a liquid can be undercooled.
Two mechanisms of melting/crystallization of solids/liquids are known: heterogeneous(surface or frontal) and homogeneous (volumetric). In the first case, it was theoretically andexperimentally shown that the melting process of crystalline materials occurs in areas ofheterogeneous nucleation, such as free surfaces [11], grain boundaries [12], interfaces [13]and in areas containing structural defects (vacancies, dislocations) [14, 15] and impurities[16], which significantly reduce energy barriers to nucleation.
Heterogeneous melting is interpreted as the movement of the surface of a strong
discontinuity in a continuous medium. On this surface, the thermodynamic, thermal andoptical properties of the substance change drastically. In this case, the movement of the phasefront is associated with superheating of the surface of the solid phase.
Homogeneous mechanisms of phase transformations are characterized by the emergence ofa new phase in a certain volume of superheated/undercooled matter. In the scientific literature,despite the many theories put forward, there is still no consensus on the exact definition ofhomogeneous melting.
The presence of a large number of interrelated physical processes with limited information
about their dynamics makes it difficult to solve these problems through purely experimental
research. For this reason, pico-femto-laser ablation continues to be an area of active researchin which mathematical modeling plays an increasingly important role.
The fundamental aspects of the theory of nonequilibrium laser melting of metals andfragmentation of melts are based on mathematical models of two levels: macroscopic –including systems of continuous equations, including nonequilibrium hydrodynamics [17],[18] and microscopic - generalizing the atomistic systems of equations [19], [20].
This work is devoted to the mathematical modeling of the processes of heterohomogeneous melting with subsequent fragmentation of the melt underlying thenonequilibrium laser ablation of metal (Al). The study of the mechanisms of melting andfragmentation is carried out within the framework of a continuum approach using anonequilibrium hydrodynamic model [18], [21].
2 MATHEMATICAL MODEL
The main features of high power laser action on metals are associated with a high speed
and volumetric method of energy release of laser radiation in the electronic component, whichleads to a strong deviation of the system from the state of local thermodynamic equilibrium(LTE). Along with the violation of thermodynamic equilibrium, it is necessary to take into
account the kinetic disequilibrium of high-speed phase transitions in the irradiated substance[22], accompanied by the appearance of highly overheated metastable states [23,24]. Powerfulflows of matter through the interphase boundaries, together with the volumetric distribution ofenergy in the electronic component (Te >> Ti), contribute to the formation of near-surface
temperature maxima in the solid phase. The model must also take into account the powerful
dynamic effects associated with the rapid propagation of phase fronts. Accounting for theseprocesses requires an explicit description of the kinetics of phase transitions and theformulation of conservation laws at phase fronts [25], which are the hydrodynamic

81

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

discontinuities.
2.1 System of equations
The mathematical description of these processes is carried out within the framework of a
two-temperature, one-dimensional multi-front nonequilibrium hydrodynamic Stefan problem
supplemented by the laser radiation transfer equation. With the emergence of spallation, themathematical formulation is applied both for the bulk target and for each detached fragment.
The system of equations is written in a domain with three kinds of movable boundaries
Γsh,s(t), Γsℓ(t), Γkυ(t)The full mathematical formulation is given below in the following form:   (u )

0
 

x
 t

  (u )  (u 2 ) P


 0


x
x
 t



me
 e  z  ama  (  )  (  u )  W
G  
 e e  e e   e  g (Te )(Te  Ti )   x
x   x t  (  )  (  u )  u W
 i i  i i   P  i  g (Te )(Te  Ti )  x t x x G  Te G  0 x k  s , , t  0,(1)[s  x  sh,s (t )]  [sh,s (t )  x  s (t )]  s (t )  x   (t )

Te
T, Wi   i (Ti ) i , P, T   Pe , Te   Pi , Ti  .
xx
Here: ρ, u, ε, T, P - are the density, gas-dynamic velocity, internal energy, temperature andpressure of the media respectively, α(Te) - coefficient of volume absorption, R(Te) – surface
reflectivity, G - density of the laser radiation, Ce,i - heat capacity, λe,i - heat conductivity, g(Te)- electron-phonon coupling factor. Indexes s,ℓ,υ correspond to solid, liquid and vapor phases
respectively, e,i - correspond to electron and ion component.

where - We   e (Te , Ti )2.2 Initial and boundary conditionst  0 : u(0, x)  0,

p  0,   0 , Te  Ti  T0  293K

(2)

On the left stationary border, we write the condition of a zero mass and heat flux:x  s : sus  0, WT  0(3)Melting front: x = Γsℓ(t). On the moving interphase melting boundary x = Γsℓ(t) we useboundary conditions according to the non-equilibrium model of melting that is a nonequilibrium kinetic analogue to the equilibrium Stephan problem, formulated for a state82

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

sufficiently far from the local thermodynamic equilibrium. The system of equations at theboundary consists of the conservation laws of mass, momentum and energy, written for astationary (laboratory) coordinate system:x  s t  : s (us  s )   (u  s )Ps  s (us  s ) 2  P   (u  s ) 2(4) Ti   Ti ne i   i  Lm  s (us  s ) x   x  s
 s   1  Tm 
LnePs  P m  Lm (Tm ( Ps ))  (CP  CPs )(Ts  Tm )  
2




  s

Here: us, uℓ, Ps, Pℓ, ρs, ρℓ, Cps,Cpℓ, are the hydrodynamic velocity, pressure, density, specificheat at constant pressure for solid and liquid phases, Lne
m is the nonequilibrium melting heat,1    is the volume thremal expansion coefficient of liquid, which value is  T  P
determined from the thermal equation of state P = Pℓ(ρℓ,Tℓ).
The law of conservation of energy is supplemented by the dependence of the equilibrium
melting temperature on pressure Tm(Ps)  

Tm ( Ps )  T0  kPs
and kinetic velocity υsℓ(∆Tsℓ) of the melting/crystallization front on the magnitude ofsuperheating/undercooling ∆Tsℓ = Tsℓ ─ Tm(Ps) [22, 23], where Tm(Ps) is the equilibrium
melting temperature, which depends on pressure, Lm is the equilibrium melting heat, υsℓ, Tsℓ is the melting front velocity and melting/crystallization temperature at the moving phaseboundary.
In the laboratory coordinate system, the temperature dependence of themelting/crystallization front velocity in the diffusion-constrained model is expressed by theequation [24]: A 
 Ts    Lm maVs (Ts )  us   (Ts )  us  D exp 
exp




  1

 kTs   Tm Ps    k BTs

(5)The equation contains 3 constants D, A, β, two of which are dimensional (D and A), one (β)
is dimensionless. Their values were determined from molecular-dynamic modeling [23], [24].
For Al target: D = 2.699∙102 [m/s], A = 2.251 10-2 [eV], β = 3.693.
The electronic component at the transition through the phase boundary is assumed to becontinuous with respect to the heat flux and electronic temperature Te: Te   Te  e   e ,
 x  s  x 

83

Te, s  Te, .

(6)

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

Evaporation front x = Γkυ(t): A surface evaporation model written in the Knudsen layer
approximation is used as boundary conditions at the moving interphase boundary Γkυ(t). The
model consists of 3 conservation laws and 3 additional parameters on the outside of theKnudsen layer (Tυ, ρυ, uυ). In general, two of these parameters (usually Tυ, ρυ) are determinedusing certain approximation ratios, [26], and the third, usually the Mach number
M = u /usound, is found from the equations of gas dynamics.
The three conservation laws on an evaporating surface, written in a fixed (laboratory)
coordinate system, have the form:
x = Γkυ(t):k (uk  k )   (u  k )Pk  k (uk  k ) 2  P   (u  k ) 2 ,

(7)

T 

ne
   a a   k (uk  k ) Lx  k
The written conservation laws are supplemented by the following conditions:   ( M )sat , T  T ( M )Tk
 L m  T 
Psat Tk   Pb exp   a  k  k BTk  Tb  sat  Psat (Tk )(8)ma, Tk  Tk  Tb ,k BTk

where ρsat, Psat are the saturated vapor density and pressure, Lυ is the equilibrium evaporation
heat, Pb, Tb are the pressure and temperature of boiling at normal conditions, αρ(M), αT(M) –are the Crout coefficients [27], k = s, ℓ - are the indexes pointing on the possibility ofevaporation from the surface of solid or liquid phases.
The nonequilibrium heat of evaporation Lne
 is determined by analogy with thedefinition of nonequilibrium heat of melting. In this case, the enthalpies of the phases aredecomposed into a Taylor series at the point (Psat(Tk),Tk), taking into account the temperature
jump at the front. As a result, Lne takes the following form:Lne  L (Tk )  C P (T  Tk ) k  
Pk  P   1  Tk  k Pk  Psat 2  kk

For the electronic component and laser radiation, the boundary conditions (the rightboundary of the condensed phase) are formulated as follows:

 e

Te
  Te4 ,x(9)G(k (t ))  (1  Rk (Te ))G0 exp( (t / ) 2 ) .

84

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

where σ is the Stefan-Boltzmann constant, Rkυ is the reflectivity of the target surface, τ is thehalf-width of the laser pulse.
Shock wave front: x = Γsh(t). On the moving front of the shock wave x =Γsh(t), threeconservation laws are written out, the so-called Rankine-Hugoniot relations [28]:x = Γsh(t):mjsh 1 (u1  sh )  0 (u0  sh ) ,ijsh P1  1 (u1  sh ) 2  P0  0 (u0  sh ) 2 ,P1 (u1   sh )emjsh
 W1  jsh a,1  
1
2


2

(10)


P0 (u0   sh )m
  W0  jsh  a,0 
0
2



2




m
i
e
where j sh
, j sh, j sh–are the mass, momentum and energy fluxes. The values with index 0
belong to the background side and with 1 to the shock wave side.

3. METHOD OF SOLUTION AND ANALYSYS OF MODELING REZULTS
3.1 Method and algorithm of solution
One of the main computational features of the model under consideration is the presence ofmoving discontinuous solutions such as shock waves, interphase fronts, contact and freeboundaries. The approximation of the initial differential model (1) – (10) was carried out by afamily of conservative finite difference schemes written on computational grids with dynamic
adaptation [29, 30].
The dynamic adaptation method is based on the idea of using an arbitrary non-stationary
coordinate system, the transition to which is carried out by a transformation of coordinates
using the sought solution. This allows the problem of constructing and adapting
computational grids to be formulated at the differential level [31]. As a result, in the extended
differential model, part of the differential equations describes physical processes, and theother describes the behavior of grid nodes. In this case, the controlled movement of the grid
nodes is carried out by the sought solution [32]. Dynamic adaptation allows calculations withan arbitrary number of discontinuous solutions, such as shock waves [32], propagating phase[33, 34] and temperature fronts [35], contact and free boundaries. It also allows automating
the problem of creating new regions [36], including in the form of detached fragments [18,21].
The construction of a computational grid with a total number of nodes equal to 200 was
carried out on a segment equal to the initial thickness (200 nm) of the target. In addition, eachof the detached fragments additionally contained an average of about 30 nodes.
3.2 Analysis of the modeling results
Heating of aluminum target with thickness 200 nm is initiated by a laser pulse withwavelength λL=0.8μm and with Gaussian temporal profile G =G0 exp(-(t/τ)2) with theduration of τ = 0.1 ps and a maximum intensity of G0 = 1.97×1012 W/cm2 , whichcorresponds to the fluence F = 0.35 J/cm2. The modeling starts from the negative time t0 = 4τ. The laser radiation propagates from right to left and, hitting the surface of the metaltarget, is partially absorbed and partially reflected. The absorbed energy is spent on heating85

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

the irradiated target, generating phase transformations and shock waves in the solid phase, andremoving material in the form of a collective release of matter from the irradiated target zone.
Hetero-homogenous melting.
As already noted, closely related mechanisms of heterogeneous and homogeneous meltingof metals, accompanied by fundamental problems of superheated/undercooled metastablestates, play a decisive role in mechanisms of disintegration and spallation.
It is known that the appearance of the heterogeneous melting front, as a rule, occurs on freesurfaces, grain boundaries, interfaces and in areas containing defects (vacancies, dislocations)and impurities, which significantly reduce energy barriers to nucleation. The fundamental
aspects of the heterogeneous mechanism of phase transformations to be studied and analyzed
include the nonequilibrium melting/crystallization temperature Tsℓ and the temperaturedependence of the propagation velocity of the heterogeneous front of the interphase boundary
υsℓ = υ(Tsℓ). The mathematical formulation of the process of appearance of the heterogeneous
melting front is made in the form of a system of equations (4), (5).
Homogeneous mechanisms of phase transformations are characterized by the emergence ofa new phase in a certain volume of superheated/undercooled matter. The fundamental
problem of homogeneous melting/solidification mechanisms in theoretical and experimental
studies is to determine the maximum degree to which a solid can be superheated and a liquid
undercooled.
In field experiments, achieving significant overheating in metals is hampered by thedominance of heterogeneous nucleation, especially in the case of low heating rates. The
solution to this problem consists of suppressing heterogeneous nucleation of melting while
ensuring a high rate of volumetric heating. The fulfillment of these conditions contributes tothe emergence of highly overheated metastable states in the initial system, which arenecessary for the development of only homogeneous nucleation inside the crystal volume. If
extremely high heating rates are not applied under external influence, then homogeneous
nucleation of melting occurs when the maximum overheating of the solid is reached. The
degree of maximum overheating of the crystal is determined before the start of calculations
using molecular dynamic modeling [37], followed by comparison with the results oftheoretical and experimental studies [38-41]. Typically for aluminum, the maximum
overheating of the solid phase is in the range Tsmax=(1.2  1.5)Tm, and can increase at higher
heating rates. The calculations used the value Tsmax = 1.4Tm.
The absorbed fraction of the energy of the fs laser pulse is released in the electronic
subsystem of the metal, causing its rapid heating to high temperatures Temax ≈ 1.9·104K withan almost cold lattice. The slow energy exchange between the metal components leads to alarge temperature difference between them, Te >> Ti (Fig.1).
Melting starts at the moment t ≈ 70 fs with the formation of a heterogeneous phase frontΓsℓ(t) at the irradiated target surface. This front propagates with the velocity υsℓ = υ(Tsℓ) into
the bulk target. Significant overheating of the melting surface and large spatial temperature
gradients ensure a high rate of propagation of the melting front υsℓ ~ 1-2 km/s. The powerfulflow of matter through the interface Γsℓ(t) and the continued volumetric heating of the lattice
due to the energy of electrons lead to the formation of a temperature maximum in the solid
phase, for which the following inequalities are fulfilled: Tsmax (P) >>Tsℓ > Tm.(P) (Fig. 2).
When the temperature of the ultimate overheating in the solid-state phase of the metal isreached, the prerequisites for the emergence of a homogeneous melting mechanism arecreated in the metastable region of a highly overheated surface layer. In this paper, the86

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskayamechanism of homogeneous melting is modeled by introducing quasi-nuclei of the liquid
phase without explicit consideration of the internal structure of the material and thecorresponding kinetic processes.

Ti

Fig.1. Time dependence of electron and phonon temperature at the target surface.

Fig. 2. Spatial profiles of temperature and location of the first quasi-nucleus.

The quasi-nucleus is introduced at the point of maximum overheating with a temperatureof Tsmax = 1.4Tm. (the position of this point in Fig. 2 is indicated by a vertical arrow) and withthe size of Δх ≈ (0.1 – 1) nm. The kinetic conditions of heterogeneous melting (4), (5) are
written out at the boundaries of a quasi-nucleus, which ensures the appearance of two new87

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskayamelting fronts propagating in opposite directions with velocities of υsℓ,2 and υsℓ,3. As a result,
we obtain a system with three boundaries (phase fronts) moving at speeds of υsℓ,1 υsℓ,2 and υsℓ,3(Fig.3). The movement of all boundaries is implemented within the framework of the
dynamic adaptation method.
In this case, the resulting fragment of a solid melts on both sides (velocities υsℓ,1, υsℓ,2). The
second front with the velocity υsℓ,2 moves towards the melting front, which moves from theside of the irradiated surface with the velocity υsℓ,1. The combined effect of the two fronts
contributes to the rapid melting of the superheated fragment. The third front with the velocity
υsℓ,3 moving from the quasi-nucleus into the bulk target forms a new superheated region in therest of the target. When the ultimate superheating is reached in this region, a new quasinucleus if formed with two fronts of heterogeneous melting.

Fig. 3. Dynamics of the first three fronts Γsℓ,1 - Γsℓ,3, moving at speeds υsℓ,1 - υsℓ,3, respectively.

The process of generation of quasi-nuclei is multiplied as long as a state of strong
thermodynamic unequilibrium is maintained in the metal, for which the ratio Te >> Ti isfulfilled. When the system approaches the state of local thermodynamic equilibrium,characterized by Te ≥ Ti, the process of homogeneous melting passes into the stage ofheterogeneous melting, the phase front of which moves due to the energy stored in the liquid
phase.
The total number of quasi-nuclei for the considered regime is about 20. The rate ofhomogeneous melting is defined as the time derivative of the total spatial size of the set of
liquid phases. Figure 4 shows the melting rate. The sections of the curve indicated by the reddashed lines correspond to the rate of heterogeneous melting, and the blue solid line
corresponds to homogeneous melting. The maximum speed value is about 45 km/s, and theaverage value during the first 5 picoseconds is about 20 km/s. After equilibrtation, thetemperatures of Te ~ Ti, melting is again determined by a slow heterogeneous mechanism at a
speed of ~ 200-500 m/s. The rate of homogeneous melting is several times higher than the88

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

speed of sound and significantly exceeds the speed of the heterogeneous melting front. The
obtained values of the homogeneous melting rate are qualitatively consistent with the valuesof the homogeneous melting rate of the gold film ( ≈ 13.3 km/s) calculated using themolecular dynamics approach [42], which is also several times higher than the speed ofsound.

Fig. 4. Velocity of heterogeneous (red dashed line) and homogeneous (blue solid line) melting.

Fragmentation. Dynamic fragmentation consists of a set of thermodynamic and mechanicalprocesses that, under certain conditions, lead to the occurrence of stress states exceeding thestrength of the substance. Unloading of these states can lead to mechanical spallation, i.e. theformation of voids in the target substance. Given the complex dynamics of the formation ofvoids in the substance, the choice of a specific type of spallation criterion plays an importantrole in the calculations.
Spallation criterion. In this work, the spallation criterion was determined taking into accountthe high strain rates 1010 – 1012 s-1 typical to femtosecond laser action and was calculated out
using molecular dynamic modeling of the stretching of the sample to its complete rupture. For
the solid phase, a linear pressure dependence passing through two points was obtained as aspallation criterion
P(T = 273 K) = - 9.2 GPa, P(T = 900 K) = - 3.92 GPa(11)

In the melt of Al, the spallation threshold is written as:P   0.3883 / kT  / 109 GPa,where   1  T / Tcr 0.659(12)The evolution of hydrodynamic processes is closely related to melting. The high velocityof the melting front and rapid thermal expansion lead to the formation of a region of strong
compression of the substance near the melting front.

89

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

Fig. 5. Introduction of spallation at the moment t = +6.5 ps, Pi → 0 at ∆x.

Fig. 6. Spatial profile of density ρ in Al target and spalled fragments at the time t = +32 ps.

The propagation of the compression wave has a loading effect on the substance, whichcorresponds to a positive half-wave of pressure (Fig. 5). It is followed by unloading, whichcorresponds to a negative half-wave of pressure, Fig. 5. If the pressure in the rarefaction wave
reaches a sufficiently large negative value necessary to fulfill the criterion of mechanicalspallation (11), (12), then its fulfillment leads to a rupture in the target material, i.e. tospallation under the influence of negative pressure. At the moment of the spallation, thepressure in the negative half-wave tends to zero. At the spatial point where P ~ 0, two new90

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya

moving boundaries with saturated vapor between them and boundary conditions (7), (8)describing the kinetics of surface evaporation and condensation are introduced.
Figure 6 shows 15 fragments of the liquid phase ranging in size from 2 to 37 nm. About 30
mechanical spallations occurred in this mode of laser action; but some of them connected toeach other, forming larger fragments. The total amount of fragmented material in terms ofsolid density is 92 nm.
4

CONCLUSIONSA mathematical model of superheated Al melting under the influence of fs laser radiationis constructed. The kinetics of superheated metal melting is analyzed on the basis of anonequilibrium continuum model, which explicitly links homogeneous and heterogeneous
melting mechanisms.
The limits of overheating of solids and melting mechanisms at very high heating rates arefundamental problems with numerous applications, one of which is the dynamic
fragmentation of metals in a state of partial or complete melt, followed by the release of
spalled material.
The mechanism of homogeneous melting is modeled by the introduction of quasi-nuclei ofthe liquid phase in the region of maximum overheating of the solid phase.
The calculations have shown that the average speed of homogeneous melting ~ 20 km/s isseveral times higher than the speed of sound and significantly exceeds the speed of theheterogeneous melting front ~ (200 - 500) m/s.
The modeling made it possible to determine the thickness of the spalled layers, toinvestigate the regimes in which the spallation occurs in the melt and to determine the amountof the fragmented material.
Acknowledgements: The results were obtained using the equipment of Shared Resource
Center of KIAM RAS (http://ckp.kiam.ru).
REFERENCES[1] V.E. Fortov, L.V. Altshuler, R.F. Trunin, A.I. Funtikov, High-Pressure Shock Compression ofSolids VII: Shock Waves and Extreme States of Matter, Springer, (2004).
[2] G.I. Kanel, “Spall fracture: methodological aspects, mechanisms and governing factors”, Int. J.
Fracture, 163, 173–191 (2010). Doi: 10.1007/s10704-009-9438-0[3] L. Davison, D.E. Grady, Shahinpoor M. High Pressure Shock Compression of Solids II: Dynamic
Fracture and Fragmentation, Springer, (1996).
[4] Y. Cao, Y.C. Shin, “Shock Wave Propagation and Spallation Study in Laser Shock Peening”, J.
Eng. Mater. Technol., 132(4), 041005 (1-8) (2010). Doi:10.1115/1.4002048.
[5] S.I. Tkachenko, K.V. Khishchenko, V.S. Vorob'ev, P. R. Levashov, I. V. Lomonosov, V. E.
Fortov, “Metastable States of Liquid Metal under Conditions of Electric Explosion”. High Temp.,
39, 674–687 (2001). Doi: 10.1023/A:1012324925983[6] G. I. Kanel, S.V. Razorenov, & V.E. Fortov, “Shock-wave compression and tension of solids atelevated temperatures: superheated crystal states, pre-melting, and anomalous growth of the yieldstrength”, J. Phys: Condensed Matter, 16(14), S1007–S1016 (2004). doi:10.1088/09538984/16/14/010.

91

V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya[7] P.N. Mayer, V.V. Pogorelko, D.S. Voronin, A.E. Mayer, “Spall Fracture of Solid and Molten
Copper: Molecular Dynamics, Mechanical Model and Strain Rate Dependence”, Metals, 12(11),1878 (2022). Doi: 10.3390/met12111878[8] S.I. Ashitkov, M. Agranat, G.I. Kanel, & V.E. Fortov, “Approaching the ultimate shear andtensile strength of aluminum in experiments with femtosecond pulse laser”, AIP Conf. Proc.,
1426, 1081 (2012). doi: 10.1063/1.3686466[9] V.I. Mazhukin, M.M. Demin, A.V. Shapranov. “High-speed laser ablation of metal with picoand subpicosecond pulses”, Appl. Surf. Sci., 302, 6–10 (2014). Doi: 10.1016/j.apsusc.2014.01.111[10] M. Xiang, S. Jiang, J. Cui, Y. Xu, & J. Chen, “Coupling of dynamic ductile damage and melting
in shock-induced micro-spalling: Modeling and applications”, Int. J. Plast., 136, 102849 (2021).
doi:10.1016/j.ijplas.2020.102849.
[11] R.W. Cahn, “Materials science: Melting and the surface”, Nature, 323(6090), 668–669 (1986).
doi:10.1038/323668a0[12] J. Maddox, “Melting is merely skin-thick”, Nature, 330, 599 (1987). Doi: 10.1038/330599a0[13] J. Fu, Ch. Zhang, T. Liu, J. Liu, “Room temperature liquid metal: its melting point, dominating
mechanism and applications. Review article”, Front. Energy, 14(1): 81‒104 (2020). Doi:
10.1007/s11708-019-0653-8[14] D.R. Nelson, B. Halperin, “Dislocation-mediated melting in two dimensions”, Phys. Rev. B,
19(5), 2457 (1979). Doi: 10.1103/PhysRevB.19.2457[15] F. Delogu, “Molecular dynamics simulations of homogeneous and heterogeneous melting
scenarios in metals: Volume scaling and concentration of defects”, Phys. Rev. B, 73, 184108(2006). Doi: 10.1103/PhysRevB.73.184108[16] J.W. Christian, The Theory of Transformations in Metals and Alloys, 1st Edition, Pergamon(2002).
[17] A.V. Mazhukin, V.I. Mazhukin, M.M. Demin, “Modeling of femtosecond laser ablation of Al
film by laser pulses”, Appl. Surf. Sci., 257, 5443–5446 (2011). Doi:10.1016/j.apsusc.2010.11.154[18] V.I. Mazhukin, M.M. Demin, A.V. Shapranov, “High-speed laser ablation of metal with picoand subpicosecond pulses”, Appl. Surf. Sci., 302, 6–10 (2014). Doi: 10.1016/j.apsusc.2014.01.111[19] L.V. Zhigilei, Z. Lin, & D.S. Ivanov, “Atomistic modeling of short pulse laser ablation of metals:
connections between melting, spallation, and phase explosion”, J. Phys. Chem. C, 113(27),11892–11906 (2009). doi: 10.1021/jp902294m[20] P.N. Mayer, & A.E. Mayer, “Late stages of high rate tension of aluminum melt: Molecular
dynamic simulation”, J. Appl. Phys., 120(7), 075901 (2016). doi:10.1063/1.4959819[21] V.I. Mazhukin, M.M. Demin, A.V. Shapranov, A.V. Mazhukin, “Role of electron pressure in the
problem of femtosecond laser action on metals”, Appl. Surf. Sci., 530, 147227 (2020). Doi:
10.1016/j.apsusc.2020.147227[22] V.I. Mazhukin, A.V. Shapranov, V.E. Perezhigin, O.N. Koroleva, A.V. Mazhukin, “Kinetic
melting and crystallization stages of strongly superheated and supercooled metals”, Math. Models
Comput. Simul., 9(4), 448-456 (2017). Doi:10.1134/S2070048217040081[23] V.I. Mazhukin, A.V. Shapranov, O.N. Koroleva, “Atomistic modeling of crystal-melt interface
mobility of fcc (Al, Cu) and bcc (Fe) metals in strong superheating/undercooling states”, Math.
Montis., 48, 70-85 (2020). doi: 10.20948/mathmontis-2020-48-7[24] V.I. Mazhukin, A.V. Shapranov, O.N. Koroleva, A.V. Mazhukin, “Modification of the Wilson–
Frankel kinetic model and atomistic simulation of the rate of melting/crystallization of metals”,
Matematicheskoe modelirovanie, 35(11), 103–121 (2023). Doi: 10.20948/mm-2023-11-08[25] V.I. Mazhukin, A.V. Shapranov, A.V. Mazhukin, O.N. Koroleva, “Mathematical formulation of akinetic version of Stefan problem for heterogeneous melting/crystallization of metals”, Math.
Montis., 36, 58 - 77 (2016). www.montis.pmf.ac.me/allissues/36/Mathematica-Montisnigri-365.pdf

92V.I. Mazhukin, M.M. Demin, A.V. Shapranov and E.N. Bykovskaya[26] V.I. Mazhukin, A.A. Samokhin, “Boundary conditions for gas-dynamical modeling of
evaporation
processes”,
Math.
Montis.,24,
8-17(2012).
www.montis.pmf.ac.me/allissues/24/Mathematica-Montisnigri-24-7.pdf[27] D.Crout, “An application of kinetic theory to the problems of evaporation and sublimation ofmonatomic gases”, J. Math. Phys., 15, 1–54 (1936).
[28] L.D. Landau, E.M. Lifshitz, Fluid Mechanics. Vol. 6 (2nd ed). Butterworth-Heinemann, (1987).
[29] P.V. Breslavskiy, V.I. Mazhukin, “Computational algorythm of a hydrodynamical version of
Stefan problem by dynamic adapting grid”, Matem. Mod., 3(10), 104–115 (1991).
[30] A.V. Mazhukin, V.I. Mazhukin, “Dynamic adaptation for parabolic equations”, Comput. Math.
and Math. Phys., 47, 1833–1855 (2007). Doi: 10.1134/S0965542507110097[31] V.I. Mazhukin, M.M. Demin, A.V. Shapranov, I. Smurov, “The method of construction
dynamically adapting grids for problems of unstable laminar combustion”, Numer. Heat Tr. BFund., 44(4), 387 – 415 (2003). Doi:10.1080/10407780390219469[32] P.V. Breslavsky, V.I. Mazhukin, “Dynamically adapted grids for interacting discontinuous
solutions”. Comp. Math. and Math. Phys., 47(4), 687- 706 (2007). Doi:
10.1134/S0965542507040124[33] V.I. Mazhukin, I. Smurov, C. Dupuy, D. Jeandel, “Simulation of laser induced melting andevaporation processes in superconducting”, Numer. Heat Transf. A, 26, 587-600 (1994).
[34] V.I. Mazhukin, M.M. Chuiko, “Solution of multi-interface Stefan problem by the method ofdynamic adaptation”, CMAM, 2(3), 283-294 (2002). Doi: 10.2478/cmam-2002-0017[35] P.V. Breslavskii, V.I. Mazhukin, “Dynamic adaptation method in gasdynamic simulations withnonlinear heat conduction”, Comp. Math. Math. Phys., 48(11), 2102–2115 (2008). Doi:
10.1134/S0965542508110158[36] V.I. Mazhukin, M.M. Chuiko A.M. Lapanik, “Dynamic adaptation method for numerical solutionof axisymmetric Stefan problems”, Math. Model. Anal., 8(4), 303-314 (2003). Doi:
10.3846/13926292.2003.9637232[37] S.-N. Luo, T.J. Ahrens, T. Çağın, A. Strachan, W. A. Goddard, & D.C. Swift, “Maximum
superheating and undercooling: Systematics, molecular dynamics simulations, and dynamic
experiments”, Phys. Rev. B, 68(13) (2003). Doi: 10.1103/physrevb.68.134206[38] J.L. Tallon, “A hierarchy of catastrophes as a succession of stability limits for the crystalline
state”, Nature, 342, 658–660 (1989). Doi: 10.1038/342658a0[39] K. Lu, Y. Li, “Homogeneous nucleation catastrophe as a kinetic stability limit for superheatedcrystal”. Phys. Rev. Lett., 80(20), 4474-4477 (1998). Doi: 10.1103/physrevlett.80.4474[40] B. Rethfeld, K. Sokolowski-Tinten, D. von der Linde, & S.I. Anisimov, “Ultrafast thermal
melting of laser-excited solids by homogeneous nucleation”, Phys. Rev. B, 65(9), 092103 (2002).
Doi:10.1103/physrevb.65.092103[41] Q.S. Mei, & K. Lu, “Melting and superheating of crystalline solids: From bulk to nanocrystals”,
Prog. Mater. Sci., 52(8), 1175–1262 (2007). Doi:10.1016/j.pmatsci.2007.01.001[42] Y. Gan, & J.K. Chen, “Nonequilibrium phase change in gold films induced by ultrafast laser
heating”, Opt. Lett., 37(13), 2691 (2012). Doi: 10.1364/ol.37.002691

Received December 10, 2023

93