Binder Jetting (BJ) additive manufacturing emerges as a promising forming process for fabricating geometrically complex ceramic components, including thin walls, internal channels, and topology-optimized features. The process achieves forming through the selective deposition of a liquid binder onto a powder bed, typically without involving a melt pool. Consequently, it is inherently compatible with high-melting-point, brittle ceramic powders and is less constrained by melt-related instabilities. Despite these advantages, achieving robust manufacturing remains challenging. Defects and variability rarely originate from an isolated processing step. They emerge from the interaction of a series of strongly coupled steps like powder spreading, binder deposition and infiltration, drying/curing and depowdering, thermal removal of organics (debinding), and sintering. Heterogeneities introduced early in the process chain are often transmitted and amplified, ultimately manifesting as non-uniform shrinkage, warpage, cracking, residual porosity, and property scatter. In this context, full-process numerical modeling is increasingly required, serving for mechanistic interpretation and for predictive process-window design.

This review consolidates recent advances on numerical modeling of ceramic Binder Jetting (BJ) and proposes an organizing principle to enhance the coherence and utility of the literature, framing the process chain through transferable state variables that quantify process inheritance. Modeling approaches are evaluated not as isolated stage-specific problems, but according to three dimensions, i.e., 1) the specific physical quantities predicted; 2) the transferability of these quantities as consistent inputs to subsequent stages; and 3) the practicality of calibrating and validating the relevant quantities through accessible experiments. Key transferable state variables highlighted include the spatial distribution of powder-bed porosity and pore-size statistics; permeability and its directional dependence; local binder saturation and organic distribution in the green body; evolving transport coefficients and gas-generation rates during debinding; and the initial relative-density and defect fields that control densification and shape evolution during sintering. Centering on these variables enables an explicit clarification of inter-stage coupling logic and charts a path toward engineering-grade predictability.

Powder spreading modeling defines the initial mesoscale structure controlling all downstream responses. The discrete element method (DEM) is widely adopted as it enables the direct representation of particle rearrangement, contact interactions, and localized compaction under blades or rollers. The existing studies demonstrate that particle size distribution, morphology, cohesion, layer thickness, and recoater kinematics jointly regulate bed density, surface roughness and structural anisotropy. For full-process simulation, key advantage of DEM converts these attributes into downstream-ready descriptors-including porosity maps, pore-network statistics and anisotropic permeability closures.

Binder deposition and infiltration determine the saturation field and organic distribution, which simultaneously control green body strength and subsequent gas evolution. At the component level, Darcy-type infiltration models remain an effective tool for rapidly exploring printing conditions. When permeability and capillary effects are expressed through measurable parameters (e.g., viscosity, surface tension, and wetting), these models can be enhanced to incorporate drying, evaporation, and curing/solidification kinetics. When pore-scale interfacial dynamics and droplet impact dominate, higher-resolution approaches are required. Interface-capturing computational fluid dynamics (CFD) and lattice Boltzmann method (LBM) simulations have been used to resolve transient interfaces and flow pathways in complex pore geometries, while coupled DEM-CFD frameworks have been introduced to capture infiltration-induced particle rearrangement in fine powders.

Debinding is typically the most failure-prone stage since gases generated by binder decomposition must vent through a pore network while the solid skeleton remains weak and thermally stressed. Consequently, predictive modeling necessitates a coupled treatment of reaction kinetics, transport and mechanical response. Decomposition kinetics is typically constrained through thermogravimetric analysis (TGA), commonly implemented with multi-step schemes that represent complex degradation pathways. Gas transport is treated through diffusion, convection, and pressure-driven Darcy flow, with permeability and effective diffusivity evolving with porosity and temperature. Coupled heat and mass transfer in actual geometries can be solved through the finite element method (FEM) implementations, with a failure risk evaluated by relating pore pressure, strength reduction, and constrained deformation.

Sintering consolidates the debound microstructure and determines densification, grain growth, and macroscopic shape change. Continuum constitutive laws embedded within FEM remain a dominant route for component-scale prediction, providing that temperature-dependent viscosities/creep parameters are calibrated and initial relative-density fields are supplied. These models are essential for predicting anisotropic shrinkage and warpage induced by inherent heterogeneity. In the microstructural scale, the Monte Carlo Potts models and phase-field approaches reveal the mechanisms of neck growth, pore evolution, grain-boundary migration, and coarsening. Despite high computational cost, they provide a basis for reduced-order kinetics and constitutive relations applicable to engineering-scale simulations, linking fundamental mechanisms to manufacturability.

A consistent thread across all stages is that predictive capability relies on effective parameter identification and validation. Accordingly, this review recommends establishing a validation matrix aligned with transferable variables. Powder spreading should be validated via bed density, surface roughness, and anisotropy metrics; infiltration via penetration depth, saturation mapping, and spatial green strength, debinding via multi-rate TGA/DTG constraints, independent characterization of pore-network evolution, and statistical measures of cracking/distortion as well as sintering via shrinkage anisotropy, warpage, final density, and defect populations. Designing validation around variables required for downstream transfer allows uncertainty to be localized and gradually reduced, thus avoiding implicit redistribution across the chain.