%~Mouliné par MaN_auto v.0.40.4 (550756fa) 2026-07-13 16:21:35
\documentclass[CRMECA,Unicode,biblatex]{cedram}

\addbibresource{CRMECA_Hure_20260346.bib}

\usepackage{bm}
\usepackage{subcaption}
\usepackage{mathtools}

\captionsetup{subrefformat=parens}

\newcommand{\eg}{e.g.}
\newcommand{\ie}{i.e.}
\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}

\newcommand{\rmin}{\mathrm{in}}
\newcommand{\rmout}{\mathrm{out}}
\newcommand{\rmlower}{\mathrm{lower}}
\newcommand{\rmmicro}{\mathrm{micro}}
\newcommand{\rmvol}{\mathrm{vol}}
\newcommand{\rmeq}{\mathrm{eq}}
\newcommand{\rmsurf}{\mathrm{surf}}
\newcommand{\rminter}{\mathrm{inter}}
\newcommand{\rmt}{\mathrm{t}}

\newcommand{\void}{\,\cdot\,}

\DeclareMathOperator{\vol}{vol}
\DeclareMathOperator{\trace}{tr}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator*{\argmin}{arg\;min}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\bracks}{[}{]}
\DeclarePairedDelimiter{\angles}{\langle}{\rangle}
\renewcommand{\llbracket}{[\mkern -2.7mu[}
\renewcommand{\rrbracket}{]\mkern -2.7mu]}
\newcommand{\bbracks}[1]{\mathopen{\llbracket} #1 \mathclose{\rrbracket}}

%%%---------------------------------------------------------------------------

%%% ENVIRONEMENTS

%%% Env. {system} (mimicks {cases} but only supports a single column)
\makeatletter
\newenvironment{system}{
  \let\@ifnextchar\new@ifnextchar
  \left\lbrace
  \def\arraystretch{2}
  \array {@{}l@{}}
}{\endarray\right.}
\makeatother

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

\title{A coalescence criterion for isotropic porous materials with inhomogeneous yield stress}
\alttitle{Un critère de coalescence pour matériaux poreux isotropes de limite d'élasticité inhomogène}

\author{\firstname{Jérémy} \lastname{Hure}\CDRorcid{0000-0003-0074-9903}\IsCorresp}
\address{Université Paris-Saclay, CEA, Service d’Étude des Matériaux Irradiés, 91191, Gif-sur-Yvette, France}
\email{jeremy.hure@cea.fr}

\author{\firstname{Léo} \lastname{Morin}\CDRorcid{0000-0002-6694-1212}}
\address{Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\address{Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France}
\email{leo.morin@u-bordeaux.fr}

\keywords{\kwd{Porous ductile solids} \kwd{void coalescence} \kwd{limit analysis} \kwd{FFT-based simulations}}
\altkeywords{\kwd{Matériaux poreux} \kwd{coalescence de cavités} \kwd{analyse limite} \kwd{simulations FFT}}

\begin{abstract}
The aim of this work is to study theoretically and numerically the effect of the inhomogeneity of the yield stress on void coalescence for isotropic porous ductile materials. Such inhomogeneities may arise due to strain hardening prior to void coalescence. Sequential limit analysis is applied to a cylindrical unit cell containing a coaxial cylindrical void in a von~Mises matrix material with a radially dependent yield stress. An estimate of the coalescence criterion is obtained for combined tensile and shear loading. The criterion relies on 1D integrals, but two approximations are also provided to obtain analytical expressions. Numerical limit analysis based on FFT simulations is performed to get exact, up to numerical errors, coalescence stresses, and to assess the theoretical expressions. A good agreement is found between the analytical coalescence criterion and the numerical results for a large set of void shape, porosities and yield stress distributions. The coalescence criterion is finally used to assess the effect of strain-hardening on the orientation of void coalescence plane, as well as to describe the full yield locus ---~accounting for both void growth and coalescence~--- of porous isotropic materials, for axisymmetric loading conditions.
\end{abstract}

\begin{altabstract}
L'objectif de ce travail est d'étudier théoriquement et numériquement l'effet de l'inhomogénéité spatiale de la limite d'élasticité, résultant par exemple de l'écrouissage, sur la coalescence de cavités dans les matériaux poreux. L'analyse limite séquentielle est appliquée à une cellule cylindrique contenant une cavité cylindrique coaxiale. La plasticité de la matrice est décrite par le critère de von~Mises en considérant une dépendance radiale de la limite d'élasticité. Un critère de coalescence est obtenu pour des chargements combinant traction et cisaillement. Le critère dépend d'intégrales unidimensionnelles, mais deux approximations sont proposées pour obtenir des expressions analytiques. Des simulations d'analyse limite numérique par méthode FFT sont réalisées pour obtenir le critère exact de coalescence et pour évaluer les expressions théoriques obtenues. Un bon accord est observé entre le critère théorique et les résultats numériques pour une large plage de valeurs de rapports d'aspects de la cavité, de porosités et de distribution spatiale de limite d'élasticité. Le critère de coalescence est finalement utilisé pour évaluer l'effet de l'écrouissage sur l'orientation du plan de coalescence, et pour décrire le critère de plasticité de matériaux poreux isotropes pour des conditions de chargement axisymétriques.
\end{altabstract}

\COI{The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.}

\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\section*{Introduction}

Ductile fracture is the main failure mechanism in metallic alloys at room temperature due to the nucleation, growth, and coalescence of microvoids. These processes are strongly influenced by the loading conditions, particularly the stress triaxiality (the ratio between the mean stress and the von~Mises equivalent stress) and whether the loading is monotonic or cyclic~\cite{benzerga_ductile_2010,pineau_failure_2016,benzerga_ductile_2016}. Void nucleation, that comes mostly from particle debonding or cracking in metallic alloys~\cite{NOELL2023101085}, is followed by growth where void volume fraction increases as a result of plastic flow~\cite{RICE1969201} without interactions between voids. Void coalescence, that corresponds to a localization of plastic strain in the ligament between neighboring voids, has been classified into three main modes which depend on microstructural factors, loading conditions and plastic flow characteristics~\cite{benzerga_micromechanics_2002,pineau_failure_2016}: coalescence by internal necking of the intervoid ligament (which is the most dominant coalescence mode), coalescence by internal shearing and necklace coalescence (or coalescence in columns). In situ experimental observations, \eg, using X-ray tomography~\cite{Maire01012014} and laminography~\cite{DU2025267}, have permitted to assess void growth and coalescence in detail. In the same time, porous unit cell simulations \cite{koplik_void_1988, lecarme_void_2011, YERRA20101016, tekoglu_representative_2014} have allowed a detailed description of void growth and coalescence.

In terms of ductile fracture modeling, the most advanced framework is the one based on the local approach to fracture, relying on a detailed and physically-based description of the local damaged process zone~\cite{mealor}. Homogenization combined with limit-analysis offers a suitable mathematical framework to derive micromechanical models~\cite{benzerga_ductile_2010}. This framework has notably been used to derive \emph{void growth models} where plasticity is diffuse, following Gurson's pioneering work~\cite{gurson_continuum_1977} based on the limit-analysis of a spherical cell containing a spherical void and made of a von~Mises material. The addition of phenomenological modeling of void nucleation and coalescence along with ad hoc parameters led to the Gurson--Tvergaard--Needleman (GTN) model~\cite{TVERGAARD1984157}. This model has known numerous extensions in order to include key features of ductile fracture such as void shape effects~\cite{gologanu_approximate_1993,madou_gurson-type_2012}, plastic anisotropy~\cite{benzerga_plastic_2001}, strain hardening effects~\cite{morin_gurson-type_2017}, crystal plasticity~\cite{paux_approximate_2015,scherer_strain_2021}, among others. The modeling of void growth was subsequently followed by the \emph{modeling of void coalescence}, for which plasticity is localized between neighboring voids. In the seminal contribution of Thomason~\cite{thomason_three-dimensional_1985}, a semi-analytical model of coalescence has been derived, corresponding to the limit-load promoting coalescence by internal necking. This work has then been revisited by the important contribution of Benzerga and Leblond~\cite{benzerga_effective_2014}, who developed the first micromechanical model of coalescence (by internal necking) using explicitly limit-analysis with velocity fields compatible with the kinematics of coalescence. Several extensions have been proposed to account for more realistic situations, including refined velocity fields for flat voids~\cite{morin_coalescence_2015,hure_theoretical_2016}, anisotropy of the matrix~\cite{keralavarma_criterion_2016,gallican_anisotropic_2017}, shear effects~\cite{torki_void_2015,torki_theoretical_2017}, distribution effects~\cite{barrioz_void_2018-1}, unified growth and coalescence~\cite{morin_unified_2016}, necklace coalescence~\cite{torki_approximate_2023}, among others, culminating to date to refined coalescence yield criterion and associated evolution equations~\cite{VIGNESHWARAN2024105804,VIGNESHWARAN2025105973}.

In the vast majority of homogenized models for porous ductile materials, strain hardening is modeled by considering an \textit{average} homogeneous yield stress in the matrix material, as proposed by Gurson~\cite{gurson_continuum_1977}. This simplifying assumption, although relevant in practice for low stress triaxiality or low strain hardening capability, is in contradiction with porous unit cell simulations showing a strong gradient of yield stress for high stress triaxiality / high strain hardening capability~\cite{koplik_void_1988,lecarme_void_2011,tekoglu_representative_2014}. Calibrating the so-called Tvergaard parameters of the GTN model as a function of the hardening modulus~\cite{faleskog,pardoentrends,KANIADAKIS2025106171} has been shown to improve the accuracy of homogenized models regarding the modeling of strain hardening, but is expected to be limited to proportional loading conditions. For void growth, the assumption of homogeneous yield stress has been overcome in~\cite{lacroix_numerical_2016,morin_gurson-type_2017,ROUBAUD2024105114,HURE2026106400} that provide analytical \textit{void growth} yield criteria for spherical and ellipsoidal voids in an inhomogeneous matrix material along with evolution equations for the hardening parameters based on sequential limit analysis~\cite{leblond_classical_2018}. These models allow for a physically-based description of isotropic and kinematic hardening. For void coalescence, the assumption of homogeneous yield stress is all the more inaccurate as significant plastic strain gradient close to cavities is observed in micromechanical simulations~\cite{koplik_void_1988,lecarme_void_2011,tekoglu_representative_2014}, leading to a highly inhomogeneous distribution of hardening parameters (isotropic and/or kinematic) prior to the coalescence regime. This inhomogeneity of hardening is therefore expected to have consequences upon the coalescence criterion. However, only very few works have been devoted to the effect of inhomogeneous yield stress on void coalescence criteria. An early contribution proposed to use the Thomason coalescence criterion, developed for homogeneous yield stress, with an estimation of the local yield stress close to the void surface~\cite{YERRA20101016}. An analytical \textit{void coalescence} yield criterion for spheroidal voids in an inhomogeneous matrix material has been proposed recently~\cite{HURE2026106400}, but limited to tensile loadings.

The aim of this work is thus to derive theoretically and assess numerically a void coalescence criterion for general loading conditions for isotropic porous materials with inhomogeneous yield stress. Sequential limit-analysis~\cite{leblond_classical_2018} is used which permits to incorporate the effects of strain hardening during the limit-analysis as done for void growth models~\cite{morin_gurson-type_2017}. The paper is organized as follows. In Section~\ref{sec:position}, the framework of sequential limit-analysis and the problem considered are presented. The macroscopic plastic dissipation and the associated coalescence criterion are computed in Section~\ref{sec:criterion}. The criterion is assessed and calibrated in Section~\ref{sec:num} with respect to numerical calculations performed using an FFT-based solver \cite{moulinecsuquet}. The implications regarding the orientation of the void coalescence plane and yield locus of porous materials are finally discussed in Section~\ref{sec:discussion}.

\section{Position of the problem}\label{sec:position}
\label{sec2}

Sequential limit analysis is first reviewed in this section to provide the key equations used to evaluate the coalescence criterion. The unit cell and boundary conditions used to describe void coalescence are then detailed.

\subsection{Sequential limit-analysis}

Limit-analysis combined with Hill--Mandel homogenization is a convenient framework to derive yield criteria for porous ductile solids~\cite{benzerga_ductile_2010}. Let us consider an elementary cell $\Omega$ containing a void~$\omega$. The material is assumed to be rigid-plastic (no elasticity) with isotropic hardening which is supposed to obey von~Mises criterion:
\begin{equation}\label{eq:micro_crit}
\phi\parens[\big]{\bm{\sigma}(\mathbf{x})} = \sigma_{\rmeq}^{2}(\mathbf{x}) - \overline{\sigma}^2(\mathbf{x}) \leq 0, \quad 
\quad \forall \mathbf{x} \in \Omega-\omega,
\end{equation}
where $\bm{\sigma}$ is the Cauchy stress tensor, $\overline{\sigma}$ the current (local) yield stress and $\sigma_{\rmeq}$ the von~Mises equivalent stress defined by
\begin{equation}
\sigma_{\rmeq}^{2} = \frac{3}{2}\bm{\sigma}':\bm{\sigma}',
\end{equation}
with $\bm{\sigma}' = \bm{\sigma} - \frac{1}{3}(\trace\bm{\sigma})\,\mathbf{I}$ (where $\mathbf{I}$ is the second-order unit tensor) the deviator of $\bm{\sigma}$. The Prandtl--Reuss flow rule associated to the criterion via the normality property reads
\begin{equation}\label{eq:micro_flowrule}
\mathbf{d} = \dot{\lambda} \frac{\partial \phi}{\partial \bm{\sigma}}(\bm{\sigma}) = 3 \dot{\lambda} \bm{\sigma}',
\end{equation}
where $\dot{\lambda}\geq 0$ is the plastic multiplier.

\emph{Classical limit-analysis} is limited to \emph{elastic-ideal-plastic} materials within a small displace\-ment-small strain (linearized) framework. The macroscopic yield locus can then be determined by fundamental inequality,
\begin{equation}\label{eq:IneqLA}
\bm{\Sigma}:\mathbf{D} \leq \Pi(\mathbf{D}),
\end{equation}
in which the macroscopic stress and strain rate tensors $\bm{\Sigma}$ and $\mathbf{D}$ are defined as the volume averages of their microscopic counterparts $\bm{\sigma}$ and $\mathbf{d}$:
\begin{equation}\label{eq:macro}
\mathbf{D} = \frac{1}{\vol(\Omega)} \int_{\Omega} \mathbf{d} \dif V,
\qquad \bm{\Sigma} = \frac{1}{\vol(\Omega)} \int_{\Omega} \bm{\sigma} \dif V.
\end{equation}
Note that Eq.~\eqref{eq:IneqLA} is valid for arbitrary macroscopic strain rate $\mathbf{D}$. The yield locus is deduced by the parametric equation
\begin{equation}\label{eq:YS}
\bm{\Sigma} = \frac{\partial \Pi}{\partial \mathbf{D}}(\mathbf{D}).
\end{equation}
The macroscopic plastic potential $\Pi(\mathbf{D})$ in Eqs.~\eqref{eq:IneqLA} and~\eqref{eq:YS} is defined by:
\begin{equation}\label{eq:macPD}
\Pi(\mathbf{D}) = \inf_{\mathbf{v} \in \mathcal K(\mathbf{D})} \angles[\big]{\pi (\mathbf{d})}_{\Omega},
\end{equation}
where the notation $\langle \void \rangle_{\Omega}$ stands for volume averaging over the volume $\Omega$. In this definition the set $\mathcal K(\mathbf{D})$ consists of velocity fields $\mathbf{v}$ which are kinematically admissible with $\mathbf{D}$ and verify the property of incompressibility, and $\pi (\mathbf{d})$ is the microscopic plastic potential defined for any traceless $\mathbf{d}$ by the formula
\begin{equation}\label{eq:pid}
\pi (\mathbf{d}) = \sup_{\bm{\sigma}^*\in \mathcal{C}_{\rmmicro}} \bm{\sigma}^*:\mathbf{d},
\end{equation}
where $\mathcal{C}_{\rmmicro}$ is the microscopic convex domain of reversibility. For von~Mises plasticity:
\begin{equation}\label{eq:pid2}
\pi (\mathbf{d}) = \overline{\sigma} d_{\rmeq}
\end{equation}
where the von~Mises equivalent strain is defined as
\begin{equation}\label{eq:deqmises}
d_{\rmeq}^2 = \frac{2}{3} \mathbf{d} : \mathbf{d}.
\end{equation}

\emph{Sequential limit-analysis} extends the methods (and results) of classical limit-analysis by incorporating the effects of strain hardening and geometric changes~\cite{leblond_classical_2018} while disregarding elasticity. The idea is to consider a hardenable material as the sequence of (different) successive rigid-ideal plastic materials. Thus, at a given instant, a hardenable material without elasticity behaves like a rigid-ideal plastic material with some initial (fixed) pre-hardening modifying its yield criterion and flow rule. The instantaneous limit-load can be determined using the theorems of classical limit-analysis. Then, in order to account for changes of the strain hardening and geometry, the local hardening parameters and present positions are then updated approximately using the trial velocity field used in the limit-analysis, integrated in a small time step. Sequential limit-analysis can then be used to describe the effect of pre-hardening at the vicinity of a cavity, due to the void growth phase, upon the coalescence limit-loads.

\subsection{Geometry}

A standard elementary volume representative of the coalescence phase~\cite{benzerga_effective_2014} is considered in this study (Figure~\ref{fig:RVE}): the random distribution of polydisperse voids (Figure~\ref{fig:RVE}\subref{fig:RVE_a}) observed in experiments is approximated as a space filling hexagonal distribution of monodisperse voids (Figure~\ref{fig:RVE}\subref{fig:RVE_b}), the latter being again approximated as a cylindrical elementary cell $\Omega$ containing a cylindrical void $\omega$ (Figure~\ref{fig:RVE}\subref{fig:RVE_c}) on which calculations are performed. The geometry is characterized by several non-dimensional parameters: the \emph{ligament parameter} $\chi \equiv R/L$, the \emph{void aspect ratio} $W \equiv h/R$, the \emph{cell aspect ratio} $\lambda \equiv H/L$, and the \emph{volume fraction of the voided band} $c \equiv h/H = W \chi/\lambda$. In the coalescence regime, plastic flow is confined within the ligaments between cavities~\cite{koplik_void_1988,benzerga_effective_2014}; the modeling of the coalescence state therefore implies that the RVE is composed of a central porous layer with a plastic ligament comprised between two rigid zones. The interface between plastic and rigid zones is denoted by $S_{\rminter}$. The local orthonormal basis associated with the cylindrical coordinates $r$, $\theta$, $z$ is denoted $(\mathbf{e}_{r},\mathbf{e}_{\theta},\mathbf{e}_{z})$ and that associated with the Cartesian coordinates $x_{1}$, $x_{2}$, $x_{3}$ is denoted $(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3})$, with $\mathbf{e}_{3}=\mathbf{e}_{z}$. The yield stress of the matrix material is assumed to have a purely radial dependence $\overline{\sigma}(r)$ in order to model isotropic hardening that occurred during the void growth phase prior to coalescence.

\begin{figure}[!h]
\begin{subcaptiongroup}
\includegraphics[width=.7\linewidth]{dessin}
\phantomcaption\label{fig:RVE_a}
\phantomcaption\label{fig:RVE_b}
\phantomcaption\label{fig:RVE_c}
\end{subcaptiongroup}
\caption{Cylindrical void in a cylindrical unit cell with inhomogeneous yield stress used as an approximation of realistic void distribution.}\label{fig:RVE}
\end{figure}

\subsection{Boundary conditions and velocity fields}

The elementary cylindrical cell shown in Figure~\ref{fig:RVE}\subref{fig:RVE_c} is subjected to the general coalescence loading conditions:
\begin{equation}\label{eq:D}
\bm{D} =
\setlength\arraycolsep{4pt}
\begin{pmatrix}
0 & 0 & D_{31} \\
0 & 0 & D_{32} \\
D_{31} & D_{32} & D_{33}
\end{pmatrix}
\end{equation}
where ${D}_{31}$ and ${D}_{32}$ are the imposed shearing rates and $D_{33}$ is the imposed axial strain rate, with $D_{11} = D_{22} = D_{12} = 0$ due to the presence of the rigid zones (Figure~\ref{fig:RVE}\subref{fig:RVE_c}) characterizing coalescence. Following~\cite{torki_theoretical_2017}, the boundary conditions for the velocity field $\mathbf{v}$ on the elementary cell, compatible with Eq.~\eqref{eq:D} should verify the following constraints:
\begin{equation}\label{eq:conditions_limites}
\begin{cases}
	v_r(L, \theta, z)\mathbf{e}_r + v_\theta(L, \theta, z)\mathbf{e}_\theta = \frac{2z}{c}({D}_{31}\mathbf{e}_1 + {D}_{32}\mathbf{e}_2)
	& (-h \leq z \leq h; \: 0 \leq \theta \leq 2\pi),
\\	v_z(r, \theta, \pm H) = \pm D_{33}H
	& (0 \leq r \leq L; \: 0 \leq \theta \leq 2\pi).
\end{cases}
\end{equation}
Eqs.~\eqref{eq:pid} and~\eqref{eq:pid2} allow using discontinuous velocity fields such that
\begin{equation}\label{eq:condition_saut}
\bbracks{\mathbf{v}} \cdot \mathbf{n} = 0 \quad \forall \mathbf{x} \in S_{\rminter},
\end{equation}
\ie{}, the discontinuity is purely tangential. In that case, the associated equivalent strain rate integrated across the discontinuity is equal to $\abs[\big]{\bbracks{v_{\rmt}(\bm{x})}} \big/ \sqrt{3}$~\cite{benzerga_effective_2014}.

Various velocity fields that are incompressible and compatible with Eqs.~\eqref{eq:conditions_limites} and~\eqref{eq:condition_saut} have been proposed in the literature~\cite{keralavarma_criterion_2016,torki_theoretical_2017}. The simplest velocity field is used here~\cite{torki_theoretical_2017}:
\begin{multline}\label{velo-deq}
\begin{system}
	\displaystyle v_r = \frac{D_{33}}{2c}\parens[\Bigg]{\frac{L^2}{r} - r} + \frac{2z}{c}(D_{31} \cos \theta + D_{32} \sin \theta)
\\	\displaystyle v_\theta = \frac{2z}{c}(-D_{31} \sin \theta + D_{32} \cos \theta)
\\	\displaystyle v_z = \frac{D_{33}}{c}z
\end{system}
\\[-8ex]	\Longrightarrow \quad
\begin{system}
	\displaystyle d_{\rmeq} = \frac{1}{\sqrt{3} c} \sqrt{ D_{33}^2 \parens[\Bigg]{3 + \frac{L^4}{r^4}} + 4 \parens{D_{31}^2 + D_{32}^2}}
\\	\displaystyle \abs[\big]{\bbracks{v_{\rmt}(r)}} = \frac{\abs{D_{33}}}{2c} \parens[\Bigg]{\frac{L^2}{r}-r}.
\end{system}
\end{multline}
Eq.~\eqref{velo-deq} can be used to compute the associated macroscopic plastic dissipation (Eq.~\eqref{eq:macPD}) and the coalescence yield criterion (Eq.~\eqref{eq:YS}). Without loss of generality, $D_{32}$ is set to zero in the following (which corresponds to a rotation of the frame around $\mathbf{e}_3$).

\section{A coalescence criterion for isotropic materials with inhomogeneous yield stress}\label{sec:criterion}
\label{sec3}

A coalescence criterion accounting for inhomogeneous yield stress is derived theoretically in this section based on sequential limit analysis and the velocity field described in Section~\ref{sec2}. The macroscopic plastic dissipation associated with Eq.~\eqref{velo-deq} is first detailed, leading to a coalescence criterion that depends on 1D integral expressions. Two different approximations are then made to obtain two analytical coalescence criteria.

\subsection{Macroscopic plastic dissipation}

As recalled in Section~\ref{sec2}, an upper-bound estimate of the macroscopic plastic dissipation reads
\begin{equation}
\Pi(\mathbf{D}) = \langle \overline{\sigma} d_{\rmeq} \rangle_{\Omega}.
\end{equation}
For velocity fields involving a tangential discontinuity across an interface $S_{\rminter}$, the macroscopic plastic dissipation is written as the sum of volumetric and surface contributions as
\begin{equation}
\Pi = \Pi^{\rmvol} + \Pi^{\rmsurf},
\end{equation}
where
\begin{equation}\label{eq:vr_vtheta_vz}
\begin{system}
	\displaystyle \Pi^{\vol} = \frac{1}{\vol(\Omega)} \int_{\Omega -\omega}\overline{\sigma}(\bm{x})d_{\rmeq}(\bm{x}) \dif V
\\	\displaystyle \Pi^{\rmsurf} = \frac{1}{\vol(\Omega)} \int_{S_{\rminter}} \frac{\overline{\sigma}(\bm{x})}{\sqrt{3}} \abs[\big]{\bbracks{v_{\rmt}(\bm{x})}} \dif S.
\end{system}
\end{equation}
For the geometry (Figure~\ref{fig:RVE}\subref{fig:RVE_c}), purely radial dependence of the yield stress $\overline{\sigma}(r)$ and velocity field (Eq.~\eqref{velo-deq}) considered, the macroscopic plastic dissipation is:
\begin{equation}\label{eq:dissip1}
\Pi = \frac{2}{\sqrt{3} L^2} \int_{R}^L \overline{\sigma}(r) \sqrt{ D_{33}^2 \parens[\Bigg]{3 + \frac{L^4}{r^4}} + 4 D_{31}^2} \, r \dif r
+ \frac{\abs{D_{33}}}{\sqrt{3}hL^2} \int_{R}^{L} \overline{\sigma}(r) \parens[\Bigg]{\frac{L^2}{r}-r} r \dif r.
\end{equation}
As shown in~\cite{torki_theoretical_2017}, both integrals from Eq.~\eqref{eq:dissip1} can be computed analytically for a constant yield stress, albeit leading to complex expressions. In the general case, the first integral is cumbersome due to the coupling of the mechanical loading $\{ D_{33}, D_{31} \}$ and yield stress inhomogeneity $\overline{\sigma}$. The following approximation of the macroscopic plastic dissipation is thus made to uncouple these effects:
\begin{equation}\label{eq:dissip2}
\Pi \approx \frac{2}{\sqrt{3} L^2} \int_{R}^L \overline{\sigma}(r) \bracks[\Bigg]{2 \sqrt{D_{33}^2 + D_{31}^2} + \abs{D_{33}} \parens[\Bigg]{\frac{L^2}{r^2} - 1}} \, r \dif r
+ \frac{\abs{D_{33}}}{\sqrt{3}hL^2} \int_{R}^{L} \overline{\sigma}(r) \parens[\Bigg]{\frac{L^2}{r}-r} \, r \dif r,
\end{equation}
where the integrand of the first integral is replaced by an expression keeping the limits ($r\rightarrow \{0,L\}$) unchanged. Moreover, this approximation is exact for $D_{33}=0$ (only shear). For $D_{31}=0$, the same approximation has been shown in~\cite{HURE2026106400} to be accurate. The change of variable $r = L u $ leads to the final expression of the macroscopic plastic dissipation:
\begin{multline}\label{eq:dissip3}
\Pi \approx \bracks[\Bigg]{\frac{4}{\sqrt{3}} \int_{\chi}^1 \overline{\sigma}(L u) u \dif u} \sqrt{D_{33}^2 + D_{31}^2}
\\	+ \bracks[\Bigg]{\frac{2}{\sqrt{3}} \int_{\chi}^1 \overline{\sigma}(L u) \parens{u^{-1} - u} \dif u + \frac{1}{\sqrt{3} W \chi} \int_{\chi}^1 \overline{\sigma}(L u) \parens{1 - u^2} \dif u} \abs{D_{33}},
\end{multline}
where the effects of mechanical loading $\{ D_{33}, D_{31} \}$ and yield stress inhomogeneity $\overline{\sigma}$ are separated. Eq.~\eqref{eq:dissip3} can be used to assess the coalescence criterion.

\subsection{Coalescence criterion}

\subsubsection{General expression}

The fundamental inequality of limit analysis (Eq.~\eqref{eq:IneqLA}) can be written as
\begin{equation}\label{eq:dissip4a}
\Sigma_{33} D_{33} + 2\Sigma_{31} D_{31} \leq \Sigma_1 \abs{D_{33}} + \Sigma_2 \sqrt{D_{33}^2 + D_{31}^2},
\end{equation}
where
\begin{equation}\label{eq:dissip4b}
\begin{system}
	\displaystyle \Sigma_1 = \bracks[\Bigg]{\frac{2}{\sqrt{3}} \int_{\chi}^1 \overline{\sigma}(L u) \parens{u^{-1} - u} \dif u + \frac{1}{\sqrt{3} W \chi} \int_{\chi}^1 \overline{\sigma}(L u) \parens{1 - u^2} \dif u}
\\	\displaystyle \Sigma_2 = \bracks[\Bigg]{\frac{4}{\sqrt{3}} \int_{\chi}^1 \overline{\sigma}(L u) u \dif u}
\end{system}
\end{equation}
that can be used to estimate the yield locus. For $\braces{ D_{33} \neq 0, \: D_{31}}$, the right-hand side of Eq.~\eqref{eq:dissip4a} is differentiable, leading to:
\begin{equation}
\begin{aligned}
	\Sigma_{33} & = \frac{\partial \Pi}{\partial D_{33}} = \parens[\Bigg]{\frac{\Sigma_2}{\sqrt{1+\xi^2}} + \Sigma_1} \sgn(D_{33}),
\\	\Sigma_{31} & = \frac{1}{2}\frac{\partial \Pi}{\partial D_{31}} = \parens[\Bigg]{\frac{\Sigma_2 \xi}{2\sqrt{1+\xi^2}}} \sgn(D_{33}),
\end{aligned}
\end{equation}
with $\xi = D_{31} / D_{33}$. Eliminating $\xi$ and $\sgn(D_{33})$ leads to:
\begin{equation}\label{eqA1}
\parens[\Bigg]{\frac{\abs{\Sigma_{33}}-\Sigma_1}{\Sigma_2}}^2 + \parens[\Bigg]{2 \frac{\Sigma_{31}}{\Sigma_2}}^2 - 1 = 0
\quad \text{for $\abs{\Sigma_{33}} \geq \Sigma_1$}.
\end{equation}
For $\braces[\big]{ D_{33} = 0, \: D_{31} \neq 0}$, the right-hand side of Eq.~\eqref{eq:dissip4a} is only left and right differentiable with respect to $D_{33}$:
\begin{equation}
\begin{gathered}
	-\Sigma_1 = \parens[\Bigg]{\frac{\partial \Pi}{\partial D_{33}}}_{D_{33}=0^-} \leq \Sigma_{33} \leq \parens[\Bigg]{\frac{\partial \Pi}{\partial D_{33}}}_{D_{33}=0^+} = \Sigma_1,
\\	\Sigma_{31} = \frac{1}{2}\frac{\partial \Pi}{\partial D_{31}} = \frac{\Sigma_2}{2} \sgn(D_{31}),
\end{gathered}
\end{equation}
leading to:
\begin{equation}\label{eqA2}
\parens[\Bigg]{2 \frac{\Sigma_{31}}{\Sigma_2}}^2 - 1 = 0
\quad \text{for $\abs{\Sigma_{33}} \leq \Sigma_1$}.
\end{equation}
Combining Eq.~\eqref{eqA1} and Eq.~\eqref{eqA2} finally leads to the yield locus associated with Eq.~\eqref{eq:dissip4a}:
\begin{equation}\label{eqA3}
\parens[\Bigg]{\frac{ \bracks[\big]{\abs{\Sigma_{33}} - \Sigma_1}^+}{\Sigma_2}}^2 + \parens[\Bigg]{2\frac{\Sigma_{31}}{\Sigma_2}}^2 - 1 = 0
\end{equation}
where $\bracks{x}^+ = \max{\parens{x,0}}$ is the positive part function.

The coalescence criterion defined by Eqs.~\eqref{eq:dissip4b} and~\eqref{eqA3} extends the criteria proposed in~\cite{torki_void_2015,keralavarma_criterion_2016,torki_theoretical_2017} for inhomogeneous porous isotropic materials, and the criterion proposed in~\cite{HURE2026106400} to account for combined tension and shear. The two parameters $\{\Sigma_1, \Sigma_2\}$ are defined by integrals that can easily be computed numerically for a known function $\overline{\sigma}$, \eg, using quadrature rules. However, it may be useful for practical applications to have analytical expressions for $\{\Sigma_1, \Sigma_2\}$ that require approximations for $\overline{\sigma}$.

\subsubsection{Piecewise constant approximation}

Assuming that $\overline{\sigma}$ is a piecewise constant function $ \overline{\sigma}(r_i \leq r < r_{i+1}) = \overline{\sigma}_i$ leads to:
\begin{equation}\label{eq:s1s2approx1}
\begin{system}
	\displaystyle \Sigma_1 = \frac{1}{\sqrt{3}} \sum_{i=1}^{n-1} \overline{\sigma}_i \parens[\Bigg]{2 \ln{\frac{\chi_{i+1}}{\chi_i}} - \chi_{i+1}^2 + \chi_i^2} + \frac{1}{3\sqrt{3} W \chi} \sum_{i=1}^{n-1} \overline{\sigma}_i \parens{3 \chi_{i+1} - 3 \chi_i - \chi_{i+1}^3 + \chi_i^3}
\\	\displaystyle \Sigma_2 = \frac{2}{\sqrt{3}} \sum_{i=1}^{n-1} \overline{\sigma}_i \parens{\chi_{i+1}^2 - \chi_i^2}
\end{system}
\end{equation}
where $\chi_i = r_i / L$. This approximation has been used in~\cite{morin_gurson-type_2017} to derive a Gurson-like criterion for inhomogeneous porous isotropic materials. Of course, for $n \rightarrow +\infty$, Eq.~\eqref{eq:s1s2approx1} tends to Eq.~\eqref{eq:dissip4b}. In practice, the choice of the number of terms of the finite sums $n$ is a trade-off between the required accuracy and the number of internal variables ($\sigma_i$) that will be needed to be stored in a numerical implementation. In order to evaluate the effect of $n$, let us consider a yield stress profile of the form:
\begin{equation}\label{eqsb}
\overline{\sigma} (r) = \sigma_{\rmin} + (\sigma_{\rmout} - \sigma_{\rmin}) \frac{1 - \exp{\bracks[\big]{-(r - R) / \ell_0}}}{1 - \exp{\bracks[\big]{-(L - R)/\ell_0}}}
\end{equation}
where $\sigma_{\rmin}$ and $\sigma_{\rmout}$ correspond to the yield stress at the void surface ($r=R$) and at the cell surface ($r=L$), respectively, and $\ell_0$ sets the typical decay lengthscale. Depending on the choice of $\sigma_{\rmin} / \sigma_{\rmout}$ and $\ell_0 / L$, Eq.~\eqref{eqsb} can describe mildly ($\sigma_{\rmin} / \sigma_{\rmout} \sim 1$, $\ell_0 / L \sim 1$) and strongly ($\sigma_{\rmin} / \sigma_{\rmout} \gg 1$, $\ell_0 / L \ll 1$) inhomogeneous yield stress profile.

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra8}
\caption{}\label{fig:effectn1_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra7}
\caption{}\label{fig:effectn1_b}
\end{subfigure}
\caption{Coalescence criterion defined by Eq.~\eqref{eqA3}, where the parameters $\Sigma_1$ and $\Sigma_2$ are computed using Eq.~\eqref{eq:s1s2approx1} for various values of $n$, \subref{fig:effectn1_a}~with or \subref{fig:effectn1_b}~without shear stress.}\label{fig:effectn1}
\end{figure}

Figure~\ref{fig:effectn1} shows the effect of $n$ on the coalescence stresses for pure tension ($\Sigma_{31}=0$) and combined tensile and shear loading, for a strongly inhomogeneous yield stress profile. A rather large value of $n$ is needed in that particular case to approach the limit $n \rightarrow +\infty$ up to few percents. The minimal value of $n$ needed to have an accuracy lower or equal to 5\% of the limit $n \rightarrow +\infty$ is shown in Figure~\ref{fig:effectn2} for various yield stress inhomogeneities (Figure~\ref{fig:effectn2}\subref{fig:effectn2_a}) and void parameters (Figure~\ref{fig:effectn2}\subref{fig:effectn2_b}). The higher the yield stress gradient / the lower the ligament parameter $\chi$, the higher the required value of $n$.

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=5cm]{gra27}
\caption{$W=1$, $\chi=0.5$}\label{fig:effectn2_a}
\end{subfigure}
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=5cm]{gra28}
\caption{$\sigma_{\rmin} / \sigma_{\rmout} = 4$, $\ell_0=L/4$}\label{fig:effectn2_b}
\end{subfigure}
\caption{Minimal value of $n$ needed to have an accuracy lower or equal to 5\% of the limit $n \rightarrow +\infty$ as a function of \subref{fig:effectn2_a}~yield stress parameters and \subref{fig:effectn2_b}~void parameters.}\label{fig:effectn2}
\end{figure}

