Semi-solid deformation of Al-Cu alloys: A quantitative comparison between real-time imaging and coupled LBM-DEM simulations

Semi-solid alloys are deformed in a wide range of casting processes; an improved understanding and modelling capability is required to minimise defect formation and optimise productivity. Here we combine thin-sample in-situ X-ray radiography of semisolid Al-Cu alloy deformation at 40e70% solid with 2D coupled lattice Boltzmann method discrete element method (LBM-DEM) simulations. The simulations quantitatively capture the key features of the in-situ experiments, including (i) the local contraction and dilation of the grain assembly during shear deformation; (ii) the heterogeneous strain fields and localisation features; (iii) increases in local liquid pressure in regions where liquid was expelled from the free surface in the experiment; and (iv) decreases in liquid pressure in regions where surface menisci are sucked-in in experiments. The verified DEM simulations provide new insights into the role of initial solid fraction on the stress-deformation response and support the hypothesis that the behaviour of semi-solid alloys can be described using critical state soil mechanics. © 2018 Acta Materialia Inc. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Pressurised casting processes such as high-pressure die-casting [1e3], squeeze casting [4,5], and twin-roll casting [6e8] are widely used in metal processing. These casting techniques often induce shear deformation with compressive stresses on semi-solid alloys containing a solid network, and the complex stress and strain localisation can lead to casting defect formation such as concentrated porosity [9e12], macrosegregation [7,13e16] and shear cracking [15]. These defects have been related to the occurrence of shear-induced dilation [9,10,16e19], and subsequent strain localisation. The behaviour is similar to the Reynolds' dilatancy phenomenon [20] in densely packed granular materials [21], whose load:deformation behaviour is well-described by critical state soil mechanics [22]. In the last decade, in-situ deformation experiments on semi-solid alloys using thin-sample radiography [23e30] and bulk-sample tomography [31e33] have enabled direct observation of semi-solid deformation mechanisms. These have directly confirmed that shear-induced dilation [24,34] leads to grain rearrangement and strain localisation, and can cause macrosegregation [27,34] and shear-cracking [34].
Modelling studies of semi-solid deformation have included thixotropic viscosity-based models [35e37]; strain localisation criteria [38,39]; finite-element microstructure-based models at high solid fraction [40,41]; models of deformation-induced macrosegregation (e.g. in the continuous casting of steel [42e45]); and models of hot tearing under tensile load [46e48]. In most of these cases, the models do not account for isothermal shear-induced volume changes including dilatancy. To address this, here we adopt the particulate discrete element method (DEM), as it is wellsuited to modelling deformation involving grain rearrangement within a solid network with dilatancy, coupled with the lattice Boltzmann method to simulate the behaviour of the interstitial liquid; this approach has not previously been used in semi-solid metals research. In DEM, each particle (grain) in an assembly is explicitly modelled, the force imposed at each grain-grain contact is determined by a specified force-displacement law [49], and grain displacements are determined explicitly by considering dynamic equilibrium. Particulate DEM has been broadly applied to simulate deformation in granular matter such as soils [49], rock-liquid mixtures [50,51], partially molten magma [52], cement pastes [53,54], and powder metallurgy [55,56]. Furthermore, coupling with the lattice Boltzmann method (LBM) allows the behaviour of the interstitial liquid phase to be captured and solid-liquid momentum interactions to be modelled [57e61]. A few studies have also applied DEM to simulate semi-solid alloy deformation [62e64]; Sistaninia et al. [62,63] developed a 3D model for high solid fraction deformation based on a combined discrete/finite element method, while Yuan et al. [64] simulated dendrite rearrangement and shear-induced dilation behaviour as a function of solid fraction during 2D equiaxed dendritic solidification. However, these studies did not apply a coupled DEM-LBM approach to account for the particle-liquid interactions.
In this paper, we begin by analyzing four synchrotron radiography datasets of shear deformation in globular semi-solid alloys at 44e71 vol % solid, each containing~1000e2000 grains in the field of view and approximately one grain thick. The relatively large grain assemblies enable a study of the development of strain heterogeneities and localisation features. The synchrotron datasets are quantified by X-ray intensity analysis and digital image correlation (DIC). A 2D liquid-solid coupled LBM-DEM model is then developed using the initial microstructures from the synchrotron experiments as inputs. A single set of DEM parameters was determined by calibrating the model across all deformation experiments from 44 to 71% solid so that only the initial solid fraction and contact stress were varied in the parametric study. The simulations and experiments are compared quantitatively, and then the simulations are used to extract deeper insights into the semi-solid deformation mechanisms occurring in the experiments.

Synchrotron radiography
Al-Cu alloys containing 8 and 15 wt% Cu and grain refined with 1 wt% Al-5Ti-e1B were cast into a steel mould to create fine equiaxed dendritic microstructures. Slices with thickness 150e200 mm were prepared from the casting. A slice of alloy was placed in the 200 mm thick cavity between two Al 2 O 3 windows, and two BN plates were attached to complete the specimen cell as described in Ref. [26].
In-situ experiments were conducted on beamline 20B2 at the SPring-8 synchrotron using in-situ heating and recording apparatus described previously [26,65]. A 16 keV X-ray beam was used. Signals were recorded at 1 frame per second in 2048 Â 2048 pixel format with 16-bit depth. In each experiment, the whole deformation system was heated in the furnace to a temperature in the solid þ liquid region, and a semi-solid alloy with globular morphology was attained by isothermally holding the equiaxed microstructure for 30 min. Isothermal shear deformation was then applied by the upward displacement of a mobile Al 2 O 3 push plate at 30 mm s À1 . Experiments were performed on samples containing 44%, 50%, 66%, and 71% solid. The experimental parameters are listed in Table 1.

