Shuitao Gu

Yu Cong

Wei Shen

Gongye Zhang
email: gyzhang@seu.edu.cn

Yu 

Cong

A Transversely Isotropic Magneto-electro-elastic

Keywords: Circular Kirchhoff plate, Magneto-electro-elastic plate, Microstructure effect, Couple stress theory, Static bending

HAL is

Introduction

Due to the piezomagnetic/piezoelectric/magnetoelectric effect, thin plates made of magneto-electro-elastic (MEE) materials are able to realize the mutual conversion between magnetic, electric and mechanical energies. Such characteristics have found significant applications in sensors, resonators, energy harvesters, semiconductors and other systems [START_REF] Ghayesh | Size-dependent performance of microgyroscopes[END_REF][START_REF] Salas | A topology optimization formulation for transient design of multi-entry laminated piezocomposite energy harvesting devices coupled with electrical circuit[END_REF][START_REF] Devoe | Piezoelectric thin film micromechanical beam resonators[END_REF][START_REF] Liang | Transient bending vibration of a piezoelectric semiconductor nanofiber under a suddenly applied shear force[END_REF][START_REF] Hu | Coupled analysis for the harvesting structure and the modulating circuit in a piezoelectric bimorph energy harvester[END_REF][START_REF] Yang | Electronic band energy of a bent ZnO piezoelectric semiconductor nanowire[END_REF]. It is observed that these thin plates exhibit microstructure-dependent size effects at the micron and nanometer scales [START_REF] Li | Static bending and free vibration of a functionally graded piezoelectric microplate based on the modified couple-stress theory[END_REF][START_REF] Lim | Size-dependent nonlinear response of thin elastic films with nano-scale thickness[END_REF][START_REF] Yin | Vibration analysis of microscale plates based on modified couple stress theory[END_REF].

Owing to the lack of any size-dependent material parameter [START_REF] Lam | Experiments and theory in strain gradient elasticity[END_REF], classical elasticity theories cannot effectively describe such size effects. Hence, higher-order elasticity theories that contain additional parameters should be introduced to develop models for thin MEE plates considering microstructure-dependent size effects.

Several studies for MEE plates using higher-order elasticity theories have been completed. For instance, a non-local constitutive relation proposed by Eringen [START_REF] Eringen | On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves[END_REF][START_REF] Eringen | Nonlocal continuum field theories[END_REF] has been used to develop a model for the MEE Mindlin plate (i.e., Li et al. [START_REF] Li | Buckling and free vibration of magnetoelectroelastic nanoplate based on nonlocal theory[END_REF]).

Vinyas et al. [START_REF] Vinyas | Numerical analysis of the vibration response of skew magneto-electro-elastic plates based on the higher-order shear deformation theory[END_REF] addressed a free vibration problem of skew MEE plates by employing Reddy's third-order shear deformation theory (e.g. [START_REF] Reddy | Mechanics of laminated composite plates and shells: theory and analysis[END_REF]). Zheng et al. [START_REF] Zheng | Nonlinear bending analysis of magnetoelectroelastic rectangular plates using higher order shear deformation theory[END_REF] also applied Reddy's high-order shear deformation theory and combined von Karman's nonlinear geometric equation to analyze the nonlinear bending of MEE plates. Wang et al. [START_REF] Wang | Two-dimensional linear elasticity theory of magneto-electro-elastic plates considering surface and nonlocal effects for nanoscale device applications[END_REF] established a model for MEE plates considering both nonlocal and surface effects. Ebrahimi and Dabbagh [START_REF] Ebrahimi | On flexural wave propagation responses of smart FG magneto-electro-elastic nanoplates via nonlocal strain gradient theory[END_REF] developed a model for MEE nanoplates via a nonlocal strain gradient theory to capture size effects. Qu et al. [START_REF] Qu | A microstructure-dependent anisotropic magneto-electro-elastic Mindlin plate model based on an extended modified couple stress theory[END_REF] developed a new anisotropic MEE Mindlin plate model by using an extended modified couple stress theory.

