Self-consistent clustering analysis for homogenisation of heterogeneous plates

TL;DR

This paper presents a novel reduced-order multiscale homogenisation method for heterogeneous plates, combining self-consistent clustering analysis with the Lippmann-Schwinger equation, achieving high accuracy with significantly reduced computational effort.

Contribution

It introduces a plate-specific SCA scheme with an offline-online strategy and self-consistent updates, enabling efficient multiscale analysis of linear and nonlinear plate problems.

Findings

Matches FFT-based simulation accuracy

Reduces computational cost by over an order of magnitude

Handles both linear and nonlinear plate problems

Abstract

This work introduces a reduced-order model for plate structures with periodic micro-structures by coupling self-consistent clustering analysis (SCA) with the Lippmann-Schwinger equation, enabling rapid multiscale homogenisation of heterogeneous plates. A plate-specific SCA scheme is derived for the first time and features two key elements: (i) an offline-online strategy that combines Green's functions with k-means data compression, and (ii) an online self-consistent update that exploits the weak sensitivity of the reference medium. The framework handles both linear and nonlinear problems in classical plate theory and first-order shear deformation theory, and its performance is verified on linear isotropic perforated plates and woven composites, as well as on non-linear elasto-plastic perforated plates and woven composites with damage. Across all cases the proposed model matches the accuracy of FFT-based direct numerical simulation while reducing computational cost by over an order of magnitude.

Tables

Table 1. Table 1 : Homogenized properties of the symmetric plate ( A i ​ j A_{ij} [MPa ⋅ \cdot m], D i ​ j D_{ij} [N ⋅ \cdot mm]).

	$H=32$		$H=64$		$H=128$		$H=256$		$H=512$		FFT
	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT
$A_{11}$	27.760	28.144	27.648	27.792	27.342	27.390	27.272	27.312	27.215	27.212	27.127	27.127
$A_{12}$	7.636	7.543	7.745	7.494	7.634	7.608	7.637	7.619	7.634	7.622	7.629	7.629
$A_{33}$	8.893	8.817	8.988	8.837	8.916	8.868	8.911	8.873	8.908	8.891	8.890	8.890
$D_{11}$	105.254	104.486	105.189	103.441	104.539	102.835	104.306	102.473	104.066	101.923	103.715	101.320
$D_{12}$	21.649	21.262	21.551	20.963	21.231	20.712	21.203	20.507	21.170	20.405	21.115	20.483
$D_{33}$	44.713	43.373	44.529	43.080	44.060	42.766	43.782	42.338	43.442	41.969	41.974	40.857

Table 2. Table 2 : Homogenized properties of the asymmetric plate ( A i ​ j A_{ij} [MPa ⋅ \cdot m], B i ​ j B_{ij} [kN], D i ​ j D_{ij} [N ⋅ \cdot mm]).

	$H=32$		$H=64$		$H=128$		$H=256$		$H=512$		FFT
	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT
$A_{11}$	37.564	37.710	37.396	37.153	36.980	36.948	36.937	36.943	36.834	36.833	36.737	36.737
$A_{12}$	10.518	10.307	10.408	10.338	10.536	10.477	10.526	10.495	10.528	10.503	10.512	10.512
$A_{33}$	11.923	11.732	12.017	11.806	11.963	11.856	11.984	11.917	11.983	11.939	11.961	11.961
$B_{11}$	1.751	1.729	1.753	1.736	1.750	1.741	1.753	1.747	1.753	1.749	1.752	1.750
$B_{12}$	0.436	0.422	0.447	0.442	0.463	0.457	0.466	0.462	0.469	0.467	0.471	0.470
$B_{33}$	0.580	0.561	0.585	0.566	0.583	0.572	0.585	0.578	0.585	0.581	0.584	0.580
$D_{11}$	556.817	543.396	556.945	540.198	553.332	537.799	551.954	536.508	550.563	533.935	548.522	531.358
$D_{12}$	124.165	121.924	125.059	120.129	122.797	118.647	122.559	118.432	122.269	117.918	122.051	117.618
$D_{33}$	225.397	214.526	225.389	214.012	222.559	213.086	221.852	212.145	220.070	210.259	212.592	204.360

Table 3. Table 3 : Elastic properties of yarn and matrix material [ 30 , 49 ] .

Material	$E_{11}$ (GPa)	$E_{22}$ (GPa)	${\nu }_{12}$	${\nu }_{23}$	$G_{12}$ (GPa)
Yarn	53.120	14.660	0.266	0.268	4.240
Matrix	3.170	-	0.350	-	-

Table 4. Table 4 : Homogenized properties of the twill weave composite material ( A i ​ j A_{ij} [MPa ⋅ \cdot m], D i ​ j D_{ij} [N ⋅ \cdot mm]).

	$H=32$		$H=64$		$H=128$		$H=256$		$H=512$		FFT
	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT
$A_{11}$	6.411	6.237	6.367	6.182	6.327	6.160	6.311	6.155	6.298	6.142	6.288	6.147
$A_{12}$	1.090	1.018	1.072	0.999	1.049	0.989	1.041	0.987	1.034	0.978	1.030	0.985
$A_{33}$	1.189	1.177	1.186	1.172	1.179	1.165	1.175	1.161	1.174	1.158	1.172	1.157
$D_{11}$	58.162	57.146	57.361	56.191	56.636	55.799	56.333	55.722	56.113	55.331	55.898	55.078
$D_{12}$	8.417	8.792	8.346	8.640	8.405	8.613	8.438	8.603	8.441	8.542	8.430	8.517
$D_{33}$	12.481	12.096	12.469	12.050	12.468	12.034	12.464	12.025	12.461	12.021	12.456	12.048

Table 5. Table 5 : Computation time of Example 1–Example 3 (Unit: [s]).

	$H=32$		$H=64$		$H=128$		$H=256$		$H=512$		FFT
	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT	CPT	FoPT
Example 1	0.002	0.004	0.007	0.011	0.055	0.083	0.231	0.341	1.249	1.675	14.304	38.781
Example 2	0.017	0.028	0.055	0.077	0.386	0.565	1.688	2.300	8.730	12.432	24.186	65.460
Example 3	0.006	0.009	0.021	0.030	0.134	0.213	0.543	0.888	2.696	3.980	20.170	30.880

Table 6. Table 6 : Computation time of elasto-plastic symmetrical perforated plate (Unit: [s]).

	$H=32$	$H=64$	$H=128$	$H=256$	$H=512$	$H=1024$	FFT
Uniaxial Tension	0.044	0.111	0.407	1.665	5.144	28.972	70.148
Pure Shear	0.035	0.094	0.348	1.515	4.473	20.821	52.828
Complex Loading	0.159	0.382	1.562	6.479	18.323	112.496	352.827

Table 7. Table 7 : Computation time and error of Example 1 for different discretizations (Error: [%]), Time: [s]).

Resolution	$H=32$		$H=64$		$H=128$		$H=256$		$H=512$		FFT
	Err	Time	Err	Time	Err	Time	Err	Time	Err	Time	Err	Time
80$\times$80	2.410	0.0027	1.924	0.0077	1.555	0.0553	1.133	0.2319	0.747	1.2113	0.129	1.3926
160$\times$160	2.081	0.0026	2.260	0.0074	1.252	0.0548	1.051	0.2298	0.718	1.2211	0.087	2.5978
240$\times$240	2.165	0.0025	2.228	0.0082	1.357	0.0554	1.097	0.2280	0.804	1.1940	0.056	6.16003
320$\times$320	2.138	0.00249	2.323	0.0074	1.215	0.0553	0.999	0.2307	0.753	1.2495	0.028	14.3043
400$\times$400	2.128	0.0026	2.276	0.0078	1.180	0.0547	1.018	0.2329	0.748	1.1895	0	24.9872

Table 8. Table 8 : Computation error comparison between A-1, A-2 and A-3 (Unit: [%]).

	CPT					FoPT
	$H=32$	$H=64$	$H=128$	$H=256$	$H=512$	$H=32$	$H=64$	$H=128$	$H=256$	$H=512$
A-1
$A_{11}$	1.666	1.265	0.829	0.241	0.122	1.656	0.797	0.424	0.113	0.055
$A_{12}$	3.776	2.871	1.423	0.596	0.294	3.575	1.275	0.265	0.303	0.694
$A_{33}$	1.323	0.918	0.666	0.237	0.124	1.675	1.138	0.658	0.365	0.176
$D_{11}$	3.039	2.237	0.873	0.475	0.241	3.258	1.655	1.054	0.494	0.234
$D_{12}$	1.764	0.757	0.288	0.137	0.074	2.913	1.525	1.129	0.606	0.303
$D_{33}$	0.104	0.062	0.048	0.024	0.008	0.477	0.091	0.212	0.270	0.288
Mean Err	1.945	1.352	0.688	0.285	0.144	2.259	1.080	0.624	0.358	0.292
A-2
$A_{11}$	1.276	1.242	0.476	0.206	0.107	1.033	0.570	0.253	0.002	0.107
$A_{12}$	4.600	4.034	1.323	0.468	0.215	2.996	1.454	0.432	0.431	0.858
$A_{33}$	1.572	1.150	0.644	0.305	0.162	1.635	1.268	0.734	0.422	0.174
$D_{11}$	3.589	2.618	1.026	0.485	0.252	3.528	2.020	0.938	0.433	0.200
$D_{12}$	0.462	0.989	0.388	0.028	0.025	2.873	1.446	0.598	0.306	0.098
$D_{33}$	0.144	0.105	0.065	0.039	0.017	0.247	0.012	0.129	0.223	0.259
Mean Err	1.941	1.690	0.654	0.255	0.130	2.052	1.128	0.514	0.303	0.283
A-3
$A_{11}$	1.943	1.242	0.615	0.358	0.155	1.473	0.570	0.225	0.134	0.081
$A_{12}$	5.836	4.034	1.874	1.047	0.357	3.427	1.454	0.418	0.286	0.668
$A_{33}$	1.402	1.150	0.528	0.254	0.108	1.701	1.268	0.650	0.315	0.106
$D_{11}$	4.052	2.618	1.322	0.779	0.385	3.755	2.020	1.308	1.169	0.458
$D_{12}$	0.154	0.989	0.288	0.100	0.141	3.229	1.446	1.126	1.013	0.294
$D_{33}$	0.199	0.105	0.095	0.065	0.045	0.393	0.012	0.122	0.198	0.231
Mean Err	2.264	1.690	0.787	0.434	0.198	2.329	1.128	0.641	0.519	0.306

Table 9. Table 9 : Material properties of resin and yarn phases.

