Molecular dynamics simulation of interface atomic diffusion in ultrasonic metal welding

The formation of a reliable joint between a large number of aluminum strands for battery applications is crucial in automotive industry, especially in the technology of autonomous vehicles. Therefore, in this study, mechanical deformations and diffusion patterns of the mating interface in ultrasonic welding of aluminum were investigated using molecular dynamics simulations. Furthermore, microscopic observations of the joints between aluminum strands from ultrasonic welding illustrating the influence of two process parameters were done. To study the nanomechanics of the joint formation, two aluminum crystallites of different orientations were built. The impact of the sliding velocity and the compression rate of the upper crystal block on the diffusion pattern at the interface of the two crystallites were quantified via the diffusion coefficient. Tensile deformations of several joint configurations were performed to investigate the load-bearing capacity of the solid state bond, taking into account the compression rate, the sliding velocity and the crystallite orientation. The atomic scale simulations revealed that the orientations of the crystallites govern the interface diffusion and the tensile strength of the joint significantly. Furthermore, interface atom diffusion increased with increasing the sliding velocity. Additionally, it was observed that a higher sliding velocity enhances the friction heat generation between the crystallites and significantly increases the interface temperature.


Introduction
The formation of a reliable joint between multiple parts of similar or dissimilar metals is crucial for developing of diverse engineering products. The vital need for a reliable connection, in particular for the technology of autonomous vehicles, drives the automotive and aerospace industries towards using more reliable methods of metal joining. Ultrasonic metal welding (USMW) is a solid state joining technique, during which the growing global interface temperature remains mainly below the melting point of the participating metals [2,7]. Formation of a Shimaalsadat Mostafavi mostafavi@iam.rwth-aachen.de 1 stated that heat and diffusion are important factors in bond formation. Plastic deformation of the surface asperities and the increased real contact area can generate a bond [57]. Langenecker [22] showed that with the application of intense ultrasonic waves, the level of static shear stress of a metal reduces, which eases its plastic deformation. Zhang et al. [64] stated that severe plastic deformations and softening of the mating materials at the interface facilitate the grain boundary (GB) migration and the interfacial diffusion. Therefore, dynamic recrystallization occurs which leads to a high strength ultrasonic joint. Mohammed et al. [33] related the ultrafine structure at the ultrasonic spot weld interface of aluminum alloys to the dynamic recrystallization. They also reported an increase in the tensile shear strength with increasing the welding energy until a threshold was reached. Many current researches applied X-ray diffraction spectra of the joint fracture surface in order to evaluate the diffusion state and recrystallization mechanism of the interface layer. To this end, Liu et al. [25] identified recrystallization near the interface of Al to pure copper at 360 • C. Macwan and Chen [28] reported an increase in the thickness of the interface diffusion layer with increasing welding energy or temperature in ultrasonic spot welding of Cu to Mg. Similar results were obtained from the observations of Mirza et al. [31] for Al to Fe welds. Haddadi [17] related the rapid growth of the intermetallic layer to the high-cycle plastic deformations followed by the interface diffusion in ultrasonic spot welding.
Bakavos and Prangnell [1] introduced three zones in the weld of aluminum sheets after ultrasonic welding, namely the close-to-the-interface, the shear and the forged (strongly deformed) zones. They stated that the weld started with the evolution of microwelds at the contact surfaces with the highest pressure and material softening was crucial for forming a strong weld. Hazlett and Ambekar [18] stated that four bonding mechanisms might occur during ultrasonic welding depending on the mating materials, such as interface melting, mechanical interlocking, interface atomic forces and interface diffusion. Through EDS analysis of the dissimilar metals joined by ultrasonic welding, Zhu et al. [67] reported evidence of diffusion across the interface. They stated that ultrasonic welding results in atomic bonding between metals. Grain boundaries have substantial effects on the mechanical properties and deformation mechanism of the polycrystalline materials. Ward et al. [56] founded clear evidence of grain growth as a result of frictional heat generation at the weld interface of nanocrystalline alloys during ultrasonic welding by transmission electron microscopy (TEM). However, in situ monitoring of the welding process and observation of the metallic bonds require precise and costly measuring techniques. Therefore, in recent years computational modeling on different scales has been deployed to further understand the underlying mechanics of bond formation in diverse joining techniques. Shen et al. [49] developed a 3D FE model to describe the deformation behavior of the material subjected to longer welding cycles in a computationally effective way. Various theoretical and numerical models have been proposed, most of which are based on fundamentals of continuum mechanics [6,13,51]. Molecular dynamics (MD) simulations can serve as an effective computational tool for studying the sliding mechanism of different materials under diverse loading conditions. Furthermore, MD simulations are suitable tools for analyzing the deformation mechanism and the GB sliding at the nanoscale [54]. Yang and Persson [61] conducted MD simulations to derive the relation between the applied pressure and the contact area as well as the interfacial separation between solid surfaces. According to Chen et al. [6], large plastic deformations do not take place at the initial stages of the atomic sliding. The stored energy of the system releases with time and causes the deformations to grow rapidly. Through a MD model of self-mating Ni, they reported plastic deformations, temperature rise and the developing of a mixing layer at the mating interface. Dong et al. [14] introduced two different theories, which indicate that macroscopic friction is independent of sliding velocity. While one theory reports an increase in atomic friction with sliding velocity, the other argues that a smaller velocity allows more time for the formation of bonds and leads to larger friction. Using MD simulations, Long et al. [27] showed that during USMW, the microwelds were deformed or even broken due to the horizontal motion of the contact surfaces. Zhang and Liu [63] showed that effective diffusion occurred due to rather GB diffusion than lattice diffusion with high activation energy. They discussed the importance of GB diffusion in intermetallic compound (IMC) layer growth at Al-Ti interface. Due to Xu and Dávila [59], grain boundary migration plays an important role in determining the tensile deformation mechanism of nanocrystalline Al bulk material for a mean grain size of a few nanometers. They also presented the complete Hall-Petch relationship for nanocrystalline Al.
Chen et al. [6] studied microstructure evolution, heat dissipation and generation of mixing layers in the region near the contact interface during nanoscale friction using MD simulations. Jiao et al. [20] presented an atomic model to study linear friction welding between Ni and Al. They reported that the pore closure occurred due to atomic diffusion and interface deformation in Al and Ni parts, respectively. Li et al. [23] studied the mechanism of diffusion bonding of Al-Cu by means of MD simulations. They concluded that the nearest neighbor hopping mechanism was the main diffusion mechanism of Cu atoms in Al atoms. Furthermore, they stated that the diffusion of Cu atoms in Al lattice occurred deeper and at lower ratios than diffusion of the Al atoms in Cu lattice.
Panteli et al. [40] pointed out that the rapid development of a thick layer of IMC at the joint interface and the fast diffusion rates through the IMC layer could be possible reasons for the poor failure behavior during lap shear test of the dissimilar Al-Mg ultrasonic weld. Robson et al. [44] proposed a model to describe the formation and growth of the IMC during ultrasonic welding. These models predicted that welding time and layer thickness have an almost proportional relationship. Shah et al. [47] stated that the integration of ultrasonic vibrations into the process of resistance spot welding helps to remove the interfacial cracks in Al-Fe joints.
In order to be able to predict the deformation behavior of the polycrystalline materials at different loading conditions, some properties of the grain boundaries, including GB migrations, GB diffusion and dislocation nucleation, must be well understood at the atomic level. Xydou et al. [60] examined the effect of crystallographic orientations of the grains on the void closure at the grain boundaries for elevated temperatures using MD simulations. They realized the dominant effect of the GB diffusion at the interface voids in Cu with a diameter larger than 3.5 nm. Cheng et al. [10] showed through MD simulations that the starting atomic configuration around the symmetric tilt grain boundary in Al influenced the GB motion under shear deformation at different temperatures. Chen et al. [4] showed by MD simulations that most of the interfacial energy in Ni-Cr alloys was related to the difference in the structure from f.c.c to b.c.c across the interface. Furthermore, the atomic-scale stick and slip phenomenon happens due to the application of normal load and sliding. In this context, the existence of a mismatch in interfacial crystal orientation of small contact surfaces influences the friction force and the stick and slip process, as stated by Sorensen et al. [53]. To this end, the interface diffusion scheme and the ultimate tensile strength were discussed. Wei et al. [58] pointed out to the surface energy and density of the interface structure as the main reasons for differences in diffusion ability of various interface orientations between Fe and tungsten.
However, in most previous investigations, the differently oriented crystallites are located inside one polycrystalline metal and their crystallographic changes under loading and the resulting mechanical behavior are discussed. In this paper, the effect of the crystallographic orientation of two single crystals on the diffusion pattern of the interface is investigated. It is assumed that the two crystals belong to two separate polycrystalline metals and are subjected to mechanical contact by the exerted sliding and compression velocities. With the application of normal load to mating surfaces in ultrasonic welding, nearly full contact may occur at the interface. From this moment on, the rates of the displacements, normal and parallel to the surfaces, play an essential role, as they cause the generation of friction heat, plastic deformations and the diffusion of the interface atoms. This new study aims to add knowledge to the limited volume of work available related to the sensitivity of the interface diffusion and the weld strength to displacement rates and the crystal orientation during USMW in the literature. The results of microsection observations of real weld samples are presented in order to express the importance of the process parameters from the microscopic point of view. To this end, MD simulations were performed to investigate the nanomechanics of the interface. The actual effect of the two welding parameters, namely amplitude of the vibrations and the sonotrode pressure on the nanoscale joint formation has not yet been completely understood. Therefore, in this study, their effects at the nanoscale are translated into more intelligible features, namely the sliding velocity and the compression rate. To this end, the influence of the latter parameters on the growing interface temperature and the diffusion behavior of the atoms at the interface of two aluminum crystallites are discussed. In addition, the lattice orientation of the model is purposely changed to investigate the influence of the lattice structure on the interface diffusion under sliding and compression. This study shows that the initial orientation of aluminum crystallites significantly affects the tensile strength of the final joint. Evidence of local plastic deformations at the interface over the welding time is detected from the change in the statistical distribution of the atoms in the disordered structure. The significant influence of increasing the sliding velocity on the disorder of the atomic structure is illustrated through radial distribution function (RDF). This paper is divided into four sections and structured as follows: Section 1 is the literature review of the previous works related to the field of macro-, micro-and nanoscale contact and joining, Section 2 illustrates the microsections of the real samples, Section 3 describes the applied simulation method, Section 4 presents and discusses the results of the current study and Section 5 concludes this work.