Synchrotron image processing and quantification
In order to calculate the initial volumetric solid fraction, g 0 S , Xray intensity processing was applied to the radiographs. For an isothermal mush with fully mixed liquid, in each transmitted image the intensity at a pixel, I SL , is mainly a function of the solid fraction in that beam path, the X-ray absorption coefficients, and the sample thickness [23,66]: where I 0 is the incident X-ray intensity, m L , m S , and m cell are the X-ray absorption coefficients of the liquid, solid and shear cell (i.e., Al 2 O 3 windows). L alloy , L S are the length of semi-solid alloy and solid in the beam path respectively. As the deformation experiments were conducted isothermally, the X-ray absorption coefficients can be regarded as constants. For the initial radiograph I 0 SL ( 0 indicates parameters prior to deformation), it is assumed that the sample thickness L 0 alloy was uniform, so the intensity of the transmitted beam through the alloy I 0 SL is only a function of L 0 S and can be linked to the solid fraction g 0 S:A in a given averaging area A [23]: where I 0 L is the transmitted intensity through 100% liquid and equals to I 0 SL ðL 0 S ¼ 0Þ, I 0 S is the transmitted intensity through 100% solid and equals to I 0 SL ðL 0 S ¼ L 0 alloy Þ, and A is the averaging area on L 0 S = L 0 alloy to derive solid fraction, g 0 S:A . The I 0 L and I 0 S fields were evaluated from the local intensity minimum/maximum images I 0 SL:min ; I 0 SL:max using a rolling-ball algorithm [67] on the radiograph I 0 SL prior to deformation. Throughout the paper, the average stress and strain measurements were determined by using a consistent normalised REV size to account for the effect of grain size on the quantification of stress and strain heterogeneity. In all cases the REV size was normalised by the grain size of the corresponding experiment/ simulation. The rolling-ball radius was the nearest integer of ffiffiffiffiffiffi 10 p d= 2 (d is the mean grain size measured by a line intercept method) in order to obtain representative minimum and maximum intensities.
The I 0 SL:min and I 0 SL:max can be simply rewritten as I 0 SL ðL 0 S ¼ L 0 S:min Þ and I 0 SL ðL 0 S ¼ L 0 S:max Þ where L 0 S:min and L 0 S:max are the local minimum and maximum thickness of solid phase defined from rolling-ball scanning, respectively. To relate the I 0 SL:min and I 0 SL:max fields to the I 0 L and I 0 S fields, two multipliers r L and r S were introduced to let r L ,I 0 SL:min converge to the centre of liquid interstices, and let r S ,I 0 SL:max capture the brightest pixel of grains. The volume fraction of solid, g 0 S , was then obtained by averaging L 0 S =L 0 alloy values within the FOV, which can be represented as g 0 with the assumption of I 0 L ¼ r L ,I 0 SL:min and I 0 S ¼ r S ,I 0 SL:max . With this approach, the measured initial volume fraction of solid was 44%, 50%, 66%, and 71% for these four samples.
The initial volumetric solid fraction, g 0 S , coupled with the bulk alloy compositions were applied to calculate the X-ray absorption coefficient of the liquid, m L , and solid, m S [68]: where r L and r S are liquid and solid densities calculated using reference [69], w Al;L , w Cu;L , w Al;S , and w Cu;S are the mass fractions of aluminium and copper in the liquid (L) and solid (S) from the phase diagram [70], and ðm=rÞ Al , ðm=rÞ Cu are the mass X-ray absorption coefficients of pure aluminium and copper at 16 keV X-ray energy [71]. The evaluated m L and m S were used to calculate L 0 alloy and L 0 S:min prior to deformation: With this approach, the initial thickness of the semi-solid was calculated to be 180 ± 30 mm for the four samples, consistent with expectation. The solid fraction field, g 0 S:REV , was obtained for images before and during deformation by volume averaging L 0 S =L 0 alloy over projected areas of ffiffiffiffiffiffi 10 p d Â ffiffiffiffiffiffi 10 p d which is a representative elementary volume (REV) that is sufficiently large to average the local microstructure and small enough that important variations in the solid fraction g 0 S:REV can be captured [72,73]. The averaging REV size and X-ray absorption coefficients are summarised in Table 1. A detailed example of finding the g 0 S and g 0 S:REV fields from radiographs is illustrated in Supplementary Information (SI) Section 1.
During deformation, the local sample thickness can change. The sample thickness field in frame n, L n alloy can be found through rearranging Eq. (1) and obtaining I n SL:min by applying the rolling-ball algorithm on each radiograph I n SL : where L n S:min is the local minimum thickness of solid phase. Alternatively, L n alloy can be evaluated from the transmitted intensity through a point of 100% solid in frame n, I n S through rearranging Eq. (1): Whether Eq. (5) or Eq. (6) is better suited to calculating L n alloy depends on whether L n S:min or I n S can be more accurately measured. To proceed, we identified regions that meet one of the following three criteria.
For regions with sufficiently small deformation, the change of L n S:min can be neglected. Therefore, it was assumed that for any REV with an average grain displacement < 0:05u n y , L n S:min ¼ L 0 S:min in Eq. (5). (u n y is the y-component of push-plate displacement in frame n): For regions with a sufficiently large liquid interstice, the I n SL:min can be regarded as the intensity of liquid, I n L . Therefore, it was assumed that for any REV with a liquid interstice ! 0:2d Â 0:2d, L n S:min ¼ 0 in Eq. (5). For regions with highly-compacted solid grains (i.e., in front of the push plate), the I n SL:max can be regarded as the intensity of the solid, I n S . Therefore, it was assumed that for any REV with average grain displacement > 0:95u n y and u n y > 0:25 mm, the I n S ¼ I n SL:max and L n alloy was found using Eq. (6).
Sufficient regions met one of these three criteria that the whole L n alloy field could be computed from interpolation of z100 scattered L n alloy ðx; yÞ points per frame. A step-by-step guide to the procedure is given in SI Section 2. The volumetric solid fraction field for radiograph n, g n S:REV , was then obtained: The volumetric strain of any representative element, ε n vol , was evaluated from the volume-averaged solid fraction g 0 S:REV and g n S:REV using Eq. (9) which is derived in SI Section 2:

