Numerical simulation of transport phenomena and its effect on the weld profile and solute distribution during laser welding of dissimilar aluminium alloys with and without beam oscillation

Remote Laser Welding (RLW) of Aluminium alloys has significant importance in lightweight manufacturing to decrease the weight of the body in white. It is critical to understand the physical process of transport phenomena during welding which is highly related to the mechanical performance of the joints. To investigate the underlying physics during welding and to understand the influence of beam oscillation on heat transfer, fluid flow and material mixing a transient three-dimensional Finite Element (FE) based Multiphysics model has been developed and validated from the experiments. The effect of welding speed, oscillation amplitude and oscillation frequency on the fusion zone dimensions, flow profile, vorticity profile, cooling rate and thermal gradient during the butt welding of Al-5754 to Al-6005, with sinusoidal beam oscillation, is analysed. It was found that one additional vortex is formed during beam oscillation welding due to the churning action of the oscillating beam. With the increase in oscillation amplitude, welds become wider and the depth of penetration decreases. An increase in oscillation frequency leads to an increase in the flow rate of the molten metal suggesting that the beam oscillation introduces a churning action that leads to an increase in mixing. It was highlighted that the material mixing depends on both diffusion and convection.


1. Introduction
Lightweight manufacturing is the prime requirement in the automotive industry due to more than ever stricter environmental regulations [1]. Stamping with tailor-welded blanks (TWB) can be used to decrease the weight of the body in white while maintaining low cost and high quality [2]. Aluminium is widely used in the automotive industry for its lightweight and corrosion resistance [3,4]. AA-5754 is frequently used in the automotive industry due to its good formability, machinability and mechanical strength [5]. AA-6005 has improved mechanical strength, but the weldability of AA-6005 is poor [6]. Due to this reason, both these alloys are welded together to utilise the advantages both alloys have to offer [7].
Recently, the demand for dissimilar welding processes has increased considering their design flexibility and economic benefits. Welding of aluminium alloys is challenging due to the high thermal conductivity of aluminium which requires a high power density to weld the aluminium alloys [8,9]. Laser welding is the most popular and widely used welding technique for the welding of aluminium alloys. It provides several advantages such as a non-contact welding process with narrow fusion and heat-affected zone (HAZ), high depth-to-width ratio, less distortion due to less heat input to the workpiece and can be easily automated for high productivity [10][11][12]. But, laser welding has limited applicability due to its small bridging ability when joining imperfect edges such as in TWB [13]. This limitation can be reduced by using the beam oscillation technique as it improves fit-up gap tolerance with wider and smoother welds [14].
The previous studies have suggested that the fusion zone geometry (i.e., weld depth, weld width, overall weld shape) in conjunction with microstructural weld attributes (i.e., cooling rate, thermal gradient) governs the weld quality and properties [5,9,15,16]. Fusion zone geometry, cooling rate and thermal gradient are determined by the history of temperature distribution and the history of melt pool development. These characteristics have an important causal relationship with transport phenomena (heat transfer, fluid flow and mass transfer) during welding [17][18][19]. Numerical simulations of transport phenomena during laser welding can be used to calculate temperature distribution, thermal cycle, cooling rate, flow profile, weld pool geometry and alloying element distribution [20,21]. This can be used for the estimation of weld geometry, solute concentration, solute intermixing layer thickness, weld compositional changes and process monitoring. These numerical calculations have enabled us to gain detailed insights into the physical processes that occur in a weld pool to assess and improve the laser welding process. A thorough understanding of the geometrical weld parameters and molten pool behaviour is critical to using beam oscillation. Several models have been proposed to develop the relationship between heat transfer, fluid flow and the shape of the weld formed [20][21][22].
Several studies were proposed to predict temperature distribution analytically during laser welding [23,24]. Rosenthal [25] first proposed the point and line heat source model for laser welding which was improved by Mazumdar et al. [26] by developing a moving Gaussian distributed surface heat source; however, it is not adequate to simulate the weld pool characteristics. Paley et al. [27] proposed a volumetric heat source and considered the keyhole surface temperature equivalent to the boiling temperature of the material. Several combined volumetric heat sources were proposed to model laser beam welding [17,22,28]. Various three-dimensional heat and mass transfer models are developed by employing the equations of energy, mass, momentum and solute transport in laser welding [29][30][31][32] to calculate temperature distribution, velocity fields and solute distribution in the weld pool. Various studies with high-frequency beam oscillation have shown that the weldability of aluminium alloys has improved with lesser defects and improved mechanical properties [33][34][35][36][37]. These studies paved a way for important development in understanding beam oscillation, but the effect of process parameters on the complex molten pool behaviour is still not well explained. Numerical simulation provides a scientific understanding of the relevant physical details of the complex molten pool that is formed for very short timescales.
In the current state of the art, the effect of beam oscillation on the solute distribution during laser welding of dissimilar aluminium alloys has not been widely reported in the literature. In the present work, a coupled three-dimensional heat transfer, fluid flow and mass transfer model are developed to investigate the effect of process parameters on the temperature distribution, flow profile, weld profile and solute distribution. The thickness of intermixing layer is estimated during laser welding with and without beam oscillation. A hybrid heat source model has been simulated to predict the fusion zone dimensions. The effects of buoyancy, gravity, surface tension and solute diffusion with and without beam oscillation are considered. The developed model is validated with the corresponding experimental measurements. It is worth mentioning that to reach optimum weld parameters, lots of experiments are required. The numerical model intends to reduce or eliminate the number of experiments required, which saves processing time and production costs.