Microscopic observations
Figures 1, 2 and 3 illustrate the microsections of the interfaces between multiple aluminum strands (EN AW-1370) of 0.43-mm diameter before and after ultrasonic welding. The ultrasonic joints were examined by backscattered electron (BSE) imaging at working distances (WD) ranging from 6 to 6.8 mm and 10 kV. Two of the main welding parameters, namely amplitude of the vibrations and sonotrode pressure, were changed and their influences on the joint formation have been studied. As shown in Fig. 3(c), high pressure of 1.5 bar results in excessive plastic deformations of the strands. Moreover, high pressure magnitudes prevent the  Figure 2(c) shows more local bonded sites for a higher amplitude of vibrations, where the boundaries of the mating parts are mostly merged into each other [34]. To better observe the joint interface, 500× magnification of some of the local bonded and non-bonded sites are further illustrated in Figs. 4, 5 and 6. Images of an obvious bonded spot and a non-boned interface line can be seen in Figs. 4(b) and 6(a), respectively. It can be observed that the shape and the average size of the aluminum grains in each strand are also affected during welding, when compared to the base metal, as shown in Fig. 1(b) and (c). The two mentioned features are essential indicators of structural characteristics for a polycrystalline material [8,19] that can change during the ultrasonic welding process. Furthermore, the atomic mismatch between two adjacent grains or crystals that have different crystallographic orientations next to a grain boundary is common and has mostly been considered as an interfacial defect [11]. Microsection observations showed that the process parameters impact the grain size, particularly at the joint interface. This can be caused due to the atomic diffusion of the mating parts owing to process parameters. Furthermore, the application of ultrasonic vibrations showed to refine the grain size and enhance the material mixing at the interface [5].

