Multiphysics Modeling and Numerical Simulation in Computer-Aided Manufacturing Processes

: The concept of Industry 4.0 is deﬁned as a common term for technology and the concept of new digital tools to optimize the manufacturing process. Within this framework of modular smart factories, cyber-physical systems monitor physical processes creating a virtual copy of the physical world and making decentralized decisions. This article presents a review of the literature on virtual methods of computer-aided manufacturing processes. Numerical modeling is used to predict stress and temperature distribution, springback, material ﬂow, and prediction of phase transformations, as well as for determining forming forces and the locations of potential wrinkling and cracking. The scope of the review has been limited to the last ten years, with an emphasis on the current state of knowledge. Intelligent production driven by the concept of Industry 4.0 and the demand for high-quality equipment in the aerospace and automotive industries forces the development of manufacturing techniques to progress towards intelligent manufacturing and ecological production. Multi-scale approaches that tend to move from macro- to micro- parameters become very important in numerical optimization programs. The software requirements for optimizing a fully coupled thermo-mechanical microstructure then increase rapidly. The highly advanced simulation programs based on our knowledge of physical and mechanical phenomena occurring in non-homogeneous materials allow a signiﬁcant acceleration of the introduction of new products and the optimization of existing processes.


Introduction
The implementation of new systems of Computer Aided Engineering (CAE) and products requires a systematic approach to the development and design of technical systems and products. Many customisable product development methodologies have been developed in recent decades, i.e., Autogenetic Design Theory (ADT) [1], Engineering Design [2], Product-Service Systems (PSS) [3], Integrated Product Development (IPD) [4], and the Verein Deutscher Ingenieure (VDI) 2221 guidelines ( Figure 1) [5]. A systematic product development methodology for developing engineering systems in smart manufacturing has been presented in the paper by Gogineni et al. [6].
The essence of numerical modeling of metal forming processes is to replace a continuous object model with a discrete model that takes the form of a system of algebraic equations. In mathematical modeling of a continuous medium, the body is treated as a geometric object in Euclidean space. During calculations, physical quantities, such as loads, displacements, stresses, and deformations, are subject to discretization and are represented by continuous functions. Numerical simulation is a complex process that takes into account The implementation of the idea of Industry 4.0 requires the computerization of traditional manufacturing industry and the development of numerical methods supporting the decision-making process and accelerating the process of designing tools and new products [7,8]. In recent years, there has been a dynamic development of two-and three-dimensional numerical modeling of metal forming processes using the Finite Element Method (FEM), Finite Volume Method (FVM), Computational Fluid Dynamics (CFD) method, Boundary Element Method (BEM), mesh-free, multi-grid, and mesh-free methods, Discrete Element Method (DEM), Extended Finite Element Method (XFEM), FEM/DEM coupling methods, Crystal Plasticity Finite Element Method (CPFEM), Cellular Automata (CA), Fast Fourier Transformation (FFT), Arbitrary Lagrangian-Eulerian (ALE) method, etc. [9,10]. The numerical methods of approximate solving of continuous boundary-initial problems include the finite difference method (FDM) [11], smoothed particle hydrodynamics (SPH) method [12,13], finite point method [10,14], and particle-in-cell (PIC) method) [10,15].
Scientific conferences on plastic working allow for the exchange of experiences and views of scientists. Comparing different views on the same topic allows us to verify the direction of further research works. The aim of the cyclically organized NUMISHEET and the Numerical Methods for Industrial Forming Processes (NUMIFORM) conferences is to integrate the scientific community and create a platform for building interdisciplinary research teams. The International conference "ESAFORM on Material Forming" is organized under the auspices of the European Scientific Association For Material Forming. The aim of the association is to promote scientific research in the field of material development, disseminate scientific information in this field, and create a platform for the exchange of experiences between research centers and industries. The NUMISHEET conference series has been established as a world-class forum through which new technologies in the area of sheet metal forming simulation. NUMIFORM international conference series aim to provide a forum where future directions in the numerical modeling of plastic working processes are discussed by engineers from industry and by scientists. The authors encourage readers to review proceedings of the above-mentioned conferences, which are a source of innovative technologies and the state-of-the-art in numerical simulations of material forming processes.
This article summarizes recent trends in the numerical contributions to metal forming developed over the last decade which fit into the development of the Industry 4.0 concept. The main numerical approaches analyzed in the context of the modeling of complex phenomena existing in metal forming ( Figure 2) include the FEM, BEM, DEM, and meshless methods.