Different from the classical couple stress theory [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Mindlin | Influence of couple-stresses on stress concentrations[END_REF][START_REF] Tang | Interpretation of bend strength increase of graphite by the couple stress theory[END_REF], the modified couple stress theory [START_REF] Yang | Couple stress based strain gradient theory for elasticity[END_REF][START_REF] Park | Variational formulation of a modified couple stress theory and its application to a simple shear problem[END_REF] and its extended version [START_REF] Qu | A microstructure-dependent anisotropic magneto-electro-elastic Mindlin plate model based on an extended modified couple stress theory[END_REF][START_REF] Zhang | A non-classical model for an orthotropic Kirchhoff plate embedded in a viscoelastic medium[END_REF][START_REF] Zhang | A transversely isotropic magneto-electro-elastic Timoshenko beam model incorporating microstructure and foundation effects[END_REF] only consider the symmetric part of the curvature tensor conjugated with the deviatoric part of the couple stress tensor, which contain fewer material length scale parameters than the classical counterparts.

Considering the difficulties in determining additional material parameters and the physical interpretation of the microstructure, these modified theories have been successfully employed to develop non-classical beam and plate models [START_REF] Qu | A microstructure-dependent anisotropic magneto-electro-elastic Mindlin plate model based on an extended modified couple stress theory[END_REF][START_REF] Zhang | A non-classical model for an orthotropic Kirchhoff plate embedded in a viscoelastic medium[END_REF][START_REF] Zhang | A transversely isotropic magneto-electro-elastic Timoshenko beam model incorporating microstructure and foundation effects[END_REF][START_REF] Zhou | A nonclassical model for circular mindlin plates based on a modified couple stress theory[END_REF][START_REF] Zhang | A non-classical model for circular Kirchhoff plates incorporating microstructure and surface energy effect[END_REF][START_REF] Ma | A non-classical Mindlin plate model based on a modified couple stress theory[END_REF].

Because of the complexity of cylindrical coordinates in mathematics, little research has been done on circular plates. However, due to the special geometry of circular and annular plates, circular MEE plates are often used in a wide range of engineering applications. Hence, a model for transversely isotropic MEE circular Kirchhoff plates at very small scale and can capture microstructure effect is very desirable.

The current work aims to develop a microstructure-dependent transversely isotropic MEE model for circular Kirchhoff plates. The rest of the paper is organized as follows. In Section 2, a new model is established through a variational method by using the extended modified couple stress theory. In Section 3, several special cases are provided. In Section 4, the static bending problem of a clamped circular plate subject to a uniformly distributed constant load is solved by applying the formulas in Section 2 through Fourier-Bessel series. The paper concludes in Section 5 with a summary.

Formulation

Figure 1 shows an annular circular thin plate of uniform thickness h, inner radius a and outer radius b. The displacement field in such a circular Kirchhoff plate can be written in the form of the cylindrical coordinate system as [START_REF] Zhang | A non-classical model for circular Kirchhoff plates incorporating microstructure and surface energy effect[END_REF] 

     , , r w r u r θ z u r z r     (1a) ( , , ) 0 θ u r θ z  (1b)   ( , , ) z u r θ z w r  (1c)
where ur, uθ and uz are, respectively, the radial, tangential and transverse components of the displacement vector u. u and w are, respectively, the radial and transverse components of the displacement vector of point on the mid-plane of the plate. According to the work of Zhang et al. [START_REF] Zhang | A transversely isotropic magneto-electro-elastic Timoshenko beam model incorporating microstructure and foundation effects[END_REF] and Wang [START_REF] Wang | On buckling of column structures with a pair of piezoelectric layers[END_REF], the electric and magnetic potentials in such a Kirchhoff plate can be written as

    0 2 , , cos π z r θ z z γ r γ h h           (2a)     0 2 , , cos π z M r θ z z ζ r ζ h h          (2b)
where Φ represents the electric potential, and M represents the magnetic potential. γ and ζ represent the spatial variations of electric potential and magnetic potential on the mid-plane, respectively. γ0 and ζ0 denote the external electric potential and magnetic potential, respectively.

The constitutive equations considering the extended version of modified couple stress theory for transversely isotropic MEE materials can be expressed as [START_REF] Qu | A microstructure-dependent anisotropic magneto-electro-elastic Mindlin plate model based on an extended modified couple stress theory[END_REF][START_REF] Zhang | A transversely isotropic magneto-electro-elastic Timoshenko beam model incorporating microstructure and foundation effects[END_REF] e q e q e q H e q

C C C C C C C C C σ ε C σ ε C σ ε C C σ ε σ ε ε                                                                   31
                                                          (3a)
A A A A A A m χ A A A m χ A m χ A m χ A A m χ χ                                                                     (3b)
                                       
z z E d H E d H E d H                                           (3c)
q ε B q ε B q ε ε ε q q                                          1 1 1 1 1 1 1 1 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 r r θ θ z z μ H d E μ H d E μ H d E                                           (3d)
where σij are the Cauchy stresses, mij are the couple stresses, Di are the electric displacements, Bi are the magnetic fluxes. Cij and Aij denote the components of the elastic stiffness tensor and the couple stress stiffness tensor, respectively. eij and qij denote the piezoelectric and piezomagnetic tensors, respectively. ϵij, dij and μij denote the components of the dielectric tensors, magneto-dielectric tensors and magnetic permeability tensors, respectively. ε ij and χij are, respectively, the components of strain and curvature tensors defined by

  T 1 2         ε u u (4a)   T 1 2         χ ψ ψ (4b) with ψ being 1 c u r l 2  ψ u (4c)
in which, u is the displacement vector, ▽ denotes the gradient, and the superscript T represents the transpose. In addition, the components of the electric field intensity Ei and the magnetic field intensity Hi are given by

  E Φ,   H M (5a
rr θθ zz rθ zr zθ u w w ε z ε u z ε ε ε ε r r r r                     (6a) 1 1 , 0 2 rθ θr rr θθ zz rz θz w χ χ r χ χ χ χ χ r r r                   (6b) 0 2 cos , 0, sin r θ z γ π γ π π E z E E z γ h r h h h                    (6c) 0 2 cos , 0, sin r θ z ζ π ζ π π H z H H z ζ h r h h h                    (6d)
The first variation of the Gibbs-type energy for the deformed MEE Kirchhoff plate is given by [START_REF] Qu | A non-classical theory of elastic dielectrics incorporating couple stress and quadrupole effects: part I -reconsideration of curvature-based flexoelectricity theory[END_REF]  

d ij ij ij ij i i i i δU σ δε m δχ D δE B δH V       ( 7
)
where Ω is the volume of the plate, and

      2 2 2 2 2 0 2 0 d 2 d 3 d d 2 ij ij ij ij
σ δε m δχ V σ δε σ δε m δχ V M Y Y rN rM N δu r δw r θ r r r r r Y rM δw rN δu M Y r δw rM rY r r r                                                                      d π θ  (8)
N σ z N σ z M σ z z M σ z z Y m z                (9a-e)     2 2 0 0 d + d d d d i i r r θ θ z z π b π b r z r a a D δE V D δE D δE D δE V r δγ r δγ r θ r δγ θ r                         (10) / 2 / 2 / 2 / 2 cos d , sin d h h r r z z h h π π π D z z D z z h h h                     (11a, b)     2 2 0 0 d d d d d i i r r θ θ z z π b π b r z r a a B δH V B δH B δH B H V r δζ r δγ r θ r δζ θ r                          (12) / 2 / 2 / 2 / 2 cos d , sin d h h r r z z h h π π π B z z B z z h h h                     (13a, b)
Note that Eqs. [START_REF] Lim | Size-dependent nonlinear response of thin elastic films with nano-scale thickness[END_REF][START_REF] Yin | Vibration analysis of microscale plates based on modified couple stress theory[END_REF] stand for the elastic part of the first variation of the Gibbs-type energy, while Eqs. [START_REF] Lam | Experiments and theory in strain gradient elasticity[END_REF][START_REF] Eringen | On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves[END_REF] and Eqs. [START_REF] Eringen | Nonlocal continuum field theories[END_REF][START_REF] Li | Buckling and free vibration of magnetoelectroelastic nanoplate based on nonlocal theory[END_REF] represent the electric potential-related and magnetic potential-related parts, respectively.

The first variation of the virtual work (with c θ = 0) can be obtained as [START_REF] Zhou | A nonclassical model for circular mindlin plates based on a modified couple stress theory[END_REF][START_REF] Zhang | A non-classical model for circular Kirchhoff plates incorporating microstructure and surface energy effect[END_REF]  

2 2 0 0 d d d b π b π r z r θ z a a δw δW rf δu rf δwδw r θ rt δu rM rt δw θ r                 ( 14
)
where fi is the body force resultant (force per unit area) of the i-direction, ti is the Cauchy traction resultant (force per unit length) of the i-direction, and θ M is the applied moment per unit arc length along the circular boundary.

Applying the Gibbs-type variational principle and the fundamental lemma of the calculus of variations [START_REF] Qu | A non-classical theory of elastic dielectrics incorporating couple stress and quadrupole effects: part I -reconsideration of curvature-based flexoelectricity theory[END_REF], the governing equations result in

  11 d 1 d 0 d d r hC ru f r r r         (15a)   3 4 2 2 11 11 12 31 31 2 2 0 12 2 z h h h h C A A w e γ q ζ f π π                 (15b)
h h h π π μ ζ d γ q w μ ζ d γ π h h         (15d)
and the boundary conditions yield

                      (
         (16d)   3
h w h w h h h C C e γ q ζ A A w r r r π π w w M r a r b r r                     (16e)
where

2 4 1 1 1 , r r r r r r r r r r r r                                         (17)

Specific Models

In this section, a few specific models reduced from the current non-classical transversely isotropic MEE circular Kirchhoff plate model developed in Section 2 are provided.

Transversely Isotropic Piezoelectric Circular Kirchhoff Plate Model

Incorporating the Microstructure Effect

When the piezomagnetic effect is neglected by setting q ij = 0, µ ij = 0 and

d ij = 0, Eqs. (15a-d) reduce to   11 1 d 0 d r d hC ru f dr r r         (18a)   3 4 2 11 11 12 31 2 0 12 2 z h h h C A A w e γ f π               (18b)
                   (
                   (19d)

Transversely Isotropic Piezomagnetic Circular Kirchhoff Plate Model

Incorporating the Microstructure Effect

When the piezoelectric effect is neglected by setting e ij = 0, ϵ ij = 0 and d ij = 0, Eqs.

(15a-d) reduce to   11 d 1 d 0 d d r hC ru f r r r         (20a)   3 4 2 11 11 12 31 2 0 12 2 z h h h C A A w q ζ f π               (
h h h ζ C A A w q t r π r w w r a r b                    (
h w h w h h C C q ζ A A w M r r r π w w r a r b r r                    (
z h h h C w e γ q ζ f π π         (22b)
h h h π π μ ζ d γ q w μ ζ d γ π h h         (
          (23a)
                 (23e)

Example: Static Bending of a Clamped Circular Plate

To study the newly developed circular Kirchhoff plate model, a static bending problem of a solid circular plate (with a = 0 and b = 20h) is solved herein. The plate is clamped at its outer edge and subjected to a uniformly distributed constant load q pointing downward.

For static bending problems, the governing equations are consistent with Eqs.

(15a-d).

In view of Eqs. (16a-e), the boundary conditions for the current circular Kirchhoff plate at its outer edge r = b are given by

( ) 0 u b  , ( ) 0 w b  , d 0 d r b w r   , ( ) 0 γ b  , ( ) 0 ζ b  (24a-e)
Besides, the following symmetry conditions need to be considered:

(0) 0 u  , 0 d d 0 r w r   , 0 d d 0 r γ r   , 0 d d 0 r ζ r   (25a-d)
Obviously, the radial displacement u is uncoupled with w, γ and ζ, so that it can be calculated by solving the boundary value problem (BVP) defined by Eqs. (15a), (24a) and (25a). According to this BVP (with fr = 0), it is not hard to get that u (r) = 0 at any point in the plate.

The transverse displacement w, electric potential γ and magnetic potential ζ can be calculated by solving the BVP defined by Eqs. (15b-d), (24b-e) and (25b-d). The coupled high-order ordinary differential equations (ODEs) are not easy to solve by direct integration because the solutions containing modified Bessel functions are usually very complicated [START_REF] Zhou | A nonclassical model for circular mindlin plates based on a modified couple stress theory[END_REF][START_REF] Ariman | On circular micropolar plates[END_REF][START_REF] Wang | An analytical solution for a multilayered magneto-electro-elastic circular plate under simply supported lateral boundary conditions[END_REF].

In order to deal with such difficulties, the Fourier-Bessel series [START_REF] Zhou | A nonclassical model for circular mindlin plates based on a modified couple stress theory[END_REF] are adopted.

Solutions for w, γ and ζ can be written as:

    1 1 1 d d n n n w r W J λ r r     (26a)   0 0 1 ( ) n n n γ r H J λ r     (26b)   0 0 1 ( ) n n n ζ r Z J λ r     (26c)
where Ji (x) is the Bessel function of the first kind of order i, λin is the nth real zero of the function Ji (x) = 0, and / When considering the boundary conditions in Eq. (24b), the expression of w can be derived from integrating Eq. (26a)

      0 1 0 1 1 1 1 n n n n n w r W J λ b J λ r λ          (27)
The uniformly distributed constant load q can be expanded as [START_REF] Zhang | A non-classical model for circular mindlin plates incorporating microstructure and surface energy effects[END_REF][START_REF] Kreyszig | Advanced engineering mathematics[END_REF]    

0 1 1 2 1 ( ) n n n q q r J λ r J λ     or     0 0 1 1 0 0 2 ( ) n n n n q q r J λ r J λ λ     (28a
n n n n K K K W Q K K K H K K K Z                                (29)
where

      3 3 11 11 11 12 1 0 1 12 2 n n h h K C A A λ J λ r           (30a)     2 12 31 0 0 0 2 n n h K e λ J λ r π  (30b)     2 13 31 0 0 0 2 n n h K q λ J λ r π  (30c)   21 31 1 0 1 2 n n h K e λ J λ r π   (30d)       2 2 22 11 0 0 0 33 0 0 2 2 n n n h π K λ J λ r J λ r h    (30e)       2 2 23 11 0 0 0 33 0 0 2 2 n n n h π K d λ J λ r d J λ r h    (30f)   31 31 1 0 1 2 n n h K q λ J λ r π   (30g)       2 2 32 11 0 0 0 33 0 0 2 2 n n n h π K d λ J λ r d J λ r h    (30h)       2 2 33 11 0 0 0 33 0 0 2 2 n n n h π K μ λ J λ r μ J λ r h    (30i)
Then, substituting Wn, Hn and Zn obtained in Eqs. [START_REF] Ma | A non-classical Mindlin plate model based on a modified couple stress theory[END_REF], respectively, into Eqs. ( 27), (26b) and (26c) will give the exact solutions of the transverse displacement w, electric potential γ and magnetic potential ζ based on the current circular Kirchhoff plate model for the clamped plate subjected to the axisymmetrically distributed load.

Note that the current Fourier-Bessel series method possesses the same limitations as those of the classical model. For an annular circular plate (with a ≠ 0), the FEM and HDQ techniques have been successfully employed to determine the static deflections of a thin annular plate [START_REF] Yuan | Comparisons of methods for solving static deflections of a thin annular plate[END_REF]. In constant to the inner radius a, the load and boundary conditions are the main factors influencing the solutions.

In order to verify the correctness of the current model, a comparative study of the deflection of the clamped circular Kirchhoff plate subjected to a uniform load between the current model and the model with couple stress effect only provided by Zhang et al. [START_REF] Zhang | A non-classical model for circular Kirchhoff plates incorporating microstructure and surface energy effect[END_REF] are plotted in Fig. 2. For a convenient comparison, the current model is degenerated to a pure elastic model by neglecting both piezoelectric and piezomagnetic effects (i.e., e 31 = q31 = μ11 = μ33 = ϵ11 = ϵ33 = d11 = d33 = 0). The plate parameters, geometric dimensions and boundary conditions of the two models are adopted from Zhang et al. [START_REF] Zhang | A non-classical model for circular Kirchhoff plates incorporating microstructure and surface energy effect[END_REF], and the Young's modulus E and Poisson's ratio v for the isotropic case are related to the current elastic stiffness constant through

11 (1 ) (1 )(1 2 ) v E C v v     (31)
From Fig. 2, it is obvious that the results of the current and classical pure elastic models are the same as those from Zhang et al. [START_REF] Zhang | A non-classical model for circular Kirchhoff plates incorporating microstructure and surface energy effect[END_REF] This validates the current model and indicates that the couple stress effect will always lead to reduced values of the deflection, as expected. Table 1 [START_REF] Zhou | A nonclassical model for circular mindlin plates based on a modified couple stress theory[END_REF][START_REF] Wang | Axisymmetric bending of functionally graded circular magneto-electro-elastic plates[END_REF], in which vf means the volume fraction of BaTiO3. When vf = 0%, the material turns to be a purely piezomagnetic (PM) material CoFe2O4, and when vf = 100%, only a purely piezoelectric (PE) material BaTiO3 is left. In addition, the uniformly distributed constant load q is fixed at 1 MPa.

2 ij ij A C l  (31)
Theoretically, the results will be more accurate as the expansion term n of Fourier Bessel series increases. After a careful study, n = 30 here is sufficient to meet the accuracy requirement. Note that the two sharp changes, as shown in Fig. 7 (with v f = 20%) and Fig. 8 (with vf = 80%), arise directly from the effective material properties of the BaTiO3-CoFe2O4 composite in Table 1. This is similar to what was found in Zhang et al. [START_REF] Zhang | A transversely isotropic magneto-electro-elastic Timoshenko beam model incorporating microstructure and foundation effects[END_REF] based on a thin magneto-electro-elastic beam model.

Summary

A new model for the microstructure-dependent transversely isotropic MEE circular Kirchhoff plates is developed by using the extended modified couple stress theory. The governing equations are derived via a variational method. The newly derived model can capture the microstructure-dependent size effect at the microscale.

To illustrate the new model, the static bending of a clamped transversely isotropic MEE circular plate subject to a uniformly distributed constant load is solved numerically via Fourier-Bessel series. The numerical results show that the values of transverse displacement, electric and magnetic potentials of the current model are always smaller than those of the classical model. However, with the increase of the plate thickness h, the differences are diminishing. In addition, the numerical results

show that the transverse displacement, electric and magnetic potentials vary significantly with the MEE couplings at all length scales.

Fig. 1 .

 1 Fig. 1. The annular circular thin plate's geometry and coordinate system

  , b) It follows from Eqs. (1a-c), (2a, b), (4a-c) and (5a, b) that in the current

  equations, and the boundary conditions in Eqs. (16a-e

21d) 3 . 3 .

 33 Classical Transversely Isotropic Magneto-electro-elastic Circular Kirchhoff Plate Model When the microstructure effect is neglected by setting (A 11-A12) = 0, Eqs. (15a-d)

  22d) as the equilibrium equations, and the boundary conditions in Eqs. (16a-e)

  Hn and Zn are the nth-order coefficients to be determined. Eqs. (26a-c) inherently satisfy the boundary conditions in Eqs. (24c-e) and the symmetry conditions in Eqs. (25b-d).

  , b) Substituting Eqs. (26b-28b) into Eqs. (15b-d

Fig. 2 .

 2 Fig. 2. Comparison of the deflection of the clamped plate subjected to a uniform load

Figures 3 ,

 3 Figures 3, 4 and 5 show the numerical results of transverse displacement w,

Fig. 3 .Fig. 4 .Fig. 5 . 7 .

 3457 Fig. 3. The transverse displacement of the clamped circular Kirchhoff plate

Table 1 .

 1 Effective properties of BaTiO3-CoFe2O4 composite Cij, in GPa, eij, in C/m 2 , ϵij, in 10 9 C 2 /(N•m 2 ), qij, in N/(A•m), dij, in 10 12 Ns/(V•C), μij, in 10 6 Ns 2 /C 2 ,ρ, in kg/m 3 , and l, in µm.

		0%						100%
	vf		20%	40%	50%	60%	80%
		(PM)						(PE)
	C11	286	262	238	226	214	190	166
	C12	173	153.8	134.6	125	115.4	96.2	77
	C13	170	151.6	133.2	124	114.8	96.4	78
	C33	269.5	248	226.5	215.75	205	183.5	162
	C44	45.3	44.84	44.38	44.15	43.92	43.46	43
	e15	0	2.32	4.64	5.8	6.96	9.28	11.6
	e31	0	-0.88	-1.76	-2.2	-2.64	-3.52	-4.4
	e33	0	3.72	7.44	9.3	11.16	14.88	18.6

* Units:

The couple stress stiffness constants can be written as

[START_REF] Qu | A microstructure-dependent anisotropic magneto-electro-elastic Mindlin plate model based on an extended modified couple stress theory[END_REF]

Acknowledgements. The work reported here is funded by the National Natural Science Foundation of China [grant numbers 12002086 and 11672099]. These supports are gratefully acknowledged.

To study the effects of the magneto-electro-elastic coupling, the maximum transverse displacement w (deflection at the midpoint of the plate) under different vf conditions is plotted in Fig. 6a. In addition, Fig. 7a and Fig. 8a present the electric potential and magnetic potential at the midpoint of the plate varying as a function of vf, respectively. The blue square lines describe the results of the current model with h = l, and the red circular lines describe the results of the classical model with h = l.

From Fig. 6a, it is observed that the transverse displacements with h = l predicted by both the current and classical models decrease with the increase of v f. However, the transverse displacement given by the classical model decreases more quickly than that predicted by the current model as vf increases. From Fig. 7a, it is shown that the electric potential of the current model increases with the increase of vf. This variation trend is similar to that of the classical model. In addition, the electric potential of the classical model is always larger than that of the current model, as shown in Fig. 7a.

From Fig. 8a, it is shown that the magnitudes of the magnetic potential of both the current and classical models gradually increase with the increase of v f before reaching their maximum values around 60%, after which the magnetic potentials decrease with vf for the remaining values. When the piezoelectric effect disappears, no electric potential is generated in the entire circular plate. Similarly, there is no magnetic potential generated for the plate when the piezomagnetic effect disappears, as expected.

To further illustrate this effect, Fig. 6b, Fig. 7b and Fig. 8b show the maximum


Shuohui Yin

Ning Zhang

Peng Liu

Jingang Liu

Tiantang Yu
email: tiantangyu@hhu.edu.cn

Shuitao Gu

Yu Cong

Dynamic fracture analysis of the linearly uncoupled and coupled physical phenomena by the variable-node multiscale XFEM

Keywords: Dynamic fracture, Piezoelectric structures, V-XFEM, Variable-node elements, Dynamic stress intensity factor piezomagnetic moduli

This paper aims to develop a variable-node multiscale extended finite element method (V-XFEM) for dynamic fracture analysis of the linearly uncoupled and coupled physical phenomena in a compact formula. The general governing equations for the linearly uncoupled and coupled physical phenomena are presented in a compact form. The local mesh refinement technique for modeling cracks is used to improve the accuracy and efficiency, in which variable-node elements without modifying the system matrix or impose additional boundary conditions are taken to connect/link different scale elements. In addition, the time-dependent equations are solved by the unconditionally stable implicit Newmark time integration method, and the dynamic intensity factors (DIFs) are derived from the domain forms of the interaction integrals. Numerical results of the linearly elastic and piezoelectric problems show that V-XFEM is an efficient numerical approach to simulate the dynamic fracture problems of the linearly uncoupled and coupled physical phenomena.

An accurate description of the singular fields near the crack-tip is often of great importance in modeling of fracture mechanics problems. However, in contrast to the static loading conditions, numerical simulation of the dynamic fracture problems remains a challenge in many practical engineering applications. Song et al. [START_REF] Song | A comparative study on finite element methods for dynamic fracture[END_REF] investigated three finite element methods for dynamic crack propagation, it shows that XFEM and interelement method show similar crack speeds and crack paths, and element deletion method is unable to predict crack branching. The extended finite element method (XFEM) [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] is a very powerful numerical tool for modeling arbitrary cracks without remeshing. Recently, it has received much attention in a wide range of engineering problems including the dynamic fracture mechanics problems, e.g., see [START_REF] Bui | Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading[END_REF][START_REF] Liu | Transient dynamic crack analysis in non-homogeneous functionally graded piezoelectric materials by the X-FEM[END_REF][START_REF] Liu | Transient thermal shock fracture analysis of functionally graded piezoelectric materials by the extended finite element method[END_REF]. A detailed overview on the application of extended finite element method in modeling of arbitrary discontinuities can be found in Refs. [START_REF] Fries | The extended/generalized finite element method: An overview of the method and its applications[END_REF][START_REF] Li | A review of the extended finite element method on macrocrack and microcrack growth simulations[END_REF].

In fracture mechanics problems, a fine mesh is often required in order to accurately describe the singular fields near the crack-tip and model crack growth path. However, the region containing cracks is usually only a small portion of the whole structure, it is unreasonable to adopt a small-scale mesh in the whole structure due to the computational efficiency. Therefore, numerous multiscale methods which adopt a locally refined mesh have been proposed. In multiscale methods, two-scale meshes are nonconforming, how to link or connect meshes at different level and scales is the main challenge. Guidault et al. [START_REF] Guidault | A two-scale approach with homogenization for the computation of cracked structures[END_REF][START_REF] Guidault | A multiscale extended finite element method for crack propagation[END_REF] decomposed the structure into two substructures with distinct scales and linked the two distinct scale meshes by transmitting interactions at multiscale interfaces. Based on a two-scale decomposition, Loehnert and Belytschko [START_REF] Loehnert | A multiscale projection method for macro/microcrack simulations[END_REF] proposed a multiscale projection method for macro/microcrack simulations in which macrocracks and microcracks are taken into account by the XFEM. Zhou and Yang [START_REF] Zhou | Multiscale numerical modeling of propagation and coalescence of multiple cracks in rock masses[END_REF] and Holl et al. [START_REF] Holl | 3D Multiscale crack propagation using the XFEM applied to a gas turbine blade[END_REF][START_REF] Holl | An adaptive multiscale method for crack propagation and crack coalescence[END_REF] adopted the multiscale projection method to simulate propagation and coalescence of multiple cracks. Murotani et al. [START_REF] Murotani | Adaptive finite elements using hierarchical mesh and its application to crack propagation analysis[END_REF] proposed an adaptive finite element using a hierarchical mesh and applied it to crack propagation analysis. Hettich et al. [START_REF] Hettich | Modeling of failure in composites by X-FEM and level sets within a multiscale framework[END_REF] proposed a variational multiscale method to model failure, in which the Lagrange multipliers was adopted to connect the large and small scales meshes. Dhia and Rateau [START_REF] Dhia | On the use of XFEM within the Arlequin framework for the simulation of crack propagation[END_REF] used the Arlequin method to link fine scale model and a large scale model for simulating crack propagation through the XFEM. For 3D crack propagation in large structures, Rannou et al. [START_REF] Rannou | A local multigrid X-FEM strategy for 3-D crack propagation[END_REF] developed a local multigrid approach to coupe the grids at different scales within the XFEM framework. Gibert et al. [START_REF] Gibert | A 3D automatic mesh refinement X-FEM approach for fatigue crack propagation[END_REF] presented a multi scale approach coupling the X-FEM method with automatic adaptive mesh refinement for 2D and 3D fatigue crack propagation.

But, a projection of the mechanical state from one mesh to another during the propagation is necessary. Fries et al. [START_REF] Fries | Hanging nodes and XFEM[END_REF] proposed a XFEM with locally refined mesh through two kinds of hanging nodes which have or have no degrees of freedom. Budarapu et al. [START_REF] Budarapu | An adaptive multiscale method for quasi-static crack growth[END_REF] developed an adaptive multiscale method in which the different scales are connected by ghost atoms. Patil et al.

[22] constructed a multiscale XFEM for modeling heterogeneous materials using the multiscale basis functions as a bridge between two scales. Teng et al. [START_REF] Teng | An adaptively refined XFEM with virtual node polygonal elements for dynamic crack problems[END_REF] integrated the virtual node polygonal elements with XFEM for modelling dynamic crack growth. Kumar et al. [START_REF] Kumar | A homogenized multigrid XFEM to predict the crack growth behavior of ductile material in the presence of microstructural defects[END_REF] developed a homogenized multigrid XFEM for the crack growth simulations in ductile materials, in which the domain contains the microstructural defects is simulated by fine mesh while other domains are simulated by coarse mesh with homogenized properties. Among those methods, extra efforts such as modification on the system matrix or imposing additional boundary conditions are unavoidable.

To overcome this numerical problem, the variable-node elements [START_REF] Lim | Variable-node element families for mesh connection and adaptive mesh computation[END_REF] with an arbitrary number of nodes on each of their edges were developed based on the generic point interpolation, which can provide a flexibility to connect different scale meshes without modify the system matrix or impose additional boundary conditions. Recently, the variable-node XFEM (V-XFEM) has been developed for modeling cracks, voids, and inclusions in both 2D and 3D elastic mediums [START_REF] Wang | Numerical modeling of 3-D inclusions and voids by a novel adaptive XFEM[END_REF][START_REF] Wang | 3-D Local mesh refinement XFEM with variable-node hexahedron elements for extraction of stress intensity factors of straight and curved planar cracks[END_REF][START_REF] Yu | Numerical simulation of 2-d weak and strong discontinuities by a novel approach based on XFEM with local mesh refinement[END_REF] by introducing the variable-node elements into the traditional XFEM codes. Due to its advantages of having no additional constraints, the variable-node elements can be conveniently embedded into the traditional XFEM codes. The V-XFEM has been proved to be a general and efficient approach in the framework of multiscale methods, and achieves a faster convergence rate in energy norm error than the standard XFEM [START_REF] Yu | Numerical simulation of 2-d weak and strong discontinuities by a novel approach based on XFEM with local mesh refinement[END_REF]. Ding et al. [START_REF] Ding | An efficient variable-node XFEM for modeling multiple crack growth: A Matlab object-oriented implementation[END_REF][START_REF] Ding | Modeling strong/weak discontinuities by local mesh refinement variable-node XFEM with object-oriented implementation[END_REF] present a Matlab object-oriented implementation of an efficient V-XFEM for the problems of multiple crack growth simulation and strong and weak discontinuities. Han et al. [START_REF] Han | Extended stochastic finite element method enhanced by local mesh refinement for random voids analysis[END_REF] proposed a triangular extended stochastic finite element method (T-XSFEM) for simulation of random void problems with the aid of variable-node elements to couple/link different mesh-scales. More recently, the V-XFEM are extended for detection of multiple complicated flaw clusters in large structures [START_REF] Ma | Detection of multiple complicated flaw clusters by dynamic variable-node XFEM with a three-step detection[END_REF]. Even much efforts in the fracture analyses, to the best of the authors' knowledge, the attention has not been paid to the dynamic fracture problems.

Piezoelectric materials have been widely used as actuators, transducers or sensors because of the inherent coupling effect between mechanical and electric fields. Piezoelectric materials possess many advantages but suffer the brittleness and low fracture toughness. As the growing demand for wide range of applications of piezoelectric materials, the dynamic fracture problem has been a serious concern nowadays for the design engineers of piezoelectric structures. Li and Mataga [START_REF] Li | Dynamic crack propagation in piezoelectric materials[END_REF] studied transient response of a semi-infinite, anti-plane crack propagating in a hexagonal piezoelectric medium by transform methods together with the Wiener-Hopf and Cagniardde Hoop techniques. However, deriving analytical solutions for time dependent problems of cracks under dynamic impact loading condition are very restricted and numerical methods are required for practical application of piezoelectric structures. Therefore, Enderlein et al. [START_REF] Enderlein | Finite element techniques for dynamic crack analysis in piezoelectrics[END_REF] applied finite element method (FEM) to model dynamic fracture problems in piezoelectric materials. Studies on fracture problems in piezoelectric materials under dynamic loading have also modeled by other numerical approaches including the BEM [START_REF] Garcia-Sanchez | 2D transient dynamic crack analysis in piezoelectric solids by BEM[END_REF][START_REF] Garcia-Sanchez | 2-D transient dynamic analysis of cracked piezoelectric solids by a time-domain BEM[END_REF], meshless methods [START_REF] Wunsche | A 2D time-domain collocation-Galerkin BEM for dynamic crack analysis in piezoelectric solids[END_REF], scaled boundary finite element method (SBFEM) [START_REF] Li | 2D dynamic analysis of cracks and interface cracks in piezoelectric composites using the SBFEM[END_REF] and the extended finite element method (XFEM) [START_REF] Bui | Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading[END_REF][START_REF] Liu | Transient dynamic crack analysis in non-homogeneous functionally graded piezoelectric materials by the X-FEM[END_REF][START_REF] Liu | Transient thermal shock fracture analysis of functionally graded piezoelectric materials by the extended finite element method[END_REF][START_REF] Yu | Interfacial dynamic impermeable cracks analysis in dissimilar piezoelectric materials under coupled electromechanical loading with the extended finite element method[END_REF]. Bechet et al. [START_REF] Bechet | Application of the X-FEM to the fracture of piezoelectric materials[END_REF] devised an enrichment scheme with six crack tip functions based on the asymptotic expansion around a crack tip in piezoelectric materials.

In the present paper, we extend the applicability of V-XFEM to the dynamic fracture analysis for the linearly uncoupled and coupled physical phenomena in a compact formula by introducing the general governing equations. The time-dependent equations to be solved by the unconditionally stable implicit Newmark time integration method and the dynamic intensity factors (DIFs) are derived from the domain forms of the interaction integrals. The numerical results of the elastic and piezoelectric problems considered as illustrated examples are discussed in order to show the efficiency and accuracy of V-XFEM. Thus, the main features of the present work consists of: (i) extension of V-XFEM to the dynamic fracture simulation, (ii) treatment of problem for all the linear uncoupled phenomena as an uncoupled problem in a unified and compact way.

The rest of the paper is structured as follows. The formulation of general governing equations is given in Section 2. Considering piezoelectric structures as illustrated examples, the V-XFEM is elaborated in detailed for dynamic crack problem including solution of dynamic fracture, variable-node elements and evaluation of generalized dynamic intensity factors in Section 3.

Section 4 the numerical illustrations are given for several elastic and piezoelectric problems, which show that the obtained solutions agree well with results available in the literature and V-XFEM has more efficient than the standard XFEM. Finally, some conclusions are summarized in Section 5.

General governing equations

In this section, we present general governing equations for the linearly uncoupled and coupled physical phenomena. Firstly, we introduce two 3D vector fields a J and a E written by ( )

, , i i div J t a a a g = = J (1) a a w = Ñ E or , . i i E a a w = (2)
Then, the constitutive relations for all the known linearly uncoupled phenomena like elasticity and the linearly coupled phenomena like piezoelectricity and magneto-electro-elasticity can be written in a unified way (see, e.g. [START_REF] Benveniste | New exact results for the effective electric, elastic, piezoelectric and other properties of composite ellipsoid assemblages[END_REF][START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF])

a ab b = J L E or i i j j J L E a a b b = (3)
where the Greek letters a and b indicate field subscripts assuming the values 1 to r, the Latin letters i and j represent space subscripts taking the values 1 and 2 for two dimensional problem or 1, 2 and 3 for three dimensional problem, and a w are the components of an r-dimensional potential vector w . It is often rewritten as

= J E ( 4
)
where J and E can be viewed as two vectors of 3r components and is a 3r×3r matrix such

é ù ê ú ê ú = ê ú ê ú ë û L L L L L L L L L , 1 2 r é ù ê ú ê ú = ê ú ê ú ë û E E E E (5)
Here, an element ab L of is a 3×3 matrix. In addition, the matrix is invertible, so that the dual expression of (4) is E= J (6) with 1 -= . It should be noted that these above formulations allow us to systematically treat the problem for all the linear uncoupled phenomena as an uncoupled problem.

Let us consider the magneto-electro-elasticity as an example for the foregoing formulation Eq. ( 4). The constitutive law of this problem is usually defined by (see, e.g. [START_REF] Alshits | On the existence of surface waves in half-infinite anisotropic media with piezoelectric and piezomagnetic properties[END_REF])

ij ijkl kl kij k kij k i ikl kl ij j ij j i ikl kl ij j ij j C h a d D h K a m K h H a s e e e ì = -Õ -L ï = Õ + + í ï = L + + î (7)
where ij s and ij e are the components of the stress and strain tensors; di and hi are the components of the electric displacement and electric field vectors; mi and ai are the components of the magnetic flux and magnetic field vectors. Here, ijkl C , ij D and ij H are the elastic, dielectric and magnetic permeability moduli; kij Õ , kij L and ij K are the piezoelectric, piezomagnetic and magnetoelectric moduli. These material moduli have the following symmetry properties:

ijkl ijlk klij C C C = = , ij ji D D = , ij ji H H = , ij ji K K = , kij kji Õ = Õ , kij kji L = L (8)
The strain, electric field, and magnetic fields can be derivable from a displacement field u, electric potential j , and magnetic potential y , respectively, as follows:

( )

, , 1
2 ij i j j i u u e = + , , i i h j = -, , i i a y = - (9)
In addition, for simplifying our analysis, inertia effects are considered and other unsteady effects are ignored. Thus, the stress tensor , the electric displacement vector d and the magnetic flux vector b satisfy:

2 , 2 i ij j u t s r ¶ = ¶ , , 0 i i d = , , 0 i i b = (10)
The magneto-electro-thermo-elastic constitutive law ( 7) can be recast into the compact form Eq. ( 5)

as a coupled field phenomena with five coupled fields ( 5 r = ). More precisely, the following definitions are adopted:

1 1 i i J s = , 2 2 i i J s = , 3 3 i i J s = , 4i i J d = , 5i i J b = ( 11
) , i i E a a w = , ( )
1 2 3 T u u u j y w = (12) 1 1 1 2 1 3 1 1 2 1 2 2 2 3 2 2 3 1 3 2 3 3 3 3 1 2 3 1 2 3 i j i j i j ji ji i j i j i j ji ji i j i j i j ji ji i j i j i j ij ij i j i j i j ji ij C C C C C C C C C L D K K H ab Õ L é ù ê ú Õ L ê ú ê ú Õ L é ù = ë û ê ú Õ Õ Õ - - ê ú ê ú L L L - - ë û = (13) with ( ) 2 2 ,0,0 . u t t g r ae ö ¶ ç ÷ ¶ è ø i = ( 14
)
The general prescribed conditions on the boundary S consisting of SJ and Sw can be expressed by

( ) ( ) i i J w J J n J on S and on S a a a a a = = = w w (15)
where n is the unit outward normal on S. And, the crack-faces in this study are assumed to be general traction-free, i.e.

0 i i J a u = ( 16
)
Here, is the unit vector normal to the crack-faces Sc .

The weak form of the general governing equations is gained from Hamilton's variational principle:

( ) ( )

1 1 0 0 2 1 1 [ ] 0 2 2 dVdt J dt a a d r d - + ò ò ò t t t V t u E E w = ( 17
)
After doing some basic calculus, the above equation finally yields the weak form as follows:

( )

1 1 0 0 2 [ ] 0 dVdt J dt a a r d d d + ò ò ò t t t V t - u u + E E w = (18)
Here, the superscript • means the derivative with respect to time.

It should be noted that the general weak form for the linearly uncoupled and coupled physical phenomena is expressed in a compact form. In the following, the elastic and piezoelectric problems will be considered as linearly uncoupled and coupled physical examples, respectively, for the dynamic fracture analysis.

V-XFEM for dynamic crack problem in linearly elastic and piezoelectric structures

Now, we proceed to elaborate a V-XFEM to simulate the dynamic fracture in the problem formulated in the previous section. Thus, we particularize the general governing equations to the linearly elasticity and piezoelectricity. It should be noted that the linearly piezoelectric problems can reduce to linearly elastic problems by ignoring the coupled piezoelectric tensor and suppressing the electric field and the electric displacements. In order not to be redundant, the corresponding formulations associated to piezoelectricity are only described in the following.

Enriched field approximation

In the X-FEM [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], the standard finite element approximation is enriched by additional functions based on the partition of unity to model the field discontinuities. For piezoelectric materials, Bechet et al. [START_REF] Bechet | Application of the X-FEM to the fracture of piezoelectric materials[END_REF] devised an enrichment scheme with six crack tip functions. A comparison of the normalized dynamic intensity factors with the six crack tip enrichments and the standard four crack tip enrichments is presented in this paper. It shows that the results are matched very well and consistent with the conclusion of Bechet et al [START_REF] Bechet | Application of the X-FEM to the fracture of piezoelectric materials[END_REF] that the four-fold enrichment is almost as efficient, concerning accuracy both in energy and in the SIFs. However, the four-fold enrichment is simpler to implement and involves less computational overhead, because it adds only four degrees of freedom (dofs) per regular dof, instead of six. Therefore, the four-fold enrichment is adopted in this paper. With general governing equations, the enriched field ( ) h a w x can be expressed in the form:

( ) ( ) ( ) ( ) ( ) ( ) 4 1
.

s cut tip h i j k i j k I I I i N j N k N N N H N F a a a a a a a = Î Î Î = + + å å å å w w x x x x a x x b ( 19
)
In the above equation, Obviously, in the piezoelectric problem, four vector fields (r = 4) involves with ( )

1 2 3 T u u u a j w = .
In this case, the enriched displacements and electrical potential approximation for the crack problem involve. That means Eq. ( 19) can be rewritten by

( ) ( ) ( ) ( ) ( ) ( ) 4 1 , s cut tip h i i j j k I kI I i N j N k N N N H N F = Î Î Î = + + å å å å u x x u x x a x x b (20a) ( ) ( ) ( ) ( ) ( ) ( ) 4 1
.

s cut tip h i i j j k I kI I i N j N k N N u N H a N F j * * * = Î Î Î = + + å å å å x x x x x x b (20b)
In addition, with reference to these works [START_REF] Bui | Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading[END_REF][START_REF] Sharma | Analysis of a subinterface crack in piezoelectric bimaterials with the extended finite element method[END_REF], the corresponding crack-tip branch enrichment functions are taken as [START_REF] Budarapu | An adaptive multiscale method for quasi-static crack growth[END_REF] where r and q are the usual crack-tip polar co-ordinates.

( ) { } 4 1 sin cos sin sin cos sin 2 2 2 2 I I F r r r r q q q q q q = ì ü í ý î þ = x
We need to emphasize that Eq. ( 19) can similarly be reduced to some prescribed formulation associated to the linearly physical phenomena under consideration.

Numerical discretization

In order to discretize the linearly piezoelectric problem, in the case of neglecting of damping, substituting the enriched approximated functions in Eq. ( 20) into the weak-form of the equilibrium equations as illustrated in [START_REF] Enderlein | Finite element techniques for dynamic crack analysis in piezoelectrics[END_REF] and after some appropriate manipulations, a system of discretized piezoelectric XFEM equations can be obtained as follows:

+ = , M K F ( 22
)
where denotes the nodal unknowns, and M , K and F are the generalized global mass matrix, global stiffness matrix and external force vector, respectively, with 19) into the weak-form of the general equilibrium equations, a system of discretized equations for all the linear uncoupled phenomena as an uncoupled problem can be obtained by introducing the corresponding number of physical fields.

0 0 0 = é ù ê ú ë û M M , = ì ü í ý î þ , = mech elec ì ü ï ï í ý ï ï î þ F F F
In the linearly piezoelectric problem, the element contribution to mass and stiffness matrices (superscript e) are as follows:

uu ua ub ij ij ij e au aa ab ij ij ij ij bu ba bb ij ij ij é ù ê ú = ê ú ê ú ë û m m m m m m m m m m ; uu ua ub ij ij ij e au aa ab ij ij ij ij bu ba bb ij ij ij é ù ê ú = ê ú ê ú ë û k k k k k k k k k k (23)
where

( ) = T uu ij i j d r W W ò m N N (24a) ( ) 2 = T aa ij i j H d r W W ò m N N (24b) ( ) ( ) 2 T bb ij I i j F d r W W ò = m N N ; ( 1, 2,3, 4 I = ) (24c) ( ) = T ua ij i j H d r W W ò m N N (24d) ( ) T ub ij I i j F d r W W ò = m N N ; ( 1, 2,3, 4 I = ) (24e) ( ) T ab ij I i j HF d r W W ò = m N N ; ( 1, 2,3, 4 I = ) (24f) and ( )
= e T rs r s ij i j d W W ò k B DB ; ( , , , = r s u a b ) (25) with , , , , , , 0 0 0
0 0 0 0 0 0 i x i y u i y i x i i x i y N N N N N N é ù ê ú ê ú ê ú = ê ú ê ú ê ú ë û B ( 26a) , , , , , , 0 0 0
0 0 0 0 0 0 i x i y a i y i x i i x i y HN HN HN HN HN HN é ù ê ú ê ú ê ú = ê ú ê ú ê ú ë û B ( 26b
) 1 2 3 4 b b b b b i i i i i é ù = ë û B B B B B (26c) ( ) ( ) ( ) ( ) ( ) ( ) , , , , , , 0 0 0
0 0 0 0 0 0 i I x i I y bI i I i I i y x i I x i I y N F N F N F N F N F N F é ù ê ú ê ú ê ú = ê ú ê ú ê ú ê ú ë û B ( 26d
)
where D is generalized elastic coefficient matrix in eq. ( 13).

Additionally, the element contribution to external force vector mech F is given by

T mech u a b i i i i é ù = ë û f f f f (27) in which , e u i i d ¶W = G ò f N t (28a) , e a i i H d ¶W = G ò f N t (28b) e bI I i i F d ¶W = G ò f N t ; ( 1, 2, 3, 4. I = ) (28c)
The external nodal charge vector elec F can be computed in a similar manner as above.

It is noted that the Gauss quadrature scheme is employed for non-enrich elements and 4×4

Gauss quadrature is employed for blending elements. To obtain an accurate integration for crack tip elements and elements cut by crack the triangular sub-domain technique [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF] is used in the same way as that of the XFEM.

In this study, the unconditionally stable implicit Newmark method is adopted to solve the discrete dynamic equilibrium equations of Eq. ( 22). At time step n, the discretized equations Eq. ( 22)

is thus rewritten as ( ) ( )

2 2 1 1 1 1 2 2 n n n n t t t b b - - - é ù D D = + D + - ê ú ë û + M K F -K ( 29
)
in which the velocity and displacement can be evaluated by ( )

1 1 1 n n n n t a a - - é ù - + D ë û = + (30a) 2 1 1 1 1 2 n n n n n t t b b - - - é ù ae ö D - + D ç ÷ ê ú è ø ë û = + + ( 30b
)
where t D denotes the time increment at the present time-step and the unconditionally stable parameter are 0.5 a = and 0.25 b = .

Variable-node elements

In order to improve computing efficiency and accuracy, multiscale mesh, i.e., mesh locally refined in the vicinity of cracks, is used and the variable-node elements [START_REF] Teng | An adaptively refined XFEM with virtual node polygonal elements for dynamic crack problems[END_REF][START_REF] Wang | Numerical modeling of 3-D inclusions and voids by a novel adaptive XFEM[END_REF] are adopted to as linking elements (or transition elements) between large scale elements and small scale elements, which is schematically sketched in Fig. 1. The advantage of variable-node elements is any number of nodes can be attached to the edges of the quadrilateral isoparametric element which is shown in Fig. 2 and the interpolation property between any two adjacent nodes remains linear as the number of nodes in an element increases.

Therefore, variable-node elements are suitable for linking between small scale mesh and large scale mesh. The approximation displacement ( )

h
x u of the variable-node element is given as:

( ) ( ) ( ) 1 = p N h T i i i= = å u N u a p ( 31
)
where p N is the number of nodes; 0 0

= i i i N N é ù ê ú ë û
N is the shape function matrix with i N is the shape function of the i-th node;

[ ]

= T i i i u v u
is the i-th nodal displacement vector; T a is the 2 p N ´ unknown coefficients matrix and ( ) p is a column vector of the polynomial basis given by [START_REF] Teng | An adaptively refined XFEM with virtual node polygonal elements for dynamic crack problems[END_REF] ( )

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 5

xh x x h h x x h h h h x x h h x x x x h h x x h h h h x

+ + + + + + + + é ê ê ê ê ê ê - + ê ê ê - + ê ê - + ê = ê ê - + - + - + - + - + ë p ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú ê ú û ( 32
)
in which x and h are the local coordinates in the isoparametric element.

Using the point interpolation, the approximated displacements is rewritten as

( ) ( ) ( ) 1 = h T T - = u a p U q p ( 33
)
where

( ) ( ) 1
- + + + + = N N q p ( 34
)
The shape functions for electrical potential field j can be computed in a similar manner as above.

Additionally, according to Eq. ( 34), the slope discontinuity of the shape functions appears at extra nodes. Herein, the integration technique applied to sub-triangles [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF] for cracked elements in the XFEM is used for variable-node elements.

In this study, we define refined zone artificially instead of based on posteriori error estimation.

The reason is that it would be cumbersome to compute posteriori error in every step of crack growth.

In each time step, mesh is refined around cracks by partitioning original large-scale elements. The crack element is easy to obtain by level set method. And the variable-node elements are generated accompanied with the mesh refinement. Once the small-scale elements are created, variable-node elements should be defined at the outermost layer of the refined area to work as transition between small and large elements.

Evaluation of generalized dynamic intensity factors

In the present work, the domain-form of the contour interaction integral is used to accurately calculate the generalized dynamic intensity factors in the piezoelectric materials, which can be found in [START_REF] Bui | Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading[END_REF] for details. The path independent electro-mechanical J-integral for a cracked piezoelectric materials is given by

( ) ( ) ,1 ,1 1 , ,1 ij i j j j i A A J u D W q dA uu qdA s j d r = + - - ò ò (35)
where

( ) / 2 = ij ij j j W D E
s eis the electric enthalpy density; A is the area inside an arbitrary contour enclosing the crack-tip and q is a smooth weight function, which has a value of unity at the crack-tip, zero along the boundary of the domain A, and a smooth linear variation in-between.

Let us now consider two independent dynamic equilibrium states of a piezoelectric cracked body. The first state

( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 1 1
, , , , ij ij j j j u D E s e corresponds to the actual state, whereas the second one

( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 2
, , , , ij ij j j j u D E s e is an auxiliary state. The interaction integral for the two states is obtained by superposition of these two states which leads to another dynamic equilibrium state.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1,2 1 2 2 1 1 2 2 1 1,2 1 2 ,1 ,1 ,1 ,1 1 , ,1 ij ij j j j j i A A I u u D D W q dA u u qdA s s j j d r = + + + - + ò ò ( 36
)
where

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1,2 1 2 2 1 2 1 1 2 1 2 = ij ij ij ij j j j j W D E D E s e s e + - - . ( 37
)
For linear piezoelectric solids, the electro-mechanical J-integral is equal to the energy release rate. Applying J-integral to the two states (1) and ( 2) and then the interaction integral ( )

1,2 I can be rewritten as ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
I K K Y K K Y K K Y K K K K Y K K K K Y K K K K Y + + + + + + + + (38)
By choosing the auxiliary state 2 as the crack opening mode with ( )

2 1 I K = ; ( ) ( ) 2 2 0 II IV K K =
= , then the interaction integral is obtained as

( ) ( ) ( ) ( ) 1, 1 1 1 22 12 23 = I I II IV I K Y K Y K Y + + (39)
Similarly, other modes can be obtained as

( ) ( ) ( ) ( ) 1, 1 1 1 12 11 13 = II I II IV I K Y K Y K Y + + (40) ( ) ( ) ( ) ( ) 1, 1 1 1 23 13 33 = IV I II IV I K Y K Y K Y + + (41)
As a consequence, the generalized dynamic intensity factors are finally obtained by solving the following linear algebraic equations

( ) ( ) ( ) ( ) ( ) ( ) 1 1
= I I II II IV IV K I Y Y Y K Y Y Y I Y Y Y I K ì ü ì ü é ù ï ï ï ï ï ï ï ï ê ú í ý í ý ê ú ï ï ï ï ê ú ë û ï ï ï ï î þ î þ (42)

Numerical results

In this section, several benchmark numerical examples for stationary dynamic cracks in linear elastic and piezoelectric structures are presented to illustrate the accuracy and efficiency of the elaborated V-XFEM. The accuracy and efficiency is numerically confirmed through the comparison of the present approach with analytical solutions and the standard XFEM. Plane-strain condition and the impermeable crack-face boundary condition are assumed throughout the study.

Dynamic fracture of edge crack in the semi-infinite domain

The first numerical example for dynamic fracture analysis deals with a semi-infinite mode-I crack loaded by a tensile stress perpendicular to the crack-face. Due to the computational limits, a limited domain is modeled for analysis, whose geometry is sketched in Fig. can be written as [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF]:

( ) ( )( ) 0 0 0, 1 2 2 1 c dyn I d c c t t K t c t t t t n s n p < ì ï = í - - ³ ï - î (43)
The relative error of normalized dynamic stress intensity factor is defined as:

100% analytic I I analytic I K K err K - = ´ (44)
where I K and analytic I K are the normalized dynamic stress intensity factor of numerical solution and analytical solution, respectively.

Validation and accuracy

In order to validate the elaborated V-XFEM and demonstrate its performance in modeling dynamic fracture problems, Fig. 4 and Fig. 5 present, respectively, the comparative study of the normalized mode-I DSIFs and their percentage errors with the standard XFEM and the analytical solution. In variable-node XFEM, the original coarse mesh use 65×25 elements as shown in Fig. As shown in Figs. 4 and5, the normalized values of the DSIFs achieved by the V-XFEM match well to the analytical solution and the standard XFEM with the small mesh of 195×75 elements (Dof is 30212). Nevertheless, the V-XFEM with 6812 Dofs has a higher accuracy than standard XFEM with 3592 Dofs (65×25 elements) and 6908 Dofs (79×41 elements). It can be concluded that the V-XFEM has a higher accuracy than standard XFEM.

Numerical efficiency

For comparison of the numerical efficiency, Table 1 presents a comparative study of the CPU-time consumed for three standard XFEM meshes and the variable-node XFEM mesh. As can be seen in Table 1 that the variable-node XFEM (see Fig. 6(b)) needs less computational cost than the standard XFEM with fine mesh (see Fig. 6(c)).

Dynamic fracture of the bevel edge crack plate

In order to prove the applicability and reliability of the V-XFEM for dynamic fracture problems, a mixed-mode crack is thus considered. A rectangular plate with an inclined edge crack as depicted in Fig. 7 , where H(t) represents the Heaviside step function. The total computation time is 3×10 -5 s and the number of iterations is 50 times. As shown in Fig. 8(a), the variable-node mesh with 15×8 coarse elements and 3×3 small scale elements in the vicinity of crack is used for dynamic fracture problems with inclined edge crack. For validation purpose, the computed result is then compared with those obtained from standard XFEM with fine mesh of 45×25 elements as shown in Fig. 8(b). Figs. 9 shows that the normalized mixed-mode DSIFs achieved by the V-XFEM agree well to XFEM with a fine mesh. are matched very well and consistent with the conclusion of Bechet et al [START_REF] Garcia-Sanchez | 2-D transient dynamic analysis of cracked piezoelectric solids by a time-domain BEM[END_REF] that the four-fold enrichment is almost as efficient, concerning accuracy both in energy and in the SIFs. However, the four-fold enrichment is simpler to implement and involves less computational overhead, because it adds only four degrees of freedom (dofs) per regular dof, instead of six. Therefore, the four-fold enrichment is adopted in the rest of this paper. The problem is solved by variable-node XFEM with 19×39 coarse elements and 3×3 small scale elements in the vicinity of crack (see in Fig. 12a), and standard XFEM with a regular fine mesh of 50×100 quadrilateral elements for comparison (see in Fig. 12b). We first consider the plate only subjected to an impact mechanical loading i.e. =0 l . Fig. 13 presents a comparison of the V-XFEM results for the normalized dynamic stress intensity factors For the fracture behavior of piezoelectric materials, the most important and interesting issue is the investigation of the effects induced by both the mechanical and electrical loading on the dynamic fracture parameters. Here, we consider the same example but subjected simultaneously to a combined mechanical and electrical impact load. First, the intensity of the electrical impact =1 l is studied and the computed results for the normalized dynamic intensity factors are shown in Fig. 15. A good agreement with each other is found, which further confirms the effective of the present variable-node XFEM. Next, we discussed the effects of the intensity of the electrical impact loading on the normalized dynamic intensity factors. The electrical loading parameter l is varied and the normalized dynamic intensity factors are evaluated individually and then depicted in Fig. 16. It can be seen that the maximum values of the normalized dynamic intensity factors are reduced as increase l and the electrical impact affects the normalized dynamic stress intensity factors significantly. It also can be concluded that the normalized dynamic electrical displacement intensity factor seems weakly dependent on the time but has a strong dependence on the load parameter l , which is also a consequence of the quasi-electrostatic assumption of the electrical field. Finally, the influence of the orientation of the material poling direction on the normalized dynamic intensity factor is discussed here. Fig. 17 shows a comparison of the normalized dynamic intensity factors derived from both methods for different poling direction. It can be observed in Fig. 17 that the normalized dynamic intensity factors have a significant dependence on the poling angle q . Beside, the electrical displacement intensity factors increase as the poling angle increase, and they are equal to zero when q = 90°, due to the piezoelectric effect vanishes for a crack parallel to the poling direction. piezoelectric material PZT-5H is taken for this study and the poling angle q = 90°. Both mechanical and electrical loadings are considered in the study. Although there are lot of crack growth criteria [START_REF] Xu | A theoretical study of branched cracks in piezoelectrics[END_REF][START_REF] Lei | A simplified evaluation of the mechanical energy release rate of kinked cracks in piezoelectric materials using the boundary element method[END_REF] available to determine the direction of crack propagation in structural materials, appropriate exact crack growth criteria for piezoelectric solids are still being developed due to their complicated inherent coupled electro-mechanical characteristics and material anisotropic features. As the main aim of this paper is to implement V-XFEM approach for modeling crack growth for piezoelectric solids, the modified hoop SIF criterion [START_REF] Mishra | Thermo-electro-mechanical fatigue crack growth simulation in piezoelectric solids using XFEM approach[END_REF] is utilized as the fracture criterion in order to predict crack growth direction. The local refined area is obtained by extending 3 layers from initial elements where crack locates and then 3×3 small scale elements are generated in the refined zone. The crack increment length is set to be 0.12 which is about twice the initial element size. The crack growth profiles for the first, second, fourth and the fifth crack increment steps which corresponds to critical crack length is presented in Fig. 19. It can be see that the V-XFEM can effectively simulate crack propagation.

Conclusions

This paper elaborates the V-XFEM to simulate dynamic fracture problems subjected to impact loads for linearly uncoupled and coupled phenomena. The linearly elastic and piezoelectric problems are considered as illustrated examples in order to show the efficiency and accuracy of V-XFEM. A dynamic XFEM integrated with variable-node elements as well as the implicit time integration scheme is developed for this purpose. An interaction integral for linear piezoelectric materials utilizing the domain-form is implement to extract the relevant dynamic intensity factors.

To illustrate the accuracy and efficiency of the proposed method, numerical results for the normalized dynamic intensity factors are presented and compared with the standard XFEM with refined mesh and analytical solution. The effects of the combined mechanical and electrical impacts, polarization direction on the normalized dynamic intensity factors is analyzed and discussed in details. Some important features are concluded as follows:

l The numerical results predicted by the V-XFEM agree well with analytical solution and standard XFEM with refined mesh, as expected.

l It can be concluded that the V-XFEM has a higher accuracy and more and efficiency than standard XFEM.

l For piezoelectric structures, numerical results show that the normalized dynamic intensity factors have a significant dependence on the poling direction and electrical impact loading.

l Because of a unified way for all single and coupled linear phenomena, the proposed V-XFEM has the potential to simulate all the linearly physical phenomena.

  nodes in the whole discretized domain, associated with crack-faces and the crack-tips

N

  a x are the standard finite element shape functions; i a w are the nodal vectors defined in standard finite elements; j a N are the set of all nodes in the whole discretized domain, the set of enriched nodes associated with crack-faces and the set of the enriched nodes associated with the crack-tips, respectively;( ) H x is a modified Heaviside step function (jump function), which takes the value of +1 above the crack while -1 below the crack; ( ) I F x is the crack-tip branch enrichment functions.

  . Analogously, substituting the enriched approximated functions in Eq. (

Fig. 1 .

 1 Fig. 1. Illustration of multiscale mesh.

Fig. 2 .

 2 Fig. 2. Illustration of an (4 + k + l + m + n)-node element.

4 ,

 4

3 . 3 .

 33 The dimensions of the limited plate are: length L =10m, half-height H=2m and crack length a=5m; The material The top of the plate is subject to a tensile stress of 0 500 s = MPa. To avoid the reflection when the tensile stress wave reaches the crack-tip, the valid time for the simulation is t≤3tc=3H/cd=1.009×10 -3 s[START_REF] Bui | Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading[END_REF], where cd is the dilatational wave speed.

Fig. 3 .

 3 Fig. 3. Geometry of a semi-infinite mode-I crack

  Fig.6(c)) is the same as variable-node XFEM local refined mesh in the vicinity of crack.

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Comparison of the normalized mode-I DSIFs among different methods.

.

  is examined. The geometry of the plate are L =44mm, H=32mm, D=16mm, crack length a=22.63mm, and the incline angle of crack =45 a °. The material parameters are set to be modulus of elasticity E=29.4GPa, Poisson's ratio =0.286 n The left, top and bottom of plate supported in normal direction, as shown in Fig.7, the right end is subjected to a

Fig. 7 .

 7 Fig. 7. Geometry of a rectangular plate with an inclined crack.

Fig. 9 .

 9 Fig. 9. Comparison of the normalized mixed-mode DSIFs for a rectangular plate with a inclined edge crack.

4. 3 .

 3 Dynamic fracture of a center crack in a piezoelectric structureIn this dynamic study, we extend the V-XFEM to piezoelectric structure. Now, we consider a PZT-5H piezoelectric plate contains a central crack of length 2a, as depicted in Fig.10. The geometry of the plate is h=20mm and a = 2.4 mm. The material parameters of PZT-5H are set as: in the study, where l is electrical impact loading parameter. of the longitudinal wave along the second principal material axis. A comparison of the normalized dynamic intensity factors with the six crack tip enrichments[START_REF] Garcia-Sanchez | 2-D transient dynamic analysis of cracked piezoelectric solids by a time-domain BEM[END_REF] and the standard four crack tip enrichments is presented inFig 11. It can be seen that the results

Fig. 10 .

 10 Fig. 10. A rectangular piezoelectric plate with a central crack under impact loading.

Fig. 11 .

 11 Fig. 11. Comparison of the normalized dynamic intensity factors for different crack-tip enrichment functions

Fig. 12 .

 12 Fig. 12. Mesh of the rectangular piezoelectric plate with a central crack: (a) variable-node local refined mesh, (b)

Fig. 13 .IK

 13 Fig. 13. A comparison of the normalized dynamic intensity factors (a) normalized dynamic stress intensity factor,

Fig. 14 .

 14 Fig. 14. A comparison of the normalized dynamic intensity factors (a) normalized dynamic stress intensity factor,

Fig. 15 .Fig. 16 .

 1516 Fig. 15. A comparison of the normalized dynamic intensity factors (a) normalized dynamic stress intensity factor,

Fig. 17 .

 17 Fig. 17. Normalized dynamic intensity factor for different poling direction (a) normalized dynamic stress intensity

4. 4 .Fig. 18 .

 418 Fig. 18. Piezoelectric Plate with an edge crack.

  (a) The first crack increment step. (b) The second crack increment step. (c) The fourth crack increment step. (d) The fifth crack increment step.

Fig. 19 .

 19 Fig. 19. The crack growth profiles for the first, second, fourth and the fifth crack increment steps

Table 1 .

 1 Comparison of numerical efficiency between the variable-node XFEM and standard XFEM.

	Methods	Nodes	Elements	Dofs	CPU Time(s)
		1716	1625	3592	36.4263
	XFEM	3360	3239	6908	48.7061
		14896	14625	30212	179.4868
	Variable-node
		3196	3025	6812	72.1416
	XFEM

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 11802261, 11672099, 52075465), the Hunan Science Foundation for Distinguished Young Scholars (Grant No.2019JJ20015), and the science and technology innovation Program of Hunan Province (Grant No. 2020RC4038). The financial supports are gratefully acknowledged. CRediT authorship contribution statement Shuohui Yin: Writing -original draft, Conceptualization, Methodology. Ning Zhang: Programing (elastic problems), Numerical results. Peng Liu: Programing (piezoelectric problems), Numerical results. Jingang Liu: Check and improved the code, review & editing. Tiantang Yu: Conceptualization, numerical simulation. Shuitao Gu: Writing -review & editing, Conceptualization, theoretical analysis. Yu Cong: Writing -review & editing.


Mohamed Sami Cherif

Djamal Habet

André Abramé

Understanding the power of Max-SAT resolution through UP-resilience

Maximum Satisfiability (Max-SAT)

Max-SAT Problem

Input: a formula Φ in Conjunctive Normal Form (CNF) Output: the maximum (resp. minimum) number of satisfied (resp. falsified) clauses in Φ over all possible variable assignments Branch & Bound (BnB) for Max-SAT Binary search algorithm which maintains and constantly updates two values : Upper Bound (UB): value of the best known solution Lower Bound (LB): estimation of the best accessible solution Cut: if LB ≥ UB then backtrack

Lower Bound Estimation

At each node with current assignment I, we compute a new estimation of LB.

LB = F C(Φ| I ) + IS(Φ| I )

F C(Φ| I ): number of falsified clauses in Φ| I IS(Φ| I ): number of disjoint Inconsistent Subsets (IS) detected in Φ| I by Simulated Unit Propagation (SUP) ψ = {{x}, {x, x 1 }, {x, x 2 }, {x 1 , x 3 }, {x 2 , x 3 }} is an IS detected through SUP represented in the form of an implication graph :

x x 1 x 3 x 2 {x} {x, x 1 } {x 1 , x 3 } { x , x 2 } { x2 , x3 } {x 2 , x 3 }
To ensure that detected ISs are disjoint, they are temporarily deleted or transformed by Max-SAT resolution in which case they can be maintained in the current subtree.

Max-SAT resolution Transformation

{x 2 , x 3 } {x 1 , x 3 } {x 1 , x 2 } {x 1 , x 2 } {x, x 2 } {x, x 1 } {x, x 1 } {x, x 1 } {x} {x} {x 1 , x 2 , x 3 } {x 1 , x 2 , x 3 } {x, x 1 , x 2 } {x, x 1 , x 2 } x 3 x 2 x 1
x

UP-resilience Definition

The transformation of an IS ψ is UP-resilient for a literal l in ψ iff : The transformation of ψ is UP-resilient iff is UP-resilient for all the literals in ψ.

The IS ψ = {{x}, {x, x 1 }, {x, x 2 }, {x 1 , x 3 }, {x 2 , x 3 }} also corresponds to another implication graph:

x x 2 x 3 x 1 {x} {x, x 2 } {x 2 , x 3 } { x , x 1 } { x 1 , x 3 } {x 1 , x 3 } pneigh(x 2 ) = {{x, }, {x, x 3 }} Θ(ψ)| {x,x 3 } = {{x 1 , x 2 }, {x 1 , x 2 }, } → x 2 cannot be propagated in Θ(ψ)| {x,x 3 } (fragmentation phenomenon) → The transformation of ψ w.r.t RPO is not UP-resilient

Properties

If the transformation of an IS ψ is UP-resilient for a set of literals L in ψ then ∀N ∈ pneigh(L) : ∈ N or ∀l ∈ L, l can be propagated in Θ(ψ)| N \{l} .

The order of application (variable sequence) of Max-SAT resolution has a direct impact on the UP-resilience of the transformations.

Efficient learning schemes (Patterns) are UP-resilient.

→ quantifies the impact of transformations on the SUP mechanism → provides a theoretical understanding of the efficiency of learning schemes

UP-resilience and Main Patterns

Since Max-SAT resolution transformations can negatively affect the efficiency of BnB solvers, they are only performed when they correspond to certain patterns.

{l 1 , l 2 }, {l 1 , l 2 } (P 1 ) {l 1 } {l 1 , l 2 }, {l 1 , l 3 }, {l 2 , l 3 } (P 2 ) {l 1 }, {l 1 , l 2 , l 3 }, {l 1 , l 2 , l 3 } {l 1 }, {l 1 , l 2 }, {l 2 , l 3 }, ..., {l k , l k+1 }, {l k+1 }} (P 3 )
, {l 1 , l 2 }, {l 2 , l 3 }, ..., {l k , l k+1 } → Patterns do not augment the size of formula and produce unit clauses.

If a subset of an IS ψ matches the premises of pattern (P1) or (P2) or (P3). Then, the Max-SAT resolution transformation described in the pattern is UP-resilient w.r.t all possible application orders.

UP-resilience and Binary UCS Patterns

A binary k-Unit Clause Subset (k b -UCS) is a set of binary clauses {c 1 , ..., c k } in an IS such that there exists a sequence of Max-SAT resolution steps on c 1 , ..., c k that produces a unit clause resolvent.

Patterns (P 1 ) and (P 2 ) are respectively equivalent to 2 b -UCS and 3 b -UCS patterns. Many UCSs can be easily detected using the First Unit Implication Point (FUIP). RPO does not necessarily ensure the UP-resilience of k b -UCSs for k ≥ 4.

The subset ψ \ {x} of the IS ψ = {{x}, {x,

x 1 }, {x, x 2 }, {x 1 , x 3 }, {x 2 , x 3 }}
is a 4 b -UCS, detected through the FUIP x, for which RPO does not ensure UP-resilience.

In the implication graph of ISs with k b -UCS detected by the FUIP, there exits exactly two disjoint paths from the FUIP to .

Empirical Study on the Relevance of UP-resilience

IRS: learns transformations whose all intermediary resolvents contain less than four literals (MiniMaxSAT learning scheme). PAT: learns the transformations when the ISs match the three main patterns. PAT+: learns the transformations when the IS matches the three main patterns or certain k-UCS patterns (k ∈ {4, 5}). UPRes: learns only UP-resilient transformations of ISs.

∀N ∈ pneigh

  ψ (l) : ∈ N or l can be propagated in Θ(ψ)| N pneigh ψ (l): possible neighborhood of l, i.e. the set of its different neighborhoods in all the implication graphs enabling the detection of the IS ψ Θ(ψ)| N : result of the Max-SAT resolution transformation of the IS ψ

→

  k b -UCSs are UP-resilient with respect to PRO and PRCO → New approach to extend patterns by UP-resilience

  weighted (partial) instances in total tested with the solver ahmaxsat, S = number of solved instances, T = average solving time, D = average number of decision, % L = percentage of learned transformations, % UPR = percentage of learned UP-resilient transformations.

  Transformation of the IS ψ = {{x}, {x, x 1 }, {x, x 2 }, {x 1 , x 3 }, {x 2 , x 3 }} with respect to the Reverse Propagation Order (RPO) x 3 , x 2 , x 1 , x :

		Max-SAT
		resolution
	Compensation clauses
	cc 1 = {x, y 1 , ..., y s , z 1 }	cc t+1 = {x, z 1 , ..., z t , y 1 }
	cc 2 = {x, y 1 , ..., y s , z 1 , z 2 }	cc t+2 = {x, z 1 , ..., z t , y 1 , y 2 }
	...	...
	cc t = {x, y 1 , ..., y s , z 1 , ..., z t-1 , z t } cc t+s = {x, z 1 , ..., z t , y 1 , ..., y s-1 , y s }

c = {x, y 1 , ..., y s } c = {x, z 1 , ..., z t } cr = {y 1 , ...,

y s , z 1 , ..., z t }, cc 1 , ..., cc t , cc t+1 , ..., cc t+s