Simulation methods
To discuss the effect of the compression rate, the sliding velocity on the nanostructure and the mechanical properties of the interface of aluminum crystallites, MD simulations have been carried out. In MD simulation methods, atomic trajectories are obtained through integrating the classical Newtonian equations of motion for a system of N atoms as particles. Appropriate potential functions are used to compute the inter-atomic forces in these equations. The velocity and the positions of the atoms are updated using the Verlet time integration schemes [42]. In this study, the embedded atom (EAM) interatomic potential for aluminum developed by Mishin et al. [32] was used, which reads where i and j describe the respective atom indices, r ij is the scalar distance between atoms i and j , φ ij is the pairwise interaction potential, ρ i is the electron density of atom i due to all its neighbors and F i (ρ i ) denotes the energy arising from embedding atom i in an electron gas of density ρ i . In MD computational model, the kinetic and the potential energy of each atom is calculated at every time step. In order   to compute the temperature as a thermodynamic quantity, it must be averaged over a certain number of atoms. To this end, both kinetic energy terms for a system of N atoms from Eq. 2 are set equal, which results in the thermodynamic temperature of the system T as in Eq. 3: where k B is the Boltzmann constant and m i and v i are the mass and the velocity of the ith atom, respectively. After defining the atomic configuration, the simulation box was allowed to relax by a maximum volume change of 0.001Å 3 in one minimization iteration. To this end, the conjugate gradient method was used to obtain the local energy minimum of the structure at the beginning of the simulation. The Velocity-Verlet algorithm was used for the numerical integration of Newtonian equations of motion with a time step of 1 fs. In order to achieve an equilibrium state of the model, it was necessary to run the simulation for 50 ps before extracting averaged quantities such as temperature or displacements from the model. The equilibrium configuration of the system was achieved by coupling the system to an external thermal bath, hence conducting a canonical (isothermal-isochoric or NVT) ensemble at 300 K. After successfully equilibrating the system, two displacements over time, namely sliding velocity and compression rate, were assigned to rigid pulling layers in xand y-directions, respectively. This step lasted 96 ps at the longest depending on the magnitude of the compression rate until a predefined value of compression in y-direction was achieved. In the following step, compression was removed from the system and the rigid atoms were subjected to only sliding in x-direction for another 50 ps. Finally, the whole system was released from loading and the simulation continued for 50 ps, to simulate the end of the welding process.
With the application of sliding and compression, aluminum crystals were bound at the interface. Before tensile loading, each sample was relaxed by applying NVT ensemble at 0.1 K for the duration of 100 ps, allowing the configurations next to the interface to stabilize and relax. The relaxed samples were then subjected to uniaxial tensile loading with a constant engineering strain rate ofε = 0.001 ps −1 perpendicular to the interface plane in y-direction. To this end, the effect of the temperature on the deformation behavior of the samples was minimized in order to investigate the pure mechanical structure-property relationship. Hereby, the isothermal-isobaric (NPT) ensemble was considered. In doing so, the other two directions of the simulation box, namely x and z, were assigned to a zero pressure during the uniaxial tensile loading to ensure that the structure is stress-free perpendicular to the loading direction. The simulation was performed for another 1500 ps and the data for the stressstrain response until fracture was recorded.