Digital image correlation (DIC) analysis
DIC analysis was applied on the synchrotron datasets using the 2D Strain Module in DaVis 8.3 software (LaVision Imaging Company, G€ oettingen, Germany). The side length of the square subset used for correlation was the nearest odd integer of ffiffiffiffiffiffi 10 p d pxl (d pxl is the mean grain size in pixels) similar to the REV for deriving g S:REV in order to offer sufficient greyscale pattern for correlation while mapping heterogeneous strain. The step size was set to the nearest integer of ffiffiffiffiffiffi 10 p d pxl =4 to reconstruct the displacement vector field with reasonable resolution. The DIC parameters used are listed in Table 2. The accumulated displacement was conducted by Table 1 A summary for in-situ observation on semi-solid deformation experiments performed at the SPring-8 synchrotron facility, Hyogo, Japan. The push-plate moving rate was measured directly from time-series images. The average semi-solid alloy grain size was measured by applying the line interception method on the initial radiographic projection image (ASTM E112-96 standard test method). The volume fraction of solid was evaluated by applying a volume averaging method to the initial image (see Supplementary Information 1). summing all displacement vectors from the correlation between the 1st/2nd frames to (n-1) th /nth frames for each calculation point. The conversion to rotation-free strain fields was performed by partial differentiation on the accumulated displacement vector field and using the polar decomposition algorithm [74]. Finally, all of the displacement vector fields and strain fields were filtered by the bilinear interpolation method [75,76] to smooth out discontinuous features.

Discrete element method (DEM) simulations
Since DEM is not commonly used in semi-solid alloy modelling, the general principles are outlined here. DEM is a dynamic modelling technique that idealises an assembly of solid particles as quasi-rigid bodies in mechanical contact. Initial and boundary conditions are specified and the main calculations take place at discrete time intervals; these calculations determine the forces and moments acting on each particle/grain, and update the particle positions and orientations. The force and moment (scalar in 2D DEM) components considered here arise from Ref. [49]: 1. Pair-wise contact: direct interaction between two grains, or between a grain and wall. These are determined via a contact detection phase that identifies contacting particles and using contact models that relate grain overlap to a contact force. 2. Pair-wise bonding: non-contact bonding that allows tensile contact forces between grains (this is used to simulate the membrane boundary here) 3. Gravitational (body) force. The DEM model used here does not include a gravitational field since the movement of grains in the experiment is observed to be majorly controlled by the physical contact between grain surfaces rather than gravity. In each DEM calculation cycle, the dynamic equilibrium of each particle is considered to determine an acceleration, a Verlet-type explicit time integration is then used to determine the particle displacements and rotations for the current time step from the acceleration value. The particle positions, orientations and velocities are updated and these new positions and velocities are used to calculate the contact and drag forces for the next time increment. DEM is algorithmically similar to molecular dynamics. 4. Drag force and moment from liquid flow.
Digital 2D grain assemblies were generated to represent each microstructure from the thin sample (~1 grain thick) experiments. To obtain the size and shape distributions, 120 grains were picked from a radiograph prior to deformation for each dataset ( Fig. 1a and b). The perimeter of each grain was identified using the "point picker" plugin in ImageJ (National Institute of Health, USA), and a least-squares elliptical function was fitted to the perimeter points, as shown in Fig. 1c. The usage of the elliptical function was because this gave a good fit to the actual globular grain shape while capturing the asphericity in a computationally efficient way. To reduce the calculation cost further, each ellipse was simplified as two equivalent overlapping circles (Fig. 1d). Grain size and shape distributions were reproduced by picking and simplifying 120 grains from the radiograph and converting them into 20 grain templates to generate a grain assembly in the DEM package PFC2D Ver. 5.0 (Itasca Consulting Group, Inc). To generate a DEM sample, these grains with normal stiffness 1 Â 10 4 N=m were constrained by a rectangular box with the same size as the thin-plate sample for each experiment, and then isotropic compression was applied by a servo-controlled mechanism until a target 2D stress of 500 N=m was reached.
Examination of experimental radiographs showed that the outside Al 2 O 3 wall did not touch the semi-solid samples as shown in SI Section 3 ( Fig. SI-6) except for the left edge of the 71% solid sample and there was only a small degree of surface deformation (with displacement < 0:15d). We can infer from this that the sample is constrained by surface tension, i.e. the motion of a grain at the free surface would increase the local curvature of the liquid-air interface, which would decrease the local liquid pressure and apply a suction that resists grain motion. In this way, the surface grains have an apparent cohesion and surface tension stabilises the free surfaces against gravity and prevents slumping. To simulate the apparent cohesion, the rigid walls were deleted and the surface grains were assigned partial cohesion using parallel bonds in PFC2D that create a movable membrane boundary that constrains the packing of grains similar to reference [77,78]. These parallel-bonded surface grains are shaded red and linked by red lines in Fig. 1eeg. With this approach, a constraint force, F N , would be installed on any membrane grain if the grain started to displace outward, and the direction of F N for the membrane circle was determined from the position of its two adjacent membrane grains ( Fig. 1h and i). The sample generation stage finished when the maximum magnitude of grain velocity in the system was less than 10 À7 m/s, and the average contact force magnitude excluding membrane circles, F 0 , was extracted to obtain the extent of the pre-existing force network for each sample. The initial 2D packing fraction, f 0 S:pk , was measured for each assembly by randomly distributing 5000 measurement circles with diameter 20d across the whole assembly (see SI Section 7 for the descriptions of measurement circle). To simulate the push-plate, a new rigid wall was introduced as illustrated in Fig. 1e. The use of measurement circles as a representative volume (RV) follows measurement approaches in the DEM literature [79e81]. The concept of a circular RV in DEM analysis here is different from the use of square/rectangular REVs for image analysis suggested by Ref. [23].
The key parameters used for the DEM simulations are summarised in Table 3. The solid phase densities differed in the four virtual samples because experiments were either carried out at different temperatures or with different bulk alloy compositions to achieve the range of solid fractions, and the density was calculated using the composition and temperature-dependent equations of ref. [69]. However, the simulated results would not be significantly altered if a single density had been used across all simulations. A low sliding friction coefficient, m, was used at contacts as the liquid phase is expected to provide lubrication, while the rolling-  Correlation subset size  pixel  101  151  165  139  Correlation step size  pixel  25  38  41  35 resistance coefficient, m r , was specified from the consideration of the non-spherical shape of particles and the existence of contact surfaces rather than point contacts [82]. High bonding strengths were assigned to the membrane contacts to prevent bond breakage. Past DEM studies have shown that, provided changes in the mean grain size do not significantly alter the number of grains in the sample, many of the deformation features scale with the mean grain size [83,84]. The contraction or dilation behaviour is sensitive to initial solid fraction and the average grain size plays a more minor role, mainly affecting the magnitude of volume change.
To account for the plastic deformation of grains during shear deformation, the contact force-overlap relationship in Fig. 2a was applied. This represents the flow curve of Al alloys deformed at constant shear rate without work hardening after yielding ðU y ; F y Þ at homologous temperature near to 1.0 [85]. The setting of an elastic unloading/reloading path between an unloading point ðU un ; F un Þ and residual overlap point ðU res ; 0Þ shown in Fig. 2a enables the grain-grain contact to follow the behaviour of cyclic loading in alloys [86]. This plastic force-overlap relationship did not apply on the grain-wall contacts or any bonded contacts between two circular membrane grains. The membrane boundary and plastic forceoverlap relationship were incorporated using the embedded scripting language FISH in PFC2D, and the detailed code descriptions are in SI Section 4.  1. The procedure to simulate semi-solid deformation adapting microstructures from deformation experiments: (a) in-situ synchrotron radiography imaging on 66% solid sample before deformation, (b) 120 particles selected for shape analysis, (c) circled particle in (b) fitted using an elliptical equation, and (d) the corresponding two-circle grain used to represent the shape of the fitted ellipse in (c). After adapting particle size and shape distribution by measuring 120 particles, a grain assembly (simulation C) was created which was constrained by the push plate and membrane as shown in (e), and (f) is the deformation microstructure during push-plate displacement with rate duy=dt ¼ 30 mm s À1 . The black lines in (f) indicate the contact force chain and the thickness of each line is proportional to the contact force magnitude. The blue solid square in (eef) indicates the corresponding field of view in X-ray imaging. (geh) The cropped region from the dash-line square in (e) and (f), respectively. (h) Showing that constraint force normal to the membrane was applied as the membrane boundary deformed, and (i) the determination of normal force direction from two adjacent circles. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