2. Numerical model
A three-dimensional heat transfer, fluid flow and solute distribution (movement of alloying elements) due to diffusion and convection have been developed using a hybrid volumetric heat source. A three-dimensional Cartesian coordinate system has been utilised with the heat source moving along the positive x-axis, the y-axis is along the weld cross-section direction and the positive z-axis is the direction of weld penetration. The oscillating laser beam is composed of two parts in the x-y plane, one is a sinusoidal motion, and the other is a linear movement in the welding direction. The following assumptions were made to reduce the computational time without affecting the accuracy of the model: (a) the fluid is considered Newtonian and incompressible, and flow is considered laminar. To account for change in density due to temperature and concentration variations Boussinesq's approximation is used [20]. (b) The weld pool surface is considered flat, and to simplify the model and decrease the computational time, keyhole dynamics are not considered, and a volumetric heat source model is used to replicate keyhole which is a common practice as evidenced by [17,20,22,28,38,39]. (c) Temperature-dependent material properties are considered. (d) No vapour and plasma flow is considered in the simulation. (e) The absorption coefficient is considered independent of temperature.

Governing equations
To determine the temperature distribution, velocity field and solute distribution, a coupled transient model was developed based on the solution of the equations of conservation of energy, mass and momentum and solute transport as given in Eqs. (1)-(4).
Mass conservation

Momentum conservation
Energy conservation

Solute transport equation
where ρ is the density, t is the time, u is the fluid velocity, η is the dynamic viscosity, P is the pressure, T is the temperature, λ is the thermal conductivity, C P is the specific heat capacity of the material, Q laser is the energy input of the laser heat source and Q vap is the energy loss by evaporation, C i is the concentration of species i, D i is the mass diffusion coefficient for species i and D Ti is the thermal diffusion for species i, F is the momentum source term which is defined as follows: The first term in Eq. (5) is according to the Carman-Kozeny equation for flow through a porous media [17,40] representing the frictional dissipation which ensures a smooth transition of velocity from zero to a large value in the mushy zone, B is a merely computational constant, very small positive number to avoid division by zero is set at 0.001 and J is a mushy zone constant related to the morphology of the porous media which is a large number (a value of 1.6 × 10 4 was used in the present study) to force the velocity of the solid zone to be zero and represents mushy zone morphology, β is the coefficient of volume expansion, g is the acceleration due to gravity, T melting is the melting temperature which is average of solidus and liquidus temperature and f l is the fraction of liquid which is defined in Eq. (6), where, T liq and T sol are liquidus and solidus temperature of the materials respectively. The second term on the righthand side of Eq. (5) accounts for the natural convection.
The phase changes from solid and liquid are considered to include temperature change due to latent heats by using the apparent heat capacity method which includes an additional term for latent heat as shown in Eq. (7) [41] where, C p,solid is the heat capacity of the solid phase, C p,liquid is the heat capacity of the liquid phase, H m is the melting latent heat. (3)