Method
The methodology of systematic review consisted of three steps, which involve a document search, paper selection, and data extraction. The latest progress in the main computational approaches has been conducted following the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines [16]. To fulfill the goal of this paper, searches have been made using international databases, including the DOAJ (Directory of Open Access Journals), Google Scholar, ScienceDirect, Springer, Web of Science, and Worldwide Science. Because computational approaches advance exponentially, it was assumed that the last ten years forms an accurate timeframe to evaluate the most important approaches. Special attention has been placed on the 2015-2020 period. Several aspects have been considered in paper selection: date of publication, subject matter, relevance, and title. The searching for articles has been restricted to the English language. In the next stage, the introduction and conclusions in the article were examined to check whether they provided the necessary evidence to answer the research questions. Finally, duplicate articles found by various databases were removed.

Finite Element Method
Engineering simulations using the latest technological solutions in Computer Aided Engineering (CAE) are one of the most important success factors for modern manufacturing companies. The most frequently used method has become the FEM, which is an increasingly popular tool for structure analysis. There are many finite element-based pieces of software, both paid-for (i.e., Adina, Ansys, Abaqus, COMSOL, Ls-Dyna, MSC Marc) and open source (Z88 Aurora, Calculix, Elmer, MFEM, etc.). FEM has been widely used in most areas of engineering, most frequently in the aviation, automotive, and machine industries for strength calculations and optimization of the forming processes. FEM makes it possible to simulate the flow of fluids [17] and heat transfer problems [18,19].
The correctness of the numerical model is primarily determined by the boundary conditions and the adopted constitutive model of the material [20]. Mesh parameters, such as element size and element type, play a crucial role in ensuring high-quality numerical calculations [21][22][23][24]. Contact models with friction developed over the years and take into account the effect of material deformation, the flattening of the surface asperities under different lubrication regimes, lubricant viscosity, and sliding speed [25][26][27]. Attempts are being made to link the physical parameters of the lubricant, such as film thickness and viscosity, the surface roughness of surfaces in contact, the process parameters (contact pressure, sliding speed), and the elasto-plastic behavior of the interacting asperities with the amount of frictional resistance [28][29][30]. Mixed and boundary lubrication models of contact allow us to determine frictional resistances by considering the effect of load, lubricant properties and changes in sliding speed [31,32]. More complex thermo-mechanical phenomena should be taken into account during friction welding processes [33,34], which are currently of wide interest in the aerospace industry.
An enhanced modeling of material behavior is demanded for the integration of experimental and numerical analysis. Constitutive models of plasticity at the micro-and macroscales are able to describe the elastic and plastic behavior of metallic materials [35][36][37][38][39], ceramic materials [40], porous media [41], fibered sheets [42][43][44][45], gradient materials [46,47], composites [48][49][50][51], and forms [52]. Novel constitutive modeling approaches take into account the change in flow stress depending on such parameters as deformation, deformation speed, the strain hardening phenomenon, the change of orientation of strain state components, temperature, deformation history, and the hardening rule [53][54][55]. The advantages of FEM include the ease of discretization of complex shapes, the ease of determining the boundary conditions, and the ease of adaptive compaction and coarsening of the mesh [56].
The FEM equations are obtained from the formulation of an integral (global) problem using two methods: the variation principle and the weighted residual method. In the first method, a weak analytical formulation is used, e.g., the equation of the principle of virtual work or defining the problem of minimizing a certain function. The weighted residues method transforms the local formulation of the boundary problem into a weak integral form [56]. The Ritz method, also known as the Rayleigh-Ritz method, is a direct method to find an approximate solution for boundary value problems. The Galerkin method is a class of method for converting a continuous operator problem into a discrete problem. It is a more general method than the Ritz method, and it can also be used when the minimization problem cannot be defined and there is no variational formulation (the weak form is then obtained from the weighted residual method) [56].
The FEM is widely used in solving the complex problems of metal forming practice. One of the most numerically demanding metal-forming processes is electromagnetic forming (EMF). It is fundamentally an electro-thermo-mechanical non-contact high-speed process. The theoretical analysis and microstructure observation of Liu et al. [57] with regard to electromagnetic dynamic tests found that the inertia force suppresses the structural instability and improves the mechanical properties by increasing the cross-slip structure. Cui et al. [58] established the 3D FE-based model to investigate the free-bulging analysis of sheet metal. The numerical model takes into account the contact between the die and the sheet and the effect of sheet deformation on magnetic analysis. It was found that failure may first arise at the center. Li et al. [59] analyzed the EMF for Ti-6Al-4V sheet with a copper driver sheet by coupled numerical simulation. In the simulation, AN-SYS/EMAG program was first used to calculate the plate and the magnetic force intensity distribution, and then the magnetic force calculated was introduced into ABAQUS/Explicit. Part et al. [60] experimentally investigated the effect of an aluminium driver sheet on the height of forming a peak of advanced high-strength steel sheet. Moreover, 3D coupled electromagnetic-mechanical FE-based simulations were conducted to analyze the effect of the aluminium driver sheet on the height at which steel sheets form a peak and to predict the proper configuration of the AA1050 driver sheet. As the driver sheet covers more area of the flat spiral copper coil, the Lorentz force and an eddy current are generated more strongly. A thinner driver sheet produces a lower eddy current. As the driver sheet becomes thicker, the force required for deforming the driver sheet also increases. Gies et al. [61] reviewed the current literature on the effect of the thickness of the driver sheet on the Lorentz force and consumer force. It was concluded that the consumer force and Lorentz force increase as the driver sheet becomes thicker.
Finite element (FE) simulations based on beam elements are most efficient for the prediction of elastic-plastic deformations of nanoporous materials. Odermatt et al. [62] predicted that elastic-plastic deformation of nanoporous gold (NPG) is an ideal model material for fundamental research on nanoporous metals. It was found that the regions that undergo plastic deformation (grey color) are located in the thin cross sections in the FEM solid model, whereas they are much more extended towards the junctions in the FEM beam model (Figure 3). Sun et al. [63] developed an Integrated Computational Materials Engineering (ICME) framework ( Figure 4) to analyze the failure behavior of carbon fiber reinforced polymer (CFRP) composites used for lightweight vehicle applications. The properties of the interphase region have been determined by a molecular dynamics analysis and the analytical gradient model. The micro-scale unidirectional representative volume element (RVE) has been implemented in the FE-based Abaqus program. This study illustrates the potential of integrated multi-scale computational modeling for the prediction of failure mechanisms and failure strengths of woven and sheet molding compact (SMC) composites. The RVE plays an important role in the homogenization of random heterogeneous microstructures, especially for composites. The homogenization methods have reached a high level of efficiency, especially in the case of thermal conductivity and mechanical properties [64]. El Moumen et al. [65] show that the RVE size of composites reinforced with the random distribution of particles is increased if the geometrical shape becomes more complex. The development of polycrystalline modeling of material behavior allows fully multiscale predictive models in which the purpose is capturing the fundamental aspects of material behavior under heterogeneous deformation [66]. In multi-scale modeling, the finite element method is coupled with discrete methods, such as the phase-field approach, Cellular Automata (CA) [67,68], Digital Material Representation (DMR) [69], Multi-resolution Models (MRM) [70], Molecular Dynamics (MD) [71], or the Monte Carlo (MC) [72,73]. The CA lattice of cells usually represents a 3D microstructure and reproduces topological relations between grain structures. The DMR approach is defined as a material description, based on measurable quantities (i.e., grain boundaries, grain size, inclusions, precipi-tates), which provides the transition between the experiment and numerical modeling ( Figure 5). The DMR approach permits a description of material properties with different micro-scale features. The multi-resolution approach improves the behavior of complex materials through the forming processes. The properties of materials are based on dislocation mechanics, fracture prediction, damage accumulation, void mechanics, dislocation density, and material localization. MD is a simulation method that is applied at the atomistic level to simulate stacking fault tetrahedra of face centered cubic (FCC) metals and alloys [74], Crystal-melt coexistence in FCC and body centered cubic (BCC) metals [75], and plastic damage [76,77].
The MC is an algorithm based on a random sampling of a solutions space for application in physical numerical models. Ou et al. [78] used the MD to quantify probabilistic characteristics of the shape and dimensional errors in net-shape metal-forming processes. This approach can be applied to numerical modeling of metal forming processes for both 2D and 3D cases. The influence of the fidelity of physical phenomena and material properties on computational performance is shown in Figure 6. The nonlinear increase in computational costs with the predictive capability of the numerical models can be revealed [66,79]. Multi-scale models that take into account the microstructure of the material are computationally demanding. The review of machine learning for materials development in metals additive manufacturing (MAM) has been presented by Johnson et al. [80]. The design space of MAM ( Figure 7) spans many engineering disciplines since the material and part are made at the same time according to the Industry 4.0 concept. The pre-build phase consists of the parts, alloys, and manufacturing processes. The final microstructure is dictated by the physics of the process, including solidification, melting, and thermal treatment. Akessa et al. [81] proposed an approach to numerical simulation of additive manufactured parts by adopting methodologies developed for composite materials. Koepf et al. [82] presented the innovation of coupling an FE-based model with a cellular automata (CA) model to predict microstructure evolution during selective electron beam melting (SEBM). The maximum CA cell size is restricted by the size of powder particles and initial grain size. Because the average cell size in the FE mesh exceeds this limit, the cell size of the CA is considerably smaller than the cell size of the FE mesh ( Figure 8). Comparison of experimental and numerical results shows good correlation between the experimentally observed and predicted microstructure of the FE-CA compound. A sequential process, as it might appear in a 3D grain growth simulation using cellular automata (CA), is shown in Figure 9.  In order to improve the formability of metallic materials, servo press machines are increasingly employed in metal forming industry. The mechanical servo-drives offer the flexibility of a hydraulic press with the accuracy, speed and reliability of a mechanical press [84]. The formability of high strength stamping parts in servo press forming can be easily improved based on lubrication conditions, tool shape design, and forming process design. Ma et al. [85] developed the method to measure the coefficient of friction under a pulse servo motion. Newly developed, the nonlinear friction model is used to simulate a deep-drawing process. The validity of the developed nonlinear friction model for servo press simulation has been successfully verified experimentally.
High strength steel is widely used in the fabrication of automotive body parts. The strength enhancement of these steels tends to increase the springback of stamped parts, which causes shape deviations of the drawpieces. To obtain the desired part shape, a CAE system JSTAMP [86] has been developed. Shindo et al. [87] used system JSTAMP to compensate computer aided design (CAD) surfaces of the stamping tools based on the springback prediction. JSTAMP consists of two major algorithms for springback compensation known as Displacement Adjustment Method (DAM) [88] and Spring Forward Method [89]. Shindo et al. [87] used the Yoshida-Uemori kinematic hardening model [90,91] which describes the work hardening stagnation and the transient Bauschinger effect. It was found that the compensation capabilities in the simulation system JSTAMP are powerful methods to compensate the springback. The effectiveness of JSTAMP in the compensation of springback of anisotropic steel sheets has been studied by Ma et al. [86].
To manufacture automotive parts of high strength steels, FE-based simulation has become an indispensable tool for the prediction of fracture behaviors and crack occurrence during metal forming. Ma et al. [92] developed a Digital Image Grid (DIG) method to measure local strains at the cracking position during a uniaxial tensile test. The method consists of a two-camera system for static tensile testing which is controlled by a computer. The local strain distribution was calculated by the tracking of the 3D reconstructed grids during tensile testing process. In the past decade, several fracture criteria have been proposed [93][94][95][96][97][98], and various approaches to the prediction of forming cracks were made. Among them, Takuda et al. [99] have shown that Cockcroft-Latham damage criterion, which integrates the plastic strain and the maximum principal stress, can accurately predict the crack occurrence during sheet metal, and the crack behavior due to tube forming was presented by Ma [100].
In order to improve the accuracy of FE-based simulation, Takada et al. [101] measured local anisotropic parameters based on the plastic strain using the DIG method and incorporated them into Hill's anisotropic yield condition. It was experimentally found that local anisotropic r-values gradually decreased due to the effects of local necking and diffused necking. Combining the measured transient displacement field with the FEM, a measurement-based FEM (M-FEM) was developed by Ma et al. [102] for the computation of ductile damage accumulation and distribution of local stresses and strains. The identified damage limit of advanced high strength steel sheets agreed very well compared with that measured by a conventional press test. Digital image correlation (DIC) technique is an optical method that employs tracking and image registration for accurate 2D and 3D measurements of complex structures. Ramazani et al. [103] and Paul [104] investigated failure modes and deformation characterization in a dual phase steels. Sato et al. [105] measured the strain rate dependent fracture behaviors of advanced high strength steel for automotive structures. To improve the strain measurement technique of the Hopkinson bar method, the DIC analysis method is applied to determine the dynamic stress-strain relations. Coppieters et al. [106] measured the large strains using the DIC technique and identified post-necking strain hardening phenomena of ductile sheet metal. An identified hardening law enables disentangling pre-and post-necking strain hardening behavior.