LBM-DEM two-phase coupling
In order to incorporate the liquid phase, the lattice Boltzmann method (LBM) was applied and LBM-DEM coupling as proposed by Cook et al. [57] was used. The LBM solves the Navier-Stokes equations for the liquid on a regular grid. Each grid node contains a packet of liquid particles that are allowed to move in horizontal, vertical and diagonal directions. The grid nodes in the vicinity of a liquid-solid surface act as a no-slip bounce-back boundary to the interstitial liquid. The detailed description of governing equations and key parameters in the LBM-DEM coupling adopted in this research is given in SI Section 5.
The rectangular LBM grid extended beyond the DEM sample boundaries (Fig. 2b). The lattice spacing was 40 mm, and the average grain size d to lattice spacing ratios were 2e3; an immersed boundary scheme was used to achieve sub-lattice resolution. The liquid density, r L , calculated from Ref. [69], dynamic viscosity, m L , from Ref. [87], and LB liquid calculation timestep are listed in Table 4. The kinematic viscosities (m L /r L ) used in the four simulations were in the narrow range 6.10e6.75 Â 10 À7 m 2 /s, showing that the Cu concentration in the liquid has a minor effect on the liquid flow behaviour. Following past work on interstitial liquid flow and permeability in 2D simulations of packed granular materials and porous media [88e90], the concept of a "hydrodynamic radius" was used and set to be 0.8. This enables the liquid and solid phases to both exist as percolating networks as in the thin sample experiments. Before deformation, the liquid velocity was initialised as zero. The LB bounce-back boundary condition [91] was applied on boundary nodes (dashed black line in Fig. 2b), and the LB nodes outside of the DEM membrane boundary were identified as an effective "solid zone" to inhibit liquid Al from flowing out of the DEM sample. The solid zone applied in the LBM simulations has no effect on the DEM simulation of the grains. When the membrane deforms, the LBM solid zone adapts to the new shape of the membrane.  During shear deformation, the push plate acts as a moving bounce-back boundary [92] and nodes overlapping with the push plate are automatically identified as solid. The liquid pressure change Dp in response to shear deformation was derived from the change of liquid density for each LBM node, and the speed of sound c s , in the LBM domain (detail in SI Section 5): where r 0 is the initial liquid density listed in Table 4, and rðx i ; tÞ is the liquid density at position x i at simulation time t. The liquid drag force and moment acting on DEM grains were updated by a cycle of the LBM-DEM coupling algorithm (see SI Section 5 for detail).