Laser heat source
Here, a hybrid heat source is adopted to simulate the heat flux on the workpiece. The general trajectory of the moving heat source is given in Eq. (8) where s is the welding speed, A is the oscillation amplitude and f is the oscillation frequency.
According to previous research [38,42], laser absorption is described as Fresnel absorption inside a keyhole and inverse Bremsstrahlung absorption above the keyhole. The hybrid volumetric heat source is a combination of a double ellipsoid heat source and a Gaussian damped heat source. The Bremsstrahlung absorption is modelled by the double ellipsoid heat source and Fresnel absorption by the modified Gaussian damped heat source The power density distribution of the double ellipsoid heat sources with origin (x 0 ,y 0 ,z 0 ) is described as follows [43]: where Q f (x,y,z,t) and Q r (x,y,z,t) are the power densities in the front and rear quadrant of the double ellipsoid heat source (Eqs. (9-10)) where P l is the power of the heat source, α is the absorption coefficient of material and (x (t) , y (t) , z (t) ) is the trajectory of the moving heat source which is defined in Eq. (8). f f and f r are the power distribution in the front and rear quadrant of an ellipsoid which follows f f + f r = 2, and in this study, it is set as 0.4 and 1.6 respectively. The shape of the weld pool depends upon heat source distribution parameters obtained by measuring the rear (a r ), front (a f ), width (b) and depth (c) of the half ellipses. The Gaussian damped heat source model suggested by Yuewei et al. [38] was adopted to describe the heat distribution, which is given in Eq. (11) as follows where m is the damping coefficient which is adjusted to minimise the simulated errors; in this study it is set as 0.1.
The effective heat absorbed by the hybrid heat source, i.e. double ellipsoid and 3D Gaussian conical heat source, is given by Eq. (12). f 1 and f 2 are the power distribution coefficient between the modified Gaussian damped and double ellipsoid heat source respectively which follows f 1 + f 2 = 1.

Boundary conditions
The initial temperature of the workpiece is assumed to be maintained at room temperature (T 0 ). At the top surface, the heat loss due to convection, radiation and evaporation is given as follows in Eq. (13).
where h is the convection coefficient, ε is the material's optical emissivity, σ is the Stefan-Boltzman constant and Q vap = WL v , W is the evaporation rate and L v is the latent heat of vaporisation. Flow condition for the free liquid surface due to the surface tension gradient due to variations in temperature and concentration is given as in Eqs. (14)- (15), and the velocity along the z-direction is zero.

Calculation domain and numerical implementation
The numerical model is solved by the finite element method (FEM). The simulations were performed using COMSOL Multiphysics software. The algebraic multigrid (AMG) solver was adopted as it provides robust solutions for large CFD simulations. For problems solved for space and time, discretisation in space is done using the finite element method, and time discretisation is done using the backwards differential formula (BDF). As the system is solved fully coupled, the damped Newton method is used for solving systems of non-linear equations. For solving fluid flow, the generalised minimal residual (GMRES) method is used.
The size of the simulation domain for a single plate is 100 mm × 50 mm × 3 mm with complete mesh divided into 750,639 elements. Butt welding configuration is utilised for this study with the beam at the junction between both plates. Tetragonal mesh is used with a minimum mesh size of r/8 mm. The mesh is fine at the weld centreline and coarser towards the base metal as shown in Fig. 1a. The thermophysical properties of AA-5754 and AA-6005 alloys are taken from [44,45].

3. Experimental details
A 10 KW Coherent ARM FL10000 laser system with an optical fibre of 150-μm diameter was used. The beam parameter product is 16 mm mrad. The laser system coupled with the WeldMaster remote welding head (Precitec GmbH, Germany) has a 150-mm collimating length, a focal length of 300 mm and a Rayleigh length of 2.8 mm. Beam oscillation is generated by using a motorised mirror and collimator integrated with the WeldMaster Scan&Track remote welding head (YW52 Precitec GmbH, Germany). Neither filler wire nor shielding gas was used throughout the experiments.
The materials used in this study were 3-mm-thick sheets of AA-5754 and AA-6005 whose chemical composition is listed in Table 1. The sheets were machined into 100 mm × 50 mm coupons. These coupons are cleaned with Welding experiments are carried out in butt joint configurations with a zero nominal gap, and the incident beam is inclined at an angle of 4°. A high-speed camera (FASTCAM Nova S6) was used to monitor the dynamic behaviour of the weld pool at a photographing speed of 30,000 frames per second. A laser illuminator was used to illuminate the weld zone. The experimental setup is shown in Fig. 1b.
After welding the welded samples were sectioned normal to welding direction and polished to a surface finish of 0.05 μm, followed by etching in caustic sodium fluoride reagent (2%NaOH + 5%NaF + 93% water). The weld dimensions were measured in a Keyence VHX7000 optical microscope. The samples were characterised using Zeiss Sigma scanning electron microscopy. The energy dispersive spectroscopy (EDS) data were acquired using the Aztec software platform. Three repetitions for each parameter set have been performed, and two optical samples were prepared from each parameter set. The sectioning for optical microscopy has been performed 25 mm away from the start and end of the weld seam.