\subsubsection{Exponential approximation}

Assuming that $\overline{\sigma}$ is an exponential function $
\overline{\sigma}(Lu) = \sigma_0 g(Lu) = \sigma_0 \bracks[\big]{1 + \alpha \exp{\parens[\big]{-\beta \parens{u \chi^{-1} - 1}}}}$ leads to:
\begin{equation}\label{eq:s1s2approx2}
\begin{system}
	\displaystyle \Sigma_1 = \frac{2}{\sqrt{3}} \sigma_0 \bracks[\big]{I_{1/2}\parens{\chi^3} - I_{1/6}\parens{\chi^3}} + \frac{1}{\sqrt{3} W \chi} \sigma_0 \bracks[\big]{I_{1/3}\parens{\chi^3} - I_{0}\parens{\chi^3}}
\\	\displaystyle \Sigma_2 = \frac{4}{\sqrt{3}} \sigma_0 I_{1/6}\parens{\chi^3}
\end{system}
\end{equation}
with
\[
I_k(f) = \int_{f^{1/3}}^1 \frac{g(Lu)}{u^{6k-2}} \dif u.
\]
The exponential approximation, proposed in~\cite{HURE2026106400}, allows recovering various yield stress functions, \eg, constant, linear, localized, by varying the two parameters $\alpha$ and $\beta$. The two parameters $\{\Sigma_1, \Sigma_2\}$ still rely on integrals $I_k$, but the key point is that these integral functions can be expressed with special functions, \ie{}, can be considered analytic in practice:
\begin{equation}
\begin{split}
I_k(f)
	& = \int_{f^{1/3}}^1 \frac{g(ub)}{u^{6k-2}} \dif u 