Geometry and boundary conditions
Following microscopic observations as a motivation of this study, two initially separated crystal blocks that have different crystallographic orientations were considered, as shown in Fig. 7. The upper block, which was subjected to sliding and compression, collided with the lower block and diffusion occurred at the interface. Based on the orientation of the crystals during contact, the degree of the atomic diffusion was affected. Furthermore, the distribution of the atomic structure at the interface was influenced by the compression rate and the sliding velocity.
Molecular modeling and simulations were performed using LAMMPS (large-scale atomic molecular massively parallel simulator) developed by Sandia National Laboratories [41]. To display the atomistic arrangement, the molecular visualization program VMD was used. The simulation box consisted of two crystal blocks filled with a total number of about 22e3 aluminum atoms. Each crystal had the dimension of 101.25 ×72.90 ×24.3Å 3 and was filled with Al atoms with an atomic mass of 26.982 g/mol and lattice constant of 4.05Å, representing a f.c.c. (face-centered cubic) atomic arrangement. The thickness of the simulation box and the initial distance between the two blocks within the box were chosen much larger than the cutoff radius of 6.287Å for the Al atoms. In doing so, any interference of atoms with their periodic images or with atoms of the other block before the beginning of the process was prevented [12]. The crystallographic orientation of the lower crystal was set to <100> for all simulations, which includes a family of identical crystallographic directions, i.e.  Figure 7 shows the geometry and boundary conditions of the 0 • crystal model.
Three types of atomic descriptions were used in the computational model, namely rigid atoms, thermostat atoms and Newtonian atoms, i.e. atoms that follow the microcanonical ensemble. The region of rigid atoms moved as a single entity. Thermostat atoms were kept at a constant temperature of 300 K based on the Nosé-Hoover thermostat [30] and were used to ensure heat dissipation during sliding.