Boundary Element Method
BEM is one of the numerical methods for solving partial differential equations. Its principle is to reduce a given boundary problem to an equivalent integral equation that is defined at the edge of a given area, which is equivalent to reducing the size of the problem. Unlike the finite element method, the BEM use does not require discretization of the interior of the area, but only its edge ( Figure 10). This is the main advantage of this method. The boundary element method is effective in the discretization of infinite areas [107]. The main disadvantage of this method is the difficulty in applying it to nonlinear problems and inhomogeneous media. BEM is oriented for linear problems in that many linear partial differential operators can be transformed into equivalent integral forms through the use of appropriate Green's functions [108]. The BEM was carried forward to include problems with material nonlinearities, such as viscoplasticity [109] and plasticity [110,111]. Recently, the BEM has been applied to analyze the stress and deformation fields in complicated metal forming problems, including both geometrical and material nonlinearities, as well as frictional interface conditions [112].
Foster et al. [113] applied three-dimensional boundary element analysis to analyze the components in aircraft structures. Two elements, i.e., a countersunk hole ( Figure 11a) and a section of an aircraft rib (Figure 11b), were analyzed. A new, efficient algorithm for accelerating the re-integration of the BEM demonstrated the viability of its application in the aerospace industry. Chintapalli et al. [114] proposed a methodology for the design optimization of a skin-stringer panel of an aircraft. The design of tension panels was conducted based on the linear elastic fracture mechanics, while numerical crack growth analysis was based on the BEM. The BEM allowed the structure of upper skin-stringer panels to be optimized for a minimum mass, while guarding against failures. Citarella et al. [115] applied an innovative dual boundary element method (DBEM) modeling approach to reinforced cracked aeronautic panels. A peculiar original arrangement of the rivet configuration in the two-dimensional DBEM model allows taking into account the transversal rivet stiffness and real in-plane panel stiffness. The dual boundary element method both correctly characterizes the singular stress fields near the crack front and simplifies the meshing process [116,117]. Furthermore, it can be easily used in combination with FEM [118][119][120][121]. Citarella et al. [122] performed three-dimensional mixed-mode crack propagation simulations by means of the DBEM code BEASY and two finite element method-based crack propagation codes: ZENCRACK (ZC) and CRACKTRACER3D (CT3D). The shape of the propagated crack, which is significantly out-of-plane for the shear and torsion loading, was matched very well. A comprehensive applications of DBEM to design of tension panels conducted based on the linear elastic fracture mechanics have been provided by Citarella et al. [115]. Citarella et al. [123] proposed coupled FEM-DBEM method to simulate the crack propagation, in order to take advantage of the main capabilities of the two methods. The previously calculated residual stresses were automatically transferred to a DBEM model of the holed plate and applied on the crack faces during the propagation.
An integrated FEM/BEM approach is an effective tool in the simulation of the cold forging process. Schongen et al. [124] applied a time-efficient method based on an integration of FEM and BEM domains in a numerical simulation designed to calculate tool loads in an extrusion process. Recently, Bäcker et al. [125] produced an integrated FEM/BEM method for the simulation of the deep drawing process and the incremental forming process [126]. The combined FEM/BEM approach enables a transient solution of integrated FEM and BEM domains using one solver. Coupled BEM formulation is an efficient tool for the mechanical modeling of non-homogeneous reinforced structural systems [127]. Neto et al. [128] applied three-dimensional nonlinear BEM formulation in the mechanical analysis of non-homogeneous complex 3D Kelvin-Voigt and Boltzmann reinforced materials. The state of the art on the coupling of FEM and BEM is presented by the authors in [129][130][131][132] and mentions the limitations and advantages of such a type of coupling. Coupled FEM/BEM approaches are commonly used for noise reduction improvements in the automotive industry [133][134][135]. In these applications, FEM models have been analyzed in order to predict the power-train bending modes, frequency response functions, and the BEM model is used to predict the power levels emitted as sound.