\\	& =
\begin{dcases}
	{\frac{1-f^{1-2k}}{3-6k} + \alpha \exp{(\beta)} \parens[\big]{f^{1-2k} E_{6k-2}(\beta) - E_{6k-2}(\beta f^{-1/3})}}
	& \text{for $k > \frac{1}{2}$},
\\	{-\frac{\log{f}}{3} + \alpha \exp{(\beta)} \parens[\big]{E_{1}(\beta) - E_{1}(\beta f^{-1/3})}}
	& \text{for $k = \frac{1}{2}$},
\\	{\frac{1-f^{1-2k}}{3-6k} + \alpha \exp{(\beta)} \beta^{6k-3} f^{1-2k} \parens[\big]{\Gamma_{3-6k}^{\rmlower}(\beta f^{-1/3}) - \Gamma_{3-6k}^{\rmlower}(\beta)}}
	& \text{for $k < \frac{1}{2}$},
\end{dcases}
\end{split}
\end{equation}
with ${E_n(x) = \int_1^{+\infty} \frac{\exp{\parens{- x u}}}{u^n} \dif u}$ the generalized exponential integral function~\cite{expint}, and $\Gamma_n^{\rmlower}(x) = {\int_0^{x} u^{n-1} \exp{\parens{- u}} \dif u}$ the incomplete lower gamma function~\cite{gamma}.