Calculation methods
At the interface of similar metals during ultrasonic welding, self-diffusion occurs, which describes the diffusion of atoms of one metal in another metal that are of the same type.
With the initiation of the contact between the crystal surfaces and during the compression and sliding processes, the interface temperature increases continuously. High interface temperature causes atoms to move in a direction, which will reduce atom concentration gradients. Due to Fick's first law, the flux of the diffused particle i, J i , is proportional to the negative gradient of the particle concentration C i and reads [24] where D i , the diffusion coefficient of component i, is independent of the magnitude of the concentration gradient [50]. The atomic coordinate r represents the direction parallel to the concentration gradient of the diffusing particle. In this context, the concentration gradient is the driving force.
To determine D, a vector is first calculated, which gives the total mean squared displacement (MSD) of the interface atoms, i.e. (dx 2 + dy 2 + dz 2 ), summed up and averaged over all the participating particles over time. Magnitudes of dx, dy and dz determine the atomic displacements in three orthogonal directions inside the simulation box. The slope of the MSD curve versus time is proportional to the diffusion coefficient of the diffusing atoms , representing the diffusion ability of the material. The diffusion coefficient is calculated as follows and is expressed in m 2 s −1 in SI unit: where t is the diffusion time and k is the number of the diffused particles. The magnitudes r(t) and r(0) represent the atomic coordinates at time t and at the initial time instant, respectively. The factor n equals 2 for a two-dimensional diffusion and 3 for a three-dimensional diffusion [50].
To study the changes in the crystal structure of aluminum during ultrasonic welding, the radial distribution function (RDF) (also called pair correlation function) g(r), was computed. This function defines the probability of finding particles of a particular type at a given distance around a center particle of a particular type [55]. Song et al. [52] calculated RDF from MD simulations to illustrate the disorder of the crystal structure during the linear friction welding of Ti and Ti-Al alloys. In solid structures, atoms fluctuate around their regular lattice sites (peaks in RDF) and the probability of finding atoms between these positions is almost zero (valleys in RDF). Aggregation or depletion of the diffusing particles at the interface will result in changes in the crystal structure and lead to disordered atomic structures. The average number of atoms that can be found at a given distance to a reference atom in all directions is described as where n(r) is the number of neighboring atoms within a distance r from a central atom, the magnitude ρ is the atomic density for the perfect crystal and g(r) is the radial distribution function [3]. In this study, g(r) function was determined with respect to the lower crystal and binned into a histogram.