Symbol	Value	Description
$E_r$	3.11 GPa	Elastic modulus of the resin
$v_r$	0.36	Poisson’s ratio of the resin
${\sigma }^{yl}$	0.0713 GPa	Yield stress of the resin
${\beta }_r$	3.8	Damage evolution coefficient of the resin
${\rho }_r$	1.3 g/cm3	Density of the resin
$E_{y,11}$	128.64 GPa	Longitudinal elastic modulus of the yarn
$E_{y,22}=E_{y,33}$	9.56 GPa	Transversal elastic modulus of the yarn
$v_{y,12}=v_{y,13}$	0.28	Poisson’s ratio of yarns on the transversal plane
$v_{y,23}$	0.33	Poisson’s ratio of yarns out of the transversal plane
$G_{y,12}=G_{y,13}$	4.47 GPa	Out-of-plane shear modulus
$G_{y,23}$	3.45 GPa	In-plane shear modulus
$X^T,X^c$	2710, 1630 GPa	Longitudinal tensile/compressive strength of yarns
$Y^T,Y^c$	0.059, 0.260 GPa	Transversal tensile/compressive strength of yarns
$S^{xy},S^{yz}$	0.094, 0.080 GPa	Shear strengths of the yarn material
${\beta }_{\tau }^f,{\beta }_c^f$	9.0, 7.6	Fibre damage evolution coefficients
${\beta }_{\tau }^m,{\beta }_c^m$	4.7, 4.7	Matrix damage evolution coefficients
${\beta }_s$	2.8	Shearing damage evolution coefficient
${\rho }_y$	1.6 g/cm3	Density of the yarn

Equations