Both approximations, piecewise constant (Eq.~\eqref{eq:s1s2approx1}) and exponential (Eq.~\eqref{eq:s1s2approx2}) may be useful in practice, the former being more general and able to handle complex yield stress profiles at the expense of a higher number of internal variables, the latter being more efficient with only three parameters but limited to smooth yield stress profiles.

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra1}
\caption{}\label{fig1_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra2}
\caption{}\label{fig1_b}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra3}
\caption{$W=1$}\label{fig1_c}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra4}
\caption{$W=3$}\label{fig1_d}
\end{subfigure}
\caption{Coalescence criterion as a function of ligament parameter $\chi$ and void shape $W$ \subref{fig1_a}--\subref{fig1_b}~without shear ($\Sigma_{31} = 0$) and \subref{fig1_c}--\subref{fig1_d}~combined tensile and shear loading. Comparisons between Eq.~\eqref{eqA3} using Eq.~\eqref{eq:s1s2constant} and models proposed in~\cite{torki_void_2015, keralavarma_criterion_2016}.}\label{fig1}
\end{figure}

\subsection{Consistency with other coalescence criteria and calibration}

The criterion proposed aims at extending coalescence criteria already available in the literature for inhomogeneous yield stress. In the limit of homogeneous yield stress ($\overline{\sigma}(r) = \overline{\sigma}$), Eqs.~\eqref{eq:dissip4b},~\eqref{eq:s1s2approx1} and~\eqref{eq:s1s2approx2} lead to:
\begin{equation}\label{eq:s1s2constant}
\begin{system}
	\displaystyle \Sigma_1\parens{\chi,W} = \frac{1}{\sqrt{3}} \overline{\sigma} \parens{2 \ln{ \chi^{-1}} - 1 + \chi^2} + \frac{1}{3\sqrt{3} W \chi} \overline{\sigma} \parens{2 - 3 \chi + \chi^3},