Results and discussion
3.1. Analysis of radiography image sequences Fig. 3aed shows processed in-situ X-ray images from four datasets prior to deformation. Each image shows the solid fraction in the beam path, L 0 S =L 0 alloy , at each pixel, calculated from the transmitted intensities using Eq. (2) and SI-Eq. (3). In this conversion, white represents solid (L 0 S =L 0 alloy ¼ 1) and black represents liquid (L 0 S =L 0 alloy ¼ 0). Pixels with L 0 S =L 0 alloy value not falling in the range [0,1] are coloured red, and are mainly regions of oxide and pores. Fig. 3eeh are the L n S =L n alloy field after a push-plate displacement of 5d for each dataset. The changes in intensity of the L n S =L n alloy shading indicate changes in the local solid fraction due to grain rearrangement and liquid flow.
The L 0 S =L 0 alloy fields were volume averaged using a REV of ffiffiffiffiffiffi 10 p d Â ffiffiffiffiffiffi 10 p d to obtain the solid fraction field prior to deformation (Fig. 3iel) where the red pixels in Fig. 3aed were excluded from the averaging process. Fig. 3iel shows that the solid phase distribution before deformation is slightly segregated with a higher solid fraction nearer the top of the FOV. This is especially evident for the 44% solid dataset and is due to grain buoyancy in Al-15Cu where the density of the solid is lower than the liquid [69]. After a 5d increment of push-plate displacement, Fig. 3mep shows the development of more heterogeneous solid fraction field. In all samples, Fig. 3mep shows accumulation/compaction of grains and expulsion of liquid (hotter colours) in front of the push plate where the local g n S:REV is increased to z0:8. Additionally, for the 66% and 71% solid samples, a reduced solid fraction and increased liquid fraction (cooler colours) with g n S:REV from 0.4 to 0.5 is present ahead of the compaction zone.
The local changes in solid fraction correspond to volumetric strains in the grain assembly, which are plotted as volumetric strain maps ε n vol in Fig. 3qet using Eq. (9) on the volume-averaged solid fraction maps. It is clear from Fig. 3qet that, for initial solid fractions of 44% and 50%, deformation is mostly contractive (blue) with only small pockets of local dilation (red) and, as the initial solid fraction increases, an increasing proportion of the sample undergoes shear-induced dilation (red) and the magnitude of local dilation increases. In Fig. 3q and r, the push-plate displacement on the low solid fraction datasets (44% and 50%) developed a compaction field with ε n vol z À 0:3, while the compaction zone in 66% and 71% solid datasets is less significant (Fig. 3s and t), where ε n vol falls in the range À0.2~À0.1. It should also be noted that the negative ε n vol regions to the right of the FOV in Fig. 3s and the rightbottom of the FOV in Fig. 3t represent local contraction corresponding to the suction of liquid away to the dilating regions elsewhere in the sample.
While there are differences between bulk 3D deformation and the current study on thin sample deformation with interactions with the confining walls, we note that the key features measured here have also been measured in bulk 3D tomographic imaging [31] and in microstructures after semi-solid deformation experiments. For example, shear-induced contraction in samples containing a loose solid network [17], shear-induced dilatation in samples containing a dense solid network [9,19], and the localisation of dilatancy [9].
The main deformation mechanisms that increase and decrease the solid fraction and cause volumetric strain in Fig. 3qet, and the resultant heterogeneous solid fraction field in Fig. 3mep are shown in Fig. 4 using the 66% solid sample as an example. In region A, the increase in solid fraction (contraction) can be seen to be due to grains being pushed closer together and expelling some interstitial liquid and, in this case, also due to compression of the individual grains. In region B, the decrease in solid fraction (dilation) is due to grains moving apart and liquid being drawn in interstitial. In region B, all grains are in mechanical contact with their neighbours and, therefore, the shear-induced dilation (dilatancy) is caused by grains pushing one another apart under the action of compressive and shear contact forces.
2D displacement fields were obtained by tracking the motion of grains using DIC. Fig. 5aed shows the incremental displacement vectors for a push-plate displacement of 0 to 5d for each dataset, and the background image is the time averaged radiographs from 0 to 5d push-plate displacement for comparison. The vectors in Fig. 5aed represent grain movement and show that, for the same push-plate displacement, grains are displaced deeper into the sample as the initial solid fraction increases. The displacement vector field in DIC can be directly converted to shear strain (ε xy ) fields as shown in Fig. 5eeh, where the background for each shear strain field is selected as the image at 5d push-plate displacement. In Fig. 5eeh, it is clear that the ε xy is highly localised above the topright corner of push plate. The red region near to the parting plane in the 44% solid sample (Fig. 5e) is relatively small, while the red bands of localised shear strain form in 50% solid (Fig. 5f). A moredeveloped region of shear strain along the parting plane can be found in Fig. 5g and h which is wider and deeper than in Fig. 5eef. From Fig. 5eeh, it can be seen that at the same push-plate displacement a higher solid fraction induces a larger area of shear strain (red) around the parting plane. The formation of blue regions in Fig. 5g and h is due to ε xy < 0 development near to the left parting plane aligning with the top-left corner of push-plate (outside of the left edge of the FOV).  Table 3.