Results and discussion
The formation of a joint in ultrasonic welding involves a large number of influencing factors and their interactions. Macroscale testing of the joint load-bearing capacity and microscopic observations of the joint from previous researches for similar [37] and dissimilar [16,26,65] material pairs point to the fact that process parameters are the crucial factors in the origin of the local bonds. Therefore, two factors, namely sliding velocity and compression rate, have been chosen and their influences on the joint formation and diffusion of Al atoms at the atomic scale have been studied. Furthermore, the influence of the crystal orientation of the aluminum crystallites on the atomic diffusion at the interface has been studied using MD simulations. After equilibrating the samples for 50 ps, the welding process starts and the rigid atoms of the upper block are subjected to constant velocities in x-and y-directions, which are referred to as sliding velocity u and compression rate v, respectively. For v = 0.2Å/ps, the process continues for 96 ps. Duration of the simulation is set depending on the value of v, such that the same magnitude of displacement in the y-direction is achieved. This condition ensures that only the effect of varying u is considered and all the samples are equally shrank in y-direction, reflecting the effect of a constant value of the welding pressure on the real samples. The simulation continues for another 50 ps, where the rigid atoms are subjected to only sliding. Compression is removed from the model during this time period. Finally, the sliding velocity u is removed and the simulation runs for another 50 ps to let the atoms find their final positions and let the structure stabilize. behavior of the atoms at the interface. In order to analyze the effect of the compression rate on the diffusion behavior of the interface, four values of v are considered, namely 0.2, 0.4, 0.6 and 0.8Å/ps, while u is set to a constant velocity of 0.5Å/ps. Temperature evolution of the interface depending on the sliding velocity for six different crystal angles is shown in Fig. 8. Once the upper block hits the lower block, a high repulsive force between atoms of the two crystals is induced and a small peak in the temperature history is observed after 100 ps. After 146 ps, compression is removed and only sliding motion remains. Interface temperature increases mainly during pure sliding, i.e. between 146 and 196 ps. During the last time period, i.e. from 196 to 246 ps, no significant rise in temperature is observed and small fluctuations at the end of the time period decrease. Comparing the temperature history of all models, it can be seen that the model with an orientation angle of 15 • for the upper crystal shows the maximum interface temperature, whereas the model with an orientation angle of 45 • roughly reaches a temperature of 600 K for the largest sliding velocity during the process. Atoms of the upper crystal with a smooth inclination angle of 15 • , move rapidly on the surface of the lower block and the interface temperature abruptly increases. The smooth inclination angle of 15 • of the upper crystal provides a relative flattened plane, where interface atoms meet. Therefore, atoms encounter fewer obstacles, which could constrain their motion. As a result, atoms move faster and their kinetic energy increases, resulting in an increased global temperature in the interface region. However, the specific collocation of atoms in the crystal structure of 45 • causes atoms to diffuse easily into the 0 • crystal of the lower block at the interface instead of sliding on the surface. To this end, the atomic motion on the contact surface and the nanoscale friction are restricted. Thus, the temperature does not rise substantially.
MD simulations show that a higher sliding velocity leads to a faster atomic migration, which increases the atomic kinetic energy. The global temperature of the atomic structure is obtained by computing the average kinetic energy per degree of freedom as in Eq. 3. Briefly speaking, a higher value of sliding velocity results in larger kinetic energy of the atoms and consequently in an increased global temperature in the interface region. Temperature acts as a driving force for atomic migrations and subsequently the diffusion of atoms at the interface. However, the higher temperature is not necessarily the reason for the higher diffusivity of the interface atoms. This phenomenon can be observed in Fig. 8, where 45 • model has the lowest value of interface temperature among all the proposed models for different sliding velocities. This model, however, shows the highest diffusion coefficients, as shown in Fig. 9(a). As a consequence, significant increase in the interface atomic diffusion is due to the dynamic of the atoms based on the Fig. 8 Temperature evolution history dependent on the sliding velocity sliding velocity for each crystallite inclination angle. This is in good agreement with the results of the work by Yang et al. [62], who introduced the diffusion process between dissimilar metals (Al and Cu) as a result of shear plastic deformations and not only the interfacial heat generation.
Although an increased temperature results in a softening of the material, the diffusion of the interface atoms is also essential in forming a reliable joint. For this purpose, the MSD of the interface atoms is calculated, as presented in Eq. 5. The slope of the MSD curve versus time is proportional to the diffusion coefficient of the diffusing atoms. This represents the diffusion ability of the material. Figure 9(a) shows an increase in the diffusion coefficient with an increase in the sliding velocity for each model. Atomic diffusion is almost negligible in the absence of the sliding velocity, i.e. u = 0.
As the velocity u increases from 0 to 0.3Å/ps for the 45 • crystal angle, the atomic diffusion slightly increases and the diffusion coefficient changes from 0.13 to 0.46Å 2 /ps. A further increase of the sliding velocity causes rapid growth of the diffusion process and results in a diffusion coefficient of 4.46Å 2 /ps for a velocity of 1.1Å/ps. This result on the dependency of the solid state diffusivity on sliding velocity is in good agreement with the results from Samanta et al. [45]. Furthermore, Cheng et al. [9] showed through MD simulations that higher frictional shear velocities mix the atoms at the interface effectively and result in the formation of an interlayer with the high bond strength in the process of vaporizing foil actuator welding of Mg to steel.
Despite the significant continuous growth of the interface atomic diffusion with the sliding velocity in Fig. 9(a),  increasing the compression rate from 0.2 to 0.4Å/ps results in a sudden drop in the diffusion process. This phenomenon may result from the repulsive forces between atoms of the same type at the interface, which prevent the diffusion of Al atoms from one crystal to the other. Further increase of the compression rate v does not show a significant decrease or increase in the diffusion coefficient, which is in good agreement with the results from Mostafavi and Markert [36].
In Fig. 9(a) and (b), it can be seen that the 45 • crystal has the largest diffusion coefficients for all the compression rates v as well as the sliding velocities u. However, from Fig. 8, the interface temperature has the lowest value at each sliding velocity among all the proposed models. One reason to explain this relation can be the effect of the mobility of the atoms at the microscopic level, which is strongly influenced by their temperature and the crystal configuration. To this end, the 15 • crystal has the largest interface temperature and therefore, the highest degree of atomic mobility, but the smallest diffusion ability of all atom configurations. The observed differences in diffusivity between the models are assigned to the arrangement of the atoms in the crystal structure and their resistance against the applied loads during the welding process.
At each sliding velocity, the diffusion coefficient increases when the crystal angle increases from 15 to 30 • and reaches its maximum value for 45 • . With a further increase of the crystal angle, the diffusivity regresses such that the diffusion coefficient of 75 • is less than the corresponding value for 60 • . A similar trend can be observed for the diffusion capability of the samples under compression rates, as shown in Fig. 9(b). Following those observations, two ranges for the interface diffusion regarding crystal orientation can be stated: (1) the diffusion ascending range from 0 to 45 • and (2) the diffusion descending range from 45 to 75 • . Figure 9 clearly shows this particular effect of the 45 • angle, which constitutes a turning point for the diffusion ability of Al crystal as described in this study.
Because of the distinct behavior of the 45 • configuration regarding the interface temperature and the diffusion, this model has been further analyzed for structural changes. Before the joining process begins, the atoms are regularly distributed in the crystal lattice. Therefore, it is most likely to find atoms at the first-neighbor distance of around 2.86Å for Al [21]. As the compression and sliding motions are applied on the upper block, atoms are forced to move and their distribution around each atom changes. With the progress of the simulation, the density of the atoms at a given distance changes continuously and the atoms are rearranged at new sites inside the crystal lattice. The peaks in the RDF graph drop and the valleys rise. These changes indicate the disorder of the atomic structure due to the compression and sliding forces. Figure 10 shows the RDF of the lower crystal, subjected to deformations after initiating the contact with the upper crystal. Although no loading is applied on the crystal structure between the times 196 ps and 246 ps (blue and red curves on Fig. 10), the crystal lattice does not return to its initial form, i.e. a perfect f.c.c structure. This phenomenon can be related to the plastic deformation of the crystal, which indicates the permanent distortion of the structure under diverse loads [52].   Figure 11(a) depicts the influence of the frictional sliding motion of the upper crystal on the atomic arrangement of the lower crystal. It can be seen that by gradually increasing the sliding velocity u from 0Å/ps to 1.1Å/ps, the degree of disorder of the crystal increases. The lattice becomes distorted and the probability of finding atoms at their regular sites inside the crystal changes. However, increasing the compression rate does not significantly affect the crystal structure, as illustrated in Fig. 11(b). The Al atoms are subjected to shear forces during sliding and, therefore, they are displaced to new locations inside the disordered region of the structure. However, repulsive forces arising from the compression of Al atoms impede such dislocations at the atomic scale.
To investigate the effect of the crystal orientation on the load-bearing capacity of the joint in ultrasonic welding, the samples are subjected to uniaxial tensile loading with a constant engineering strain rate ofε = 0.001 ps −1 perpendicular to the interface plane in the y-direction. In order to minimize the influence of the temperature and to investigate the purely mechanical structure-property relationship, the temperature is set to 0.1 K, applying the NVT ensemble. The results are presented in Table 1. The mean value of the maximum tensile stresses after the joint formation is calculated for each model configuration.
The second and fourth columns of the table represent the mean value for different sliding velocities u and different compression rates v. It can be seen that the 30 • configuration has the minimum value of the mean σ , whereas the 45 • configuration has the maximum value of the mean σ . The stress-strain relations under uniaxial tensile deformation of these two models are illustrated in Fig. 12 for two rather moderate values of u and v. The results show that the initial orientation of the aluminum crystallites has a significant effect on the ultimate tensile strength of the joint. As shown in Fig. 12, both crystals of 30 • and 45 • inclination angles show a linear elastic deformation at the beginning, which is followed by discontinuous stress drops. These successive stress drops correspond to local atomicscale rearrangements, also referred to as stress avalanches or strain bursts and govern the deformation behavior of disordered solids in engineering problems [39]. The stress rise and drop points may also indicate the number of local joint spots inside the microstructure, which are detached successively. Richeton et al. [43] statistically characterized the intermittent plastic strain bursts for a single metallic crystal by acoustic emission experiments. In metal alloys, the interacting dynamic dislocations are the sources of deformation avalanches [29]. With the uniaxial tension, the material yields and plastic deformations occur [15]. After reaching the maximum value, the stress decreases abruptly and finally approaches zero. Depending on the inclination angle and the atomic arrangement of the upper crystal, aluminum atoms diffuse at the interface in the atoms of the lower crystal and a joint with a defined mechanical strength is made. Thus, applying a tensile load to the joints between crystals of different orientations results in different stress-strain relations with varying ultimate stresses. In this context, the 45 • orientation reaches a higher magnitude of the tensile stress σ before total failure of the joint occurs, as shown in Fig. 12. The 30 • angle joint experiences a higher tensile strain before the sample ruptures, while the 45 • sample represents a slightly lower tensile strain for the same sliding velocity and compression rate. The difference in the