\\	\displaystyle \Sigma_2\parens{\chi} = \frac{2}{\sqrt{3}} \overline{\sigma} \parens{1 - \chi^2}.
\end{system}
\end{equation}
The predictions of Eq.~\eqref{eqA3} using Eq.~\eqref{eq:s1s2constant} are compared to the predictions of the other available coalescence criteria~\cite{torki_void_2015, keralavarma_criterion_2016} in Figure~\ref{fig1}. Except for very small flat voids ($\chi \ll 1$ and $W < 1$), the coalescence criterion defined by Eq.~\eqref{eq:s1s2constant} is outside the ones defined by the two other models. Without shear ($\Sigma_{31} = 0$), the three models give close coalescence stress whatever the values of ligament parameter $\chi$ and void shape $W$ (Figure~\ref{fig1}\subref{fig1_a}--\subref{fig1_b}). For pure shear loading ($\Sigma_{33} = 0$), the three models lead to the same coalescence stress. The main discrepancies are observed for combined tension and shear loading (Figure~\ref{fig1}\subref{fig1_c}--\subref{fig1_d}), as a result of the approximation done in Eq.~\eqref{eq:dissip2}.

The comparisons shown in~\cite{torki_void_2015, keralavarma_criterion_2016} indicate that, while the theoretical coalescence criteria obtained by limit analysis capture qualitatively the trends of numerical results, they do not lead to quantitative predictions. The root cause lies into the choice of the trial velocity field and the various approximations made to compute the plastic dissipation analytically. Calibration of a few heuristic parameters is therefore required. The following form of the coalescence criterion is used:
\begin{equation}\label{eq:ylapprox}
\parens[\Bigg]{\frac{ \bracks[\big]{\abs{\Sigma_{33}} - a \Sigma_1\parens{\chi,W}}^+}{ b\Sigma_2\parens{\chi}}}^2 + \parens[\Bigg]{2\frac{\Sigma_{31}}{b \Sigma_2 \parens{\chi}}}^2 - 1 = 0
\end{equation}
where $\{a,b\}$ are parameters, that may be dependent on $\chi$ and $W$, to be calibrated.

\section{Numerical evaluation of the coalescence criterion}\label{sec:num}