4. Results and discussions
The process parameter combinations used for welding trials and corresponding experimental values for the weld width and weld depth calculated experimentally using optical microscopy are given in Table 2. The setup number is used to refer to the process parameter combination in the rest of the paper. The setups #1 to #5 are selected to understand the welding speed, setups #6 to #9 are selected to understand the effect of oscillation amplitude and setups #6 and #10 to #12 are selected to understand the effect of frequency of oscillation on fusion zone dimensions and flow profile. Figure 2 provides the main numerical results at the three different welding speeds (setup #1 to #3) at a constant laser power and the two different welding speeds at a constant specific point energy density (SPED) [46] (setup #4 and #5). Figure 2a-e shows the shape and dimension of the weld on the transverse section, Fig. 2f-j demonstrates the fluid flow fields in the x-z plane and Fig. 2k-o depicts the vorticity formed at the top surface of the weld pool formed. As expected, the higher the welding speed, the narrower and shallower the weld pool whose dimensional values are given in Table 3. The length of the weld pool decreases with an increase in welding speed (setup #1-3); this is due to a decrease in net heat input to the workpiece which leads to less formation of liquid metal. But at constant SPED (setup #4-5), the length of the weld pool increases with an increase in welding speed while width and depth almost remain constant. The velocity field of the weld pool is indicated by the arrows whose  magnitude can be estimated by the colour of the arrows. The peak velocity increases with an increase in welding speed even at a constant SPED due to the flow caused by the increase in movement of the laser beam. The velocity profile in Fig. 2f-j shows a distinct vortex formed along the vertical mainly moving up and down [47]. The liquid metal tends to flow first through the length section as compared to the depth, and this flow tendency is a result of a change in the temperature gradient. Figure 2k-o shows the obvious vortexes formed near the keyhole opening position on the horizontal plane. A higher fluid flow rate magnitude indicates higher vorticity. The vorticity profile confirms that the liquid metal is flowing from the centre to the outward periphery. This outward flow of the liquid metal is due to the Marangoni convection where flow is from a higher temperature area to a lower temperature area [34]. The values of peak flow rate and vorticity magnitude are given in Table 3. In the core region of the weld pool, especially during the cooling cycle, the temperature gradient increases as welding speed increases. When welding speed increases from 4 to 6 m/min, the thermal gradient increases from 2367 to 2607 K/mm. The thermal gradient during the cooling cycle is important as it governs the grain morphology during solidification [48]. The cooling rate calculated in the core region shows an increasing trend with an increase in welding speed. This is due to the decline in heat input, and the values of thermal gradient and cooling rate are given in Table 3. The cooling rate governs the scale of the microstructure; the higher the cooling rate, the finer the structure will be. oscillation amplitude, the penetration depth and length of the weld pool decrease while weld width increases. This is because the laser beam covers a larger distance in the same amount of time so the energy input per unit length decreases though the net heat input to the workpiece remains constant. It was found that with the increase in oscillation amplitude the mode of welding changes from the keyhole to the conduction mode of welding which can also be confirmed from the peak temperature calculated from the model as the temperature drops down below the boiling point of the material. As compared to no oscillation condition, beam oscillation leads to the formation of three distinct vortices: one vortex along the vertical plane as shown in Fig. 3e-h and two vortices on the horizontal lane as shown in Fig. 3i-l. The two vortices formed in the horizontal plane are distinguishable for the higher oscillation amplitude as shown in Fig. 3k-l. The two vortices formed on the horizontal plane are because of the stirring action due to the application of beam oscillation. The peak flow rate (values given in Table 4) decreases with the increase in oscillation amplitude due to the decrease in the Marangoni driven convection because of the decrease in the temperature. The thermal gradient at the core calculated during the cooling cycle decreases with an increase in oscillation radius while the cooling rate increases. This is due to the less retention of heat with increasing oscillation amplitude. The values of thermal gradient and cooling rate are given in Table 4 for different setups. When compared with no oscillation and beam oscillation conditions at the same welding speed and laser power, penetration depth decreases due to a decrease in the energy input per unit length, whereas there is an increase in weld width and length of the weld pool which suggests that a larger area is in molten form. The thermal gradient of beam oscillation is lower than the no oscillation due to the churning action of the oscillating beam in the weld pool. This is evident from the formation of one more vortex on the horizontal plane. The beam oscillation has a higher cooling rate which suggests that the heat retention for each point is less as compared to no oscillation welding. The peak flow rate is higher for beam oscillation when compared to no oscillation due to the stirring effect caused by the oscillating beam. The decrease in peak flow rate at higher oscillation amplitude, when compared to the no oscillation condition, is due to the decrease in peak temperature which leads to a decrease in Marangoni driven flow.   Figure 4 shows the weld shape in the transverse direction, fluid flow profile and vorticity profile formed with different oscillation frequencies while other process parameters are kept constant. There are not many differences found in terms of weld width and length of the weld pool with increasing frequency. The peak flow rate increases with increased oscillation frequency due to the stirring produced at a higher frequency [49]. The re-heating and re-melting due to sinusoidal beam oscillation are less as compared to circular or infinity beam oscillation as the movement of the beam in case of sinusoidal oscillation is always in the forward direction. The thermal gradient decreases with increasing frequency, and at higher oscillation frequency it almost becomes constant. This is due to the decrease in peak energy accumulation which leads to a decrease in the thermal gradient of the weld pool. The increasing oscillation frequency accelerates the flow and promotes heat transfer which leads to an increase in the cooling rate. The values of the cooling rate and the thermal gradient are given in Table 5. There is a gradual Fig. 4 Weld shape comparison between experiment and simulation (a-d), velocity profile in the x-z plane (e-h) and vorticity profile in the x-y plane (i-l) increase in high-speed flow area with the increase in the frequency of oscillation which is depicted in the flow profile and vorticity profile in Fig. 4. The numerical model is validated by comparing the fusion zone dimensions measured experimentally using optical microscopy and estimated from the numerical model. The comparison between the shape of the fusion zone obtained from the numerical model and optical microscopy is given in Figs. 2, 3 and 4. The measured value obtained experimentally is given in Table 2, and the estimated from the numerical model is given in Tables 3, 4 and 5. Also, the length and width of the weld pool formed on the top surface are used for the validation of the model. The comparison of the weld pool top surface morphology between the experiment and the numerical model is shown in Fig. 5. A high-speed camera is used to obtain the top surface morphologies experimentally. Due to the low contrast between the solid metal and molten pool, the weld pool boundary is recognised due to the fluctuation of the molten metal between the two successive frames. In Fig. 5a-c, the upper image depicts the top surface of the weld at time t and the bottom image at time t + 10 μs. Figure 5d-f shows the simulated weld pool morphologies which are indicated by isothermal contours from liquidus temperature to boiling point. The temperature above the boiling point where the keyhole is formed is shown in white colour in Fig. 5d-f. It can be seen that the keyhole profile is roughly circular (Fig. 5a-c) which matches the simulated contour map (Fig. 5d-f). The simulated contour maps also show that a greater thermal gradient is formed in front of the keyhole. The comparison of the length of the weld pool and the width of the weld pool between the experiment and simulation manifests the validity of the model to investigate heat transfer and fluid flow in laser welding.