Discrete Element Method
The DEM is a particle-scale numerical method for modeling the bulk behavior of granular materials. Discrete models take into account material discontinuities or its fragmentation, treating it as a medium composed of discrete objects. Currently, the structure of the material and the functionality of the elements of this structure, possibly on a small scale, are more and more frequently the object of interest. The ability to design materials and atomic structures increases interest in those mechanics that lead to the modeling of materials at the nanostructural or atomic level, such as nanomechanics. In the discrete element method, a material is modeled as a set of rigid elements (particles) interacting with each other through contact forces. Particles (discrete elements) can be of any shape-(i) in two-dimensional problems: cylinders with a circular base or prisms; and (ii) in threedimensional problems: ellipsoids, polyhedrons, or spheres ( Figure 12). It is also possible to use discrete elements of any irregular shape, defined mathematically or formed by the connection of spherical particles (Figure 13) [136]. DEM is based on the modeling of the motion of granular material as discrete particles of various shapes. Using the soft particle approach, each contact is modeled with a dumper in the contact normal direction and a dumper in the contact tangential direction, as is shown schematically in Figure 14.   In discrete models, the macroscopic behavior of a material is obtained by adopting appropriate models of interaction between specific elements of a discrete model. Modeling at a lower level allows for the development of a constitutive model for the immediately higher level. This is used in multi-scale modeling, which describes the behavior of materials at multiple scales. The constitutive model of a material in the DEM can be treated as a discrete micromechanical model, which is defined by the constitutive relations for the contact interaction between discrete elements [56].
Several contributions on the numerical modeling of Solid Granule Medium Forming (SGMF) using the discrete element method can been found in the literature. Most researchers use LS-DYNA and Abaqus package software, which include the DEM algorithm [137]. In discrete approaches, DEM is a common numerical algorithm that has been widely applied to simulations of granular media. Flow regimes of granular flows can be clearly identified by combining granular temperature and the Savage number [138]. Combination of the discrete element method with other numerical approaches have also been carried out, such as CFD-DEM, FEM-DEM. MFIX is an open-source code developed at the National Energy Technology Laboratory, which is used for multi-phase CFD simulations.
FEM-DEM coupling is becoming the hot-spot for contemporary research. In 2014, Su et al. [139] developed a new algorithm RIGID-II, which is able to detect the interactions of spherical particles with arbitrarily defined shapes. The advantage of the algorithm that has been developed is that it has a simple procedure for the identification of contact states and classification of contact types with a simultaneous reduction in computation time. Wu et al. [140] presented a clump model which is closer to the real particle than a spherical particle. Calculations using DEM showed good agreement with experimental results. Du et al. [141] conducted an FEM-DEM coupling analysis of SGMF tube formed by the outer pressure compression method. The authors developed a conservation principle for systems for solving numerical models of DEM and FEM which involves the discrete element (the particles) and the continuum solid element (the workpiece). The FEM-DEM simulation of coupling is applied to model outer pressure compression forming in a solid granular medium tube ( Figure 15). The non-metallic granules (NMG) which developed by themselves were selected as the granule medium. Experimental validation of the numerical results showed that the algorithm developed was an efficient way to solve the coupling elastic-plastic deformation problems in the contact zone between the continuum and discontinuum. Lattice models permit an analysis of the fracture process at the microstructure level [142], including material anisotropy [143], and they are also used for modeling composites [144].