$$\{\begin{array}{rl}& \nabla \cdot \mathbf{N}(\mathbf{x})=0,\mathbf{x}\in \Omega \\ & \nabla \cdot \left[\nabla \cdot \mathbf{M}(\mathbf{x})\right]=0,\mathbf{x}\in \Omega \\ & \mathbf{N}(\mathbf{x})=\mathbb{A}(\mathbf{x}):\mathbf{\epsilon }(\mathbf{x})+\mathbb{B}(x):\mathbf{\varphi }(\mathbf{x}),\mathbf{x}\in \Omega \\ & \mathbf{M}(\mathbf{x})=\mathbb{B}(\mathbf{x}):\mathbf{\epsilon }(\mathbf{x})+\mathbb{D}(\mathbf{x}):\mathbf{\varphi }(\mathbf{x}),\mathbf{x}\in \Omega \\ & \mathbf{\epsilon }(\mathbf{x})=\nabla \mathbf{u}(\mathbf{x}),\mathbf{x}\in \Omega \\ & \mathbf{\varphi }(\mathbf{x})=\mathbf{\Delta }\mathbf{w}(\mathbf{x}),\mathbf{x}\in \Omega \end{array}$$

$$\{\begin{array}{rl}& \nabla \cdot \mathbf{N}\left(\mathbf{x}\right)=0,\mathbf{x}\in \Omega \\ & \nabla \cdot \mathbf{M}\left(\mathbf{x}\right)-\mathbf{Q}\left(\mathbf{x}\right)=0,\mathbf{x}\in \Omega \\ & \nabla \cdot \mathbf{Q}\left(\mathbf{x}\right)=0,\mathbf{x}\in \Omega \\ & \mathbf{N}\left(\mathbf{x}\right)=\mathbb{A}\left(\mathbf{x}\right):\mathbf{\epsilon }\left(\mathbf{x}\right)+\mathbb{B}\left(\mathbf{x}\right):\mathbf{\varphi }\left(\mathbf{x}\right),\mathbf{x}\in \Omega \\ & \mathbf{M}\left(\mathbf{x}\right)=\mathbb{B}\left(\mathbf{x}\right):\mathbf{\epsilon }\left(\mathbf{x}\right)+\mathbb{D}\left(\mathbf{x}\right):\mathbf{\varphi }\left(\mathbf{x}\right),\mathbf{x}\in \Omega \\ & \mathbf{Q}\left(\mathbf{x}\right)=\mathbf{S}\left(\mathbf{x}\right)\cdot \mathbf{\gamma }\left(\mathbf{x}\right),\mathbf{x}\in \Omega \\ & \mathbf{\epsilon }\left(\mathbf{x}\right)=\nabla \mathbf{u}\left(\mathbf{x}\right),\mathbf{x}\in \Omega \\ & \mathbf{\varphi }\left(\mathbf{x}\right)=\nabla \mathbf{\theta }\left(\mathbf{x}\right),\mathbf{x}\in \Omega \\ & \mathbf{\gamma }\left(\mathbf{x}\right)=\nabla \mathbf{w}\left(\mathbf{x}\right)+\mathbf{\theta }\left(\mathbf{x}\right),\mathbf{x}\in \Omega \end{array}$$

$$\{\begin{array}{rl}& u_i=u_i^{\ast }+{\bar{\epsilon }}_{ij}{x}_j,\text{and}i\in 1,2,3,j\in 1,2,3\\ & {\theta }_i={\theta }_i^{\ast }+{\bar{\varphi }}_{ij}{\theta }_j,\text{and}i\in 1,2,j\in 1,2\end{array}$$

$$\{\begin{array}{rl}& u_i=u_i^{\ast }+{\bar{\epsilon }}_{ij}{x}_j,i\in 1,2,j\in 1,2,u_i^{\ast }\#\\ & u_3=u_3^{\ast }+{\bar{\gamma }}_{3j}{x}_j,j\in 1,2,u_3^{\ast }\#\\ & {\theta }_i={\theta }_i^{\ast }+{\bar{\varphi }}_{ij}{\theta }_j,i\in 1,2,j\in 1,2,{\theta }_i^{\ast }\#\end{array}$$

$$\{\begin{array}{rl}& {\mathbf{\tau }}^{(N)}\left(\mathbf{x}\right)=\mathbf{N}\left(\mathbf{x}\right)-\left[{\mathbb{A}}^0:\mathbf{\epsilon }\left(\mathbf{x}\right)+{\mathbb{B}}^0:\mathbf{\varphi }\left(\mathbf{x}\right)\right],\mathbf{x}\in \Omega \\ & {\mathbf{\tau }}^{(M)}\left(\mathbf{x}\right)=\mathbf{M}\left(\mathbf{x}\right)-\left[{\mathbb{B}}^0:\mathbf{\epsilon }\left(\mathbf{x}\right)+{\mathbb{D}}^0:\mathbf{\varphi }\left(\mathbf{x}\right)\right],\mathbf{x}\in \Omega \\ & {\mathbf{\tau }}^{(Q)}\left(\mathbf{x}\right)=\mathbf{Q}\left(\mathbf{x}\right)-{\mathbf{S}}^0\cdot \mathbf{\gamma }\left(\mathbf{x}\right),\mathbf{x}\in \Omega \end{array}$$

$$\{\begin{array}{rl}& \mathbf{\epsilon }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathbb{G}}^{0\left(NN\right)}\left(\mathbf{x}-\mathbf{y}\right):{\mathbf{\tau }}^{\left(N\right)}\left(\mathbf{y}\right)d\mathbf{y}-\bar{\mathbf{\epsilon }}=0\\ & \mathbf{\varphi }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathbb{G}}^{0\left(MM\right)}\left(\mathbf{x}-\mathbf{y}\right):{\mathbf{\tau }}^{\left(M\right)}\left(\mathbf{y}\right)d\mathbf{y}-\bar{\mathbf{\varphi }}=0\end{array}$$

$$\{\begin{array}{rl}& \mathbf{\epsilon }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathbb{G}}^{0\left(NN\right)}\left(\mathbf{x}-\mathbf{y}\right):{\mathbf{\tau }}^{\left(N\right)}\left(\mathbf{y}\right)d\mathbf{y}-\bar{\mathbf{\epsilon }}=0\\ & \mathbf{\varphi }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathbb{G}}^{0\left(MM\right)}\left(\mathbf{x}-\mathbf{y}\right):{\mathbf{\tau }}^{\left(M\right)}\left(\mathbf{y}\right)d\mathbf{y}+{\int }_{\Omega }{\mathcal{G}}^{0\left(MQ\right)}\left(\mathbf{x}-\mathbf{y}\right):{\mathbf{\tau }}^{\left(Q\right)}\left(\mathbf{y}\right)d\mathbf{y}-\bar{\mathbf{\varphi }}=0\\ & \mathbf{\gamma }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathcal{G}}^{0\left(QM\right)}\left(\mathbf{x}-\mathbf{y}\right):{\mathbf{\tau }}^{\left(M\right)}\left(\mathbf{y}\right)d\mathbf{y}+{\int }_{\Omega }{\mathbf{G}}^{0\left(QQ\right)}\left(\mathbf{x}-\mathbf{y}\right):{\mathbf{\tau }}^{\left(Q\right)}\left(\mathbf{y}\right)d\mathbf{y}-\bar{\mathbf{\gamma }}=0\end{array}$$

$$\{\begin{array}{rl}& \Delta \mathbf{\epsilon }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathbb{G}}^{0\left(NN\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{N}\left(\mathbf{y}\right)-{\mathbb{A}}^0:\Delta \mathbf{\epsilon }\left(\mathbf{y}\right)\right]d\mathbf{y}-\Delta \bar{\mathbf{\epsilon }}=0\\ & \Delta \mathbf{\varphi }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathbb{G}}^{0\left(MM\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{M}\left(\mathbf{y}\right)-{\mathbb{D}}^0:\Delta \mathbf{\varphi }\left(\mathbf{y}\right)\right]d\mathbf{y}\\ & +{\int }_{\Omega }{\mathcal{G}}^{0\left(MQ\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{Q}\left(\mathbf{y}\right)-{\mathbb{S}}^0\cdot \Delta \mathbf{\gamma }\left(\mathbf{y}\right)\right]d\mathbf{y}-\Delta \bar{\mathbf{\varphi }}=0\\ & \Delta \mathbf{\gamma }\left(\mathbf{x}\right)+{\int }_{\Omega }{\mathcal{G}}^{0\left(QM\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{M}\left(\mathbf{y}\right)-{\mathbb{D}}^0:\Delta \mathbf{\varphi }\left(\mathbf{y}\right)\right]d\mathbf{y}\\ & +{\int }_{\Omega }{G}^{0\left(QQ\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{Q}\left(\mathbf{y}\right)-{\mathbb{S}}^0\cdot \Delta \mathbf{\gamma }\left(\mathbf{y}\right)\right]d\mathbf{y}-\Delta \bar{\mathbf{\gamma }}=0\end{array}$$

$$\{\begin{array}{rl}& \frac{1}{c^I|\Omega |}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\Delta \mathbf{\epsilon }\left(\mathbf{x}\right)d\mathbf{x}\\ & +\frac{1}{c^I|\Omega |}{\iint }_{\Omega }{\chi }^I\left(\mathbf{x}\right){\mathbb{G}}^{0\left(NN\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{N}\left(\mathbf{y}\right)-{\mathbb{A}}^0:\Delta \mathbf{\epsilon }\left(\mathbf{y}\right)\right]d\mathbf{y}d\mathbf{x}-\Delta \bar{\mathbf{\epsilon }}=0\\ & \frac{1}{c^I|\Omega |}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\Delta \mathbf{\varphi }\left(\mathbf{x}\right)d\mathbf{x}\\ & +\frac{1}{c^I|\Omega |}{\iint }_{\Omega }{\chi }^I\left(\mathbf{x}\right){\mathbb{G}}^{0\left(MM\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{M}\left(\mathbf{y}\right)-{\mathbb{D}}^0:\Delta \mathbf{\varphi }\left(\mathbf{y}\right)\right]d\mathbf{y}d\mathbf{x}\\ & +\frac{1}{c^I|\Omega |}{\iint }_{\Omega }{\chi }^I\left(\mathbf{x}\right){\mathcal{G}}^{0\left(MQ\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{Q}\left(\mathbf{y}\right)-{\mathbb{S}}^0\cdot \Delta \mathbf{\gamma }\left(\mathbf{y}\right)\right]d\mathbf{y}d\mathbf{x}-\Delta \bar{\mathbf{\varphi }}=0\\ & \frac{1}{c^I|\Omega |}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\Delta \mathbf{\gamma }\left(\mathbf{x}\right)d\mathbf{x}\\ & +\frac{1}{c^I|\Omega |}{\iint }_{\Omega }{\chi }^I\left(\mathbf{x}\right){\mathcal{G}}^{0\left(QM\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{M}\left(\mathbf{y}\right)-{\mathbb{D}}^0:\Delta \mathbf{\varphi }\left(\mathbf{y}\right)\right]d\mathbf{y}d\mathbf{x}\\ & +\frac{1}{c^I|\Omega |}{\iint }_{\Omega }{\chi }^I\left(\mathbf{x}\right)G^{0\left(QQ\right)}\left(\mathbf{x}-\mathbf{y}\right):\left[\Delta \mathbf{Q}\left(\mathbf{y}\right)-{\mathbb{S}}^0\cdot \Delta \mathbf{\gamma }\left(\mathbf{y}\right)\right]d\mathbf{y}d\mathbf{x}-\Delta \bar{\mathbf{\gamma }}=0\end{array}$$

$${\chi }^I\left(\mathbf{x}\right)=\{\begin{array}{rl}& 1,\mathbf{x}\in {\Omega }^I\\ & 0,\mathbf{x}\notin {\Omega }^I\end{array}$$

$$\{\begin{array}{rl}& \Delta \mathbf{\epsilon }\left(\mathbf{x}\right)=\sum _{I=1}^H{\chi }^I\left(\mathbf{x}\right)\Delta {\mathbf{\epsilon }}^I\left(\mathbf{x}\right)\\ & \Delta \mathbf{\varphi }\left(\mathbf{x}\right)=\sum _{I=1}^H{\chi }^I\left(\mathbf{x}\right)\Delta {\mathbf{\varphi }}^I\left(\mathbf{x}\right)\\ & \Delta \mathbf{\gamma }\left(\mathbf{x}\right)=\sum _{I=1}^H{\chi }^I\left(\mathbf{x}\right)\Delta {\mathbf{\gamma }}^I\left(\mathbf{x}\right)\end{array}$$

$$\{\begin{array}{rl}& \Delta \mathbf{N}\left(\mathbf{x}\right)=\sum _{I=1}^H{\chi }^I\left(\mathbf{x}\right)\Delta {\mathbf{N}}^I\left(\mathbf{x}\right)\\ & \Delta \mathbf{M}\left(\mathbf{x}\right)=\sum _{I=1}^H{\chi }^I\left(\mathbf{x}\right)\Delta {\mathbf{M}}^I\left(\mathbf{x}\right)\\ & \Delta \mathbf{Q}\left(\mathbf{x}\right)=\sum _{I=1}^H{\chi }^I\left(\mathbf{x}\right)\Delta {\mathbf{Q}}^I\left(\mathbf{x}\right)\end{array}$$

$$\{\begin{array}{rl}& \Delta {\mathbf{\epsilon }}^I+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(NN\right)}\left[\Delta {\mathbf{N}}^J-{\mathbb{A}}^0:\Delta {\mathbf{\epsilon }}^J\right]-\Delta \bar{\mathbf{\epsilon }}=0\\ & \Delta {\mathbf{\varphi }}^I+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(MM\right)}\left[\Delta {\mathbf{M}}^J-{\mathbb{D}}^0:\Delta {\mathbf{\varphi }}^J\right]+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(MQ\right)}\left[\Delta {\mathbf{Q}}^J-{\mathbf{S}}^0\cdot \Delta {\mathbf{\gamma }}^J\right]-\Delta \bar{\mathbf{\varphi }}=0\\ & \Delta {\mathbf{\gamma }}^I+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(QM\right)}\left[\Delta {\mathbf{M}}^J-{\mathbb{D}}^0:\Delta {\mathbf{\varphi }}^J\right]+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(QQ\right)}\left[\Delta {\mathbf{Q}}^J-{\mathbf{S}}^0\cdot \Delta {\mathbf{\gamma }}^J\right]-\Delta \bar{\mathbf{\gamma }}=0\end{array}$$

$$\{\begin{array}{rl}& {\mathbf{D}}^{IJ\left(NN\right)}=\frac{1}{c^I\left|\Omega \right|}{\iint }_{\Omega }{\chi }^I(\mathbf{x}){\chi }^J(\mathbf{x}){\mathbb{G}}^{0\left(NN\right)}(\mathbf{x}-\mathbf{y})dydx\\ & {\mathbf{D}}^{IJ\left(MM\right)}=\frac{1}{c^I\left|\Omega \right|}{\iint }_{\Omega }{\chi }^I(\mathbf{x}){\chi }^J(\mathbf{x}){\mathbb{G}}^{0\left(MM\right)}(\mathbf{x}-\mathbf{y})dydx\\ & {\mathbf{D}}^{IJ\left(MQ\right)}=\frac{1}{c^I\left|\Omega \right|}{\iint }_{\Omega }{\chi }^I(\mathbf{x}){\chi }^J(\mathbf{x}){\mathcal{G}}^{0\left(MQ\right)}(\mathbf{x}-\mathbf{y})dydx\\ & {\mathbf{D}}^{IJ\left(QM\right)}=\frac{1}{c^I\left|\Omega \right|}{\iint }_{\Omega }{\chi }^I(\mathbf{x}){\chi }^J(\mathbf{x}){\mathcal{G}}^{0\left(QM\right)}(\mathbf{x}-\mathbf{y})dydx\\ & {\mathbf{D}}^{IJ\left(QQ\right)}=\frac{1}{c^I\left|\Omega \right|}{\iint }_{\Omega }{\chi }^I(x\mathbf{x}){\chi }^J(\mathbf{x})G^{0\left(QQ\right)}(\mathbf{x}-\mathbf{y})dydx\end{array}$$

$$\{\begin{array}{rl}& \mathbb{A}={\int }_z{\bar{\mathbb{Q}}}_{m}dz=\sum _{l=1}^{n_l}\left[{\bar{\mathbb{Q}}}_m^{(l)}\left({\bar{z}}_l-{\underline{z}}_l\right)\right]\\ & \mathbb{B}={\int }_{z}z{\bar{\mathbb{Q}}}_{m}dz=\frac{1}{2}\sum _{l=1}^{n_l}\left[{\bar{\mathbb{Q}}}_m^{(l)}\left({\bar{z}}_l^2-{\underline{z}}_l^2\right)\right]\\ & \mathbb{D}={\int }_z{z}^2{\bar{\mathbb{Q}}}_{m}dz=\frac{1}{3}\sum _{l=1}^{n_l}\left[{\bar{\mathbb{Q}}}_m^{(l)}\left({\bar{z}}_l^3-{\underline{z}}_l^3\right)\right]\\ & \mathbf{S}={\int }_z{K}_s{\bar{\mathbf{Q}}}_{s}dz=\sum _{l=1}^{n_l}\left[K_s{\bar{\mathbf{Q}}}_s^{(l)}\left({\bar{z}}_l-{\underline{z}}_l\right)\right]\end{array}$$

$$\bar{\mathbf{S}}=\mathbf{T}\mathbf{S}{\mathbf{T}}^T$$

$${\bar{\mathbf{Q}}}_m={\left[\begin{array}{ccc}{\bar{S}}_{11}& {\bar{S}}_{12}& {\bar{S}}_{16}\\ {\bar{S}}_{21}& {\bar{S}}_{22}& {\bar{S}}_{26}\\ {\bar{S}}_{61}& {\bar{S}}_{62}& {\bar{S}}_{66}\end{array}\right]}^{-1},{\bar{\mathbf{Q}}}_s={\left[\begin{array}{cc}{\bar{S}}_{44}& {\bar{S}}_{45}\\ {\bar{S}}_{54}& {\bar{S}}_{55}\end{array}\right]}^{-1}$$

$$x_i=\left(n_i-\frac{1}{2}\right)\frac{L_i}{N_i},i=1,2,3\text{and}{n}_i\in \left[1,N_i\right]$$

$$\{\begin{array}{rl}& {\xi }_1=\frac{N_1}{L_1}\cdot \sin \left({\bar{\xi }}_1\right)\cdot \cos \left({\bar{\xi }}_2\right)\\ & {\xi }_2=\frac{N_2}{L_2}\cdot \cos \left({\bar{\xi }}_1\right)\cdot \sin \left({\bar{\xi }}_2\right)\end{array}$$

$${\bar{\xi }}_i=\frac{\pi }{N_i}\{\begin{array}{rl}& n_i-\frac{N_i-1}{2}-1,\text{if }{N}_i\text{ is odd}\\ & n_i-\frac{N_i}{2}-1,\text{if }{N}_i\text{ is even}\end{array},i=1,2\text{and}{n}_i\in \left[1,N_i\right]$$

$$\epsilon \left(x\right)=\mathbb{O}\left(\mathbf{x}\right):\bar{\epsilon }$$

$$\{\begin{array}{rl}& {\int }_{\Omega }{\chi }^J\left(\mathbf{x}\right){\mathbb{G}}^{0\left(NN\right)}\left(\mathbf{x}-\mathbf{y}\right)dy={\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathbb{G}}}^{0\left(NN\right)}\left(\xi \right)\right]\\ & {\int }_{\Omega }{\chi }^J\left(\mathbf{x}\right){\mathbb{G}}^{0\left(MM\right)}\left(\mathbf{x}-\mathbf{y}\right)dy={\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathbb{G}}}^{0\left(MM\right)}\left(\xi \right)\right]\\ & {\int }_{\Omega }{\chi }^J\left(\mathbf{x}\right){\mathcal{G}}^{0\left(MQ\right)}\left(\mathbf{x}-\mathbf{y}\right)dy={\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathcal{G}}}^{0\left(MQ\right)}\left(\xi \right)\right]\\ & {\int }_{\Omega }{\chi }^J\left(\mathbf{x}\right){\mathcal{G}}^{0\left(QM\right)}\left(\mathbf{x}-\mathbf{y}\right)dy={\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathcal{G}}}^{0\left(QM\right)}\left(\xi \right)\right]\\ & {\int }_{\Omega }{\chi }^J\left(\mathbf{x}\right)G^{0\left(QQ\right)}\left(\mathbf{x}-\mathbf{y}\right)dy={\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathbf{G}}}^{0\left(QQ\right)}\left(\xi \right)\right]\end{array}$$

$$\{\begin{array}{rl}& {\mathbf{D}}^{IJ\left(NN\right)}=\frac{1}{c^I\left|\Omega \right|}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\left\{{\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathbb{G}}}^{0\left(NN\right)}\left(\xi \right)\right]\right\}dx\\ & {\mathbf{D}}^{IJ\left(MM\right)}=\frac{1}{c^I\left|\Omega \right|}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\left\{{\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathbb{G}}}^{0\left(MM\right)}\left(\xi \right)\right]\right\}dx\\ & {\mathbf{D}}^{IJ\left(MQ\right)}=\frac{1}{c^I\left|\Omega \right|}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\left\{{\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathcal{G}}}^{0\left(MQ\right)}\left(\xi \right)\right]\right\}dx\\ & {\mathbf{D}}^{IJ\left(QM\right)}=\frac{1}{c^I\left|\Omega \right|}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\left\{{\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathcal{G}}}^{0\left(QM\right)}\left(\xi \right)\right]\right\}dx\\ & {\mathbf{D}}^{IJ\left(QQ\right)}=\frac{1}{c^I\left|\Omega \right|}{\int }_{\Omega }{\chi }^I\left(\mathbf{x}\right)\left\{{\mathcal{F}}^{-1}\left[{\hat{\chi }}^J\left(\xi \right):{\hat{\mathbf{G}}}^{0\left(QQ\right)}\left(\xi \right)\right]\right\}dx\end{array}$$