Solute distribution for no oscillation and beam oscillation conditions
Transport of solute during laser welding is governed by two mechanisms, namely convection due to fluid flow and diffusion due to temperature gradient. Figure 6 shows the comparison between the experimental and simulated solute distribution. The line scan length and position are shown in Fig. 6a, c. Figure 6b, d shows the variation of the composition of Mg and Si with the distance. The fluctuation in the experimental data of the element concentration is due to the presence of pores, voids and cracks. There is no observed fluctuation in the numerical model as these attributes have not been simulated in the numerical model. Other possible reasons for the fluctuation of the experimental data can be types of the detector, dead time of the detector, oxygen content and software settings [50]. The simulated result is in good agreement with the experimental results which manifests the validity of the model to predict solute distribution during laser welding. Figure 7 shows the concentration of Mg and Si for no oscillation and beam oscillation conditions. AA-5754 is Mg-rich and AA-6005 is a Si-rich alloy. During welding, there is a mixing of Mg and Si from the alloy having a higher concentration to a lower concentration around the fusion zone. This intermixing is due to convection and diffusion. The thickness of the solute intermixing layer depends upon the process parameters. The thickness of the Si-layer for setups #1, #6 and #7 are 4.80 mm, 4.36 mm and 4.30 mm respectively. Similarly, for Mg-layer for setups #1, #6 and #7, the thickness are 3.80 mm, 3.60 mm and 3.55 mm. The thickness of this layer decreases when beam oscillation is applied due to a decrease in thermal gradient and an increase in cooling rate. Solute movement is dominant at the fusion zone boundary due to the diffusion mechanism. It should be noted that the weld width is ≤ 3.1 mm for all the cases and the width of the intermixing layer is more than the weld width. This suggests that the solute distribution occurs due to both diffusion and convection during welding. The thickness of this layer further decreases with the increase in oscillation amplitude due to a decrease in thermal gradient and an increase in cooling rate which diminishes the effect of diffusion at the boundary. Figure 7 also depicts the start and position of the heat source. The start point is 1 mm from the end of the coupons. Therefore, the mixing starts when the welding starts till the position of the heat source which logically follows the results depicted. Also, it was found that the use of mass transfer improves the prediction of fusion zone shape as shown in Fig. 8. The improvement is not very much significant between the red and yellow lines where the red line is when the solute distribution is included in the numerical model and the yellow line is when no solute distribution is included. The maximum difference occurs at the top and concavity of fusion where most of the imperfections occur. Therefore, it is reasonable to say that diffusion coefficients taken for the pure alloy did give consistent results. Future work will be conducted on the development of a process capability space framework [51] to provide a guideline for the welding of aluminium alloys to avoid some common weld defects such as porosity and cracking.

