Systematic design optimization of grabs considering bulk cargo variability

Ship unloader grabs are usually designed using the manufacturer’s in-house knowledge based on a tra- ditional physical prototyping approach. The grab performance depends greatly on the properties of the bulk material being handled. By considering the bulk cargo variability in the design process, the grab per- formance can be improved signiﬁcantly. A multi-objective simulation-based optimization framework is therefore established to include bulk cargo variability in the design process of grabs. The primary objec- tive is to reach a maximized and consistent performance in handling a variety of iron ore cargoes. First, a range of bulk materials is created by varying levels of cohesive forces and plasticity in the elasto-plastic adhesive DEM contact model. The sensitivity analysis of the grabbing process to the bulk variability allowed three classes of iron ore materials to be selected that have signiﬁcant inﬂuence on the product performance. Second, 25 different grab designs are generated using a random sampling method, Latin Hypercube Design, to be assessed as to their handling of the three classes of iron ore materials. Of this range of grab designs, optimal solutions are found using surrogate modelling-based optimization and the NSGA-II genetic algorithm. The optimization outcome is veriﬁed by comparing predictions of the optimization algorithm and results of DEM-MBD co-simulation. The established optimization framework offers a straightforward and reliable tool for designing grabs and other similar equipment. (cid:1) 2021 The Society of Powder Technology Japan. Published by Elsevier B.V. and The Society of Powder Technology Japan. This is an open access article under the CC BY-NC-ND license (http://creativecommons. org/licenses/by-nc-nd/4.0/).


Introduction
Iron ore products are transported from origin mines to customers, mainly steel manufacturers. Due to the high demand for iron ore products in the steel industry, iron ore products have the largest dry bulk trading per year, more than that of coal and grains [1]. The top 5 importing countries are located in Asia and Europe, while the majority of iron ore mines are located in Australia and Brazil. This fact necessitates the intercontinental shipping of iron ore products. In 2018, a total of 3210 million tonnes of dry bulk solids were shipped, of which 46% was iron ore [1].
Grabs are extensively used to unload iron ore cargoes from bulk carriers at destination. The general model of the grabbing process (i.e. of an iron ore cargo) is illustrated in Fig. 1 schematically. A grab is connected to a crane and a crane operator controls the grab using the motion of winches and cables. In a grab cycle, the crane operator transfers the grab to a cargo hold to collect the bulk material. Next, the grab is lifted and transferred back to the quayside to release the collected bulk solid. To unload a bulk carrier with a capacity of 150 000 tonnes (DWT), approximately 5000 to 7500 grab cycles are required.
To minimize the costs of terminal operators and other stakeholders the waiting time of bulk carriers needs to be as short as possible. A time-efficient and reliable unloading process is therefore required at destination ports. Considering the increasing global demand for iron ore [2], the unloading process could be improved in terms of productivity to use available facilities, such as cranes, in a sustainable way. In general, cranes have a significantly higher total cost of ownership than grabs. Therefore, a promising solution for improving the unloading process is to enhance the design of grabs to allow bulk carriers to be unloaded with a minimum number of grab cycles.
In addition to the grab design itself, dimensions of the ship's hold, properties of bulk cargo, crane operator, winches and cables are the main contributing elements in the grabbing process. A crane operator controls the grab by using winches and cables. Predicting performance of new design concepts is a challenge, as it requires that interactions between multiple contributing elements be considered. Assessment of a design concept involves modelling its dynamics and kinematics as the grab interacts with a bulk solid cargo. Controllable inputs of the grabbing process are operation and design parameters. Operation parameters (e.g. velocity of cables), are controlled during the handling process. It is known that in general, there is a positive correlation between skills of a crane operator and productivity. Design parameters (e.g. bucket dimensions and length of opening span) can be adjusted when the grab is manufactured. Design parameters are major controllable input prior to manufacturing a grab as they eventually determine the performance during operation.
Schott et al. [3] developed a novel design method for grabs based on the Discrete Element Method (DEM) and MultiBody Dynamics (MBD). They demonstrated how a grab design can be improved systematically by creating various virtual prototypes in interaction with a specific bulk cargo. In practice, however, there is a variety of iron ore products that need to be unloaded at a destination port. Mohajeri et al. [4] have used validated DEM-MBD cosimulation setups to demonstrate how the grabbing process differs between free-flowing and cohesive iron ore. The variability of bulk solid properties influences the grabbing process considerably [16]; this is a major possible source of deviation from the desired grab performance [4][5][6]. To design a sustainable product, a minimized performance deviation should be achieved. There are numerous DEM-based research articles on the design of equipment, including the interaction with bulk material, such as [7][8][9][10][11][12][13][14][15]. However, bulk cargo variability, as one of the main factors in the performance of bulk handling equipment, has not been explicitly incorporated in the design process. A grab is often used to handle a broad variety of iron ore cargoes that differ as to their properties, such as moisture content, shear strength and bulk density.
The grab performance can be quantified by using the Key Performance Indicators (KPIs) defined in [3]. The mass indicator (W mass ) is the most important performance indicator that determines the efficiency of a grab cycle. W mass is quantified by comparing the collected mass with the weight of the grab using Eq. (1).
where M DWT is the collected mass inside the grab, M e is the mass of the grab, and M spillage is the mass of the bulk material spilled during the grabbing process. Including the spilled material in the equations allows for capturing the effectiveness of the grabbing process [17].
For a specific hoisting capacity of a crane, maximized W mass values are desired for handling a variety of iron ore products. Fig. 2 shows an overview of how the bulk cargo properties contribute to the grabbing process as an uncontrollable input variable. Mohajeri et al. [5] measured grab-relevant bulk properties of a broad range of iron ore fines. Cohesive forces (i.e. liquid bridge) between iron ore particles are typically created when moisture is introduced [18]. Cohesive forces may influence the bulk properties of iron ore fines, such as shear strength and flowability [5]. Bulk compressibility and moisture content are also correlated for cohesive iron ore [5,19]. Pre-consolidation stress is another grabrelevant bulk property of cohesive iron ore that varies over the cargo depth during the unloading process [4]. Due to the increasing overburden pressure, a more consolidated cargo is stored at greater depths.
In this study, a novel optimization framework is developed to incorporate the bulk cargo variability in the design process of grabs. Discrete Element Method (DEM) is employed to first investigate how a virtual grab prototype can be tested considering the bulk cargo variability, including various levels of cohesive forces and bulk plastic compressibility. Such a sensitivity analysis allows for selecting bulk material classes that create significant deviation in the grab performance. This follows from optimizing a virtual prototype to reach a maximized W mass in handling a variety of significant bulk cargo classes while the deviation of grab performance is minimized. Multiple surrogate models are created to find optimal design settings, which are evaluated in a verification step.

Multi-objective optimization framework for including bulk cargo variability
This section describes the multi-optimization framework developed to incorporate the bulk cargo variability into the grab design process. Both controllable and uncontrollable types of input are included in the framework. A matrix, [W], containing performance indicator values can be quantified for a combination of uncontrollable and controllable inputs, as shown in Fig. 3. The primary aim is to minimize the undesirable effect of variability of X on Y. Thus, for an optimal design configuration, Y opt , a maximized performance, W mean , is reached on average, while its standard deviation, W SD, is minimized.
The optimization framework is designed in four sequential steps where the output of each step is used in the next step as illustrated in Fig. 4. This allows grab designers to follow a straightforward procedure when a new concept is being developed. To model grabs in interaction with bulk solids, the framework of [20] is used to make a two-way coupling between a DEM solver and a MBD solver. A coupling server communicates between two solvers at each time interval; the geometry motion is calculated using the MBD solver, and the reaction forces on the geometry are calculated using the DEM solver. This allows for creating a real-scale co-simulation between grab and bulk material [21], which requires virtual crane operator, CAD model of grab, and bulk material model as inputs.

2.1.
Step 1. Sensitivity of the grab performance to bulk cargo variability 2.1.1. Reference material model of the cohesive iron ore -X ref A DEM material model of a cohesive iron ore cargo, named Carajas Sinter Feed (CSF), has been validated for the grabbing process by Mohajeri et al. [4]. In the current study, that validated material model is used as a reference material model, X ref , to create a bulk cargo variability. Main parameters of the reference material model of the cohesive iron ore is presented in Table 1. To model interaction between particles in the reference material model, an elasto-plastic adhesive contact spring, named EEPA [22], is used. Fig. 5 shows a schematic diagram of the non-linear mode of the EEPA contact spring in the normal direction. For details of the EEPA contact spring we refer to [18,22]. In the EEPA contact spring, the cohesive forces can be adjusted by varying the constant pull-off force (f 0 ) and surface energy (Dc). Sensitivity studies on the dependency of bulk behaviour (e.g. angle of repose, shear strength, bulk density) on the variation of f 0 and Dc have been documented in [23][24][25]. For f 0 and Dc, reference values of À0.2 N and 100 J/m 2 are used respectively.
The plasticity ratio (k P ) controls the contact stiffness during unloading and reloading of the spring and, thus, this parameter controls the bulk compressibility. The plasticity ratio, k P , controls the ratio between stiffness in branch II (k 2 ) and stiffness in branch I (k 1 ), which shows the influence of plasticity ratio at contact scale. This means that by increasing the plasticity ratio, a higher level of plastic overlap occurs during contact and, thus, a higher level of bulk compressibility. For k P , a reference value of 0.2 is used.
The EEPA contact spring is able to capture a stress-historydependent behaviour [4,26,27] and, therefore, no input parameters need to be adjusted in the material model for this purpose. A preconsolidated situation can be simulated by applying a specific amount of pressure on the bulk surface and then releasing that pressure, as described in [25]. The reference material model has been validated in operational conditions for two different levels of pre-consolidation: 65 and 300 kPa. The grabbing process of the cohesive iron ore, for various levels of pre-consolidation, has been investigated in [4], which showed the negative effect of pre-consolidation on the grab performance.

Bulk cargo variability -[x]
This sensitivity analysis evaluates whether the variability of cohesive forces and bulk compressibility influences the grabbing process or not. The effect of a variable is considered significant if it creates ±5% deviation in the mass indicator. As displayed in Table 2, a bulk variability, [X], based on the reference material model is created.
The relative cohesion term, as defined in [23], is used to vary the level of cohesive forces when the EEPA contact spring is applied. The relative cohesion, C bulk , distinguishes between the expected levels of bulk cohesion in a qualitative way. For example, when using a lower relative cohesion, a lower angle of repose is expected compared to the reference material model. To create a low relative cohesion, f 0 and Dc are decreased by 50% compared to the reference material model. An increase of 100% is also applied to create  a high relative cohesion. C bulk is set to ''non" in bulk materials 1, 2, and 3, by setting both f 0 and Dc to 0. The EEPA model behaves like an elastic spring if the plasticity ratio is set to 0, while using values close to 1 the model behaves like a plastic spring. In an elastic spring there is no residual overlap once the force drops to zero. Any values between 0 and 1 result in a certain level of plastic compressibility in the contact spring. The reference material has a plasticity ratio of 0.2, which we correspond to low relative plastic compressibility, k bulk . Medium and high levels of relative plastic compressibility are defined by using 0.55 and 0.9 for k P respectively. If k P , f 0 and Dc are all set to zero, then the material model behaves like a non-cohesive elastic bulk solid. The grabbing process of non-cohesive elastic iron ore has  been already investigated in [3,4] and is, therefore, excluded from the current sensitivity analysis.

Simulation setups
Grab-relevant behaviour of all 12 bulk materials are evaluated in the following simulation setups: Angle of repose Uni-axial consolidation Penetration test These preliminary simulations are executed to verify that the created virtual bulk variability represents various states of bulk cohesion and compressibility. Next, the grabbing process is simulated at full-scale, which allows the influence of bulk variability on the process to be investigated. The particle diameter of the validated material model is relatively large (55 mm in diameter), compared to particle sizes used in typical laboratory scale DEM simulations. Therefore, relatively large domains are also created to fit enough numbers of particles without undesirable boundary effects.
The angle of repose is simulated by pouring particles from a specific height. The simulation setup is shown in Fig. 6. Particles are created in a factory that is located 1.5 m above the bottom plate; due to the force of gravity, particles drop on the bottom plate to form an angle of repose over time. 2500 particles are created with a total mass of around 800 kg. Once the simulation is finished, a stable angle of repose is formed, and the position of particles that are on the slope is analyzed. A linear regression is then fit on the data points to determine the angle of repose. The angle of repose, a M, is therefore the measurement objective of the simulation.
The uni-axial consolidation process, including loading and unloading, is simulated in four stages to evaluate the bulk compressibility as well as bulk density. The simulation domain is 1x1x2 meter. A block of material is created using a particle factory that moves upward. This kind of technique minimizes the impact force during the particle generation (Fig. 7a). Next, the particles are allowed to settle for 2 s, and a low kinetic energy in the bulk material ( 1eÀ4 J) is reached (Fig. 7b). Next, the bulk material is consolidated by applying a uniform pressure (i.e. 65 kPa) on its surface by means of a geometry plate for 2 s (Fig. 7c). The pressure is unloaded by moving the geometry plate upward with a velocity of 1 cm/s (Fig. 7d).
The initial bulk density, q b,0 is quantified when the particles relax in the second stage. The compressed bulk density, q b,c is measured at the end of loading in the third stage. The final bulk density, q b,end is quantified when the unloading is finished and a preconsolidated situation is created.
The penetration resistance of the bulk material is the third grabrelevant property that is investigated here. The penetration resistance is an influential bulk property in the grabbing process [5], as a lower resistance to penetration of grabs into the bulk solids results in a higher payload generally. The penetration process is simulated for a material block that is pre-consolidated with a vertical pressure of 65 kPa, as shown in Fig. 8a. A cube-shaped geometry with the volume of 8 m 3 is used to contain the material block.
In general, ship unloader grabs have wedge-shaped knives with a blunt tip to trade-off between the penetration resistance and amount of wear. The wedge-shaped penetration tool has a width of 40 mm and its tip is 20 mm wide. That makes the crosssection of the penetration tool similar to the setups used in [5,25,27] that focused on the grab application too. This tool (I) is driven into the pre-consolidated bulk material (II) with a constant velocity of 0.1 m/s. A plane contact 2000 mm in length is created during the penetration, which replicates the grab dimensions adequately. The reaction force on the penetration tool is quantified as the measurement objective.
Once the outcome of preliminary simulations confirms that an adequate bulk variability is achieved, the grabbing process can be simulated for the 12 bulk materials. The DEM simulation of the grabbing process is run on a combination of CPU and GPU. This allows for reducing the computation time of a MBD-DEM cosimulation by around 6 times, compared to a CPU-based cosimulation. NVIDIA Quadro GP100 is used as the graphics card in this study.
Once the co-simulation is finished, the grab performance is quantified for the 12 different bulk materials. The mass indicator,

Step 2. Random sampling from design space (LHD)
Once the significant bulk material classes are created, a parametric variation of the grab design can be investigated. In step 2 of the optimization framework, design space is searched effectively to create randomized variations of grab configurations. If all the possible combinations of variables with the design space are considered, a full factorial design is thus created. For each parameter, a series of levels, or values, N s , is defined. When every possible combination is tested, the total number of samples, N', is given by Eq. (2). Fractional factorial designs can offer more effective sampling methods compared to a full design, in terms of offering an affordable computation time [25]. The Latin Hypercube Design (LHD) method is selected in this study, as it allows for searching a parameter space effectively using a minimum number of sampling points [28]. A set of sampling points is constructed in such a way that each of the parameters is divided into p equal levels, where p is the number of samples. This is illustrated in Fig. 9 using two examples and for two parameters. In example 1, the samples are constructed with an extremely poor space filling quality, while example 2 has a better filling quality with a fine filling of the design space.
The LHD is constructed according to the algorithm developed in [29]. The U P criterion was defined, as shown in Eq. (3), to measure the performance of a LHD-based sampling.
where p is a positive integer, d ij is the inter-point distance. In the current study, p = 50 is used following the recommendation of Jin et al. [30]. By minimizing the U P criterion [31], the location of levels for each parameter is randomly, simultaneously, and evenly distributed over the parameter space. Maintaining a maximized distance between each two points allows for satisfying the U P criterion.
In total, five design variables are included in the optimization, referred here as D1, D2, D3, D4, D5. The variables and the range   of variations are selected based on a previous parametric study [3] as well as in consultation with grab designers. For example, Schott et al. [3] demonstrated that the length of the grab bucket, D1, plays an important role in the grab performance. Also, the radius of a bucket, D2, is a significant design parameter as it influences the bucket shape, and, thus, its volume. The construction of the LHD-based samples for D1 and D2 is visualized in Fig. 10. Samples for three other design variables are randomly created in a similar way, thus, minimizing U P for five variables. A range of 1650-2000 mm is considered for D1, as it is a typical range for such a grab prototype. For the same reason, D2 is also varied between 1200 mm and 2000 mm. Therefore, 25 different grab designs, N 0 = 25, are created, including 5 variables, N g .

Step 3. Multi-objective optimization using surrogate modelling
A surrogate model is a computationally affordable mathematical model that can replace the actual simulation or experiment. Surrogate models approximate a function based on a set of available data points and can then predict the function at new points [13]. A surrogate model offers a faster computation time, compared to the actual DEM-MBD co-simulation, to predict performance of a new grab configuration. Surrogate models can be also used to obtain trends and identify the influence of specific parameters on the grab performance. Three different types of regression-based surrogate models are tested in the current study:

Linear regression Linear Support Vector Machine Kernel Polynomial Support Vector Machine Kernel
The linear regression model is the most widely used regression model. In general, this type of regression model is a linear function between variables, response of the system, and constant coefficients [32], as shown in Eq. (4).
where f is the regression model, b is a constant coefficient, and y is a (design) variable. A surrogate model can be created by fitting a regression function, f k , for each bulk material. Therefore, f(Y) maps the relationship between the grab design variables, bulk materials, and the selected response of the system, which is the mass indicator, W mass , in the current study.
Kernel models transform variables using kernels. The transformed variables are measures for similarity or correlation between the data points. Multiplying the transformed variables with weights (constant coefficients), as with the linear models, gives an estimate of the output. The predictor function for a kernel-based regression is given by: where b i is the weight factor corresponding to data point i, and y' indicates a vector of variable (at its new location) and y i is an available data point. b is a constant to minimize the fitting error, e. One difference between a linear regression model and a kernel model is that the latter has a number of coefficients, bi, corresponding to data points rather than variables. U is the kernel function that transforms data points into another space to handle the nonlinearity. Linear and polynomial functions are used for U in the current study. The support vector machine (SVM) regression uses a kernel function to first estimate the correlation between data points before fitting coefficients (Fig. 11). The advantage of SVM is that it allows for an error between observations and predictions [33]. The cost function is not increased until the specified amount of error, e, between observation and prediction is reached, which forms an e-tube around the prediction function. Outside the tube, the cost function increases and forces the prediction function to a specific range of data points.
As discussed earlier, two objectives are considered in the current optimization: a maximized average mass indicator (W mass, mean ) and a minimized standard deviation for mass indicator (W mass,SD ) measured in different bulk materials. The unloading frequency of a certain cargo is also considered in the optimization. For example, if a grab unloads a specific cargo 20% of time, and another cargo 80% of time, the second cargo should have a higher weighting factor in the optimization for maintaining an adequate productivity. The distance between origin mine and customer, production capacity of mine, and technical demands of customer are among the influencing factors on the frequency of receiving a specific bulk cargo at destination. The unloading frequency can usually be obtained by analyzing available databases of customers. Therefore, to consider the frequency distribution of bulk variability, weighting factors with P N Ã k¼1 w k ¼ 1 are defined. w k is the weighting factor of material k in the optimization.
Once different grab samples, Y, are simulated, one can select a design configuration that may jointly satisfy the optimization objectives. However, a response optimizer can find better design configurations, compared to the simulated samples, by using the surrogate models. Creating surrogate models allows for predicting the response of the system without the necessity of running a DEM-MBD co-simulation. Once a surrogate model is created, the optimal design can be found by selecting a combination of design  variables that jointly satisfy the optimization objectives [25]. The NSGA-II genetic algorithm [35] is a proper tool to solve DEMbased optimization problems [36][37][38][39], and is therefore used in the current study to search for the optimal solution within the design range.

Step 4. Verifying the optimal design
The selected optimal design is verified by running conforming simulations. This allows for quantifying the error of surrogate models as well. The prediction error is quantified using Eq. (6).
where |e| mean is the mean of absolute relative differences for the grabbing process in N* different bulk materials. f k is the prediction for system response, W mass in the current study. f k ' is the simulated response of the optimal design solution for bulk material k. The acceptable error of |e| mean is considered to be 5% multiplied by N* . In other words, on average a prediction error of 5% for each bulk material is considered to be adequate. If the prediction error is not acceptable, the number of data points in step 2 can be increased to improve the accuracy. Additionally, based on the prediction error, the performance of the different surrogate models can be compared. The optimization ends with a verified optimal design configuration, Y opt .

Results and discussion
This section presents and discusses the outcome of four steps of the optimization in a sequential manner.

Results of step 1, sensitivity to cargo variability
Step 1 aims to identify a bulk material variability that has a significant influence on the grab performance. First, results of preliminary simulations are discussed. Second, the grab performance in handling the 12 different material models is analyzed. Third, a matrix, [X*], containing the significant bulk material classes is created. Fig. 12 shows the angle of repose results including two variables; relative cohesion (C bulk ) and relative plastic compressibility (k bulk ). The angle of repose depends on the relative cohesion significantly. Increasing cohesive force values, f 0 and Dc, results in a higher angle of repose. The relative plastic compressibility, also influences a M . When a non-cohesive material is used, the relative plastic compressibility has a positive correlation with a M . However, when the cohesive forces are present, the relative plastic compressibility has a negative correlation with a M .

Angle of repose
In case of non-cohesive materials, a higher contact plasticity results in a larger contact area upon unloading, thus increasing the required sliding distance of particles relative to each other. However, when cohesive forces are active, a higher relative plastic compressibility results in a denser pile of material. Since the particle density is constant, a denser packing of material results in a heavier failure wedge in the slope, thus a lower angle of repose could be expected with increasing the contact plasticity. The effect of contact plasticity on the packing is discussed further in the uniaxial consolidation simulation setup. Fig. 13 displays initial, compressed, and final bulk density values that are quantified for the 12 different bulk materials under 65 kPa pre-consolidation pressure. Results are presented in three separate graphs, each showing the outcome for a certain level of k bulk . All bulk density parameters decrease when cohesive forces increase, independent of the contact plasticity value.

Uni-axial consolidation
The higher the cohesive forces, the larger the restrictive forces between particles to fill the voids; consequently, a lower bulk density is created. Furthermore, by increasing the contact plasticity, the residual overlap in contact spring increases [18], thus a smaller difference between q b,c and q b,end might be expected. Therefore, Fig. 11. A tube with the radius of e is fitted to data points in the SVM regression model [34]. both variables, k bulk and C bulk , have significant influence on the bulk compressibility and the bulk density.
3.1.3. Penetration resistance W 500 , the accumulative reaction force (in Joules) on the wedgeshape tool is quantified at the penetration depth of 500 mm. That is similar to the penetration depth that occurs in the grabbing process of the CSF cargo under 65 kPa pre-consolidation pressure [4]. The outcome of the penetration test simulations is shown in Fig. 14, including two variables: k bulk and C bulk .
There is a positive correlation between the relative plastic compressibility, k bulk , and W 500 . A higher contact plasticity results in a denser packing, thus, a higher resistance against the penetration of the wedge-shaped tool. There is no clear relationship between the relative cohesion and the penetration resistance. Therefore, only k bulk is a significant bulk variable influencing the penetration resistance.
The influence of each variable on the grab-relevant bulk properties is shown above. The relative cohesion has a significant influence on the angle of repose and bulk density, while the relative plastic compressibility plays a significant role in the angle of repose, bulk compressibility, and the penetration resistance. Fig. 15 displays the influence of k bulk and C bulk on the grab performance. The influence of relative plastic compressibility on the grab performance is significant. That could be expected, based on the penetration simulations. The relative cohesion also plays a role in the grabbing process, especially when a low k bulk is used.

Grabbing process
Although the effect of the pre-consolidation pressure is not investigated in the current analysis, it is known that the preconsolidation plays a significant role in the grabbing process of cohesive iron ore [4]. Three different bulk materials are selected for further optimization of the grab design, as presented in Table 3. Material 1* is a non-cohesive iron ore with no relative plastic compressibility, the material model of which is developed in [21]. Due to lack of compressibility, no pre-consolidation is applied on material 1*. Material 2* is a cohesive iron ore with a low C bulk and a high k bulk that is pre-consolidated with a relative high pressure of 200 kPa. By contrast, material 3* has a high C bulk and a low k bulk , that is pre-consolidated with a relative low pressure of 40 kPa. Such pressure is expected at a cargo depth of around 1.5 m to 2 m. By analysing an available database of a grab customer, the weighting factors are selected for each bulk material. Summarising, three different bulk material classes with significant variability for the grabbing process are selected as the outcome of step 1.    W mass,mean values of samples are distributed between 2.05 and 2.54, thus a variation of around 24% in the performance of different design samples is captured. That confirms the adequacy of the random sampling approach that is based on LHD. Next, surrogate models are fitted on the 25 data points.

Optimal solutions
Three different surrogate models are fitted on the available data points: linear regression, linear SVM kernel, and polynomial SVM kernel. Next, optimal solutions using the NSGA-II genetic algorithm are found for each surrogate model. The outcome is illustrated in Fig. 17, indicating that different optimal solutions (red line) are found using different surrogate models. The available data points are shown in blue. The polynomial SVM kernel predicts optimal solutions that are better than the predictions of two other surrogate models. The non-linear relationships between optimization objectives and design variable could be captured well using the polynomial SVM kernel.

Verified optimal design
It needs to be verified whether the predications of the surrogate models and the optimization algorithm are sufficiently accurate. Therefore, the ''knee-point" in the line presenting optimal solutions is selected, as recommended in [40]. Co-simulations are executed for each optimal solution in three significant bulk material classes, [X*]; the outcome is shown in Table 4.
All three surrogate models have a prediction error smaller than 5% for W mass,mean . The polynomial SVM kernel shows the highest grab performance as well as the lowest prediction error, while the linear regression model shows the opposite. Therefore, the polynomial SVM kernel can be recommended as a surrogate model to find design configurations of an optimal grab, including the bulk cargo variability.

Conclusion
In this study, a sequential multi-objective optimization framework was established to include multiple grab design variables as well as a variety of bulk material properties in the design process. A wide range of bulk material properties was used in the optimization, from a non-cohesive incompressible iron ore to a pre-consolidated cohesive compressible cargo. By executing the optimization analysis, a maximized grab performance, W mass,mean , was achieved, while a minimized value for the performance deviation was maintained. The established optimization framework offers a straightforward and reliable tool for designing grabs and other similar equipment, including the bulk cargo variability.
A virtual bulk variability was created to consider various levels of cohesion and compressibility of iron ore products. Three preliminary simulations were performed to verify that a realistic bulk variability is replicated using DEM. The outcome of simulations show that the relative cohesion has a significant influence on the angle of repose and bulk density, while the relative plastic compressibility plays an important role in the angle of repose, bulk compressibility, and the penetration resistance. The simulations of the grabbing process, in which a range of virtual bulk materials is used, showed that the relative plastic compressibility has a larger influence on the product performance, compared to the relative cohesion. 25 different random grab designs were created using the Latin Hypercube Design sampling method, including 5 different geometrical variables. A variation of 24% in the grab performance was captured using the random design samples, indicating the adequacy of the sampling method. Comparing the average mass indicator values, W mass,mean , as well the corresponding standard deviation values allows for assessing performance of different grab designs. Three different surrogate models were created using linear regression, linear support vector machine kernel and polynomial supper vector machine kernel models. The outcome of the optimization was most promising and accurate when the surrogate model is constructed using a polynomial SVM kernel model, as it captures the non-linear relationships between variables and objectives.
By following steps II, III and IV of the optimization framework, design concepts can be enhanced by including the bulk cargo variability. If a design concept needs to be optimized for a different type of cargo (e.g. coal), it is recommended that all steps of the optimization framework are followed in a sequential manner. In the current study, the virtual bulk variability was created by changing levels of cohesion and compressibility, which are dominant sources of bulk variability for fine and moist iron ore cargoes. In future studies, other possible sources of bulk variability can also be investigated, such as adhesion between geometry and particles, particle shape and size distribution.
Future research should focus on finding a design solution to create an impact-less grab that is entirely insensitive to uncontrollable factors. Although the grab operation is a controllable input of the process, an ideal design solution would also be insensitive to operational parameters. This can be achieved by creating a robust control system for grabs.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.