$$\{\begin{array}{rl}& {\hat{G}}_{khij}^{0(NN)}=\frac{1}{4{\mu }_m^0}{\hat{G}}_{khij}^{1(NN)}+\frac{{\lambda }_m^0+{\mu }_m^0}{{\mu }_m^0\left({\lambda }_m^0+2{\mu }_m^0\right)}{\hat{G}}_{khij}^{2(NN)}\\ & {\hat{G}}_{khij}^{0(MM)}=\frac{1}{{\mu }_b^0{\left|\xi \right|}^2+G^0}{\hat{G}}_{khij}^{1(MM)}+\frac{1}{\left({\lambda }_b^0+2{\mu }_b^0\right)}{\hat{G}}_{khij}^{2(MM)}\\ & {\hat{G}}_{khi}^{0(MQ)}=\frac{1}{{\mu }_b^0{\left|\xi \right|}^2+G^0}{\hat{G}}_{khi}^{1(MQ)}\\ & {\hat{G}}_{kij}^{0(QM)}=\frac{1}{{\mu }_b^0{\left|\xi \right|}^2+G^0}{\hat{G}}_{kij}^{1(QM)}\\ & {\hat{G}}_{ki}^{0(QQ)}=\frac{1}{{\mu }_b^0{\left|\xi \right|}^2+G^0}{\hat{G}}_{ki}^{1(QQ)}+\frac{1}{G^0}{\hat{G}}_{ki}^{2(QQ)}\end{array}$$

$$\{\begin{array}{rl}& {\hat{G}}_{khij}^{1(NN)}=\frac{{\delta }_{ik}{\xi }_j{\xi }_l+{\delta }_{il}{\xi }_j{\xi }_k+{\delta }_{jl}{\xi }_i{\xi }_k+{\delta }_{jk}{\xi }_i{\xi }_l}{{\left|\xi \right|}^2}{\hat{G}}_{khij}^{2(NN)}=-\frac{{\xi }_i{\xi }_j{\xi }_k{\xi }_l}{{\left|\xi \right|}^4}\\ & {\hat{G}}_{khij}^{1(MM)}=\frac{{\delta }_{ih}{\xi }_j{\xi }_k+{\delta }_{ik}{\xi }_h{\xi }_j}{2}-\frac{{\xi }_h{\xi }_i{\xi }_j{\xi }_k}{{\left|\xi \right|}^2}{\hat{G}}_{khij}^{2(MM)}=\frac{{\xi }_h{\xi }_i{\xi }_j{\xi }_k}{{\left|\xi \right|}^4}\\ & {\hat{G}}_{khi}^{1(MQ)}=\frac{i{\delta }_{ih}{\xi }_k+{\delta }_{ik}{\xi }_h}{2}-\frac{i{\xi }_h{\xi }_i{\xi }_k}{{\left|\xi \right|}^2}\\ & {\hat{G}}_{kij}^{1(QM)}=\frac{i{\xi }_i{\xi }_j{\xi }_k}{{\left|\xi \right|}^2}-i{\delta }_{ik}{\xi }_j\\ & {\hat{G}}_{ki}^{1(QQ)}={\delta }_{ik}-\frac{{\xi }_i{\xi }_k}{{\left|\xi \right|}^2},{\hat{G}}_{ki}^{2(QQ)}=\frac{{\xi }_i{\xi }_k}{{\left|\xi \right|}^2}\end{array}$$

$$\frac{{\rho }^{i+1}}{{\rho }^i}=C$$

$$\frac{{\left({\mu }_b^0{\left|\xi \right|}^2+G^0\right)}^i}{{\left({\mu }_b^0{\left|\xi \right|}^2+G^0\right)}^{i+1}}=C\iff \frac{{\left({\mu }_b^0\right)}^i}{{\left({\mu }_b^0\right)}^{i+1}}=\frac{{\left(G^0\right)}^i}{{\left(G^0\right)}^{i+1}}=C$$

$$\{\begin{array}{rl}& {\mathbf{r}}^{I\left(N\right)}=\Delta {\mathbf{\epsilon }}^I+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(NN\right)}\left[\Delta {\mathbf{N}}^J-{\mathbb{A}}^0:\Delta {\mathbf{\epsilon }}^J\right]-\Delta \bar{\mathbf{\epsilon }}\\ & {\mathbf{r}}^{I\left(M\right)}=\Delta {\mathbf{\varphi }}^I+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(MM\right)}\left[\Delta {\mathbf{M}}^J-{\mathbb{D}}^0:\Delta {\mathbf{\varphi }}^J\right]+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(MQ\right)}\left[\Delta {\mathbf{Q}}^J-{\mathbf{S}}^0\cdot \Delta {\mathbf{\gamma }}^J\right]-\Delta \bar{\mathbf{\varphi }}\\ & {\mathbf{r}}^{I\left(Q\right)}=\Delta {\mathbf{\gamma }}^I+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(QM\right)}\left[\Delta {\mathbf{M}}^J-{\mathbb{D}}^0:\Delta {\mathbf{\varphi }}^J\right]+\sum _{J=1}^H{\mathbf{D}}^{IJ\left(QQ\right)}\left[\Delta {\mathbf{Q}}^J-{\mathbf{S}}^0\cdot \Delta {\mathbf{\gamma }}^J\right]-\Delta \bar{\mathbf{\gamma }}\end{array}$$