Meshless Methods
One of the inconveniences of applying the FEM is the need to generate an appropriate mesh of elements which, in the case of complex three-dimensional geometries, requires a high-density mesh in order to obtain suitable computational accuracy. In the analysis of large deformations, the mesh of elements undergoes excessive distortion, which results in the necessity to sometimes undertake multiple adaptive remeshing of the domain. These problems can be overcome by meshless methods [145,146].
Meshless methods include SPH, the mesh-free method [147], the finite point method [148,149], and the Galerkin meshless method [150]. The common feature of the meshless methods is that the local approximation of the function that is being searched for is based only on the values of this function at individual selected points in the area. This approximation does not require the existence of any rigid structure of nodes, nor the definition of any connections between nodes. Micromechanics-based second gradient continuum theory for cohesive granular materials following damage elasticity has been proposed by Yang and Misra [151]. Numerical calculations show that the proposed second gradient continuum model can produce accurate solutions without a priori assumption of the shear band path. Systems derived using granular mechanics are often encountered in many applications of engineering practice [152,153].
In the non-meshed methods, it is also relatively easy to introduce material discontinuities. For these reasons, the free-points method and SPH are used to simulate material failure under shock load. The SPH method was originally developed for discrete problems in astrophysics [154]. Recently, it has been used as a discretization method for continuous problems of high-velocity impacts [155], pseudo-spring SPH simulations on the perforation of metal targets [156], and the metal machining process with Eulerian and Total Lagrangian (TL) SPH [157]. It is often included in a broad class of particle methods. Some of the material particle methods are based on the discrete model and meet the assumptions that allow them to be treated as a special variant of the discrete element method [56].
Three main numerical approaches are used to simulate and study the electrohydraulic forming (EHF) process: the FEM, ALE, and SPH methods. The main differences between the methods of simulating the EHF process used by researchers mainly turn on the method of applying the effect of the electrical discharge into the numerical analysis. FE-based numerical analyses in different software packages have been conducted by, e.g., Melander et al. [158] and Rohatgi et al. [159]. Zahoor and Mousavi [160] simulated the EHF process using the SPH method which benefits from advantageous aspects of both the Lagrangian method and the meshless method. It was concluded that the SPH method overestimates the pressure profile of the shock that is generated and yields larger dome heights in comparison to the ALE method. It also takes a bit longer for the SPH method to complete. Lowering the difference between the numerical and experimental is possible by applying electrical discharge energy into the analysis [160]. As is observed (Figure 16), the overall shape of the shockwave that is generated at exact and identical time-steps for the two different methods is in good agreement. Kumar et al. [161] developed a meshless multi-scale approach to modeling severe plastic deformation of copper in equal-channel angular extrusion. They used a vertical homogenization approach based on a maximum-entropy meshless formulation on the macroscale. The macroscale approach has been combined with enhanced stability to describe quasi finite-strain deformations ( Figure 17). The framework applied allows one to analyze grain refinement and the evolution of both strain and texture. Mesh-free methods may be suitable in a situation where excessive mesh distortion may cause a negative elemental volume leading to inaccurate simulations. Islam and Shaw [156] investigated the behavior of Weldox plate under projectile impact by a pseudospring SPH framework. It was found that the approach used is efficient in capturing any arbitrary crack propagation. The SPH framework equipped with a suitable damage model incorporated in a particle-based impact code can satisfactorily predict plate perforation for plates of different thicknesses at different velocities ( Figure 18). Damage models deals with microcrack formation and corrosion, which imply a loss of material stiffness under loading conditions. Continuum modeling is considered the most feasible approach to deal with complex phenomena of material damaging [162,163]. Material damaging can be viewed at various scales, ranging from the atomic to many intervening scales in the middle, including crystal structures and grain boundaries [164]. Quasi-brittle materials, including high-strength steels, play a crucial role in many aerospace and aircraft applications. A reliable prediction of the failure for quasi-brittle materials has become a subject of interest in experimental studies and computational mechanics [165,166]. A novel coupling between the Total Lagrangian (TL) and Eulerian formulations of the SPH method has been presented by Young et al. [167]. The proposed coupling makes use of the mixed kernel-and-gradient correction scheme, which is applied in the analysis of high-velocity impact in order to improve the consistency of the framework. The following four combinations of particle kernel types were investigated: (i) both the debris and target were modeled using Eulerian particles, TL particles (Figure 19a), and TL particles with a central patch (Figure 19b). In the fourth configuration, both the debris and target were modeled using an approach in which the TL formulation was adaptively converted to a Eulerian kernel type (AdapTLE) (Figure 19c). The AdapTLE novel coupling method demonstrates superior results to either a TL or Eulerian formulation independently. Meshless methods that are based on a direct nodal integration scheme exhibit poor performance in modeling solid mechanics problems [168]. On the other hand, the characteristics of discretization flexibility have made mesh-free methods favorable for modeling the discrete cracks of brittle materials [169]. Wu et al. [170] developed a combined mesh-free continuous-discontinuous technique that supports arbitrary crack initiation and propagation in ductile fracture simulation. The developed new technique concerns an interactive particle enrichment algorithm that describes the formation of macro-cracks and strainmorphed method that models the degradation of ductile metals. The advantage of the interactive particle enrichment algorithm is a particle insertion-deletion scheme that leads to a better presentation of the ductile fracture process.