Numerical simulation results and comparative study
It was found that a single set of optimised LBM-DEM parameters could reasonably reproduce the displacement and strain fields in all four experiments while varying few parameters: the initial packing fraction (i.e. the solid fraction), the initial average contact stress, solid/liquid densities, and liquid viscosity directly linking to corresponding experiments. Note that, while the solid/liquid densities during deformation. (qet) The evaluation of volumetric strain ε n vol using Eq. (9). Images in the leftmost column correspond to the 44% solid sample, followed by 50%, 66%, and 71% solid samples, respectively. White-line scale bars read 2mm. and liquid viscosity were calculated for the relevant temperature and phase compositions, the variation in these parameters is small and the simulated results would not be significantly altered if a single density and liquid viscosity value had been used across all simulations. Thus, effectively two parameters were varied. The optimised set of parameters including contact stiffness values and friction coefficients is shown in Table 3, and was found by an iterative approach to reproduce the deformation microstructure in each experiment across the full range of applied displacements. This was done by starting with literature values for stiffness [64,93] and friction coefficient [64,82,93,94], and then varying these to obtain a good correlation (e norm in Eq. (11) within 10%) of the 1D displacement profile between simulation and experiment.
To enable direct comparison DIC analysis was performed on DEM simulated data, using the same DIC techniques and parameters as in the experimental analysis. To do this, DEM grains were coloured by random greyscale to construct a speckle pattern, the pixel size and FOV of DEM microstructure plots were adjusted to be the same as in radiography images, and the same DIC parameters were used (Table 2). Fig. 5iel shows the result for the displacement fields, and a good correlation of vector fields between experiment (Fig. 5aed) and simulation (Fig. 5iel) can be seen. It can also be seen that the DIC shear strain fields from the simulations (Fig. 5mep) correctly capture many of the shear strain localisation characteristics in the experiments (Fig. 5eeh). For example, the shape and extent of the red localised shear zones, and the influence of initial solid fraction on the shear strain fields are well reproduced. Also note from Fig. 5g and h, that the shear strain is positive just in front of the push plate for the 66% solid sample while it is negative for the 71% solid sample. This is because the boundary conditions are slightly different for these samples. In the 66% solid sample (Fig. 5g), the sample sides did not touch the Al 2 O 3 side walls and, in the 71% solid sample, the left edge of the sample is in contact with the Al 2 O 3 side wall (as introduced in section 2.4). The simulations in Fig. 5o and p used the boundary conditions measured for each experiment and gave a similar result. This is because constraint from the left side wall prevents grain displacement to the left, resulting in a change in shear strain direction.
The displacement vector fields in Fig. 5iel were then interpolated and scanned to construct grain displacement profiles as shown in Fig. 6a, which plots the grain displacement relative to the push-plate displacement at various distances ahead of the pushplate along a 1D line starting at the middle of the push-plate as indicated by the vertical dashed line in Fig. 5m. Comparison of the 1D displacement profiles between experiments and simulations after a 5d increment of push-plate displacement confirms that the simulation reproduces the gradient of the displacement as well as the higher relative displacement in higher solid fraction datasets. To quantify this, Fig. 6b is a bar chart of the normalised error of simulated displacement profiles compared with the experiments, e norm , with increments of push-plate displacement from 1d to 5d. The normalised error is defined as: From Fig. 6b, it can be seen that the error value e norm is <10% for all experiment/simulation pairs, indicating good agreement between simulation and experiment.
We next examine how closely the simulations can reproduce the heterogeneous strain patterns observed in the radiography datasets. Comparison between localised strains can be simplified by the conversion from two-dimensional strain fields to averaged onedimensional line profiles by averaging along x or y within a ROI, as shown from Fig. 7 to Fig. 10. In  transitions of contraction-to-dilation response ahead of the push plate. The strain profiles acquired from triangulation are "jumpy" because the strains are calculated from the relative position changes of groups of three adjacent grains without any averaging by a REV.
It can be seen in Figs. 7e10 that, in general, there is good agreement between experiments and simulations. Specifically, the strain profiles from the 44% solid fraction experiment and simulation A (Fig. 7) show net-contraction in front of the push plate from the ε n vol profile in X-ray intensity processing on the experiment and divergence (ε xx þ ε yy ) profiles by DIC and triangulation on simulation A in Fig. 7c and e, and the distinctively negative ( < À 0:1) volumetric strain or divergence to yz12d in the blue ROI (Fig. 7e) and at x < 0 in the red ROI (Fig. 7c) is captured well. The average shear strain profiles show minimal strain localisation in both the experiment and simulation (Fig. 7e and f).
For deformation at a slightly higher initial solid fraction (50% solid sample and simulation B) shown in Fig. 8, the volumetric strain ε n vol in the experiment and divergence profiles in the simulation are contractive in Fig. 8c and e, but less negative than in Fig. 7c and e. The shear strain profiles (Fig. 8d) go higher than in the 44% solid/simulation A of Fig. 7, indicating more strain-localisation at higher solid fraction. Shear strain development occurs ahead of the push plate to yz15d in both the experiment and corresponding simulation B, as shown in Fig. 8f. Fig. 9 shows that deformation at higher initial solid fraction (66%) gives rise to localisation of both positive ε n vol and shear strain ε xy in a similar location (xz À 5 $ À3d in the red ROI of Fig. 9c and  d), which indicates a dilatant region [26]. The broad positive (dilative) ε n vol and divergence zone to xz12d in Fig. 9c is a distinctive difference of strain response compared with the contractive low solid fraction datasets in Figs. 7c and 8c. In Fig. 9d, peaks of the shear strain profile are approximately located at the parting plane (xz0). The co-location of the positive zone in divergence/volumetric strain profiles and ε xy peaks in Fig. 9c and d gives a dilatant shear band width at 66% solid within the range 10d to 15d, which is consistent with the measurements in past bulk rheology tests [9,17,95]. Analysis of the blue ROI, captures the contraction field immediately ahead of push-plate (Fig. 9e), followed by relatively smooth and positive profiles (normalised y axis from 2d to 21d), in both the experiment and corresponding simulation C. The contraction field at the push plate front is due to local compaction at this location. Fig. 9c captures another negative (contractive) volumetric strain region to the far-right of the red ROI starting at x ¼ 12d, both in the ε n vol profile of the experiment and the divergence profiles from LBM-DEM simulation B, due to liquid drainage at this location (i.e. the suction of liquid from the far-right of the FOV to feed dilation near the parting-plane). Fig. 10 is the response of the highest solid fraction dataset studied here (71%) which had the most net-dilative response and simulation D. Fig. 10c shows a liquid-enriched dilation zone at x < 7d, followed by a contraction zone to the right of the red ROI (negative ε n vol or divergence at x > 10d) both in the experiment and simulation D strain profiles. In the red ROI, both the experimental   ε n vol profile and simulated ε xx þ ε yy profile by triangulation show more positive dilative strains ahead of the push plate (x < 0d) than in Fig. 9c, but the peak width of the ε xy profile in Fig. 10d is similar to Fig. 9d (from 10d to 15d). The blue ROI contains a very small contractive region immediately ahead of the push plate, both in the 71% solid experiment and simulation D (Fig. 10e), and the ε n vol profile in the experiment then seems to have a periodic pattern with multiple peaks and valleys. The explanation of this phenomenon is incomplete liquid feeding during deformation at high solid fraction. As a result, the ε n vol field of the 71% solid sample (Fig. 3t) shows some island-like contractive areas where liquid is sucked to adjacent dilating zones. Those locally dilative zones correspond to peaks at yz10:5d and 17:5d in Fig. 10e. The complex heterogeneous liquid flow phenomena were not completely captured by the LBM-DEM simulation, but a reasonable correlation with experiment in Fig. 10e can still be observed.
While the simulations are 2D, the thin-sample experiments can be considered to be 2.5D with interactions from the confining walls. For the contractive region ahead of the push plate, the LBM-DEM simulation captures grain deformation by overlap between DEM grains which is tuned in the simple approach to plasticity in Fig. 2a to match the overlap in the through-thickness averaged radiographs. Therefore, although the 2D simulation model does not capture the change in thickness directly, it does capture the 'effective grain overlap' in regions of thickening due to compaction. For other regions, the change of sample thickness was low with overlap between grains of <0.1 d both in experiments and simulations, where the assumption of in-plane deformation is reasonable.
Comparison of all experimental and simulated profiles in Figs. 7e10 demonstrates that the coupled LBM-DEM simulations quantitatively captured many of the key features of the heterogeneous volumetric and shear strain fields at solid fractions in the range 44e71 vol% solid, by only altering the initial packing fraction and initial average contact stress in the simulations. An important feature of the LBM-DEM simulations is that the local changes in packing density arise naturally as emergent phenomena. For example, in Fig. 4, Region B undergoes shear-induced dilation due to grains pushing each other apart in response to contact forces in both the experiment and simulation. The LBM-DEM simulations also reproduce further important phenomena observed in the experiments. For example, Fig. 11a shows that as the push plate continued to 11d displacement in the lowest solid fraction sample (44%), liquid was expelled from the sample-air interface at the topleft of the sample, indicating that an excess liquid pressure developed in this region. The simulated liquid pressure field in Fig. 11b correctly predicts a positive (red) change in liquid pressure in this region at 11d displacement and agrees with the location of the liquid expulsion. In Fig. 11b, note that the liquid pressure is highest on the lower-left-side of the Figure. However, the lower-left-side of the Figure is not a free surface; it is the edge of the experimental FOV. In SI Section 9, the full simulation domain is shown (which equals the whole experimental sample), where it can be seen that the true left-side free surface of the sample has a liquid pressure that is lower than the liquid pressure at the top of the FOV in Fig. 11b. It also shows that the region where liquid is expelled in Fig. 11a is the region with the highest liquid pressure near a surface in Fig. 11b.
On the other hand, Fig. 11c shows that, in the highest solid fraction sample, menisci are sucked in to the liquid at the free surface near the right edge of the push plate during deformation. This is in agreement with the simulated liquid pressure field in Fig. 11d where a decreased liquid pressure (blue) region exists at this location. The different liquid pressure behaviour reflects the net-contractive response at low solid fraction (44%, Fig. 3q) and net-dilative shear response at high solid fraction (71%, Fig. 3t), and is a result of the coupling between local volumetric strains in the grain assembly and the liquid pressure field. This behaviour is expected to depend on the permeability and the strain rate and it was found that, when similar samples with initial solid fractions >73% were deformed, shear cracking occurred (not shown here).
The simulated push plate penetrates the compacted grain assembly in the DEM domain, and simultaneously acts as a moving solid boundary on the LBM liquid domain. Fig. 12 illustrates the simulated liquid flow field after a push plate displacement of 5d. This field was constructed by colouring the LBM nodes by the relative velocity magnitude normalised by the push plate displacement rate. In Fig. 12, the localised high liquid flow rate in response to the penetration is very clear among all deformation datasets, but more prominent liquid flow across the whole sample in Figs. 12c and d is observed. The liquid flow is induced by the drag (momentum exchange) from adjacent DEM grains. Also in Fig. 12c and d, there is a subtle liquid flow towards the specimen centre from the membrane boundary with relative velocity magnitude <0.3. This can be related to a pressure differential that draws liquid to regions undergoing high shear-induced dilation. This is also seen in the experiments as a local contraction (grains moving closer together) to the right and bottom-right of Fig. 3s and t.