The coalescence criterion is evaluated numerically in this section by performing FFT simulations of porous elastic perfectly plastic materials with inhomogeneous yield stress under combined tension and shear. The geometry used in the simulations is equivalent to the one used in the theoretical analysis in order to evaluate the assumptions made (velocity field and approximation of the macroscopic dissipation).

\begin{figure}[!h]
\begin{subfigure}{.58\linewidth}
\centering
\includegraphics[trim = 80 240 150 300,clip,height=4.6cm]{mesh1}
\end{subfigure}
\hfill
\begin{subfigure}{.38\linewidth}
\centering
\includegraphics[height=4.6cm]{mesh2}
\end{subfigure}
\caption{Representative Volume Element of a hexagonal lattice of cylindrical voids discretized in voxels used in the numerical simulations: blue voxels correspond to voids while red voxels to matrix material.}\label{fig:FFTmesh}
\end{figure}

\subsection{Numerical limit analysis}

Eq.~\eqref{eq:ylapprox} along with Eqs.~\eqref{eq:s1s2approx1} or~\eqref{eq:s1s2approx2} stand as an analytical estimate of the yield locus associated with localization of plasticity in a plane containing a hexagonal distribution of cylindrical voids (Figure~\ref{fig:RVE}\subref{fig:RVE_b}). In order to assess the accuracy of this yield locus / coalescence criterion, numerical simulations are performed considering the corresponding Representative Volume Element (Figure~\ref{fig:FFTmesh}). Only one layer of voids is modeled. The RVE is discretized with voxels (\ie{}, cubic elements). The voxel size is chosen to be one-fifth of the smallest characteristic length among void radius $R$, ligament $L-R$ (Figure~\ref{fig:FFTmesh}) and yield stress profile length $\ell_0$ (Eq.~\eqref{eqsb}). The matrix follows elasticity (with Young's modulus $E$ and Poisson ratio $\nu$) and von~Mises plasticity without hardening, with an inhomogeneous yield stress defined according to Eq.~\eqref{eqsb} where $r$ is the distance from the axis of the closest void. In addition, a large value of yield stress is used for voxels such that $\abs{x_3} > \max{(h,L_1)}$ to ensure that localization occurs in the $\{ \mathbf{e}_{1}, \mathbf{e}_{2} \}$ plane. A typical yield stress profile is shown in Figure~\ref{fig:FFTys}.

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[trim = 150 220 150 280,clip,height=4.5cm]{sigmaY_1}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.5cm]{sigmaY_2}
\end{subfigure}
\caption{Typical yield stress $\overline{\sigma}$ field according to Eq.~\eqref{eqsb}, with $\sigma_{\rmin}=4$, $\sigma_{\rmout}=1$ and $\ell_0 = 0.2 L$.}\label{fig:FFTys}
\end{figure}

Periodic boundary conditions are used with macroscopic strain rate defined by Eq.~\eqref{eq:D}. Small strain simulations are performed to assess only the yield criterion. As the geometry is fixed and no hardening is considered, stress saturates upon loading, which corresponds to the limit load / numerical coalescence criterion characterized by the macroscopic stress $\bm{\Sigma}$ (Eq.~\eqref{eq:macro}). Simulations are performed with the \texttt{AMITEX\_FFTP} solver~\cite{amitex} for various values of the ligament parameter $\chi = R/L = [0.2 - 0.8]$, void shape $W= h/L = [0.5;1;2]$, stress inhomogeneity $\sigma_{\rmin} / \sigma_{\rmout} = [1;2;4;8]$ and $\ell_0 /L= [0.1;0.2;1]$, and loading $D_{31}/D_{33} = [0;1;2;4;10;3;100]$, leading to a total of $1764$ simulations. For all simulations, $D_{32} = 0$ and the elastic parameters are set to $E / \sigma_{\rmout} = 2000$ and $\nu = 0.49$. Note that the elastic constants do not affect the limit load. The coalescence stresses $\{\Sigma_{33}, \Sigma_{31}\} $ obtained numerically are compared to Eq.~\eqref{eq:ylapprox} (with Eq.~\eqref{eq:s1s2approx1} using $n=+\infty$ or with Eq.~\eqref{eq:s1s2approx2} that give the same results) considering $\chi \leftarrow 2 \chi / \sqrt{3}$ due to the cylindrical approximation of the model.

\subsection{Comparisons between analytical predictions and numerical results}

The two parameters $\{a,b\}$ involved in Eq.~\eqref{eq:ylapprox} are first calibrated based on the numerical results for a homogeneous matrix material ($\sigma_{\rmin} / \sigma_{\rmout} = 1$) for $W=1$. Figure~\ref{fig:resu1} shows that an overall good agreement is obtained between the numerical results and the theoretical predictions for $a=0.75$ and $b=0.9$. Note that even with calibration, the model fails to reproduce the shape of the coalescence criterion for shear dominated loading. This has been discussed in~\cite{torki_void_2015} and is related to the discontinuous trial velocity field used in limit analysis. Refined velocity fields allow to get a better agreement~\cite{keralavarma_criterion_2016} but with more complex coalescence criteria. The parameters $\{a,b\}$ are considered constant in the following.

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra5}
\caption{}\label{fig:resu1_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra6}
\caption{}\label{fig:resu1_b}
\end{subfigure}
\caption{Coalescence stress for \subref{fig:resu1_a}~pure tensile loading $D_{31} = 0$ and \subref{fig:resu1_b}~combined tension and shear. Comparison between numerical results and predictions from Eq.~\eqref{eq:ylapprox}.}\label{fig:resu1}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra9}
\caption{$W=0.5$, $\ell_0/L=1$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra10}
\caption{$W=0.5$, $\ell_0/L=0.2$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra11}
\caption{$W=0.5$, $\ell_0/L=0.1$}
\end{subfigure}
\caption{Coalescence stress for pure tensile loading for void shape $W = 0.5$, ligament parameter $\chi \in [0.2:0.8]$, inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ (corresponding to purple, green, blue and yellow lines, respectively) and $\ell_0 / L \in [0.1:1]$. Comparison between numerical results and predictions from Eq.~\eqref{eq:ylapprox}.}\label{fig:resu2a}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra12}
\caption{$W=1$, $\ell_0/L=1$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra13}
\caption{$W=1$, $\ell_0/L=0.2$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra14}
\caption{$W=1$, $\ell_0/L=0.1$}
\end{subfigure}
\caption{Coalescence stress for pure tensile loading for void shape $W = 1$, ligament parameter $\chi \in [0.2:0.8]$, inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ (corresponding to purple, green, blue and yellow lines, respectively) and $\ell_0 / L \in [0.1:1]$. Comparison between numerical results and predictions from Eq.~\eqref{eq:ylapprox}.}\label{fig:resu2b}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra15}
\caption{$W=2$, $\ell_0/L=1$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra16}
\caption{$W=2$, $\ell_0/L=0.2$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra17}
\caption{$W=2$, $\ell_0/L=0.1$}
\end{subfigure}
\caption{Coalescence stress for pure tensile loading for void shape $W = 2$, ligament parameter $\chi \in [0.2:0.8]$, inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ (corresponding to purple, green, blue and yellow lines, respectively) and $\ell_0 / L \in [0.1:1]$. Comparison between numerical results and predictions from Eq.~\eqref{eq:ylapprox}.}\label{fig:resu2c}
\end{figure}

The comparisons between the numerical results and the analytical predictions are shown in Figures~\ref{fig:resu2a}--\ref{fig:resu2c} for tensile loading, for the whole range of void shape $W \in [0.5:2]$, ligament parameter $\chi \in [0.2:0.8]$, inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ and $\ell_0 / L \in [0.1:1]$. Overall, a good agreement is observed for all parameter values, where the largest differences appear for the strongest inhomogeneity ($\sigma_{\rmin} / \sigma_{\rmout} = 8$). These comparisons indicate that considering constant calibration parameters $\{a,b\}$ is enough to recover quantitative predictions. This means that the approximations made in the analytical limit analysis are weakly dependent on the input parameters. In fact, for pure tensile loading ($D_{31}=0$), the true velocity field is almost fully imposed by the geometrical constraints of the \textit{coalescence} deformation mode (Eq.~\eqref{eq:D}). This is shown in Figure~\ref{fig:resu3} where the strain rate fields obtained in the numerical simulations for $\sigma_{\rmin}/\sigma_{\rmout} = 1$ (Figure~\ref{fig:resu3}\subref{fig:resu3_a}) and $\sigma_{\rmin}/\sigma_{\rmout} = 4$ (Figure~\ref{fig:resu3}\subref{fig:resu3_b}) are found to be almost identical.