Summary
Along with the continuous development of manufacturing technology, as well as the rapid development of Information Technology, optimization is also developing rapidly, allowing the development of innovative new processes and improving the quality of products. Intelligent production driven by the idea of Industry 4.0 and the demand for high-quality equipment in the aerospace and automotive industries forces the development of manufacturing techniques evolving towards intelligent manufacturing and ecological production. Due to the key role of design optimization throughout the production cycle, great emphasis was placed on supporting the product design process using virtual tools, such as FEM, Variable Frequency Microwave (VFM), CFD, BEM, SPH, XFEM, ALE, and mesh-free methods.
Large automotive companies have more interest than other enterprises in modeling and simulating all aspects of the production process. These companies have to simulate the operation of production installations in detail because thousands of vehicles are produced each year on such lines. Hence, it is very advisable to have very good optimization in the digital environment. The technology of integrating software with manufacturing processes helps industry and plays an important role in research purposes.
Supporting the manufacturing processes with CAE tools, although quite an expensive process, may contribute to increasing the efficiency of implementing new products in small companies. Large automotive companies are most interested in the possibilities for numerical modeling and simulation of manufacturing processes. The production of components taking place on lines producing thousands of products requires very good optimization of the design process in the digital environment. The virtual model allows one to modernize existing production systems or to implement the manufacturing process without the need to manufacture real tools.
The use of simulation programs to analyze plastic working processes in the Industry 4.0 concept requires a thorough knowledge of the mechanics of continuous media. The development of macroscopic constitutive models of a material must be based on all the physical mechanisms occurring in the material. The development of constitutive models of materials is now focused on a proper understanding of the anisotropy of mechanical properties of polycrystalline materials and multi-layer gradient composite materials. Multiscale approaches that tend to move from macro-to micro-and even nano-parameters become very important in numerical optimization of forming tools. Under these conditions, the dimensions of the optimization task of a fully coupled thermomechanical and multi-scale/cross-scale microstructure increase rapidly. A goal of the physics of strongly heterogeneous materials is to derive their effective properties from knowledge of the distribution of their components and constitutive laws of components. This process is called homogenization. In order to determine the representative material properties in the form of a stiffness matrix, the behavior of the microstructure of the material must be determined under the influence of macroscopic loads. A very good example of such a material is RVE, which is a statistical representation of its surroundings from which one can infer the entire structure of the material of an object.
During the last decade, the coupling of individual computational methods to provide both multi-scale and multi-physics responses have proved to have enormous predictive capabilities. There has been a dynamic development of two-and three-dimensional numerical modeling of the industrial processes using the finite element method, boundary element method, finite difference method, finite volume method, neural networks, multi-grid and mesh-free methods, crystal plasticity finite element, discrete element method, extended finite element, an arbitrary Lagrangian-Eulerian, and cellular automata (CA). The latest achievements presented in this manuscript cover several industrial application scenarios, leading to the so-called Industry 4.0 and the latest research related to the computational methods for a wide range of industrial applications.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.