Stress-deformation response
A well-calibrated set of LBM-DEM simulation results can offer extra information that is not available in X-ray imaging. Fig. 13a is the evolution of normal push-plate stress tracked to 10d push-plate displacement, using the sum of grain-plate normal contact forces divided by the length of the push plate front similar to reference [96]. The normal stress gradually increases in lower solid fraction simulations A and B, while a peak push-plate stress can be found at 2d displacement for higher solid fraction simulation C, and at 3d displacement for the highest solid fraction simulation D. It is noted that the predicted push-plate stress evolution cannot be compared with deformation experiments since the friction between the pushplate and Al 2 O 3 window compromises our ability to accurately measure the push-plate load. Still, both the trend of increased peak stress in higher solid fraction simulations and the strain softening phenomenon observed in simulation C and D are consistent with previous deformation tests on granular/particulate materials including semi-solid alloys [10,19,31,97], experimental soils [96,98,99], and DEM-simulated granular soils [100,101]. Fig. 13b shows the evolution of total volumetric strain, which was calculated from: where ε DEM vol is the volumetric strain directly evaluated from the change of area in the simulation samples, A and A def are the area of the initial and deformed sample, which were obtained from the area of the polygon that joins the centroids of all surface bonding grains minus the push-plate penetration area. Fig. 13b clearly shows the transition from net-contractive simulation sets A and B with low initial solid fraction, to net-dilative simulation sets C and D with higher initial solid fraction. The reduction of sample area in A and B occurs because the initially loosely packed solid network in a low-solid fraction LBM-DEM simulation is prone to develop a highly localised contractive field, resulting in the closure of liquid interstices. Simulation C, on the contrary, deforms with slightly net positive ε DEM vol behaviour because the contact force can be transmitted to the membrane to push the boundary outward and compensate for the local contractive response ahead of the push plate. Moreover, deformation in simulation D causes many of the liquid interstices to expand during grain rearrangement outside of the local compacted zone as well as the increased ability to push the enclosing membrane away, resulting in the net dilative deformation behaviour. Fig. 13b also shows the slight decrease of ε DEM vol after 6d penetration on simulation C and 7:5d penetration on simulation D. This is because the total solid fraction of the samples is decreased during net dilation, which increases the compressibility of the whole assembly.
The complex liquid behaviour at the free surface of semi-solid samples in the experiments can be related to the localised liquid pressure change Dp as shown in Fig. 11. We can then investigate the relationship between the change of volume of liquid-saturated samples and the overall liquid pressure changes. To get a representative trend, the mean liquid pressure change Dp was derived across the whole DEM assembly: (13) where N node is the number of LBM nodes inside the DEM assembly. Fig. 13c is the evolution of Dp for each simulation set and comparison with Fig. 13b shows the interplay between the global liquid pressure change and the volumetric strain, The net-contractive simulations A and B have a net increase of liquid pressure in response to net contraction. The high solid fraction simulations C and D have a slightly positive pressure change up to 1.5d deformation, followed by a decrease in liquid pressure in response to the global dilation, and the liquid pressure is near-constant once the peak dilatational volumetric strain is reached (Fig. 13b). It is useful to investigate the local stress evolution in a circular (d) Schematic illustration of the measurement regions for deriving local parameters. Diameter of measurement circles is 20d in order to obtain representative stresses, and there is a 1:5d gap between measurement circles and the push plate front. The horizontal position of the grey circle is in the middle of push plate, and the horizontal position of the blue circle is at the top-right corner of the push plate. The corresponding (e-f) q À p0, (geh) e À uy=d, and (iej) e À log p0 plots are shown for black and blue measurement circles. Arrows in (iej) indicate the initial state for each simulation. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) RV by explicit tracking of contact force vectors in the DEM system especially for the region in front of the push plate and the parting plane. Fig. 13d illustrates the setting of measurement circles in the DEM system to derive the local parameters such as the effective mean stress p 0 , deviatoric stress q, and void ratio e. In twodimensional DEM, the relationship between void ratio e and grain packing fraction (i.e. solid fraction) f S:pk is e ¼ ð1 À f S:pk Þ=f S:pk . These local parameters were tracked to examine whether the DEM system behaves as a granular material such as a sand [22]. The diameter of the measurement region was chosen to be 20d and includes~250 grains; this is sufficiently large to obtain the effectively averaged behaviour of the regions of interest. Fig. 13e and f illustrate the evolution of the q versus p 0 relationship to 10d displacement. Before shear deformation, the effective stress p 0 reflects the required mean stress on the DEM assembly to reach the desired packing fraction, and the mean stress values are similar between black and blue circles for every simulation set. The q-p 0 curves all start from q ¼ 0 and the stresses are increasing for simulation A and B. It is noted that the low solid fraction simulation sets (red and blue in Fig. 13e and f) soon have an approximately constant q= p 0 ratio and follow the regression line q ¼ 0:428p 0 , which can be related to the well-defined critical state in soil mechanics [22] stating that a frictional fluid should reach a macroscopic frictional constant M ¼ q= p 0 . In contrast, the maximum deviatoric stress occurs in simulations C and D at 0.98 Â 10 À4 N/m and 2.35 Â 10 À4 N/m for the black circle, and 0.64 Â 10 À4 N/m and 2.07 Â 10 À4 N/m in the blue circle, respectively. After reaching the maximum deviatoric stress, both p 0 and q stresses decrease. The simulation C then reaches the q ¼ 0:428p 0 regression line, while the q-p 0 curves in simulation D show a nearly À50% decrease of stresses from the peak deviatoric stress point and moves towards the regression line (but does not reach it by 10d push-plate displacement).
The evolution of void ratio is shown in Fig. 13g and h. The high contraction field in front of the push plate developed at an early stage of deformation (3d) reduces the void ratio by about À0.03 in the black circle and about À0.01 in the blue circle. The void ratio decrease corresponds to the increase of local mean effective stresses as shown in Fig. 13i and j. There is a log-linear relationship between the final void ratio and the final value of p 0 ; this corresponds to the critical state line concept accepted for granular soil behaviour [22]). The increase of void ratio in high solid fraction deformation is due to shear-induced dilation in both the black and blue circles.
Whilst this research used globular grains to prevent significant shape change during isothermal experiments and enable an efficient model, casting processes often create more complex morphologies such as equiaxed dendrites. Experiments [102] and DEM studies [64] have shown that dendrites form a solid network at lower solid fraction than globular grains due to the high liquid fraction within dendrite envelopes and due to the more complex shape of the envelopes. These factors will also alter the permeability [103], give more contacts per grain, and affect the local stresses acting on dendrites [64]. Future research can extend the LBM-DEM approach to dendritic microstructures in the future.