$$\{\begin{array}{rl}& {\mathbf{J}}^{\left(N\right)}=\frac{\partial {\mathbf{r}}^{\left(N\right)}\left(\Delta \mathbf{\epsilon }\right)}{\partial \Delta \mathbf{\epsilon }}\\ & {\mathbf{J}}^{\left(M\right)}=\frac{\partial {\mathbf{r}}^{\left(M\right)}\left(\Delta \mathbf{\varphi }\right)}{\partial \Delta \mathbf{\varphi }}\\ & {\mathbf{J}}^{\left(Q\right)}=\frac{\partial {\mathbf{r}}^{\left(Q\right)}\left(\Delta \mathbf{\gamma }\right)}{\partial \Delta \mathbf{\gamma }}\end{array}$$

$$\{\begin{array}{rl}& {\mathbf{J}}^{\left(N\right)}={\mathbf{D}}^{IJ\left(NN\right)}\frac{\partial \Delta \mathbf{N}}{\partial \Delta \mathbf{\epsilon }}+\mathbf{I}-{\mathbf{D}}^{IJ\left(NN\right)}\mathbb{A}\\ & {\mathbf{J}}^{\left(M\right)}={\mathbf{D}}^{IJ\left(MM\right)}\frac{\partial \Delta \mathbf{M}}{\partial \Delta \mathbf{\varphi }}+\mathbf{I}-{\mathbf{D}}^{IJ\left(MM\right)}\mathbb{D}\\ & {\mathbf{J}}^{\left(Q\right)}={\mathbf{D}}^{IJ\left(QQ\right)}\frac{\partial \Delta \mathbf{Q}}{\partial \Delta \mathbf{\gamma }}+\mathbf{I}-{\mathbf{D}}^{IJ\left(QQ\right)}\mathbf{S}\end{array}$$

Full text

Self-consistent clustering analysis for homogenisation of heterogeneous plates

Menglei Li*

National Key Laboratory of Science and Technology on Advanced Composites in Special Environments, Harbin Institute of Technology, Harbin 150001, PR China

Haolin Li*

Department of Aeronautics, Imperial College London, London, SW7 2AZ, UK

Bing Wang†

National Key Laboratory of Science and Technology on Advanced Composites in Special Environments, Harbin Institute of Technology, Harbin 150001, PR China

Bing Wang‡

National Key Laboratory of Science and Technology on Advanced Composites in Special Environments, Harbin Institute of Technology, Harbin 150001, PR China

11footnotetext: Menglei Li and Haolin Li contributed equally to this work.22footnotetext: Third author: wangbing91@hit.edu.cn33footnotetext: Corresponding author, wangbing86@hit.edu.cn

This work introduces a reduced-order model for plate structures with periodic micro-structures by coupling self-consistent clustering analysis (SCA) with the Lippmann–Schwinger equation, enabling rapid multiscale homogenisation of heterogeneous plates. A plate-specific SCA scheme is derived for the first time and features two key elements: (i) an offline–online strategy that combines Green’s functions with k-means data compression, and (ii) an online self-consistent update that exploits the weak sensitivity of the reference medium. The framework handles both linear and non-linear problems in classical plate theory and first-order shear deformation theory, and its performance is verified on linear isotropic perforated plates and woven composites, as well as on non-linear elasto-plastic perforated plates and woven composites with damage. Across all cases the proposed model matches the accuracy of FFT-based direct numerical simulation while reducing computational cost by over an order of magnitude.

Keywords:

multiscale, self-consistent clustering analysis, reduced order model, homogenization, heterogeneous plates

1 Introduction

The growing demand for lightweight and high-performance materials in aerospace, automotive, and biomedical industries has led to major developments in multiscale heterogeneous materials [1, 2]. Engineered materials such as woven composites, porous metamaterials, and thin-walled plates with periodic microstructures achieve high stiffness-to-weight ratios and multiple functions through carefully designed microscopic periodic structures [3, 4]. However, their complex multiscale structures make mechanical characterisation difficult. Classical homogenisation theory, which builds equivalent mechanical models for representative volume elements (RVEs) at the microscale, has become a key approach for connecting the micro- and macro-scale mechanical behaviours in heterogeneous materials [5, 6].

With progress in computational mechanics, numerical solutions to RVE boundary value problems have become the main approach. The finite element method (FEM), as the most established and commonly used numerical tool, remains the standard for homogenisation analysis [7, 8]. However, FEM has limitations when dealing with complex geometries and multiscale problems, as it faces rapidly increasing meshing difficulty and computational cost. The fast Fourier transform (FFT)-based homogenisation method, first proposed by Moulinec and Suquet [9], brought a major change by reformulating the RVE problem into periodic Lippmann-Schwinger equations using the convolution theorem and Green’s operator. This method allows efficient solutions for heterogeneous material unit cells with subgrid-level accuracy, reduces computational cost to $O(N\log N)$, and avoids the need for complex meshing. Over the years, FFT-based methods have been applied successfully to many areas, including composite damage modelling [10], crystal plasticity [11], and multiphysics coupling analysis [12]. A full review of these developments and improvements can be found in [13, 14].

Despite significant progress, FFT-based direct numerical simulation (DNS) remains constrained by voxel discretisation scale, posing challenges for ultra-large-scale and high-fidelity multiscale computations. To alleviate this issue, reduced-order models (ROMs) such as TFA [15], NTFA [16], POD [17], PGD [18], and ECM [19] have been developed to enhance computational efficiency. However, these methods typically rely on empirical assumptions or require extensive prior datasets to ensure predictive accuracy. The self-consistent clustering analysis (SCA) framework proposed by Liu et al. [20] represents a paradigm shift in data-driven reduced-order modelling. By incorporating mechanical response clustering into the FFT framework, SCA constructs cluster-based Lippmann-Schwinger equations that retain the computational accuracy of full-field simulations while achieving the efficiency of traditional ROMs. This cluster-based reduced-order model (CROM) has demonstrated remarkable success in critical applications, including progressive failure analysis [21], multiphysics coupling [22], and multiscale computation [23, 24, 25, 26]. Further advancements, such as virtual cluster analysis (VCA) [27] and finite cluster analysis (FCA) [28], supported by rigorous mathematical convergence proofs [29], underscore the transformative potential of CROM in the field of computational homogenisation [30, 31, 32].

It is worth noting that existing homogenisation techniques predominantly address cell problems formulated within the framework of solid mechanics. However, lightweight engineering structures, such as aircraft skins [33] and automotive body panels [34], exhibit distinct thin-walled plate or shell characteristics that necessitate homogenisation descriptions based on plate theory [35]. Classical plate theory provides an efficient theoretical foundation for plate homogenisation by establishing the mapping between in-plane stress/moment resultants and macroscopic strain/curvature through ABD stiffness matrices. While FEM remains the dominant computational tool in this domain, recent efforts have explored FFT-based alternatives. Notably, Li et al. [36] introduced a Lippmann-Schwinger equation-based plate homogenisation model, extending the FFT method to plate structures. However, this DNS-type approach remains computationally prohibitive for complex geometric configurations and multiscale analysis scenarios. Therefore, the development of CROM-based computational methodologies for plate models represents a crucial breakthrough, offering substantial theoretical and engineering benefits in multiscale computation [37], topology optimisation [38], and material database construction [39, 40, 41].

In this study, we propose an innovative approach that systematically integrates the CROM framework into plate theory to establish a self-consistent clustering analysis methodology for thin-walled structures. A voxel-based finite prism discretisation scheme is employed to enable efficient numerical solutions for complex heterogeneous plate unit cells. The methodological scope is further expanded through the combined use of classical plate theory (CPT) and first-order shear deformation theory (FoPT). To address the coupling challenges between the reference medium and discrete wave vectors in FoPT, we introduce a physically justified decoupling hypothesis that removes redundant computations of interaction tensors, thereby ensuring the robustness and computational efficiency of SCA in moderately thick plate analyses.

The paper is structured as follows: Sec. 2 presents the proposed methodology, including the homogenisation equations based on plate theory, the construction of the cluster-based reduced-order model, and the numerical implementation. Sec. 3 verifies the accuracy of the method using benchmark cases covering linear elasticity, elasto-plasticity, and progressive failure. Sec. 4 analyses the convergence behaviour of reference medium parameters and evaluates the computational efficiency of the approach. Finally, Sec. 5 summarises the key findings and outlines possible directions for future research.

2 Methodology

2.1 Lippmann-Schwinger homogenization for heterogeneous plates

In this study, we assume that the heterogeneous plate is thin and consider only formulations of the Classical Plate Theory (CPT) and the First-Order Shear Deformation Theory (FoPT) [42]. The definitions of the degrees of freedom and the deformation assumptions for CPT and FoPT are shown in Fig. 1. By defining $\Omega$ as the domain of the representative volume element (RVE) and $\bm{x}$ as a material point, the heterogeneous plate satisfies the following governing equation for the CPT problem:

[TABLE]

The heterogeneous plate satisfies the following governing equation for the FoPT problem:

[TABLE]

where $\bm{N}$, $\bm{M}$, and $\bm{Q}$ denote the stress, moment resultant and shear stress, respectively. $\bm{\varepsilon}$, $\bm{\phi}$, and $\bm{\gamma}$ represent the strain, curvature and shear strain, respectively. $\mathbb{A}$, $\mathbb{B}$, $\mathbb{D}$, and $\bm{S}$ are the membrane stiffness, membrane-bending coupling stiffness, bending stiffness, and shear stiffness, respectively. $\bm{u}$, $\bm{\theta}$, and $\bm{w}$ denote the displacement, rotation angle and out-plane displacement, respectively.

The local displacement and rotation angle of the RVE model are decomposed into average and perturbation components:

[TABLE]

Here, ${\bar{\varepsilon}}_{ij}$, ${\bar{\phi}}_{ij}$, ${u}_{i}^{*}$, and ${\theta}_{i}^{*}$ represent the average strain, average curvature, perturbation displacement, and perturbation rotation, respectively. The RVE model for heterogeneous plates satisfies periodic boundary conditions, which differ from those used in solid formulations, as they must also account for the periodicity of the rotational degrees of freedom. The periodic boundary conditions are expressed as follows:

[TABLE]

where $\#$ denotes periodicity.

By introducing a homogeneous isotropic linear elastic reference medium (with stiffness $\mathbb{A}^{0}$, $\mathbb{B}^{0}$, $\mathbb{D}^{0}$, and $\bm{S}^{0}$), the boundary value problems (BVPs) of the RVE model for the heterogeneous plate can be solved using Green’s function. First, the polarisation stress $\bm{\tau}^{(N)}$, polarisation moment $\bm{\tau}^{(M)}$, and polarisation shear stress $\bm{\tau}^{(Q)}$ are introduced as follows:

[TABLE]

The polarisation stress is treated as an external stress on the reference medium, and the original BVPs of CPT can be rewritten as the following Lippmann–Schwinger integral equations:

[TABLE]

The original BVPs of FoPT can be rewritten as the following Lippmann–Schwinger integral equations:

[TABLE]

$\mathbb{G}^{0\left(NN\right)}$, $\mathbb{G}^{0\left(MM\right)}$, $\mathcal{G}^{0\left(MQ\right)}$, $\mathcal{G}^{0\left(QM\right)}$, and $\bm{G}^{0\left(QQ\right)}$ represent different Green’s operators. The Green’s operators for CPT and FoPT have explicit formulations in Fourier space, as detailed in Appendix A.1 .

In addition, the Lippmann–Schwinger equations can be solved iteratively using the fixed-point iteration method [36], which follows the basic framework of the Fourier spectral (FFT) method [9]. This approach is used as the DNS scheme in the present study.

2.2 Solution framework based on clustered reduced-order modeling

Although the solution of Eq.(7) using FFT can effectively improve computational speed, its efficiency still depends directly on the level of discretisation and thus remains within the domain of DNS. This approach continues to face computational challenges when applied to models with complex heterogeneity. Following the concept of cluster-based reduced-order models (CROMs), the RVE model can be decomposed into a finite number of material clusters. By deriving the cluster-based Lippmann–Schwinger equations, the computational degrees of freedom can be significantly reduced, resulting in improved efficiency.

Due to the complex coupling of Green’s operators involved in FoPT, both the discretised Lippmann–Schwinger equation and the self-consistent update scheme differ significantly from the classical self-consistent clustering analysis (SCA) [20]. In contrast, the Lippmann–Schwinger equation in CPT involves only two simple Green’s operators, and its derivation closely follows that of classical SCA. For convenience, the derivation for FoPT is presented below, while the derivation for CPT is provided in Appendix A.2. First, the Lippmann–Schwinger equation, Eq. (7), is written in an incremental form:

[TABLE]

For the continuous form incremental equation, Eq. (8), discretisation is required for numerical solution. It is assumed that the equation has been discretised into regular pixels, and the RVE model is divided into $H$ material clusters. The volume averaging of Eq. (8) is performed within the I-th cluster as follows:

[TABLE]

where $c^{I}$ denotes the volume fraction in the I-th cluster. $\chi^{I}$ is the characteristic function, defined as:

[TABLE]

Assuming that the variable field is uniform in each region, the strain field of the RVE can be regarded as a piecewise linear interpolation of the strain fields within the individual material clusters:

[TABLE]

Similarly, the stress field can be expressed as:

[TABLE]

Substituting the strain and stress fields into Eq. ( 9 ):

[TABLE]

Here, $\bm{D}^{IJ\left(\sim,\sim\right)}$ represents the interaction tensor, which characterizes the interaction between different clusters. The corresponding expression is as follows:

[TABLE]

Therefore, in the cluster-based ROM, the focus is on solving the interaction tensor and the incremental Eq. (13). The computation of Eq. (14) relies on the material cluster information from the initial partition, corresponding to the "offline" stage. The solution of Eq. (13) depends on real-time loading and local material information, corresponding to the "online" stage.

2.3 Numerical implementation

2.3.1 Discretization

The simulation domain consists of a periodic RVE of the microstructure embedded in a cuboidal domain $\Omega$ with dimensions $L_{1}$, $L_{2}$, and $L_{3}$. In solid formulations, the RVE is discretised with a regular array of $N_{1}\times N_{2}\times N_{3}$ voxels, where each voxel belongs to one of the represented phases. When applying plate theory, the 3D solid structure is simplified into a 2D plate. In this study, a voxel-based finite prism approach is employed, where the discretised 3D architecture of the regular voxel array is represented as finite prisms [42], as shown in Fig. 2. In the thickness direction, the mid-plane of the plate is chosen as the origin of the coordinate system. The stiffnesses are calculated along the plate thickness based on classical laminate theory:

[TABLE]

Here, $\bar{\mathbb{Q}}_{m}$ and $\bm{\bar{Q}}_{s}$ represent the in-plane and out-of-plane stiffness in the global coordinate system, respectively. $K_{s}$ denotes the shear correction factor, while $\bar{z}_{l}$ and $\underline{z}_{l}$ represent the coordinates of the upper and lower ends of the l-th layer prism, respectively. The transformation from 3D stiffness in the local coordinate system to 2D stiffness involves a coordinate transformation. Let the 3D stiffness and flexibility matrices in the local coordinate system be denoted by $\bm{Q}$ and $\bm{S}$, respectively. The flexibility matrix in the global coordinate system is obtained by transforming the local flexibility matrix as follows:

[TABLE]

where the transformation matrix $\bm{T}$ is provided in [43]. Then, the 2D stiffness matrix is derived from the 3D flexibility matrix using the following formulation:

[TABLE]

The plate model shares the same in-plane discretisation as the initial 3D model, i.e., a regular $N_{1}\times N_{2}$ array of pixels. The spatial coordinates of the pixels are given by:

[TABLE]

Since the Green’s operator is calculated in Fourier space, discrete frequencies need to be defined. To reduce the Gibbs phenomenon [44], a rotation-based sequence is used to define the discrete frequencies in this work [45]:

[TABLE]

with:

[TABLE]

2.3.2 Offline stage

In the offline stage, the simulation domain is first partitioned into a finite number of material clusters using the K-means clustering algorithm (Fig. 2). The goal of the clustering process is to ensure that the material points within each cluster exhibit similar mechanical behaviour. A simple and effective way to identify material points with similar mechanical behaviour is to perform multiple DNS calculations on the RVE model to obtain the strain concentration tensor $\mathbb{O}$ for each material point:

[TABLE]

In this study, the strain concentration tensor is calculated using FFT simulations under six orthogonal loading conditions. After determining the strain concentration tensor, all material points are clustered into H material clusters using the K-means algorithm [46].

The static clustering-based reduced-order model assumes “similar responses once, similar responses always.” This means that material points with initially similar mechanical behaviour are assumed to maintain similar behaviour throughout different loading stages. Consequently, the material clusters defined based on the initial strain concentration tensor remain unchanged. In addition, according to the formula for calculating the interaction tensor, the interaction tensor depends only on the spatial distribution of the material clusters and the properties of the reference medium. Therefore, the interaction tensor can be pre-calculated after clustering and remains constant during the application of specific loads, i.e., the so-called “offline stage.”

Since the Green’s operators have explicit expressions in Fourier space, the integration by parts for the interaction tensor distribution can be performed efficiently in Fourier space:

[TABLE]

The interaction tensor can be calculated using the following equation:

[TABLE]

Similar to the FoPT, the calculation of the interaction tensor for the CPT is detailed in Appendix A.2.

Additionally, in the case of damage or nonlinear behaviour at local material points within the model, the properties of the reference medium can be adjusted to improve convergence, e.g., by using a self-consistent update scheme. However, the properties of the reference medium directly affect the calculation of the Green’s operators, which in turn influences the calculation of the interaction tensor. According to Eq. (23), computing the interaction tensor in its discrete form is computationally expensive. Frequent re-calculation of the interaction tensor would result in significant computational cost. A feasible approach is to decouple the parameters associated with the reference medium from the discrete frequencies. This allows the computationally intensive tasks related to the discrete frequency terms to be performed in the offline stage, while the online computations focus on updating the material coefficients corresponding to the reference medium. As a result, the Green’s operators and the interaction tensor can be updated efficiently. The separable formulation for the CPT is provided in Appendix A.2.

The separation form of the Green’s operators for the FoPT can be expressed as:

[TABLE]

where:

[TABLE]

For the operators $\hat{G}_{khij}^{0(NN)}$, $\hat{G}_{khij}^{0(MM)}$, $\hat{G}_{khi}^{0(MQ)}$, $\hat{G}_{kij}^{0(QM)}$, and $\hat{G}_{ki}^{0(QQ)}$, the coefficient $\rho=\frac{1}{\mu_{b}^{0}{{\left|\xi\right|}^{2}}+{{G}^{0}}}$ cannot be directly decoupled into discrete frequencies and the material parameters of the reference medium. To enable efficient re-calculation of the Green’s operators when the material parameters of the reference medium are varied, the coefficient $\rho$ is assumed to satisfy the following condition:

[TABLE]

where ${\rho}^{i}$ and ${\rho}^{i+1}$ represent the coefficients before and after the update, respectively, and $C$ is a constant. Under this assumption, the updated coefficient can be obtained by multiplying the previous coefficient by $C$. If this assumption is satisfied, the Lamé constants of the reference medium must meet the following condition:

[TABLE]

This means that the ratio of the Lamé constant ($\mu_{b}^{0}$) before and after the update is equal to the ratio of ${G}^{0}$ before and after the update. However, this assumption is not universally valid when the material parameters of the reference medium change. As discussed in Section 6.1 of Li et al.’s [36] work on solving the FoPT using FFT-based methods, changes in ${G}^{0}$ show strong insensitivity to the convergence of FFT methods. If ${G}^{0}$ exhibits similar insensitivity in cluster-based computations, then the Green’s operators can be efficiently recalculated based on the above assumptions. Further discussion on the sensitivity of ${G}^{0}$ is provided in Sec. 4.1.

The numerical implementation of the offline stage is summarised in Algorithm 1.

2.3.3 Online stage

When the average strain increment is specified, the discrete incremental Lippmann–Schwinger equations, Eq. (13), must be solved online. This process involves updating the constitutive behaviour at local material points based on the calculated response, i.e., the so-called "online stage". Equation (13) can be solved using the Newton–Raphson method, where the residual of the equations is expressed as:

[TABLE]

When the Newton–Raphson iterative method is used to solve the system of equations, the Jacobian matrix must be computed to improve convergence. For the FoPT, the problem is divided into three systems of equations, requiring the calculation of three residuals and their corresponding Jacobian matrices:

[TABLE]

Substituting Eq. ( 28 ) into Eq. ( 29 ) yields:

[TABLE]

where,

[TABLE]

Although the Newton–Raphson iterative method can be employed to solve the system of equations, for continuous Lippmann–Schwinger equations the solution is independent of the choice of reference material, which only affects the convergence speed. However, for cluster-based reduced-order models, the discretisation can be interpreted as a piecewise uniform assumption of material clusters, where the material points in the RVE model cannot strictly satisfy equilibrium conditions. Therefore, the solution of such reduced-order models depends on the choice of the reference medium. Although this error can be mitigated by increasing the number of material clusters, the increase in degrees of freedom consequently raises the computational cost. To balance computational efficiency and accuracy, a self-consistent scheme can be used to update the reference medium, which is also the purpose of separating the reference material parameters and the discrete frequency terms in the Green’s operator, as discussed in Sec. 2.3.2.

In the self-consistent scheme, the stiffness tensor of the reference material should approximate the tangent stiffness of the RVE model:

[TABLE]

Li et al. [36] determined the Lamé constants of the stiffness based on the eigenvalues of the stiffness tensor. However, in the clustering-based model, this approach results in lower accuracy. To achieve a more accurate approximation of the isotropic stiffness, this study adopts a projection-based self-consistent scheme with improved accuracy [28]. The implementation details of this scheme are provided in Appendix A.3. The update of ${{\bar{\bm{S}}}^{0}}$ is governed by ${G}^{0}$ and is only defined for the FoPT. The value of ${G}^{0}$ is determined according to the assumption described in Eq. (27). Additionally, during the offline stage of the FFT calculations, the Lamé constants of the reference material are determined following the procedure described by [36].

The numerical implementation of the online stage is summarised in Algorithm 2.

2.3.4 Non-linear analysis

Two critical challenges emerge in the implementation of the two-stage "offline–online" framework for non-linear analysis:

(i) Error propagation from linear clustering inheritance: Directly transferring the initial linear clustering configuration to the non-linear regime introduces cumulative approximation errors. However, previous studies have shown that the initial clustering distribution retains reasonable accuracy, and for strongly non-linear behaviour, refinement strategies such as sub-clustering can improve accuracy — an aspect beyond the scope of this study.

(ii) Constitutive modelling limitations in plate formulations: Traditional plate models inherently compromise out-of-plane material characterisation through kinematic assumptions, rendering them inadequate for capturing true non-linear constitutive responses in heterogeneous laminates. To address this issue, we implement a nested model order reduction scheme in which:

• •

Material nonlinearities (plasticity/progressive damage) are resolved at 3D material points.

• •

Updated stress states are consistently projected back to the plate model.

The strain ($\bm{\varepsilon}_{s},\bm{\gamma}_{s}$) at a given position in the solid model can be determined from the strain $\bm{\varepsilon}$ and curvature $\bm{\phi}$ of the plate model, along with the distance $z$ from the target position to the mid-plane [42]:

[TABLE]

Furthermore, we employ an incremental solution strategy for non-linear analysis, where macroscopic loading is applied as a homogeneous strain to the RVE. The RVE stress field is calculated through the SCA, followed by homogenisation to obtain the effective stress resultants ($\bar{\bm{N}},\bar{\bm{M}},\bar{\bm{Q}}$). These homogenised quantities are derived through volume averaging of the stress and moment resultant distributions of the RVE:

[TABLE]

In such cases, the constitutive laws may depend on both strain and strain rates in specific non-linear models. The computational framework for non-linear analysis is summarised in Algorithm 3.

3 Results

Linear analyses of periodic perforated plates (both symmetric and asymmetric) and a twill woven composite are performed, comparing the proposed method with the DNS method (FFT-based approach). Subsequently, non-linear analyses are carried out on an elasto-plastic periodic perforated plate and a practical plain woven composite. For the latter, the predictions obtained using the proposed method are compared with experimental results. The computation time comparison for all examples is discussed in Sec. 4.2.

3.1 Linear analysis of periodic perforated plate

The periodic perforated plates considered in this example feature “+”-shaped pores, with the symmetric plate composed of steel and the asymmetric plate consisting of an aluminium–steel laminate. The dimensions and geometric details of the RVE models are shown in Fig. 3. The Young’s modulus and Poisson’s ratio of the steel plate are 200 GPa and 0.3, respectively [47], while those of the aluminium plate are 70 GPa and 0.33, respectively [48]. Both models have an in-plane resolution of $320\times 320$, and a convergence study of the selected resolution is presented in [36].

The linear analysis is performed to determine the effective properties of the RVE in the elastic state and requires six orthogonal loading conditions. Average in-plane strain is applied to the RVE model to calculate the effective membrane properties of the periodic plate. The three cases considered are as follows:

• •

Case. 1: $\bar{\bm{\varepsilon}}=\left[\begin{matrix}0.01&0\\ 0&0\end{matrix}\right],\bar{\bm{\phi}}=\bm{0}$

• •

Case. 2: $\bar{\bm{\varepsilon}}=\left[\begin{matrix}0&0\\ 0&0.01\end{matrix}\right],\bar{\bm{\phi}}=\bm{0}$

• •

Case. 3: $\bar{\bm{\varepsilon}}=\left[\begin{matrix}0&0.01\\ 0.01&0\end{matrix}\right],\bar{\bm{\phi}}=\bm{0}$

The effective bending properties are determined from the following three loading cases:

• •

Case. 4: $\bar{\bm{\varepsilon}}=\bm{0},\bar{\bm{\phi}}=\left[\begin{matrix}0.01&0\\ 0&0\end{matrix}\right]$

• •

Case. 5: $\bar{\bm{\varepsilon}}=\bm{0},\bar{\bm{\phi}}=\left[\begin{matrix}0&0\\ 0&0.01\end{matrix}\right]$

• •

Case. 6: $\bar{\bm{\varepsilon}}=\bm{0},\bar{\bm{\phi}}=\left[\begin{matrix}0&0.01\\ 0.01&0\end{matrix}\right]$

The six load cases are analysed using both the FFT method and the SCA approach. In the SCA calculations, five different cluster numbers $H$ (i.e., 32, 64, 128, 256, and 512) are considered. In addition, both the CPT and FoPT problems are analysed. The homogenised properties of symmetric and asymmetric perforated plates are presented in Tables 1 and 2, respectively. Homogenised properties with zero values (e.g., $B_{ij}=0$ for the symmetric perforated plate) and symmetric properties (e.g., $A_{11}=A_{22}$) are omitted. Since the thickness-to-length ratio of all plates does not exceed 0.1, the homogenised properties obtained from CPT and FoPT exhibit minimal differences. However, as the thickness-to-length ratio increases, the discrepancy in $D_{ij}$ also increases. Furthermore, a comparison between the SCA and FFT results indicates that as $H$ increases, the homogenised properties progressively converge to those obtained by FFT. When $H=128$, the error remains within 1%. This demonstrates that the SCA can achieve high accuracy in predicting homogenised properties using a limited number of clusters.

Figs. 4–7 illustrate the stress distributions for four loading cases based on CPT and FoPT, along with the corresponding stress values extracted along the central vertical line at $x=2\,\mathrm{mm}$. For smaller $H$, the interpolation of the stress field, constrained by the coarse cluster resolution, lacks the accuracy needed to precisely capture stress concentration regions. In this case, the calculated stress range remains within the extreme values obtained from FFT, with SCA overestimating stresses in low-stress regions and underestimating them in high-stress regions. This discrepancy gradually diminishes as $H$ increases, and at $H=512$, the stress distributions exhibit close agreement. However, the extracted stress values indicate that SCA produces piecewise-constant stress distributions, reflecting its inherent cluster-wise uniform stress assignment. Furthermore, while the distributions of $N$ remain largely consistent between CPT and FoPT, notable differences are observed in $M$. Compared to CPT, FoPT demonstrates superior accuracy in capturing $M$ in SCA calculations, as it incorporates more detailed bending features when constructing the offline clustering database.

3.2 Linear analysis of twill woven composites

A twill weave composite material is investigated in this numerical example, as shown in Fig. 8(a), where the macroscopic configuration and the microscopic model are presented. The computational domain is discretised as $320\times 320\times 40$. The elastic properties of the yarn and matrix are listed in Table 3. Notably, due to the strongly anisotropic characteristics of woven composites, clustering analysis based solely on the strain concentration tensor may result in excessive variations in mechanical properties within each cluster, potentially compromising computational accuracy. To address this issue, a two-step clustering strategy that incorporates both the material concentration tensor and the strain concentration tensor is commonly employed. Specifically, the initial clustering is performed based on the material concentration tensor $\bm{C}$, followed by secondary clustering of the strain concentration tensor $\bm{O}$ within each sub-cluster. However, existing works have not explicitly clarified the relationship between the clustering configurations of the material concentration tensor and the strain concentration tensor. Based on this context, three progressive clustering arrangements (A-1, A-2, and A-3) are designed for this example, as shown in Fig. 8(b). A-1 directly uses the strain concentration tensor for one-step clustering; A-2 fixes the clustering number of the material concentration tensor at 8 while progressively increasing the sub-cluster number of the strain concentration tensor; conversely, A-3 fixes the sub-cluster number of the strain concentration tensor at 8 while gradually increasing the clustering number of the material concentration tensor. By comparing the computational accuracy of the different schemes (detailed data are provided in Appendix A.4), A-2 is finally adopted as the clustering approach for this example.

The homogenised properties of the twill weave composite material are presented in Table 4, based on FFT and SCA with five different cluster numbers $H$. Since the overall structure of the twill weave composite is symmetric with respect to the mid-plane at the macroscopic level, the components $B_{ij}=0$. The results indicate a significant discrepancy between the homogenised properties calculated using CPT and FoPT, mainly due to the relatively large thickness-to-length ratio of the model ($0.125>0.1$). Furthermore, similar to the periodically perforated plate model, as $H$ increases, the homogenised properties gradually converge to those obtained via FFT. For $H=128$, the error remains within 1%.

Fig. 9 and Fig. 10 present the stress distributions for four loading cases, along with the corresponding stress values extracted along the diagonal line at $x=y$, based on CPT and FoPT, respectively. The results indicate that even for $H=32$, the SCA method is able to capture the overall stress distribution, as further supported by the extracted stress values along the diagonal. In addition, the agreement between the SCA and FFT results improves as $H$ increases. However, even at $H=512$, significant local differences remain between the stress field computed by SCA and that obtained using the FFT method, as illustrated by the stress distributions in Case 1. This difference arises from the fact that clustering in SCA is performed based on the strain concentration tensors under six orthogonal load cases, aiming to generalise to arbitrary loading conditions. Consequently, the material points within each cluster do not exhibit strict similarity under a given load case. A rigorous theoretical discussion on the convergence of cluster number in SCA is available in [29]. Furthermore, the extracted stress values show that as $H$ increases, the SCA results gradually converge to those of the FFT method. Regardless of whether CPT or FoPT is used, the stress values extracted from SCA at $H=128$ are largely consistent with those of FFT, and at $H=512$ they are nearly identical.

3.3 Non-linear elasto-plastic analysis of periodic perforated plate

In the linear elastic analysis of fixed models, the performance advantages of the SCA method are not fully demonstrated due to its reliance on an "offline–online" two-stage computational framework. Its primary advantage lies in using a single offline calculation to enhance computational efficiency during the online stage. This benefit is particularly important in scenarios requiring repeated computations, such as multi-incremental and multiscale analyses. However, as discussed in Sec. 2.3.4, when SCA is used to simulate non-linear behaviour, its prediction accuracy may be affected. Therefore, it is necessary to evaluate the predictive capability of SCA in non-linear simulations. To this end, we present the elasto-plastic (weakly non-linear) analysis of the symmetric perforated plate described in Sec. 3.1, under different loading conditions. Since the thickness-to-length ratio of the perforated plate is $0.05<0.1$, the following analysis is based solely on the CPT. In this analysis, the elastic stiffness properties remain the same as in Sec. 3.1, while the plastic behaviour follows the von Mises yield criterion with the following linear hardening law [47]:

[TABLE]

where $\sigma_{y}$ and $p$ denote the yield stress and equivalent plastic strain, respectively.

Three loading conditions are applied in this example: uniaxial tension, pure shear ($\gamma_{12}=0.01$), and a complex loading path (as shown in Fig. 11). The first two loading conditions are applied incrementally with 100 uniformly distributed load steps, while the complex loading path is discretised into 400 increments to ensure an accurate representation of the loading process. The boundary conditions for the uniaxial tension case are not explicitly defined by a homogeneous strain loading. Therefore, an a posteriori correction procedure is employed to ensure that the following boundary conditions are satisfied:

[TABLE]

Fig. 12(a) and Fig. 12(b) show the stress–strain curves obtained from SCA and FFT calculations under uniaxial tension and pure shear loading, respectively. The results indicate that while the SCA method effectively captures the non-linear elasto-plastic behaviour, its accuracy decreases when $H$ is relatively small. To further illustrate this trend, we include the case with $H=1024$ and show that increasing the number of clusters consistently improves computational accuracy. Furthermore, Fig. 12(c) shows the normal and tangential stress–strain responses under a complex loading path, where again a high level of agreement between the two methods is observed.

Fig. 13 shows the final equivalent plastic strain distributions obtained using the SCA and FFT methods under three different loading conditions, along with the corresponding equivalent plastic strain values extracted along the central vertical line. Similar to the linear elastic analysis, SCA is also capable of predicting the equivalent plastic strain field. However, the local distribution appears more diffuse, resulting in a lower peak equivalent plastic strain compared to the FFT results. This is due to the fundamental design principle of SCA, which aims to achieve accurate prediction of globally homogenised mechanical responses while minimising computational information. Within the coarse clustering framework of SCA, the equivalent plastic strain in each cluster is treated as the average of the equivalent plastic strains of all pixels within that cluster. Consequently, high equivalent plastic strain pixels near the pore are averaged with low equivalent plastic strain pixels farther from the pore, thereby mitigating the local strain concentration effect. As a result, the SCA method provides a more conservative estimate, leading to a predicted stress–strain curve that is higher than that obtained using the FFT method. To increase the resolution in local regions, global accuracy can be improved by increasing the total number of clusters. Alternatively, a more efficient approach is to apply a non-linear K-means clustering strategy, where a larger number of clusters are dynamically assigned to critical regions with high strain gradients, such as areas around the pore.

3.4 Non-linear progressive failure analysis of a practical plain woven composite

To demonstrate the potential of the proposed method for solving non-linear problems and handling complex structures, a progressive failure analysis of a multilayer plain woven composite is presented in this section. In addition, experimental validation is performed to further verify the effectiveness of the proposed method. The numerical model of the plain woven composite is constructed based on real material data obtained from microscopic images. A modelling strategy known as the projected element method (PEM) is used to reconstruct the RVE of the composite, accurately capturing the deformation of the yarns introduced during the manufacturing process. The reconstructed composite structure and corresponding RVE, including the plate formulation, are shown in Fig. 14.

The non-linear behaviour of the composite is simulated using the damage models described in [36]. The resin phase follows a rate-dependent damage model, while the yarn phase follows the Hashin damage model as the failure initiation criterion. The post-damage behaviour of the composite is described by a Weibull distribution-based model; the specific material parameters can be found in Appendix A.5. The in-plane discretisation resolution of the model is set to $128\times 128$. Since the thickness-to-length ratio of the model is $0.222>0.1$, only the FoPT is used for the simulations. Uniaxial tension and bending loads are applied, both with a loading rate of $10^{-4}\ \mathrm{s}^{-1}$.

Fig. 15 shows the stress–strain curves obtained using the SCA and FFT methods under the two different loading conditions, along with the experimental results. The results show consistency across all methods, with the calculated modulus closely matching the experimental data. Under tensile loading, the calculated tensile strength is approximately 8% higher than the experimental average, while under bending loading, the calculated strength is approximately 5% lower than the experimental average.

4 Discussion

4.1 Sensitivity analysis of $G^{0}$

The accuracy of SCA strongly depends on the self-consistent update of the reference medium. In this study, the projection-based scheme described in Appendix A.3 is used to update the reference medium. However, for the FoPT problem, as formulated in Eq. (24), the material parameters ($\mu_{b}^{0}$ and $G^{0}$) of the reference medium in the Green’s operator are coupled to the discrete frequencies. This coupling requires the Green’s operator—and therefore the interaction tensor—to be recalculated at each update step. As discussed in Sec. 2.3.2, if $\mu_{b}^{0}$ and $G^{0}$ satisfy the assumption in Eq. (26), the coupled term ($\rho$) can be treated as a coefficient, eliminating the need to recalculate the interaction tensor.

Existing sensitivity analyses in the context of the FFT method indicate that $\mu_{b}^{0}$ exhibits significant sensitivity, and the term $2\mu_{b}^{0}+\lambda_{b}^{0}$ influences the overall convergence behaviour of the computation. Therefore, when solving the FoPT problem with SCA, we employ a self-consistent update strategy for $\mu_{b}^{0}$ to ensure its accuracy, while $G^{0}$ is updated accordingly based on the condition imposed by Eq. (26). Although $G^{0}$ has been shown to be relatively insensitive in the FFT method, its sensitivity in SCA remains unclear. The following is the selection strategy for $G^{0}$ in the FFT method:

[TABLE]

Here, $\text{eig}(\cdot)$ denotes the eigenvalues of a matrix, and $\beta(G^{0})$ represents a coefficient associated with the selection of the reference medium. By adjusting this coefficient, $G^{0}$ can be modified. To comprehensively investigate the sensitivity of $G^{0}$, the twill weave composite material introduced in Section 3.2 is employed. Several values of $\beta(G^{0})$ are selected, ranging from $10^{-2}$ to $10^{2}$, while all other parameters of the reference medium are kept fixed.

Fig. 16 illustrates the homogenised properties obtained from the SCA calculations for two different cluster numbers ($H=128$ and $H=512$), along with the total number of iterations required under six orthogonal loading conditions. The dashed lines in the figure represent the reference homogenised properties computed using the FFT method. All stiffness values are in close agreement with the reference values, indicating that variations in $G^{0}$ do not affect the computational accuracy. In addition, the total number of Newton–Raphson iterations required to solve the six load cases remains nearly constant for different $G^{0}$. This conclusion holds for different clustering numbers, demonstrating that the computation of membrane and bending responses is insensitive to $G^{0}$. Consequently, the assumption used in Eq. (26) to calculate the FoPT is valid.

4.2 Computation time

The computational efficiency of the developed SCA, as a reduced-order model (ROM) for heterogeneous plates, is systematically evaluated through comparative studies with direct numerical simulation (DNS) using the FFT-based homogenisation method. All numerical tests were conducted on a workstation equipped with an Intel(R) Core(TM) i7-11700 CPU.

Table 5 summarises the cumulative computational time for six loading cases across three examples: Example 1 (symmetrical perforated plate in Sec. 3.1), Example 2 (asymmetrical perforated plate in Sec. 3.1), and Example 3 (twill-woven composite in Sec. 3.2). SCA demonstrates superior computational efficiency, achieving over 10× speed-up compared to FFT even at $H=512$. The speed-up exceeds two orders of magnitude when $H=128$ is used, with acceptable solution accuracy (error $<1\%$). Similar performance advantages are observed for elasto-plastic analyses of perforated plates under three loading conditions (Table 6), where SCA achieves more than twice the computational efficiency relative to FFT at $H=1024$. Notably, the computational time exhibits quadratic scaling with respect to the cluster number $H$ ($T\propto H^{2}$), consistent with numerical observations in classical SCA frameworks [30]. These results highlight SCA’s exceptional ability to balance computational accuracy and efficiency—especially considering that the FFT-based DNS method itself is widely recognised as a computationally efficient homogenisation technique within the computational mechanics community. The speed-up achieved stems from the intrinsic dimensionality reduction of SCA’s data-driven architecture, which bypasses intensive local field computations through effective clustering-based homogenisation.

The computation time required for DNS generally increases with model resolution. In contrast, the SCA features a computational cost intrinsically governed by the number of clusters, as its degrees of freedom correspond directly to the discretised cluster quantity. To systematically validate this hypothesis, five different discretisations are implemented for the symmetrical perforated plate model in Example 1. Both SCA and FFT methods are employed to compute six loading cases under CPT. The computation time and error are summarised in Table 7, where the error represents the mean absolute error of the homogenised properties relative to the highest-resolution FFT results. The computational results show that FFT computation time increases progressively with higher resolution, while the SCA computation time remains nearly constant across all discretisations. This confirms that the computational cost of SCA is independent of model resolution, enabling refined discretisation to capture intricate geometric features without incurring additional computational cost. Furthermore, the FFT error decreases systematically with increasing resolution, while the SCA error remains essentially unchanged across different discretisations. The SCA error exhibits a decreasing trend with increasing cluster quantity ($H$), which conclusively verifies that its computational degrees of freedom depend solely on the number of clusters rather than the model resolution.

5 Conclusions

This study successfully develops a self-consistent clustering analysis (SCA) framework that integrates classical plate theory (CPT) and first-order shear deformation plate theory (FoPT), providing an efficient unit cell solution for heterogeneous plate structures. By extending the theoretical framework, the proposed reduced-order model significantly enhances computational efficiency and applicability for heterogeneous plates. Numerical validations demonstrate that the method accurately captures the non-linear behaviour of both isotropic and anisotropic materials, including elastoplastic and progressive failure mechanisms. Compared to direct numerical simulation (DNS), the proposed approach achieves a computational error below 1% with a limited number of clusters while delivering over a 100-fold speed-up.

Although the present work focuses on mechanical boundary value problems, the modular design of this framework allows seamless extension to complex multiphysics scenarios such as thermo-mechanical coupling and micropolar continuum mechanics. The findings of this study provide an innovative methodological foundation for various cutting-edge research domains, particularly in: (1) multiscale mechanical property prediction of composite materials, (2) topology optimisation of porous materials [50], (3) intelligent design of mechanical metamaterials, and (4) data-driven modelling. These directions will serve as key focal points for future research.

Acknowledgements

This work was financially supported by the National Key Research and Development Program of China (2022YFB3707800).

Appendix A Appendix

A.1 Green’s operator of CPT and FoPT

The Green’s operator of CPT consists of two components, and its explicit expression in Fourier space is as follows:

[TABLE]

The Green’s operator of FoPT consists of five components, and its explicit expression in Fourier space is as follows:

[TABLE]

Here, $\lambda_{m}^{0}$, $\mu_{m}^{0}$, $\lambda_{b}^{0}$, $\mu_{b}^{0}$, and ${G}^{0}$ represent the Lame constants and shear modulus of the reference medium. It should be noted that the $\hat{G}_{khij}^{0(NN)}$ terms are identical for both CPT and FoPT, whereas the $\hat{G}_{khij}^{0(MM)}$ terms are different.

A.2 Cluster-based reduced-order modeling for CPT

First, the Lippmann-Schwinger Eq. ( 6 ) is written in an incremental form:

[TABLE]

The RVE model is divided into H material clusters. The volume averaging of Eq. ( 40 ) is performed within the I-th cluster as follows:

[TABLE]

Substituting the discretized stress and strain fields under the uniform assumption into the above equation:

[TABLE]

where the interaction tensor $\bm{D}^{IJ\left(\sim,\sim\right)}$ is defined as:

[TABLE]

It can be calculated in the Fourier space as follows:

[TABLE]

To facilitate the self-consistent update of the reference medium parameters during the online stage, we separate the discrete frequency term and the reference medium parameter term in Green’s operator as follows:

[TABLE]

where,

[TABLE]

The system of equations is solved using the Newton-Raphson method, where the residual of the system of equations is expressed as:

[TABLE]

The Jacobian matrices are as follows:

[TABLE]

A.3 The projection-based self-consistent scheme

The effective tangential stiffness of the RVE model for the partitioned H material clusters is given by:

[TABLE]

where, $\mathbb{A}_{\text{alg}}^{I}$ and $\mathbb{D}_{\text{alg}}^{I}$ represent the local constitutive laws for the incremental strain of the i-th material cluster:

[TABLE]

Then, the effective tangential stiffness is projected onto a fourth-order isotropic tensor as follows:

[TABLE]

where,

[TABLE]

A.4 Clustering arrangement comparison of twill woven composites

For the twill woven composites in Sec. 3.2, computational analyses are performed on three clustering arrangements shown in Fig.8 . Table 8) shows the errors in the homogenized properties calculated by SCA relative to the FFT results. All three arrangements show decreasing errors with increasing $H$, while maintaining relatively small error magnitudes overall. However, in terms of both mean error and errors in the homogenized properties along the principal diagonal, the A-2 has the smallest error for both the CPT and FoPT.

A.5 Material properties of the practical plain woven composite

All the material parameters of the practical plain woven composite are summarized in Table 9.