Summary and conclusion
Although the mechanics of joint formation in ultrasonic welding and similar joining techniques has been studied for decades, the observed dispersion in the joint strength under the same verified process parameters has been remarkable. This study aims to point out that although process parameters are critical influencing factors in the physics of joining, the local configuration of the faying surfaces at the point of bond formation plays a crucial role. This implies that the crystal orientations of the grains at each contact point during the joint formation manipulate the atomic diffusion and, therefore, influence the effect of the process parameters. This phenomenon results in a low or large spread of the joint strength values [37]. Local bonded sites in microscopic observations of Figs. 2 to 6 provide evidence from the atomic point of view.
In this study, MD simulations have been deployed to elucidate the nanomechanics of the complex coupled thermomechanical process of ultrasonic welding within the picosecond timescale. To this end, the influence of the lattice orientation of two aluminum crystallites on the structure and mechanical behavior of the interface has been investigated. Moreover, this study suggests that the sliding velocity significantly affects the temperature and the diffusion pattern at the interface. Finally, the stress response of the joint crystals under tensile loading as a function of crystal orientation was studied.
The main conclusions and achievements of this study are summarized as follows: 1. With increasing the sliding velocity, the diffusion coefficient at the interface of two aluminum crystals increases as well. However, compression does not effectively influence the interface diffusion. 2. For a larger sliding velocity of the faying surfaces, a higher interface temperature is achieved. 3. The initial crystallite structure becomes disordered as the process continues. This reveals the dependency of the structure on the welding time. Plastic deformations and disorder are larger for higher magnitudes of the sliding velocity. However, the influence of the compression rate on the level of structural disorder is marginal. 4. Orientations of the two crystals significantly influence the interface diffusion and interface temperature during ultrasonic welding at the atomic scale. The outcomes of this study showed that the lattice angle of 45 • for the upper crystal results in maximum mean values of tensile stress for six values of the sliding velocity due to uniaxial loading. The same result was achieved for four values of the compression rate.
The atomic scale simulations of the nanostructure can be replicated at the microscale, leading to a bottom-up approach of the bond formation in USMW. Therefore, the current study provides valuable insights into the joint formation at the atomic scale in USMW. Molecular dynamics simulations can help delve into the origin of the bonding method and the influencing factors and improve the process at the macro scale.
Funding Open Access funding enabled and organized by Projekt DEAL. Simulations were performed with computing resources granted by RWTH Aachen University under project ID <rwth0567>. The authors gratefully appreciate this support.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.