\begin{figure}[!h]
\begin{subfigure}{.32\linewidth}
\centering
\scalebox{.7}{\raisebox{0.5cm}{\includegraphics[trim = 1018 893 670 500,clip,height=4.2cm]{p_1_5_0b}}}
\caption{}\label{fig:resu3_a}
\end{subfigure}
\hfill
\begin{subfigure}{.32\linewidth}
\centering
\scalebox{.7}{\raisebox{0.5cm}{\includegraphics[trim = 1018 893 670 500,clip,height=4.2cm]{p_4_5_0b}}}
\caption{}\label{fig:resu3_b}
\end{subfigure}
\hfill
\begin{subfigure}{.32\linewidth}
\centering
\scalebox{.7}{\includegraphics[height=5.3cm]{gra29}}
\caption{}\label{fig:resu3_c}
\end{subfigure}
\caption{\subref{fig:resu3_a}--\subref{fig:resu3_b} Strain rate fields obtained in the numerical simulations for $D_{31}=0$, $W=1$ and $\chi=0.5$ for \subref{fig:resu3_a}~$\sigma_{\rmin}/\sigma_{\rmout} = 1$ and \subref{fig:resu3_b}~$\sigma_{\rmin}/\sigma_{\rmout} = 4$. \subref{fig:resu3_c}~Strain rate field considered in limit analysis.}\label{fig:resu3}
\end{figure}

Figures~\ref{fig:resu4a}--\ref{fig:resu4c} show the comparisons between the numerical results and the analytical predictions for combined tension and shear, again for the whole range of inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ and $\ell_0 / L \in [0.1:1]$ and selected values of void shape / ligament parameter $(W,\chi) = \braces[\big]{(0.5,0.3),(1,0.5),(2,0.7)}$. Overall, the model, despite its simplicity, provides predictions in good agreement with the reference results. As expected, the discrepancies increase as the parameters deviate significantly from the ones considered to calibrate the model ($W=1$, $\chi=0.5$, $\sigma_{\rmin}/\sigma_{\rmout}=1$). A better agreement could have been obtained by calibrating the parameters $\{a,b\}$ with respect to $W$, $\chi$, $\sigma_{\rmin}/\sigma_{\rmout}$ and $\ell_0/L$, which has not been attempted here. Note finally that the discrepancy regarding the shape of the coalescence criterion for flat voids ($W=0.5$) observed in Figure~\ref{fig:resu4a} is rooted into the discontinuous velocity field used in limit analysis as already explained in~\cite{torki_void_2015}. Finally, all numerical results are compared to the model predictions in Figure~\ref{fig:resu5} where the coalescence stress obtained numerically is plotted against Eq.~\eqref{eq:ylapprox} with (Figure~\ref{fig:resu5}\subref{fig:resu5_a}) or without shear (Figure~\ref{fig:resu5}\subref{fig:resu5_b}). In both cases, most of the model predictions are within $\pm 10\%$ of the reference results.

\begin{figure}[h!]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.8cm]{gra35}
\caption{}\label{fig:resu5_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.8cm]{gra34}
\caption{}\label{fig:resu5_b}
\end{subfigure}
\caption{Coalescence stresses predicted by Eq.~\eqref{eq:ylapprox} as a function of the coalescence stresses obtained numerically \subref{fig:resu5_a}~with or \subref{fig:resu5_b}~without shear. The solid line corresponds to $y=x$, the dotted lines to $y = (1 \pm 0.10)x$.}\label{fig:resu5}
\end{figure}


\begin{figure}[!h]
\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra18}
\caption{$W=0.5$, $\chi=0.3$, $\ell_0/L=1$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra19}
\caption{[$W=0.5$, $\chi=0.3$, $\ell_0/L=0.2$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra20}
\caption{$W=0.5$, $\chi=0.3$, $\ell_0/L=0.1$}
\end{subfigure}
\caption{Coalescence stresses under combined tension and shear for the whole range of inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ (corresponding to purple, green, blue and yellow lines, respectively) and $\ell_0 / L \in [0.1:1]$ and selected values of void shape / ligament parameter $(W,\chi) = \braces[\big]{(0.5,0.3)}$. Comparison between numerical results and predictions from Eq.~\eqref{eq:ylapprox}.}\label{fig:resu4a}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra21}
\caption{$W=1$, $\chi=0.5$, $\ell_0/L=1$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra22}
\caption{$W=1$, $\chi=0.5$, $\ell_0/L=0.2$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra23}
\caption{$W=1$, $\chi=0.5$, $\ell_0/L=0.1$}
\end{subfigure}
\caption{Coalescence stresses under combined tension and shear for the whole range of inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ (corresponding to purple, green, blue and yellow lines, respectively) and $\ell_0 / L \in [0.1:1]$ and selected values of void shape / ligament parameter $(W,\chi) = \braces[\big]{(1,0.5)}$. Comparison between numerical results and predictions from Eq.~\eqref{eq:ylapprox}.}\label{fig:resu4b}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra24}
\caption{$W=2$, $\chi=0.7$, $\ell_0/L=1$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra25}
\caption{$W=2$, $\chi=0.7$, $\ell_0/L=0.2$}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\linewidth}
\centering
\includegraphics[height=5.6cm]{gra26}
\caption{$W=2$, $\chi=0.7$, $\ell_0/L=0.1$}
\end{subfigure}
\caption{Coalescence stresses under combined tension and shear for the whole range of inhomogeneity of the matrix material $\sigma_{\rmin}/\sigma_{\rmout} \in [1:8]$ (corresponding to purple, green, blue and yellow lines, respectively) and $\ell_0 / L \in [0.1:1]$ and selected values of void shape / ligament parameter $(W,\chi) = \braces[\big]{(2,0.7)}$. Comparison between numerical results and predictions from Eq.~\eqref{eq:ylapprox}.}\label{fig:resu4c}
\end{figure}

\FloatBarrier

The results presented in this section show that the coalescence criterion defined by Eq.~\eqref{eq:ylapprox} along with Eq.~\eqref{eq:s1s2approx1} or~\eqref{eq:s1s2approx2} leads to quantitative predictions with respect to numerical reference results obtained on the same geometry used for the modeling, once calibrated ($a=0.75$ and $b = 0.9$). As shown in previous studies~\cite{benzerga_effective_2014,torki_void_2015}, slight recalibration of the parameters would be needed to handle different void microstructures (\eg, cubic or random), void shape (\eg, ellipsoidal or penny-shaped). This has not been attempted in this study as none of these geometrical models can be considered more representative of the void morphologies observed in real metal alloys during ductile fracture.

In the next section, the consequences of Eq.~\eqref{eq:ylapprox} regarding the orientation of the coalescence plane and how to incorporate Eq.~\eqref{eq:ylapprox} in a general homogenized model for porous materials accounting for strain hardening are discussed.

\section{Discussion}\label{sec:discussion}
\subsection{Effect of strain hardening on coalescence plane}

Obtaining coalescence criteria accounting for general loading conditions, \ie{}, that depends on both normal and shear stresses~\cite{BENZERGA2023105344}, enabled predicting the orientation of the coalescence plane. For isotropic void distributions, \ie{}, where the void microstructure does not impose some preferential directions, the orientation of the coalescence plane corresponds to the one for which coalescence is first attained over all possible angles. This has been done in~\cite{keralavarma_multi-surface_2017} using coalescence criteria assuming homogeneous yield stress. The criterion obtained in this study allows assessing the effect of inhomogeneous yield stress distribution on coalescence plane orientation.

\begin{figure}[!h]
\begin{subcaptiongroup}
\includegraphics[width=.7\linewidth]{dessin2}
\phantomsubcaption\label{fig:local_a}
\phantomsubcaption\label{fig:local_b}
\end{subcaptiongroup}
\caption{\subref{fig:local_a}~Random distribution of voids where the anisotropic yield stress distribution around voids is shown in shaded gray. \subref{fig:local_b}~Modeling where all coalescence planes defined by an angle $\theta$ are geometrically equivalent.}\label{fig:local}
\end{figure}

As an example, let us consider a random distribution (Figure~\ref{fig:local}\subref{fig:local_a}) of voids under axisymmetric stress condition:
\begin{equation}\label{eq:axi}
\bm{\Sigma} = \Sigma_{11}
\setlength\arraycolsep{4pt}
\begin{pmatrix}
1 & 0 & 0 \\
0 & \alpha & 0 \\
0 & 0 & \alpha
\end{pmatrix}
\quad \text{with}
\quad T = \frac{\Sigma_m}{\Sigma_{\rmeq}} = \frac{1+2\alpha}{3(1-\alpha)},
\end{equation}
uniquely defined by the stress triaxiality $T$. As shown in~\cite{KANIADAKIS2025106171,HURE2026106400}, strain hardening may lead to an inhomogeneous yield stress around the voids, with both a radial and angular dependence, as sketched in Figure~\ref{fig:local}. The local yield stress is typically maximal for $\theta = 0$ and minimal for $\theta=\pi/2$, and decays as a function of the radial distance from the void, which is modeled using Eq.~\eqref{eqsb} with:
\begin{equation}\label{eq:sintheta}
\sigma_{\rmin}(\theta) = (\sigma_{\rmin} - \sigma_{\rmout}) \parens[\Bigg]{\frac{1 + \cos{ 2\theta}}{2}} + \sigma_{\rmout}.
\end{equation}
For a coalescence band corresponding to an angle $\theta$, the normal and shear stresses with respect to the band, needed to compute the coalescence criterion, are:
\begin{align}
	\Sigma_{nn}(\theta) & = \Sigma_{11} \parens{\cos^2{\theta} + \alpha \sin^2{\theta}},