Conclusions
Time-resolved synchrotron X-ray radiography and coupled LBM-DEM simulations have been applied to investigate shear deformation mechanisms in equiaxed globular AleCu alloys at 44e71% solid. Microstructural and macroscopic responses to load have been quantified by intensity processing and digital image correlation (DIC) for radiographic imaging with >1000 grains in the field of view (FOV), and also reproduced by LBM-DEM after finding a single set of DEM parameters calibrated across all deformation experiments from 44 to 71% solid. From the results obtained by the various analysis approaches, the following conclusions can be drawn: The volumetric strain during semi-solid deformation due to local changes of solid fraction is strongly influenced by the initial solid fraction. A contractive strain field developed in shear or compressive deformation in alloys initially containing 44% and 50% solid. In contrast, significant net-dilation occurred during shear on the 66% and 71% solid samples. The mechanism of solid fraction increase during deformation is grains being pushed together by compaction, or by the liquid being sucked away. The solid fraction decrease, on the other hand, is due to shear-induced dilation by grains pushing one another apart during grain rearrangement in a solid network, similar to phenomena in compacted granular/particulate materials. Deformation was highly heterogeneous and the strain was strongly localised. Both the experiments and simulations agree with the trend of net-contraction ahead of the push plate and, in high solid fraction deformation, the co-location of dilative volumetric strain and shear strain in a region with width of 10e15 mean grains wide. Additionally, a contraction field formed outside the locally dilating regions in high solid fraction samples due to liquid being sucked away to feed dilation. The coupled LBM-DEM model naturally simulates the change of packing density in primary Al grains as an emergent phenomenon, and captures the complex liquid behaviour recorded in radiographs from the simulated liquid pressure change. The calibrated LBM-DEM model provided information on both the local stresses and void ratios; these data show that the material exhibits a load deformation response that can be described by the critical state framework for soil behaviour.