Modelling of paste ram extrusion subject to liquid phase migration and wall friction

Extrusion of solid-liquid particulate pastes is a well-established process in industry for continuously forming products of defined cross-sectional shape. At low extrusion velocities, the solids and liquid phases can separate due to drainage of liquid through the interparticle pores, termed liquid phase migration (LPM). The effect of wall friction, die shape and extrusion speed on LPM in a cylindrically axisymmetric ram extruder is investigated using a two-dimensional finite element model of paste extrusion based on soil mechanics principles (modified Cam-Clay). This extends the smooth walled model reported by Patel et al. (2007) to incorporate a simplified Tresca wall friction condition. Three die entry angles (90 , 60 and 45 ) and two extrusion speeds are considered. The extrusion pressure is predicted to increase with the Tresca friction factor and the extent of LPM is predicted to increase with decreasing ram speed (both as expected). The effects of wall friction on LPM are shown to be dictated by the die shape and ram displacement: there are few general rules relating extruder design and operating conditions to extent of LPM, so that finite element-based simulation is likely to be needed to predict the onset


Introduction
Particulate pastes are used widely to manufacture products such as agrochemicals, pharmaceuticals, ceramic parts, mortar and solder pastes using techniques including ram extrusion, screw extrusion and injection moulding (Wilson and Rough, 2012). These materials often feature a high volume fraction of particulate solids mixed with a liquid binder. Their rheology exhibits complex yield stress behaviour and hardening (Götz et al., 2002;Mascia and Wilson, 2008) as well as wall slip (Meeker et al., 2004). Some pastes are viscoplastic, reflecting the use of a highly viscous binder or (less frequently) a rate-dependent particulate matrix (Mascia and Wilson, 2007). Others, such as mortar pastes, reflect 'frictional' rheology more reminiscent of dry particulate assemblies (Perrot et al., 2006b). In the latter case, the rheology is dominated by friction at the interparticle contacts as opposed to viscous shear in the liquid binder. These frictional pastes are the focus of the current work.
Several types of flaw may develop during paste extrusion that are specific to their multiphase nature (Benbow and Bridgwater, 1993a). Liquid phase migration (LPM) is one such flaw. When a force is applied to the paste to move it through the die of an extruder, the stress (extrusion force/contact area) is split between the load-bearing particulate network (the solids skeleton) and the liquid binder present in the inter-particle pore network. During extrusion, a pore pressure gradient develops within the pore network and this promotes LPM relative to the solids skeleton. The rheology of the paste is highly dependent on the local liquid volume fraction, and LPM therefore changes the flow patterns in the extruder (Tomer and Newton, 1999;Chen et al., 2000). The product (extrudate) exhibits a time-varying liquid fraction as a result of LPM, which is highly detrimental as the extrudate liquid volume fraction (ELVF) is usually a key product specification. An additional concern is damage to the extruder if LPM becomes excessive and the pressure required to extrude the compacted paste left in the barrel becomes excessive. However, a minority of paste processes exist in which LPM is intended, e.g. sugar cane juice extraction (Loughran and Kannapiran, 2005).
Detailed and reliable simulation of paste extrusion processes is required to improve the design of dies and to limit flaws such as LPM. The formulation of the paste (particle size and size distribution, solids volume fraction, binder rheology, etc.) is highly relevant to the extent of LPM as this largely decides the permeability of the solids skeleton. The formulation also determines the rheology of the material and its interaction with the walls of the forming device. As a result, formulation design is complex and incorporates a large number of independent variables that are best investigated by separate, dedicated studies (Blackburn and Bohm, 1994). This is actively being pursued in applications such as bone cements (O'Neill et al., 2016), where LPM causes problems in the delivery of such materials when they are injected from a syringe through a catheter during surgery. The aim of the current work is to provide a quantitative framework for the phenomenon that can be used to interpret results from laboratory studies and then predict the performance of such pastes in different configurations. To this extent,

Nomenclature
Roman 1-D one-dimensional 2-D two-dimensional CED conical entry die geometry CED45 45°conical entry die geometry CED60 60°conical entry die geometry d particle diameter [ the influences of ram speed, die entry angle, and the magnitude of wall friction on the extent of LPM are investigated here. Research into LPM during paste flow has tended to focus on the simpler, 2-D geometry of upsetting/squeeze flow (Poitou and Racineux, 2001;Sherwood, 2002;Kolenda et al., 2003;Roussel and Lanos, 2004). Here, the interfacial condition at the plates approaches the limiting cases of perfect lubrication (frictionless slip) or full sticking, for which fully or partially analytical solutions are available (Steffe, 1996). Ram extrusion does not permit onedimensional analysis as the flow contracts significantly as it approaches the die entry, adding significant radial flow and extensional deformation. However, almost all previous numerical investigations of LPM during ram extrusion have been conducted in 1D (Wroth and Houlsby, 1983;Rough et al., 2002;Mascia et al., 2006;Perrot et al., 2006aPerrot et al., , 2009Khelifi et al., 2013) due to the numerical complexity of simulating extrusion in higher dimensions. One exception is our prior study (Patel et al., 2007), which described a model for ram extrusion in concentric cylindrical systems with smooth walls. In this case, the coaxiality of the barrel and die allow the circumferential coordinate to be neglected, i.e. the model is axisymmetric. That work did not consider the influence of friction between the paste and the extruder wall, for simplicity, and the approach is extended here to include a simple (Tresca) wall friction model. The effects of friction on process variables such as the extrusion pressure and the extent of LPM are highlighted. This is the first time to the authors' knowledge that wall friction has been incorporated into two-dimensional simulations of ram extrusion alongside an analysis of LPM.
There are three challenges in modelling these systems. The first lies in how to describe the rheology of the paste, as the behaviour lies somewhere between classical fluid approaches (based on strain rate) and solids (based on strain). Much of the work on modelling LPM has employed fluid constitutive relationships, which are more appropriate for 'softer' materials. The work reported here employs a plasticity approach, which is arguably more appropriate for stiffer materials, and builds on the work of Horrobin and Nedderman (1998). The second is the treatment of wall friction, i.e. tribology, and this is not fully understood for dense particulate systems. A noteworthy addition to the field is the I-rheology approach for dense particulate systems by Gray and Edwards (2014). Incorporating such methods into fluid mechanics coding is challenging. For example, Bryan et al. (2015) demonstrated that including non-linear wall slip into computational fluid dynamics (CFD) simulations of viscoplastic fluids can lead to computational difficulties and the results are very sensitive to the values of the parameters used, as also reported by Ardakani et al. (2013). Finally, liquid flow through a matrix undergoing shear is likely to differ from the stationary case owing to the evolution of pore shape and size. In the absence of experimental data on this topic, classical Darcy's law approaches are used to describe the local liquid flux.

Numerical model
A detailed description of the model is given by Patel et al. (2007). For brevity, only essential features and new developments (wall friction condition; Section 2.4.1) are presented here.

Model assumptions
Key assumptions are that: i. the process is isothermal, ii. gravitational and inertial forces are negligible, iii. the individual particles are incompressible, iv. the liquid binder is Newtonian, incompressible and saturates the pores, i.e. there is no entrained air such that Darcy's law may be used to model LPM.

Domain behaviour -solids skeleton
Following Terzaghi's principle (Terzaghi, 1936), normal stresses applied to the paste are opposed by the sum of the effective stress in the solids skeleton and the pore liquid pressure, while shear stresses are sustained by the solids skeleton only. The effective stress state is represented by the tensor r ij , which is modelled as symmetric due to moment equilibrium. r ij can be resolved into deviatoric (s ij ) and pressure (p d ij ) contributions: As explained in our prior study (Patel et al., 2007), the modified Cam-Clay (MCC) constitutive model (Schofield and Wroth, 1968) was chosen to relate the effective stress to the strain in the solids skeleton.

Volumetric constitutive behaviour
Under the action of an increasing effective pressure (p), MCC materials reduce in volume (compact) elastically via Eq. (4) and then plastically via Eq. (5): The empirical lumped parameters j and k are the logarithmic elastic and logarithmic plastic bulk modulus, respectively. The term p 0 is a reference pressure that is set at 1 atm. For MCC, elastic compaction is assumed to reflect reversible deformation at the interparticle contact points, while plastic compaction is assumed to reflect particle re-arrangement and crushing. However, the material from which the particles are made does not change in volume, i.e. the particles themselves are incompressible. The voids ratio, e, is defined as: e ¼ local volume of voids local volume of solids A boundary condition is required for Eqs. (4) and (5). This is the (0, e 0 ) state in Fig. 1, which reflects a solids skeleton that has been compacted to a pressure no larger than p 0 during its stress history.
The movement of the stress state along the j-lines and the k-line, and the physical interpretation of these stress paths, are described in the prior study (Patel et al., 2007).
The ratio of the current effective pressure (p) to the effective pressure at yield during hydrostatic compression, i.e. the pressure at which the j-line intersects the k-line, is termed the overconsolidation ratio (OCR). This is described further in Section 2.2.2.

Deviatoric constitutive behaviour
Additive strain rate decomposition is assumed in this study, as in all prior modelling studies of paste extrusion. The elastic strain rate in the solids skeleton is thus assumed small relative to the plastic strain rate. The total strain rate then simplifies to the sum of the two components (ABAQUS, 2007a). This simplification eases modelling considerably and is reasonable for extrusion processes in which the plastic strain is dominant. A small increment in the deviatoric stress, ds ij , drives an increment in the (engineering) elastic deviatoric strain, dc e : Here G and E are the local shear modulus and Young's modulus, respectively. The Poisson ratio of the solids skeleton (fully drained limit) is represented by m and is assumed constant for simplicity.
When the stress state reaches the yield surface, plastic flow is initiated. The MCC yield surface is illustrated in the meridional plane ( r vs. p), the deviatoric p-plane and in principal stress space in Fig. 2.
The surface shrinks (softening) and grows (hardening) when the plastic strain increment incorporates positive volumetric strain (dilation) and negative volumetric strain (compaction), respectively. The plastic strain increment may also incorporate deviatoric strain and the ratio of deviatoric to volumetric strain is assumed to follow the associated flow rule. As a result, the plastic strain increment can be plotted on the same axes as the yield surface and is normal to it. Regardless of the extent of hardening, all MCC yield surfaces intersect at the origin of Fig. 2(a), and thus all predict zero yield stress at zero effective pressure (the solids skeleton is cohesionless). The points on the yield surfaces at which plastic deformation results in purely deviatoric strain (zero volumetric strain; see top of semicircle in Fig. 2(a)) all lie on a common line termed the critical state line. The gradient of this line, M, is modelled here as constant (unity).
As stated in Section 2.2.1, the ratio of the effective pressure at yield under purely hydrostatic loading (p ⁄ in Fig. 2) to the current effective pressure is termed the overconsolidation ratio. It is clear from Fig. 2 that (for M = 1) the OCR exceeds 2 (p/p ⁄ < 0.5) during softening, is equal to 2 (p/p ⁄ = 0.5) during critical state deformation, and is smaller than 2 (p/p ⁄ > 0.5) during hardening.
In summary, the parameters required to implement the simplest form of the MCC constitutive model are j, k, e 0 , m and M.
The values used here are those used in our prior study (Patel et al., 2007) and are presented in Table 1.

Domain behaviour -liquid phase
LPM is described using Darcy's law, i.e. the binder undergoes laminar, incompressible flow at a superficial velocity (U) that is directly proportional to the local gradient in pore liquid pressure Cross section through the MCC yield surface at positive p. All cross sections describe a circle (von Mises yield criterion). As a result, the undrained uniaxial yield stress, σ 0 , is equal to the shear yield stress, τ 0 , multiplied by √3.  3)). This plane is scaled by p (3/2) ( r-axis) and p 3 (p-axis) relative to the principal stress space used to construct (b). The material is modelled as cohesionless such that the yield surface is undefined at negative p. The term p * denotes the maximum p reached in the materials stress history. A, B and C represent dilatant, zero volume strain (critical state) and compactive plastic strain increments, respectively. These drive the yield surface to shrink in diameter along the p-axis (while still passing through the origin), to remain unchanged and to expand along the p-axis while still passing through the origin, respectively. The directions of the strain increment vectors are perpendicular to the yield surface as prescribed by the associated flow rule, although the different scaling of the axes in (a) and (b) requires multiplicative constants to be incorporated. dc p is the increment in engineering plastic deviatoric strain and must be halved to give the increment in plastic deviatoric strain. de v p is the increment in plastic volumetric strain. of elastic compaction and the sole plastic compaction path, respectively. The voids ratio e 0 is exhibited by a virgin sample of the solids skeleton at the reference effective pressure, p 0 . The states (e 1 , p 1 ) and (e 2 , p 2 ) denote two further states at which the material is at incipient yield.
(rP). The constant of proportionality is the local absolute permeability of the solids skeleton (K) divided by the binder viscosity (l). K is estimated using the Carman-Kozeny (CK) model for permeability: The term k t denotes the local tortuosity of the pores and is defined as the ratio of the actual length of flow channels to the straight-line channel length. The (common) value of 5 has been used (Kay and Nedderman, 1985). The viscosity of the liquid binder is 300 Pa s, which reflects the properties of the binder used in related experimental investigations (Bradley et al., 2004;Patel, 2008). S is the wetted surface area per unit volume of the solids skeleton, which is defined by the CK model for the case of spherical particles of uniform size (d = 175 mm; Eq. (11)).
The continuity equation for the liquid binder is given by Eq. (12), and relates U to the volumetric strain rate in the solids skeleton (de v /dt): where s s is the shear yield stress at the wall (the slip stress), m T is the Tresca friction factor and s 0 is the local shear yield stress of the solids skeleton. The friction factor is modelled as constant for simplicity and values of 0.33 and 0.67 are used here, spanning the range of likely values. A value of zero (frictionless slip) was used by Patel et al. (2007) and their results therefore constitute a parallel dataset. A value of unity indicates the extruder wall is perfectly rough: this value was not used as the balance of experimental evidence points to some wall slip (Wildman et al., 1999;West et al., 2002;Barnes et al., 2004).
Eq. (13) requires the value of s 0 at the wall. For MCC materials, s 0 is a function of the local values of both p and p ⁄ ( Fig. 2(a)). Thus, a complex series of calculations is required to obtain s s . For simplicity, a constant slip stress is used here and is based on the initial values of p and p ⁄ . This approximation is accurate only in the limit of zero LPM (high ram speed) as LPM induces changes in p and p ⁄ . The validity of this assumption is discussed alongside the results in Section 3.

Extrudate surface condition
The portion of the extrudate surface not in contact with the ram or the wall (the extrudate free surface) is subject to capillary pressure for which the following model is used. The capillary suction force, F, on a single cylindrical pore of diameter d p that is caused by surface tension between the solid and liquid phases, s t , is a standard result and is given by, e.g. Douglas et al. (2001, p. 15): In Eq. (14), h c is the solid-liquid-air contact angle: s t and h c are reported in Table 1. As with the CK permeability model, the pore exits at the exposed surface are assumed to be cylindrical and to feature a uniform diameter. S (see Eq. (11)) is assumed to give the wetted pore perimeter per unit area at the extrudate free surface as well. The local capillary pressure can be estimated by multiplying the capillary force per unit wetted pore perimeter, F/p d p , by S. For the data in Table 1, this yields an initial capillary pressure of À1.2 kPa. As this is the only stress before extrusion commences, this is therefore the initial pore pressure and the initial (hydrostatic) effective stress in the solids skeleton.
If the pore pressure at the extrudate free surface exceeds the capillary pressure, the liquid seeps from the surface until the pore pressure falls below this limit. This boundary condition is available within the finite element solver used here (ABAQUS, 2007b) and has been observed for sewage pastes (Chevalier et al., 1997;Chaari et al., 2003) and a food paste (Cheyne et al., 2005). If the pore pressure falls below the capillary pressure, the difference in pressure is transferred to the solids skeleton as a supplementary (compressive) effective stress.
This capillary pressure condition is slightly different to that used previously (Patel et al., 2007). There, atmospheric pressure was incorrectly imposed at the free surface in addition to capillary pressure. The proportion sustained by each phase was equal to the local volume fraction of the phase (initial / s = 0.5). This is incorrect Table 1 Values of material constitutive parameters used in all simulations -for definitions, see Nomenclature. This table is reproduced with permission from the authors' prior study (Patel et al., 2007).

lm
Reflects the near-monodisperse glass ballotini used in later experimental studies by the authors (Bradley et al., 2004;Patel, 2008) l 300 Pa s Reflects the high viscosity liquid binder used in later experimental studies by the authors (Bradley et al., 2004;Patel, 2008), which prevented sedimentation of the ballotini during extrusion experiments s t 72.3 mN m À1 Viscous aqueous liquid -for simplicity, the value for pure water at 23°C was used (Lide and Frederikse, 1995, Section 6-8) h c 0°Viscous aqueous liquid -perfect wetting assumed for simplicity and the extrusion pressures are then 1 atm too large. This was corrected in those results presented here.

Initial conditions
The initial voids ratio e initial = 1 (initial / s = 0.5). The initial effective stress and pore pressure are À1.2 kPa and 1.2 kPa, respectively (Section 2.4.2). The initial overconsolidation ratio is set at 2 (Section 2.2.2).
2.5. Finite element model 2.5.1. Basic construction All simulations were run using the ABAQUS/Standard v6.7-1 implicit finite element solver, on the Windows desktop PC used in the prior study (Patel et al., 2007). Fig. 3 presents the model geometry and a sample initial mesh. The mesh follows the material in all simulations, i.e. the simulations are Lagrangian, and contains $800 quadrilateral elements of type 'CAX4P', which are first order with respect to solids displacement and pore pressure: this choice is explained in the prior study.
The nodes adjacent to the notional ram move at a fixed speed V towards the die face. The nodes at the axis cannot move radially and the nodes in contact with the extruder wall cannot penetrate it. Finally, a length of extrudate, l cut , is removed from the mesh when its length reaches a value of l cut + l tol (clearance length; Table 2). The average extrudate liquid volume fraction in these cut sections (ELVF) is logged as described in the prior study.

Adaptive remeshing
The elements in the die entry become highly distorted during extrusion. Adaptive remeshing is employed to overcome this problem; specifically, a pre-existing code (Horrobin, 1999, pp. 90-101;Patel et al., 2007) is used to pause the ABAQUS simulation periodically, draw a new mesh within the confines of the old, distorted mesh, map the distributions of solution variables from the old mesh to the new mesh, and then restart the simulation. The term 'adaptive' refers to each new mesh, which features small elements wherever p varies strongly in the old mesh.
The mapping process involves (i) averaging the values of the variables in the old mesh to the nodes of the old mesh (performed by ABAQUS), (ii) interpolating the mapped variables from the nodes in the old mesh to the nodes in the new mesh (for pore pressure) and the integration points in the new mesh (for effective stress and voids ratio). Interpolation was performed using the standard bilinear shape functions for first order quadrilateral finite elements. The mass of each phase (liquid binder and particulate solids) was calculated before and after the mapping process. This was found to change by ±10 À5 (relative) after each mapping stage. The voids ratios throughout were adjusted to negate this change before recommencing the simulations as described in the prior study (Patel et al., 2007). Finally, the simulations were run on the Fig. 3. Geometry of the axisymmetric ram extrusion model and a sample initial mesh. R b and R d are the barrel and die land radii, respectively, and R d /R b (20%) is termed the reduction ratio. L is the axial distance between the barrel-die face corner and the die land exit corner, not the axial distance between the die land entry corner and the die land exit corner. The extruder dimensions are listed in Table 2, and match those used in the authors' prior study (Patel et al., 2007). The initial height of the paste billet h 0 = 3 R b for all simulations. 'A' denotes the extrudate free surface, 'B' denotes the portion of the mesh surface subject to friction -node C is not included to avoid a local numerical overconstraint (Horrobin, 1999), and D denotes the ram nodes. same PC described in the prior study and the runtimes were similar to those quoted previously.

Results and discussion
Simulations considered three die entry angles: 90°(square entry die geometry; 'SED'), 60°(conical entry die geometry; 'CED60') and 45°('CED45'), for a reduction ratio, R b /R d , of 20%. The ram nodes moved axially towards the die face at two (scaled) velocities, specifically V/R b = 0.1 s À1 and 0.002 s À1 . The values of Tresca friction factor tested were taken to be 0.33 and 0.67, giving a total of twelve simulations. Patel et al. (2007) reported parallel results for the smooth-walled case (m T = 0) and these are referred to here for comparison. Data are presented in three forms: as plots of the evolution of the extrusion, pore and solids pressure with ram displacement, termed extrusion profiles; plots of the liquid content of the extrudate (ELVF) versus ram displacement; and as contour plots of voids ratio. The extrusion pressure (P e vs. z) profiles show an initial increase as the billet enters and fills the die region. LPM is manifested thereafter as variations in extrusion pressure and ELVF. LPM was only observed at the smaller ram speed, as was observed in the absence of wall friction. The largest change in voids ratio from the initial value (e initial = 1) was ±8% (D/ s = ±2%).

Smooth walls
The main findings of the Patel et al. (2007) study with smooth walls are summarised here. At the higher ram speed (V/R b = 0.1s À1 ), the frictionless simulations are essentially LPM-free. As a result, the extrusion pressure after the initial transient is nearly constant and the extrudate liquid volume fraction (ELVF) remains close to its initial value (0.5). The extrusion pressures are summarised in Table 3 and show a slight reduction in P e for the conical dies compared to the SED. Very slow/static zones are not observed at the barrel-die face corner for any die shape. P e is smaller at the lower ram speed for all die entry angles tested. This is due to either a reduction in the pore pressure at the ram (P eL ) or the effective stress in the solids skeleton (P eS ), or both. The former is caused by an increase in permeability in the die entry zone due to dilation, which allows for more rapid dissipation of the excess P eL developed during extrusion. The latter is due to the softening (weakening) of the solids skeleton in the die entry zone. The formation of a compacted static zone at the barrel-die face corner is promoted at the lower ram speed for SED geometry.
The ELVF does increase with displacement (indicative of greater LPM) at the lower ram speed for all die shapes. The ELVF data, however, show little variation with h d .
The above results reflecting the extent of LPM may be compared with the threshold criteria described by Wroth and Houlsby (1983) and Perrot et al. (2009). These two criteria essentially reflect ratios of the estimated timescale of LPM to the timescale of extrusion (t process ), i.e. the time taken for an element of paste to traverse the die region. That of Wroth and Houlsby (1983) yields a minimum value for t process that restricts compaction of the paste billet to 5% of the maximum possible: In Eq. (15), K denotes the permeability of the paste (Eq. (10)) and p c denotes the effective pressure of the solids skeleton at t = 0 if it were to shear at critical state. The latter is given by the initial value of p (0.5 atm + 1.2/101.325 atm; Section 2.4.3) multiplied by the initial value of the OCR (2), i.e. 1.02 atm. The value of t process is given by the ram speed (V) divided by the eventual ram displacement (1 R b ; Table 2 Extruder geometry and initial conditions for ram extrusion simulations. These data match those used in the authors' prior study (Patel et al., 2007), except for the shear yield stress at the wall which is assumed to be non-zero in this study.  Table 3 Total extrusion pressure, P e = P eL + P eS predicted when extrudate first becomes visible at the die land exit. This value of ram displacement corresponds to the beginning of the final segment of the P e profiles when the die entry zone and the die land are just full; see Section 3.1. Results at m T = 0 were first presented in the authors' prior study (Patel et al., 2007).  Table 2). The logarithmic plastic bulk modulus k = 0.05 and the reference voids ratio used by the k-line (Fig. 1), e 0 , is 0.978 (Table 1). L process denotes the axial distance between the ram and the die exit. The values of K and L process are set at their initial values (8.5 Â 10 À11 m 2 & 4 R b , respectively) for simplicity, and the values of t process are 10 s and 500 s at the highest and lowest ram speeds simulated, respectively. The LHS values of Eq. (15) are then 0.0012 and 0.06 at the highest and lowest ram speeds, respectively. Eq. (15) therefore predicts negligible LPM at the higher ram speed and compaction of the paste in the barrel at the lower ram speed by approximately 5% Â 0.06/0.1 = 3% of ultimate. The corresponding prediction of the finite element model at the lower ram speed is a 4% decrease in voids ratio near the ram for SED geometry (as was modelled by Wroth and Houlsby (1983)), which represents good agreement given the assumptions inherent in the criterion. Perrot et al. (2009) describe two further criteria for t process for the case of ram extrusion of a saturated, highly frictional mortar paste -the simpler criterion is used here as it assumes the properties of the paste are essentially constant (initial values), and the criterion can therefore be evaluated analytically rather than numerically. Strictly speaking, these criteria are not applicable to the (frictionless wall) results described in this section as Perrot et al. (2009) assume the extruder wall to be perfectly rough, i.e. the slip stress at the wall (s s ) is equal to the bulk shear yield stress. However, there is a lack of simplified criteria for paste formulators to use to avoid LPM during ram extrusion, and we suggest therefore that there is still value in testing this criterion with the current dataset. Thus, we incorporate (a posteriori) the assumption that the extruder wall is perfectly rough. Perrot et al. (2009) proposed that LPM becomes significant when the timescale of the process approaches the time required for the paste at the ram to compact (due to LPM) sufficiently to double its shear yield stress (t 2s0 ): In Eq. (16), De represents the reduction in voids ratio at the ram. K, h 0 and s 0 denote the initial permeability of the paste (8.5 Â 10 À11 m 2 ; Eq. (10)), the initial height of the billet (h 0 = 3 R b = 38 mm; Table 2) and the initial shear yield stress in the paste, respectively. For the von Mises yield criterion ( Fig. 2 (b)), the value of s 0 at undrained conditions (assumed here for simplicity) is given by r 0 / p 3 = 30 kPa cf. 20 kPa for the mortar pastes studied by Perrot et al. (2009). The values of n ext are therefore 0.012 and 0.6 at the higher and lower ram speeds, respectively. These values compare well with the results of the finite element model for SED geometry (as used by Perrot et al. (2009)), which were essentially zero LPM at the higher ram speed and a 4% reduction in voids ratio at the ram at the lower ram speed, which corresponds to a (57-30) kPa/30 kPa = 90% increase in the maximum, i.e.
critical state value of s 0 (calculation omitted for brevity).
The above results demonstrate that the 2-D finite element model described here and in our prior study produces results broadly consistent with two independent, lower-order models of LPM during ram extrusion. This implies that our finite element model is constructed appropriately. We now present results obtained in the presence of wall friction. Fig. 4 presents extrusion pressure profiles (P e and P eS ) at the higher ram speed (V/R b = 0.1 s À1 ). The ram displacement is dimensionless, i.e. z ⁄ z ram /R b , where z ram denotes the instantaneous ram displacement. The P eL profile (i.e. P e À P eS ) can be estimated visually from the P e and P eS profiles. The corresponding ELVF profiles are presented in Fig. 5.

Wall friction, high ram speed
The difference in the initial profile shape between SED and CEDs arises from the filling of the conical entry region, which is not present in the former and is longer for CED45. The maximum P e also increases with h d : this is partly due to a design decision made when drawing the extruder within ABAQUS. The extruder dimension 'L' in Fig. 3, i.e. the axial distance between the barrel-die face corner and the die land exit, is identical for all three h d . This results in an increase in the die land length, and thus the contact area between the paste and the die land wall, with h d . The total contact area between the paste and the extruder is then essentially unaffected by h d , i.e. the area of the die face decreases and that of the die land increases as h d increases and these cancel out nearperfectly. It is known that P e increases significantly with the contact area in the die land (Benbow and Bridgwater, 1993b), and the increase in maximum P e with h d predicted at V/R b = 0.1 s À1 (Table 3) is therefore physically consistent. Table 3 also shows that P e increased with m T (as expected) for all three geometries in an almost linear fashion.
The six simulations at the high ram speed featured less LPM than observed at the low ram speed, as reported for a smooth wall. This is evident from (i) the proximity of the corresponding ELVF profiles to 0.5, (ii) the uniformity in the final voids ratio distributions (see Section 3.3), and (iii) the linearity of the P e profiles after the initial transients. The extrusion behaviour of a two-phase soillike paste in the presence of wall friction, and in the absence of LPM, therefore reduces to that of an elastic-perfectly plastic solid. This was also reported for a frictionless extruder wall (Patel et al., 2007). Therefore, the recommendation that extrusion should be conducted at high ram speed to avoid LPM also holds when wall friction is present.
The shapes of the P e profiles are now discussed. Fig. 4 presents four P e profiles for CED geometry at the high ram speed. Each profile features four distinct segments, which are separated by the dashed, vertical lines in Fig. 4(a-b)(ii-iii). For all four cases, the first segment of the profiles (lowest z ⁄ ) is linear. This region relates the increase in P e due to the increasing mobilisation of friction (increasing s w ) at the barrel wall as z ⁄ increases. Supporting evidence for this interpretation is that the gradient, dP e /dz ⁄ , in the first segment is similar at h d = 45°and 60°.
The second segment is curved and corresponds to the filling of the conical die entry zone. The area of the free surface of the paste decreases with increasing z at higher-than-linear order. By conservation of mass, the velocity of the free surface of the paste during filling of the die entry zone increases with ram displacement at similar order. Therefore, the average extensional strain rate in the paste increases with ram displacement at similar order, giving the second segment of the P e profiles its curved shape.
The third segment of the P e profiles is narrow and linear. This describes filling of the die land, i.e. P e increases with z ⁄ due to the increase in the total frictional force mobilised at the die land wall. The range of z ⁄ incorporated by the third region is narrower at h d = 45°that at 60°as the die land length decreases with decreasing h d . The value of dP e /dz ⁄ at each m T is the same for both values of h d , which is physically consistent with the conditions of (i) limited LPM, (ii) fixed wall shear yield stress (same m T ), and (iii) common diameter reduction ratio (R b /R d ).
The beginning of the fourth segment corresponds to the emergence of extrudate from the die exit, i.e. the onset of visible flow. This segment features a mild linear decrease in P e for all the high ram speed cases in Fig. 4. This occurs as the paste-wall interfacial contact area in the barrel (A bw ) decreases linearly with increasing ram displacement. The gradient of these profiles is now estimated. If slip is assumed to occur along the entire barrel wall, the frictional force at the barrel wall is s s A bw . As s s is fixed in these simulations, dP e /dz ⁄ in the fourth region may be estimated from dP e dz Ã % À This result is independent of h d , and gives À0.2 atm R b À1 and À0.4 atm R b À1 at m T = 0.33 and 0.67, respectively. At the high ram speed (little LPM) and m T = 0.33, dP e /dz ⁄ in the fourth region is À0.12, À0.32, and À0.39 atm R b À1 at h d = 90°, 60°and 45°, respectively. The corresponding predictions at m T = 0.67 are À0.36, À0.30, and À0.44 atm R b

À1
. The agreement between the two sets of values is poor. This occurs because friction at the barrel wall strongly affects the flow field in the extruder, as can be seen from the distributions of the dimensionless nodal speed, V n in Fig. 6, where V n is given by: In Eq. (18), V r and V z are the radial and axial velocities of the solids skeleton, respectively. For the SED, near-static zones are evident at the barrel-die face corner at the high ram speed; see Fig. 6(a)(i-ii). Thus, slip does not occur over the entire barrel wall and Eq. (17) overpredicts the magnitude of dP e /dz ⁄ (|dP e /dz ⁄ |). At h d = 60°and m T = 0.67, |dP e /dz ⁄ | at the high ram speed is again overpredicted by Eq. (17) due to a static zone covering part of the barrel wall (see Fig. 6(b)(ii)). By contrast, |dP e /dz ⁄ | is underpredicted by Eq.
(17) at m T = 0.33. This is partly because no static zone is predicted at m T = 0.33 (Fig. 5(b)(i)), and partly because more dilation occurs in the die entry zone at the smaller m T (Fig. 6). As explained previously (Patel et al., 2007), dilation in the die entry zone causes a decrease in both P eL and P eS with increasing ram displacement. The former effect is due to the increase in permeability in the die entry zone caused by dilation, which allows excess pore pressure (high P eL ) to dissipate more rapidly. The latter effect is due to dilation (softening) in the die entry zone, which reduces the stress carried by the solids skeleton (P eS ). At h d = 45°, there is little LPM and static zones do not develop at the barreldie face corner at the high ram speed for either m T ; see Fig. 6(c) (i-ii). Therefore, Eq. (17) should provide reliable predictions of |dP e /dz ⁄ |. At m T = 0.33, the prediction (À0.4 atm R b À1 ) matches the gradient of the final segment of the P e profile very well (À0.39 atm R b À1 ). A small amount of dilation in the die entry zone occurs at m T = 0.67, which explains why dP e /dz ⁄ (À0.44 atm R b

À1
) is larger in magnitude than the value predicted by Eq. (17).
Comparing P e with P eL in Fig. 4 shows that the pore pressure is the dominant contribution to the extrusion pressure for all six cases run at the higher ram speed. This matches the result for the smooth wall and indicates that ram extrusion through a die of moderate reduction ratio is similar to one-dimensional confined compaction: high values of P eL are therefore consistent and reason-able. For the SED geometry, m T has little influence on P eS and is in fact similar to the corresponding profile for frictionless walls (Patel et al., 2007). The predominant effect of wall friction is to increase P e via an increase in P eL . This occurs because P eS is controlled by the voids ratio distribution, which is essentially constant at the high ram speed for the m T values tested; see Fig. 7(a)(i-ii) and Patel et al. (2007). Similarly, the P eS profiles and voids ratio distributions for the CED case at the high ram speed do not vary significantly with m T ; see Fig. 7 Fig. 7(c)(i-ii) and Patel et al. (2007). These results are consistent with the ELVF profiles in Fig. 5, which indicate little LPM in all cases at the higher ram speed. The V n distributions in Fig. 6 demonstrate that at high ram speed, the increase in m T from 0.33 to 0.67 for h d = 90° (Fig. 6(a) (i-ii)) and h d = 60° (Fig. 6(b)(i-ii)) promotes the formation of a static zone at the barrel-die face corner. For h d = 45°, the paste near the barrel-die face corner slows as m T is increased (Fig. 6(c)(i-ii)) but does not become static. Sticking at the die land wall did not occur for any case at the high ram speed.
The voids ratio distributions presented in Fig. 7 demonstrate that at the high ram speed, there is a slight increase in LPM (greater variation in voids ratio) with m T for SED geometry (Fig. 7(a)(i-ii)). This pattern is reversed at h d = 60° (Fig. 7(b)(i-ii)), and little change occurs in voids ratio for h d = 45° (Fig. 7(c)(i-ii)). Similarly inconsistent variations with m T occur in the die land. However, these variations are all small (<0.5% of e initial ) and are likely to lie within the resolution of the numerical method. Therefore, these results are not analysed further.
To summarise, the extent of wall friction has a significant effect on the flow field at the high ram speed for all die shapes. This is reflected by (small) changes in extrusion pressure as the ram moves towards the die face which, as demonstrated here for the first time, cannot be explained other than via the use of twodimensional, two-phase simulations. However, the overall extent of LPM at the high ram speed remains small for all cases (nearundrained conditions), and is reflected by the nearly constant extrusion pressure and ELVF, as well as almost uniform voids ratio distributions at z ⁄ = 1 (close to e initial ).

Wall friction, low ram speed
The P e profiles at a given die entry angle in Fig. 4 exhibit the same basic shape regardless of ram speed or friction factor. The P e profiles again comprise two (SED) or four (CED) segments. However, for all three h d values there are noticeable differences in P e in the final segment, which begins with the peak P e in five of the six cases. The peak P e decreases with die entry angle for m T = 0 and 0.33 (as at the high ram speed), but shows no overall trend with h d at m T = 0.67 due to the coupled influences of wall friction and LPM on the components of P e at the low ram speed.
Inspection of the individual components of P e at h d = 90°and 60°s hows that P e is smaller at the low ram speed at all m T : this is due to a decrease in P eL with ram speed that occurs due to dilation in the die entry zone (Section 3.2). At h d = 45°, P e at the low ram speed exceeds its value at the high ram speed over two ranges of ram displacement, namely 0.7 < z ⁄ < 1 at m T = 0.33 and z ⁄ < 0.32 at m T = 0.67; see Fig. 4(c). The first range occurs due to P eS being much larger at the low ram speed, which more than negates the accompanying decrease in P eL . This increase in P eS is (indirectly) due to the dilation in the die entry zone, which promotes compaction (hardening) of the solids skeleton in the remainder of the barrel, thereby increasing P eS . The second range is also due to hardening of the solids skeleton in (most of) the barrel. However, the decrease in P eL with ram speed eventually exceeds the increase in P eS at z ⁄ ! 0.32. These results demonstrate that LPM causes significant variation in P e via competing effects on P eL and P eS , and that these changes are coupled to the die shape, friction factor and current ram displacement. Prediction of these effects is therefore beyond the capabilities of modelling tools that do not consider the twophase behaviour of the paste and the precise shape of the paste billet.
P eL is again the dominant component in P e for all cases, although the time-mean P eL and P eS increase and decrease with increasing die entry angle, respectively; see Fig. 4. The former result (P eL /P e ! 1) also occurs at the high ram speed and the same explanation is proposed, i.e. that ram extrusion resembles one dimensional confined compaction for which P eL is large. The decrease in P eL with die entry angle (at fixed m T ) is consistent with this analogy as the resemblance between processes decreases with die entry angle. The increase in P eS with h d then reflects how the force imposed by the ram is increasingly transferred onto the solids skeleton (as effective stress) instead of onto the liquid binder. This variation of mean P eL and P eS with die shape increases in magnitude with m T (see Fig. 4). This is consistent with an increasing proportion of the force applied at the ram being sustained as shear (effective) stress along the whole length of the extruder wall (i.e. wall friction) rather than the sum of normal effective stress and pore pressure at the die plate, which is largest for SED geometry as the die plate is parallel to the ram. These results imply that in the presence of moderate-to-substantial wall friction, the SED geometry is superior to a conical one for minimising LPM as the solids skeleton throughout the upper barrel (above the die entry zone) is less prone to compaction due to shear stress at the barrel wall. This is supported by the voids ratio distributions at the end of extrusion, which show dilation being progressively localised away from the barrel as die entry angle decreases ( Fig. 7(a-c)(iii-iv)).
With regards to its variation with ram speed, P eS is larger at the low ram speed for all h d at m T = 0.67 due to compaction in the barrel. At m T = 0.33, P eS varies with ram speed in a less simple manner. For SED geometry, P eS is smaller at the low ram speed for z ⁄ < 0.05 and z ⁄ > 0.8. At h d = 60°, P eS is smaller at the low ram speed for intermediate ram displacements (0.64 < z ⁄ < 0.9). At h d = 45°, P eS is larger at the low ram speed for all z ⁄ . Patel et al. (2007) reported similarly complex behaviour in the smooth walled case. Thus, ram speed has competing influences on P eS that are coupled to die shape, wall friction and the extent of extrusion. Further elucidation of this feature requires (i) further simulations, (ii) extending simulations to higher z ⁄ to see if clearer trends emerge, and (iii) replacement of the simplified friction condition used here (constant s s ) with a condition allowing s s to reflect the instantaneous shear yield stress (s 0 ) at the wall.
With regards to ELVF, Fig. 5 shows little influence of m T on ELVF at the low ram speed, as reported for m T = 0 by Patel et al. (2007) and noted for the high ram speed. This appears counterintuitive as the paste in the die entry zone, which eventually exits the die land as extrudate, increases slightly in voids ratio with m T at the low ram speed for h d = 90°and 45°. However, the voids ratio also decreases slightly as the paste flows through the die land at high m T , which counters the increase in dilation in the die entry zone and renders ELVF independent of m T . This result could change if the die land length were kept constant at all h d instead of dimension 'L' in Fig. 3.
The V n distributions at the low ram speed in Fig. 6 demonstrate that the formation of static zones at the barrel-die face corner is promoted at m T = 0.33 relative to m T = 0 at h d = 90°and 60°(see also Patel et al. (2007)), with little change occurring at h d = 45°. Increasing m T to 0.67 results in little further change in the flow field for h d = 90°and 60°, although static zone formation is initiated for h d = 45°. These results imply that the maximum die entry angle that can be used without promoting static zone formation decreases with increasing m T .
Comparing the flow fields in Fig. 6 between the two ram speeds indicates that static zone formation is promoted by reducing ram speed for all die entry angles at m T = 0 (Patel et al., 2007) and m T = 0.33, but not at m T = 0.67. At the largest m T , static zone formation is inhibited at h d = 60°by decreasing the ram speed: this (counterintuitive) result is due to the simplicity of the friction model used here, i.e. constant shear yield stress at the wall instead of s s with voids ratio (liquid binder content). The paste at the barrel-die face corner at h d = 60°(m T = 0.67) is more compacted at the low ram speed. Thus, the shear yield stress in the paste (s 0 ) is locally increased. If s s is modelled as constant, then by the definition of the Tresca friction condition (Section 2.4.1; Eq. (13)) the 'effective' value of m T at this corner is decreased. This inhibits static zone formation at the low ram speed, as observed. This effect also explains why decreasing the ram speed causes sticking at the portion of the die face nearest the die entry corner at m T = 0.67 and h d = 45°; compare Fig. 6(c)(ii) vs. Fig. 6(c)(iv). At the low ram speed, the voids ratio near the die entry corner is significantly increased (Fig. 7(c)(iv)). The local, effective value of m T is therefore increased promoting sticking. Both these results at m T = 0.67 are purely numerical in origin and are therefore not discussed further. However, they highlight the need for a sufficiently sophisticated wall friction model to enable these systems to be modelled accurately.
The voids ratio distributions in Fig. 7 predict that increasing m T from 0 to 0.33 at the low ram speed causes a reduction in the voids ratio (compaction) in the upper portion of the barrel at all h d ; see Patel et al. (2007) (m T = 0) and Fig. 7(a-c)(iii) (m T = 0.33). Compaction occurs simultaneously at the barrel-die face corner at h d = 60°and 45°but not at 90°. Increasing m T further to 0.67 has two effects. Firstly, the region of the paste exhibiting the most dilation increases in size at h d = 90°and 45°although not at h d = 60°. Secondly, the paste at the barrel-die face corner compacts more at the higher m T for h d = 90°and 60°, but not for h d = 45°, at which compaction is inhibited. The latter result is due to the sticking of paste to the portion of the die face nearest the die entry corner (Fig. 6(c)(iv)), which does not occur at m T = 0.33 (Fig. 7(c)(iii)). The sticking of paste at the die face promotes a pressure difference between the paste at the barrel-die face corner (high p) and the paste above the die entry (low p). This promotes the flow of paste radially inwards from the corner to the die entry, and thus inhibits compaction and static zone formation. Collectively, these results demonstrate that increasing m T from 0 to 0.33 or 0.33 to 0.67 promotes different aspects of LPM at different h d .
Comparing the voids ratio distributions across ram speeds demonstrates that at m T = 0 (Patel et al. (2007), reducing the ram speed promotes compaction everywhere in the barrel except for the die entry zone. At m T = 0.33, reducing the ram speed again causes compaction throughout the barrel except for the die entry zone, but compaction is particularly pronounced at the barrel-die face corner for h d = 45°; see Fig. 7(c)(iii). At m T = 0.66, reducing the ram speed once again drives compaction everywhere in the barrel except the die entry zone, but compaction is particularly pronounced at the barrel wall for h d = 90°. These results demonstrate that decreasing ram speed promotes different aspects of LPM at different h d , and that these changes are coupled to m T .

Discussion
The simulations have established that complex relationships exist between the extent of LPM and the ram speed (pre-eminent factor), the die shape and the friction factor (secondary, equally important factors). These relationships are coupled such that few general trends exist. For the (strongly plastic) pastes modelled here, the use of high ram speed ensures both a consistent extrusion pressure and avoids LPM, thus maintaining consistent extrudate quality (liquid binder content) and avoiding loss of the paste in the barrel due to excessive compaction. The extrusion pressure is dominated by the pore liquid pressure, which is high within the barrel and decreases through the die entry zone and die land. Maintaining a binder-tight seal at the ram and die plate is therefore critical to preventing binder loss via the ram and die plate seals. Extrusion pressure increases with wall friction factor (as expected), which affects the flow field and distribution of LPM but not its overall magnitude (i.e. the extrudate binder content). Static zone formation is promoted by wall friction (as expected) but accurate predictions for CED geometries require a more sophisticated wall friction model than used here.
At low ram speed, significant LPM occurs and this gives rise to transients in the extrusion pressure profile via competing effects on the pore liquid pressure and the effective stress at the ram. These effects are coupled to the die shape, wall friction factor and ram displacement. The only general trends are an increase in the pore liquid pressure and a decrease in the effective stress at the ram with increasing die entry angle. These trends reflect the resemblance of the SED case to one-dimensional confined compaction of a soil, a process during which the pore liquid pressure increases by much more than the effective stress. Both trends increase in magnitude with wall friction. This reflects the fact that as the die entry angle decreases, slip at the barrel wall and die plate is promoted. Therefore, an increase in the magnitude of wall friction (which is sustained by the solids skeleton) has a more dramatic effect on the flow field and voids ratio distribution (LPM) for a highly conical die than for SED geometry.
For the short duration tests considered here, the extrudate binder content is higher at the low ram speed (increased LPM), but is not a function of die shape or friction factor. The former independency is due to the extruder dimensions in the model (L = f(h d )), and the results imply that the extrusion pressure would increase significantly with die land length (increasing LPM by increasing the pore pressure differential across the paste) but would also promote greater recompaction in the die land. Thus, the role of the die land length on LPM requires further elucidation. Static zone formation is increasingly minimised as the die entry becomes more conical (smaller h d ), as expected. However, in the presence of moderate-to-significant wall friction, LPM (compaction) occurs along the entire barrel length for highly conical dies. Collectively, these results indicate that a 'static', one-dimensional model of ram extrusion is unlikely to suffice for use in an industrial setting.
The current model requires several material parameters as data inputs. The uniaxial yield stress of the paste at a range of strain rates (« 1 s À1 in the barrel, » 1 s À1 in the die entry zone) is required for accurate predictions of both extrusion pressure and extent of LPM. A sophisticated wall friction model that is calibrated across a range of slip rates (order 1 mm/s at the barrel wall, » 1 mm/s at the die land wall) is also required. The accuracy of this model is paramount as wall friction affects all calculated parameters (extrusion pressure, flow field and extent of LPM). Finally, the permeability of the solids skeleton is required to determine the onset velocity of LPM. Collectively, these data allow the estimation of the extrusion pressure and the extent of LPM at small ram displacements, i.e. before significant variations develop in the voids ratio distribution. If LPM is estimated to be significant at small ram displacements, the operating conditions (ram speed), extruder geometry (die design) and paste formulation (permeability, yield stress) can be adjusted (subject to any external constraints) to minimise LPM.
There are two major aspects where the simulations could be improved. Firstly, the list of assumptions in Section 2.1 should be reduced to widen the range of paste formulations that can be simulated, e.g. by accounting for viscoplasticity, binder compressibility and the presence of entrained air in the paste. Secondly, a robust (and rapid) experimental protocol is required for measurement of the model parameters. These would ease the adoption of finite element simulations by industry to support extruder and die design.

Conclusions
Two-dimensional, axisymmetric simulations of ram extrusion of a stiff paste, modelled as a soil, have been performed using finite element modelling and adaptive remeshing. This model is derived from an earlier simulator that assumes the extruder wall to be frictionless (Patel et al., 2007). Wall friction is now incorporated into this model, which represents the first time that friction has been considered in a ram extrusion model of this sophistication.
The extent of LPM is predicted to be coupled to the ram speed (most important factor), die shape and friction factor (of roughly equal importance) and the ram displacement. Whereas the likelihood of LPM occurring is predicted reasonably well by the two 1-D criteria considered, the extent of LPM and the effect of die geometry are not predicted by these approaches. Dilation plays an important role and this requires two-dimensional and two-phase simulations. While the material parameters used here were taken from the literature rather than determined for a specific paste, these results still demonstrate that extrusion models must incorporate at least two material phases and at least two spatial dimensions to attain high predictive accuracy. Future development of the simulator will include generalisation of the rheological model for the paste and a more sophisticated model for the frictional interactions between the paste and extruder wall.