\\	\Sigma_{nt}(\theta) & = \Sigma_{11} \parens{1 - \alpha} \cos{\theta} \sin{\theta}.
\end{align}
The orientation of the coalescence plane is finally defined as:
\begin{equation}\label{eq:theta}
\theta\parens{T, \sigma_{\rmin} / \sigma_{\rmout}, \chi,W} = \argmin_{\theta} \, \Sigma_{11}
\quad \text{with}
\quad \parens[\Bigg]{\frac{ \bracks[\big]{\abs{\Sigma_{nn}} - a \Sigma_1\parens{\chi,W}}^+}{ b\Sigma_2\parens{\chi}}}^2 + \parens[\Bigg]{2\frac{\Sigma_{nt}}{b \Sigma_2 \parens{\chi}}}^2 -1 =0.
\end{equation}

The dependence of the orientation angle $\theta$ to the inhomogeneity of the yield stress is shown in Figure~\ref{fig:discu1} for a given void microstructure ($W=1$ and $\chi=0.5$) for two stress triaxiality values. For the reference cases of homogeneous matrix material ($\sigma_{\rmin} / \sigma_{\rmout} = 1$), the angle goes from $\theta = 35^{\circ}$ for $T=1$ to $\theta = 10^{\circ}$ for $T=2$. The angle increases strongly as yield stress inhomogeneity increases, especially with higher values of $\sigma_{\rmin} / \sigma_{\rmout}$. Interestingly, Figure~\ref{fig:discu1}\subref{fig:discu1_b} shows that anisotropic inhomogeneous yield stress can change coalescence mode from internal necking ($\theta \approx 0^{\circ}$) to coalescence in columns ($\theta \approx 90^{\circ}$). Although these results depend on the choice of the angular dependence of the yield stress (Eq.~\eqref{eq:sintheta}), Figure~\ref{fig:discu1} makes clear that strain hardening can have a strong influence on void coalescence, not only by affecting quantitatively the coalescence stresses (Figures~\ref{fig:resu4a}--\ref{fig:resu4c}) but also by changing the coalescence mode, even for an isotropic void distribution and matrix material. The interplay between void microstructure, anisotropy of the matrix material, strain hardening and mechanical loading remains however to be studied to predict quantitatively coalescence plane orientation.

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=5cm]{gra30}
\caption{$W=1$, $\chi=0.5$, $T=1$}\label{fig:discu1_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=5cm]{gra31}
\caption{$W=1$, $\chi=0.5$, $T=2$}\label{fig:discu1_b}
\end{subfigure}
\caption{Orientation of the coalescence plane defined by Eq.~\eqref{eq:theta} as a function of yield stress inhomogeneity (Eqs.~\eqref{eqsb} and~\eqref{eq:sintheta}).}\label{fig:discu1}
\end{figure}

\subsection{Towards a homogenized model for porous materials accounting for strain hardening}

The coalescence criterion proposed in this study can be coupled to the growth criteria already proposed in the literature to get a yield criterion for isotropic porous materials with inhomogeneous yield stress by considering the yield criterion
\begin{equation}\label{eq:ylfinal}
\Phi(\bm{\Sigma},\overline{\sigma},W,\chi) = \max{\parens{\phi_{g}, \phi_{c}}} \leq 0,
\end{equation}
where $\phi_c$ is given by Eq.~\eqref{eq:ylapprox}, and $\phi_g$ is given in~\cite{morin_gurson-type_2017,ROUBAUD2024105114,HURE2026106400}. As an example, Figure~\ref{fig:discu2} shows the yield locus predicted by Eq.~\eqref{eq:ylfinal} for axisymmetric loading conditions (Eq.~\eqref{eq:axi}) for a cubic array of spherical voids for different stress inhomogeneity ($\sigma_{\rmin} / \sigma_{\rmout}$). The different parts of the yield locus, growth ($\Phi = \phi_g$, \cite{leblond_improved_1995, morin_gurson-type_2017}) and coalescence ($\Phi = \phi_c$) are indicated. For the latter, shear assisted coalescence ($\Sigma_{nt} \neq 0$ in Eq.~\eqref{eq:ylapprox}) is also spotted. Figure~\ref{fig:discu2} shows that the yield locus expands as the stress inhomogeneity increases and that, for a given stress triaxiality, stress inhomogeneity can trigger a change of deformation mode (from growth to coalescence; Figure~\ref{fig:discu2}\subref{fig:discu2_a}). One should remark that the yield surfaces in Figure~\ref{fig:discu2} show corners, which are due to the intersection of yield loci obtained using different RVEs and velocity fields.

\begin{figure}[t]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra33}
\caption{$f=0.001$}\label{fig:discu2_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[height=4.7cm]{gra32}
\caption{$f=0.01$}\label{fig:discu2_b}
\end{subfigure}
\caption{Yield loci obtained from Eq.~\eqref{eq:ylfinal} for axisymmetric loading conditions (Eq.~\eqref{eq:axi}) for a cubic array of spherical voids for different stress inhomogeneity ($\sigma_{\rmin} / \sigma_{\rmout})$. Solid lines correspond to $\Phi=\phi_g$ \cite{leblond_improved_1995, morin_gurson-type_2017}, red lines to $\Phi=\phi_c$ (Eq.~\eqref{eq:ylapprox}) with $\Sigma_{nt} \neq 0$ and blue lines to $\Phi=\phi_c$ with $\Sigma_{nt} = 0$.}\label{fig:discu2}
\end{figure}

Eq.~\eqref{eq:ylfinal} should be supplemented by equations for the evolutions of the plastic strain $\bm{E}_p$, yield stress $\overline{\sigma}$, porosity $f$, void shape $W$ and microstructure $\chi$ to get a fully featured homogenized model. Plastic strain rate is:
\begin{equation}
\dot{\bm{E}}_p = \dot{\lambda} \frac{\partial \Phi}{\partial \bm{\Sigma}}
\end{equation}
where the property of normality is preserved during the homogenization procedure~\cite{leblond_classical_2018}. The evolution of porosity is classically obtained following mass conservation:
\begin{equation}
\dot{f} = (1 - f) \trace\dot{\bm{E}}_p.
\end{equation}
Following sequential limit-analysis, yield stress profile $\overline{\sigma}$ can be updated by using the trial velocity field used to derive either the growth or coalescence criterion. Knowing the macroscopic plastic strain rate $\dot{\bm{E}}_p$, the velocity field estimates the microscopic plastic strain rate $\dot{\bm{\varepsilon}}_p$, hence the cumulated microscopic plastic strain rate $\dot{p}$ which is used to update the local yield stress:
\begin{equation}
\overline{\sigma}_i = R\parens[\big]{p_i + \dot{p}_i \parens{\dot{\bm{E}}_p} \dif t}.
\end{equation}
Finally, equations are needed for the evolution of void shape and intervoid spacing. These evolution equations are available for growth and coalescence but applicability for inhomogeneous yield stress remains to be assessed.

\section{Conclusion and perspectives}

A coalescence criterion for isotropic porous materials with inhomogeneous yield stress has been derived theoretically based on limit analysis. The criterion depends on normal and shear stresses with respect to the coalescence plane as well as on void shape, intervoid spacing and yield stress profile. Different versions of this criterion have been proposed, relying on 1D integrals, finite sums or special functions. The criterion has thoroughly been compared to numerical reference results and calibrated to provide predictions within 10\% in most cases. This coalescence criterion is expected to be particularly useful to describe the mechanical behavior of porous materials with strongly hardening matrix material, where strong gradient of local yield stress may be observed close to the void surface. This model allows estimating the strong effect of strain hardening on coalescence plane orientation. The modeling of strain hardening within a coalescence criterion thus appears of prime importance as it modifies both the overall yield stress but also the orientation of the coalescence plane, thus having important consequences upon ductile failure. The coalescence criterion can be combined with growth criterion already available in the literature accounting also for inhomogeneous yield stress, providing a fully featured homogenized model for porous materials suitable to describe strongly hardening materials.

Various perspectives can be foreseen. First, extending the proposed criterion to account for kinematic hardening, as done in~\cite{ROUBAUD2024105114}, is at hand. Second, the model proposed assumes that the voids, coalescence plane and yield stress profile share the same principal axes. Obtaining a fully general coalescence criterion for arbitrary voids and yield stress distribution is therefore a path to pursue, however expected to be very challenging as already intricated for homogeneous yield stress~\cite{VIGNESHWARAN2024105804,VIGNESHWARAN2025105973}. Finally, the full homogenized model for porous materials with strongly hardening matrix material sketched in Section~\ref{sec:discussion} should be implemented numerically and compared to reference porous unit cell results.

\printCOI

\printbibliography

\end{document}