5. Conclusions
This paper studied the effect of weld process parameters such as welding speed, oscillation amplitude and oscillation frequency on the transport phenomena during the laser welding of Al-5754 with Al-6005 in butt joint configuration. A remote laser welding process was employed with transverse beam oscillation, and the effect of process parameters is studied sequentially. A combination of experiments and FEM modelling has been presented to study weld profile and solute mixing within and around the weld. The FEM model predicts the weld profile and provides information about temperature gradient, cooling rate, flow profile and vorticity profile, which are difficult to measure directly using experiments. The simulated model was verified by comparing the cross-section of the weld, top surface morphologies of the weld pool and solute distribution measured experimentally, and the verification was good to predict the welding process. The following conclusions are made: (1) The length of the weld pool decreases with an increase in welding speed at constant power while it increases with an increase in weld speed at constant SPED. The fluid flow velocity increases with an increase in welding speed suggesting an increase in convection. The flow of fluid tends to flow first through the length as compared to the depth, and the fluid flow direction at the top surface is from the centre to the outward periphery due to Marangoni convection. The depth of penetration decreases, and the weld become thinner with increasing welding speed. (2) Two vortices are formed in no oscillation welding one at the horizontal plane and one along the vertical, while in the case of beam oscillation three vortices are formed one along the vertical and two on the horizontal plane. The third vortex formed is due to the churning action of the oscillating beam.  kept constant. This leads to more material mixing by convection. The thermal gradient at the core during the cooling cycle decreases with the increase in oscillation amplitude and oscillation frequency, while the cooling rate increases with the increase in oscillation amplitude and oscillation frequency. (4) Introducing solute distribution in the mass transfer model improves the prediction of fusion zone boundary as compared to the fluid flow model. The solute intermixing layer thickness decreases when beam oscillation is applied as diffusion of solute at the boundary decreases due to a decrease in thermal gradient and an increase in cooling rate. The thickness of the intermixing layer is larger than the weld width which illustrates that material mixing is governed by both diffusion and convection.
Author contribution AM contributed to the conceptualization, defined methodology, developed the numerical model, performed experiments for validation, and formal analysis and manuscript preparation. PF contributed to the development of the experimental setup and manuscript reviewing. DC contributed to supervision, and MA contributed to supervision, and manuscript writing, review and editing.
Funding This work is supported by the WMG, University of Warwick (UK), and the Indian Institute of Technology Kharagpur (India) as the financial support for the PhD studentship via the WMG-IIT partnership program.

Declarations
Ethics approval Not applicable.

Consent to participate Not applicable.
Consent for publication Not applicable.

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:// creat iveco mmons. org/ licen ses/ by/4. 0/.