The Response of Conductive-Fiber Reinforced Composites to Electric Field

TL;DR

This paper develops an analytical model coupling electric, magnetic, thermal, and mechanical effects to predict the response of fiber-reinforced composites under electric fields, with applications in sensors and actuators.

Contribution

It introduces a coupled analytical approach for predicting the multi-physics response of conductive-fiber composites subjected to electric fields, including computational methods and practical applications.

Findings

Magnetic, thermal, and mechanical field distributions are quantitatively analyzed.

The model predicts the composite's macroscopic response under electric excitation.

Application to iron fiber-epoxy composites demonstrates the approach's effectiveness.

Abstract

An analytical procedure, which couples electric, magnetic, thermal and mechanical effects is presented for the prediction of the response of unidirectional fiber-reinforced composites that are subjected to electric field, applied in the fibers' direction. It is assumed that at least one of the phases of the composite (e.g., the fibers) is electrically conductive, and that all phases are thermally conductive. The composite is assumed to occupy a finite, symmetric domain which is discretized into a double array of subcells. The governing equations, with the interfacial and boundary conditions, are satisfied in the integral sense. The externally applied field generates electric current, which induces a magnetic field as well as temperature increase. The mechanical deformation of the composite results from the combined effect of the ponderomotive force, which is created by the magnetic field, and the temperature distribution within the constituents. The purpose of the present paper is three-fold. (a) To perform quantitative analysis of the model for the ponderoromotive force in deformable media. (b) To present computational strong-form treatment of the magneto-mechanical boundary-value problem in a composite. (c) To suggest a computational apparatus for deriving the response of a sensor/actuator excited by an applied electric field or electric field gradient. Application is given, presenting the magnetic, thermal and mechanical field distributions, as well as the macroscopic (global) response of a composite, which consists, for simplicity, of two iron fibers embedded in an epoxy matrix.

Equations

$$\nabla \cdot \mathbf{D}=0,\nabla \times \mathbf{E}=0$$

$$\nabla \cdot \mathbf{B}=0,\nabla \times \mathbf{H}=\mathbf{J}$$

$$\nabla \cdot \mathbf{q}=\mathbf{J}\cdot \mathbf{E}$$

$$\nabla \cdot \mathbf{\Sigma }+\mathbf{F}=0$$

$$\mathbf{B}={\mu }_0{\mu }_r\mathbf{H}$$

$$\mathbf{J}=\sigma \mathbf{E},\mathbf{D}={ϵ}_0{ϵ}_r\mathbf{E}$$

$$\mathbf{q}=-\kappa \nabla \theta$$

$$\mathbf{\Sigma }=\mathbf{C}\mathbf{ϵ}-\mathbf{\Gamma }\theta$$

$$\mathbf{F}=\mathbf{J}\times \mathbf{B}+{\mu }_0({\mu }_r-1)(\nabla {\mathbf{H}}^{\top })\mathbf{H}+(\nabla {\mathbf{E}}^{\top })\mathbf{P}$$

$$\mathbf{E}={\left\{E_0,0,0\right\}}^{\top }$$

$$\nabla \cdot {\mathbf{H}}^h$$

$$\nabla \cdot {\mathbf{H}}^p$$

$${\mathbf{H}}^{p(\beta \gamma )}=\frac{{\sigma }^{(\beta \gamma )}{E}_0}{2}{\left\{0,-y_3^{(\gamma )},y_2^{(\beta )}\right\}}^T$$

$${\varphi }^{h(\beta \gamma )}$$

$$\displaystyle-2\left[\begin{array}[]{cccc}a&0&b&0\\ 0&c&0&d\\ c&0&d&0\\ 0&a&0&b\\ \end{array}\right]^{(\beta\gamma)}\ \left[\begin{array}[]{c}\hat{h}_{1}\\ \hat{h}_{2}\\ \hat{h}_{3}\\ \hat{h}_{4}\end{array}\right]^{(\beta\gamma)}+\frac{\sigma^{(\beta\gamma)}E_{0}}{4}\left[\begin{array}[]{c}-h_{\beta}l_{\gamma}\\ h_{\beta}l_{\gamma}\\ 0\\ 0\\ \end{array}\right]$$

$${\bar{\mathbf{H}}}^{\pm (3)(\beta \gamma )}$$

$${\bar{\mathbf{H}}}^{\pm (2)(\beta \gamma )}$$

$$a^{(\beta \gamma )}=\sin \left(\sqrt{c_1^{(\beta \gamma )}}{h}_{\beta }/2\right)\cosh \left(\sqrt{c_1^{(\beta \gamma )}}{l}_{\gamma }/2\right),b^{(\beta \gamma )}=\sinh \left(\sqrt{c_2^{(\beta \gamma )}}{h}_{\beta }/2\right)\cos \left(\sqrt{c_2^{(\beta \gamma )}}{l}_{\gamma }/2\right)$$

$$c^{(\beta \gamma )}=\cos \left(\sqrt{c_1^{(\beta \gamma )}}{h}_{\beta }/2\right)\sinh \left(\sqrt{c_1^{(\beta \gamma )}}{l}_{\gamma }/2\right),d^{(\beta \gamma )}=\cosh \left(\sqrt{c_2^{(\beta \gamma )}}{h}_{\beta }/2\right)\sin \left(\sqrt{c_2^{(\beta \gamma )}}{l}_{\gamma }/2\right)$$

$$\displaystyle\left[\begin{array}[]{c}\hat{h}_{1}\\ \hat{h}_{2}\\ \hat{h}_{3}\\ \hat{h}_{4}\end{array}\right]^{(\beta\gamma)}=\frac{1}{2(bc-ad)^{(\beta\gamma)}}\left[\begin{array}[]{cccc}d&0&-b&0\\ 0&-b&0&d\\ -c&0&a&0\\ 0&a&0&-c\\ \end{array}\right]^{(\beta\gamma)}\ \left[\begin{array}[]{c}{\bar{H}}^{+(3)}_{2}+\frac{\sigma E_{0}}{4}h_{\beta}l_{\gamma}\\ {\bar{H}}^{+(2)}_{3}-\frac{\sigma E_{0}}{4}h_{\beta}l_{\gamma}\\ {\bar{H}}^{+(2)}_{2}\\ {\bar{H}}^{+(3)}_{3}\end{array}\right]^{(\beta\gamma)}$$

$$\displaystyle\left[\begin{array}[]{c}{\bar{H}}^{-(3)}_{2}\\ {\bar{H}}^{-(2)}_{3}\\ {\bar{H}}^{-(2)}_{2}\\ {\bar{H}}^{-(3)}_{3}\end{array}\right]^{(\beta\gamma)}=\left[\begin{array}[]{c}{\bar{H}}^{+(3)}_{2}+\frac{\sigma E_{0}}{2}h_{\beta}l_{\gamma}\\ {\bar{H}}^{+(2)}_{3}-\frac{\sigma E_{0}}{2}h_{\beta}l_{\gamma}\\ {\bar{H}}^{+(2)}_{2}\\ {\bar{H}}^{+(3)}_{3}\end{array}\right]^{(\beta\gamma)}$$

$$\displaystyle\left[\begin{array}[]{c}{\bar{H}}^{+(3)(\beta,\gamma+1)}_{2}\\ {\bar{H}}^{+(2)(\beta+1,\gamma)}_{3}\\ \mu^{(\beta+1,\gamma)}_{r}{\bar{H}}^{+(2)(\beta+1,\gamma)}_{2}\\ \mu^{(\beta,\gamma+1)}_{r}{\bar{H}}^{+(3)(\beta,\gamma+1)}_{3}\end{array}\right]=\left[\begin{array}[]{c}{\bar{H}}^{+(3)(\beta\gamma)}_{2}\\ {\bar{H}}^{+(2)(\beta\gamma)}_{3}\\ \mu^{(\beta\gamma)}_{r}{\bar{H}}^{+(2)(\beta\gamma)}_{2}\\ \mu^{(\beta\gamma)}_{r}{\bar{H}}^{+(3)(\beta\gamma)}_{3}\end{array}\right]+\left[\begin{array}[]{c}-\frac{E_{0}}{2}h_{\beta}l_{\gamma+1}\sigma^{(\beta,\gamma+1)}\\ \frac{E_{0}}{2}h_{\beta+1}l_{\gamma}\sigma^{(\beta+1,\gamma)}\\ 0\\ 0\end{array}\right]$$

$$\displaystyle\left[\begin{array}[]{c}{\bar{H}}^{+(3)(\beta\gamma)}_{2}\\ {\bar{H}}^{+(2)(\beta\gamma)}_{3}\\ {\bar{H}}^{+(2)(\beta\gamma)}_{2}\\ {\bar{H}}^{+(3)(\beta\gamma)}_{3}\end{array}\right]=\left[\begin{array}[]{c}{\hat{H}}^{(\beta)}_{23}-\frac{E_{0}}{2}\sum^{\gamma}_{\gamma^{\prime}=1}h_{\beta}l_{\gamma^{\prime}}\sigma^{(\beta\gamma^{\prime})}\\ {\hat{H}}^{(\gamma)}_{32}+\frac{E_{0}}{2}\sum^{\beta}_{\beta^{\prime}=1}h_{\beta^{\prime}}l_{\gamma}\sigma^{(\beta^{\prime}\gamma)}\\ \frac{1}{\mu^{(\beta\gamma)}_{r}}{\tilde{H}}^{-(2)(1,\gamma)}_{2}\\ \frac{1}{\mu^{(\beta\gamma)}_{r}}{\tilde{H}}^{-(3)(\beta,1)}_{3}\end{array}\right]=\left[\begin{array}[]{c}{\hat{H}}^{(\beta)}_{23}-\frac{E_{0}}{2}\sum^{\gamma}_{\gamma^{\prime}=1}h_{\beta}l_{\gamma^{\prime}}\sigma^{(\beta\gamma^{\prime})}\\ {\hat{H}}^{(\gamma)}_{32}+\frac{E_{0}}{2}\sum^{\beta}_{\beta^{\prime}=1}h_{\beta^{\prime}}l_{\gamma}\sigma^{(\beta^{\prime}\gamma)}\\ \frac{1}{\mu^{(\beta\gamma)}_{r}}{\tilde{H}}^{+(2)(N_{\beta},\gamma)}_{2}\\ \frac{1}{\mu^{(\beta\gamma)}_{r}}{\tilde{H}}^{+(3)(\beta,N_{\gamma})}_{3}\end{array}\right]$$

$${\hat{H}}_{32}^{(\gamma )}=-\frac{E_0}{4}{l}_{\gamma }\sum _{\beta =1}^{N_{\beta }}{h}_{\beta }{\sigma }^{(\beta \gamma )},{\hat{H}}_{23}^{(\beta )}=\frac{E_0}{4}{h}_{\beta }\sum _{\gamma =1}^{N_{\gamma }}{l}_{\gamma }{\sigma }^{(\beta \gamma )}$$

$$\displaystyle\left[\begin{array}[]{c}{\bar{H}}^{+(3)}_{2}\\ {\bar{H}}^{+(2)}_{3}\\ {\bar{H}}^{+(2)}_{2}\\ {\bar{H}}^{+(3)}_{3}\end{array}\right]^{(\beta\gamma)}=\left[\begin{array}[]{c}\frac{E_{0}}{4}\left(\sum^{N_{\gamma}}_{\gamma=1}h_{\beta}l_{\gamma}\sigma^{(\beta\gamma)}-2\sum^{\gamma}_{\gamma^{\prime}=1}h_{\beta}l_{\gamma^{\prime}}\sigma^{(\beta\gamma^{\prime})}\right)\\ \frac{E_{0}}{4}\left(2\sum^{\beta}_{\beta^{\prime}=1}h_{\beta^{\prime}}l_{\gamma}\sigma^{(\beta^{\prime}\gamma)}-\sum^{N_{\beta}}_{\beta=1}h_{\beta}l_{\gamma}\sigma^{(\beta\gamma)}\right)\\ \frac{1}{\mu^{(\beta\gamma)}_{r}}{\tilde{H}}^{-(2)(1,\gamma)}_{2}\\ \frac{1}{\mu^{(\beta\gamma)}_{r}}{\tilde{H}}^{-(3)(\beta,1)}_{3}\end{array}\right]$$

$${\hat{h}}_1^{(\beta \gamma )}$$

$${\hat{h}}_2^{(\beta \gamma )}$$

$${\hat{h}}_3^{(\beta \gamma )}$$

$${\hat{h}}_4^{(\beta \gamma )}$$

$${\tilde{H}}_2^{-(2)(1,\gamma )}$$

$${\tilde{H}}_3^{-(3)(\beta ,1)}$$

$${\bar{\xi }}_2=\sum _{n=1}^{{\beta }^{\prime }-1}{h}_n-\sum _{n=1}^{\beta -1}{h}_n,{\bar{\xi }}_3=\sum _{n=1}^{{\gamma }^{\prime }-1}{l}_n-\sum _{n=1}^{\gamma -1}{l}_n$$

$$I_1(a,x_1,x_2)$$

Full text

Key words : Fiber-reinforced composite, Magnetic field, Ponderomotive force, Electric conduction, Heat conduction, Mechanical deformation.

1 Introduction

The micromechanical response of composite materials is usually investigated when the application of mechanical or thermo-mechanical loading is considered. Additional macroscopic physical effects, such as electric and magnetic field distribution, are commonly treated for materials with constitutive coupling, such as in the case of piezoelectric and electrostrictive or magnetostrictive materials. Fewer work deals with the response to electric loading of electrically conductive composites, as, for example, in the case of a composite in which one of the phases is electrically conductive (e.g., carbon fibers). In such a case, as a result of the application of electric field, electric current flows in the conductive phase, inducing magnetic field distribution, as well as temperature increase. The magnetic field gives rise to the so-called ponderomotive, or generalized Lorentz force, which, along with the corresponding temperature distribution, generates mechanical deformation.

The main challenge in addressing the aforementioned type of problems is the proper description of Maxwell’s equations and the associated Maxwell tensor, or alternatively, the ponderomotive force and moment couple density, for the case of deforming continua. Fundamental investigations can be found in the work of Tiersten (1971), Pao (1978) and in the recent reviews of Hutter et al. (2006) and Santapuri et al. (2013).

Beyond the establishing of the constitutive and governing equations of the aforementioned ’multi-physics’ problem, there arise mathematical issues pertinent to the specifics of the field form of Maxwell’s equations. In the case of slow electromagnetic dynamics, the main issue is related to the fact that the electric and magnetic fields each have to satisfy multiple differential equations. In the context of numerical solution, this difficulty is usually treated by introducing auxiliary gauges, e.g. the Coulomb gauge, see Zienkiewicz et al. (1977), Simkin and Trowbridge (1979) and Barton and Cendes (1987), for example. Another approach is taken in Kudryavtsev and Trashkeev (2013), where a hyperbolic system of evolution equations is derived with no auxiliary differential equations to impose.

Computationally, the coupled multiphysics response of composite materials for general three-dimensional domains is normally derived using the finite element formalism, see Sridhar et al. (2015) for a recent example. In the case of plates, the coupled problem is solved by integrating the plate equations using finite difference discretization, in conjunction with the Love-Kirchhoff, as in Zhupanska and Sierakowski (2007), Zhupanska and Sierakowski (2011), Barakati and Zhupanska (2014) or the von Kármán approximation, as in Librescu et al. (2003).

In Librescu et al. (2003), the dynamic response of a plate reinforced by electrically conducting fibers is considered. As a result of the nonlinear dependence of the mechanics on the magnetic field, linearization is performed, based on the assumption of small disturbances. In Barakati and Zhupanska (2014) on the other hand, a similar setting is assumed, also for plate dynamics, however without the small disturbances assumption. A different model for the ponderomotive force is used, and much as in Librescu et al. (2003), uniform material, rather than a composite, is considered. The solution is obtained by a finite difference method.

The purpose of the present article is three-fold. First, to perform quantitative investigation of the effect of the correction for the ponderomotive force, relatively to the standard expression of the Lorentz force, which is established for the particle-in-vacuum case. To this end we follow Santapuri et al. (2013) in choosing the Maxwell-Minkowski model, reduced for the linear quasistatic regime. The framework for the examination of the aforementioned effect is that of unidirectional deformation of a two-constituent composite specimen. The second perspective in which the present work is performed is computational. To this end, the Higher Order Theory (HOT), as in Aboudi et al. (2013), which has been employed for the analysis of functionally graded materials, is generalized to account for electric-current-induced mechanical deformations. In the framework of the HOT, the governing equations together with the interfacial and boundary conditions are imposed in the average (integral) sense. It should be noted that since the HOT is a strong-form approach to the solution of boundary-value problems, the present work opts to construct a strong-form computational apparatus for the analysis of the multi-physics response of composites in the linear quasistatic regime. In the third perspective, the present paper in effect proposes a sensor/actuator with the accompanying computational apparatus, which will have the effect of nonlinear response to either electric field or electric field gradient, in line with the view presented in Santapuri et al. (2015). The details of these perspectives are further discussed in the concluding section.

To elaborate on the third perspective, the present article offers a coupled electro-magneto-thermo-mechanical analysis for the prediction of the mechanical deformation of partially conductive composites, caused by the application of electric field. The engineering motivation for this analysis can be illustrated by the following discussion. Consider a setting (e.g. a large capacitor), in which a uniform electric field exists. It is then assumed that a local effect (perturbation – a small relative change in the field but with a high local gradient, as created by a small conducting object) has been introduced, as a result of which the electric field loses its uniformity. The present analysis shows that a device (sensor) in the form of a ring that consists of, say, two conductive wires embedded in a polymeric (insulation) matrix, will reveal the effect by showing the induced mechanical deformation. This deformation, caused jointly by the Joule heating and the ponderomotive force, can be calibrated to estimate the intensity of the perturbation-induced local electric field gradient. Such a (partially conducting) ring would be advantageous over, say, a piezoelectric sensor in that it would sense the finite electric field gradient rather than the nearly uniform local field itself. High accuracy of the device-calibration would be obtained due to the offered analysis, owing to it being based on strong-form solution of the equations (which is possible for linear constitutive relations). It should be noted that a sensor in the form of a ring is a three-dimensional object, and a 3D imperfectly conducting ring in a non-uniform electric field develops a current, see Assis et al. (1999). It is also shown in that article that the expressions for the 3D corrections to the electric field in the ring become negligible for a slender wire. Consequently, in the limiting case of a slender ring-wire the applicability of the present 2D analysis is justified.

In the present article, the analysis of mechanical deformation generated by the application of externally applied electric field on electrically and thermally conducting fiber-reinforced composites is considered. The composite material is assumed to consist of continuous fibers that are distributed in the matrix, forming a symmetric array. The electric field is applied in the fibers’ direction, and the boundaries of the composite are assumed to be traction-free and thermally convective. The composite region is divided into a double array of subcells. The magnetic field induced by the applied electric field is derived, and the resulting ponderomotive force is established. Next, the equations that govern the thermal problem related to the heat generated by the electric current are established and solved numerically, to obtain the temperature field distribution and its gradient. Finally, the mechanical deformations of the composite phases are obtained from the solution of a system of sparse algebraic equations, following the HOT approach. The presented derivation is illustrated for the case in which uniform (throughout the composite domain) electric field is applied to a thermally conductive polymeric matrix (epoxy) reinforced by two iron fibers, which are both electrically and thermally conductive. The solution demonstrates the distributions of the magnetic, thermal and mechanical fields created by the interaction between the two fibers. It should be emphasized, however, that the present theory can be generalized to establish the response of any number of fibers of identical or different properties and cross-section. Such generalization would affect only the computational complexity caused by the increase in the number of subcells.

The structure of this paper is as follows. In Section 2, the governing equations are presented. In Sections 3, 4, 5 and 6 the electric, magnetic, thermal and mechanical fields, respectively, are considered, and the procedures for their computation are derived. In Section 7, the application of the derived computational method is demonstrated. The obtained numerical results are discussed in Section 8. The conclusions section, Section 9, provides overall discussion of the derived theory.

2 Governing Equations

Consider a unidirectional fiber-reinforced composite, in which the fibers are as shown in Fig. 11(a). The composite occupies the region $-H/2\leq x_{2}\leq H/2$, $-L/2\leq x_{3}\leq L/2$, and it is assumed that the domain is symmetrically occupied with respect to the axes. The fibers are assumed to be electrically conductive. A uniform electric field, namely, ${\bm{E}_{0}}=\left\{E_{0},0,0\right\}^{\top}$ is acting at the boundary in the out-of-plane x1-direction (fibers’ direction), formally ’at $\pm\infty$’. As a result of the existence of the conductive phase, electric current is generated and magnetic field is induced. The corresponding temperature increase and ponderomotive force produce mechanical deformation of the composite.

The governing equations are given by Maxwell’s equations, along with a model for the ponderomotive force, energy balance and mechanical equilibrium (momentum balance), as follows. Starting with Faraday’s law for the static regime and Gauss’s law for the case of uncharged material:

[TABLE]

where $\bm{D}$ is electric displacement vector and ${\bm{E}}$ is the emergent electric field intensity in the composite, and following by magnetic flux conservation and Ampère’s law:

[TABLE]

where ${\bm{B}}$, ${\bm{H}}$ and ${\bm{J}}$ are the magnetic flux density, magnetic field intensity and electric current density, respectively. The energy balance in the linear elastic quasistatic limit reads:

[TABLE]

where ${\bm{q}}$ is the heat flux. Next, the mechanical equilibrium equation reads:

[TABLE]

where ${\bm{\Sigma}}$ and ${\bm{F}}$ are the Cauchy stress tensor and body force-per-unit-volume vector, respectively.

The corresponding constitutive equations are as follows. The magnetization law is assumed to be linear, even though in the application section a ferromagnetic material is chosen in order to enhance the effect of the ponderomotive force on the deformation. The resolution of this seeming contradiction is simply the assumption that the magnetic field intensity remains small, corresponding to the assumption of an applied weak electric field excitation at the boundary. In addition, quasistatic increase of the electric field at the boundary is assumed, which corresponds to linear response with no appearance of magnetic hysteresis:

[TABLE]

where $\mu_{0}$ is the free space permeability and $\mu_{r}$ is the constant relative material permeability. Ohm’s law in the vanishing material velocity limit and the electric polarization equation in the weak field approximation in terms of the electric displacement vector reads:

[TABLE]

where $\sigma$ is the constant coefficient of electric conductivity and $\epsilon_{r}$ is the relative constant electric permittivity coefficient of the material, and where $\epsilon_{0}$ is vacuum permittivity. Fourier’s law of heat conduction for thermally isotropic materials in the quasistatic limit and small temperature gradient/change assumption:

[TABLE]

where $\kappa$ is the constant scalar heat conductivity of the solid and $\theta$ is the temperature field measured as deviation from a uniform reference temperature. Finally, the Cauchy stress tensor ${\bm{\Sigma}}$ is related to the small strain tensor ${\bm{\epsilon}}$ and the temperature deviation $\theta$ as follows:

[TABLE]

where ${\bm{C}}$ and ${\bm{\Gamma}}$ are the constant material stiffness fourth-order and thermal stress second-order tensors, respectively. The body force-per-unit-volume vector in Eq. (4) is the ponderomotive force, representing the macroscopic effect of long-range magnetic dipole interactions. Here the Maxwell-Minkowski model is taken, following Santapuri et al. (2013), which in the linear quasistatic regime takes the following form:

[TABLE]

where $\bm{P}$ is electric polarization. The distributed couple $\bm{L}$ appearing in the Maxwell-Minkowski model in the dynamic case vanishes for the quasistatic regime assumed herein, rendering the Cauchy stress tensor symmetric (see Santapuri (2016)).

In the framework of the present two-dimensional analysis, the external electric field is assumed to be tangential to the global boundaries. Consequently, the electric boundary condition is that of tangential field continuity. The magnetic boundary conditions are those of tangential field intensity and normal flux density continuity at the external boundaries. The mechanical boundary conditions are given by zero tractions imposed at $x_{2}=\pm H/2$ and $x_{3}=\pm L/2$, together with a suitable displacement pinning to prevent singularity. The thermal boundary conditions are assumed to be convective. The interfacial conditions require the continuity of the tractions and normal components of magnetic flux density and heat flux. In addition, the temperature, displacements and tangential electric and magnetic field intensities must be continuous at the interfaces between the phases. According to the present method of solution, the composite region $-H/2\leq x_{2}\leq H/2$, $-L/2\leq x_{3}\leq L/2$ is divided into $N_{\beta}$ and $N_{\gamma}$ subcells in the x2 and x3 directions, respectively. Each subcell is labeled by the indices $(\beta,\gamma)$ with $\beta=1,...,N_{\beta}$ and $\gamma=1,...,N_{\gamma}$. The dimensions of subcell $(\beta\gamma)$ are denoted by $h_{\beta}$ and $l_{\gamma}$, and local coordinates $(y^{(\beta)}_{2},y^{(\gamma)}_{3})$ are introduced within each subcell whose origin is located at its center, see Fig. 1(b).

It should be mentioned that the composite materials which are dealt with in the present article are not assumed to possess periodic microstructure. Thus a homogenization technique and its associated periodic boundary conditions are not employed in the proposed analyses.

3 The Electric Field

The electric problem consists of the combination of Eqs. (1) and (6b), and suitable boundary conditions. The approach followed in the present paper is to divide the spatial domain occupied by the composite to subcells containing each only one material.

In every such subcell, the material properties are uniform and there is no discontinuity in the fields.

Consequently, the differential form of Maxwell’s equations is valid.

The corresponding boundary conditions are then continuity of the fields along the interfaces between the subcells. In the present electric problem, a possible solution of the aforementioned equations is a field $\bm{E}$ uniform across a given subcell.

Since the global boundary conditions require the electric field to point in the $x_{1}$ direction only, and since the interfaces of a subcell in the $x_{2}-x_{3}$ plane are oriented perpendicular to $x_{1}$, clearly if one assumes the electric field for every value of $x_{1}$ in the composite to point in the $x_{1}$ direction, just as the applied field, then the continuity conditions of the tangential components of the electric field are satisfied.

It is known that tangential continuity in electric problems should be enforced on the electric field. For materials with different relative permittivities, electric displacement continuity along the interface may result in an electric field that is non-uniform across material boundaries.

However, since a uniform electric field pointing in the $x_{1}$ direction satisfies the differential equations in both the fibers and the matrix, and the only interface continuity condition is that this electric field should be continuous across the boundaries as well, evidently the emergent electric field in the composite should be just as if the composite were uniform monolithic material: a uniform field in the $x_{1}$ direction, insensitive to material boundaries, namely:

[TABLE]

much like the global boundary condition itself. Then, clearly, if Eq. (10) is a solution, then due to linearity it is also the unique solution.

As for the interpretation of the boundary conditions, two possibilities exist. (a) To understand the problem as two-dimensional, straight and extending between minus to plus infinity, with prescribed values there. (b) Alternatively, to view the system as a ring with a large radius, relatively to the cross-sectional dimensions. In this scenario, an external field would be tangential to the ring and the material field would be equal to the external field, due to tangential component electric field continuity.

From the engineering perspective this would actually mean that the external field would have a large-scale gradient, such that a net current would be created in the composite fiber-reinforced ring.

The external electric field gradient would have to be large enough to create finite current in the ring, yet small enough to assume approximately uniform field distribution within the cross-section of the ring. It is in this perspective that one may suggest an electric field gradient sensing device, with the accompanying computational apparatus, which would be based on the derivation of the magneto-thermo-elastic response, as obtained in the following.

4 The Magnetic Field

The magnetic field intensity ${\bm{H}}$ is governed by Eq. (2) in conjunction with Eq. (5). Let ${\bm{H}}$ be decomposed into homogeneous ${\bm{H}}^{h}$ and particular ${\bm{H}}^{p}$ parts. These parts satisfy:

[TABLE]

It follows that ${\bm{H}}^{h}$ can be derived from a harmonic magnetic potential ${\phi}^{h}$ such that ${\bm{H}}^{h}=-\nabla\phi^{h}$. It can be verified that ${\bm{H}}^{p}$ within subcell $(\beta\gamma)$ can be uniquely satisfied by:

[TABLE]

This is the Ampère solution inside an electric current-bearing wire. The harmonic magnetic potential can be generated as a four-terms expansion, consistent of harmonic eigenfunctions of the Laplace equation for $\phi^{h}$:

[TABLE]

The reason for which the number of free coefficients in the expansion in (13) is four is as follows. Since the aim of the present work is to enhance the micromechanical HOT method to include ponderomotive force effects, and the continuity at the faces of a subcell is enforced in the average sense, clearly, for a 2D subcell, 4 independent coefficients are needed to enforce the 4 interface conditions on the faces/edges. Since the magnetic equations are first order in $\bm{H}$, only one condition is needed at one edge of the subcell for each of the two vector components of $\bm{H}$, and that is for every one of the two spatial directions. This gives a total of 4 conditions, which due to linearity correspond to 4 necessary free coefficients in a feasible expansion for the harmonic potential. In order for the coefficients to be extractable from the interface conditions, the 4 eigenfunctions have to be orthogonal. This can only be achieved if 2D Cartesian harmonic functions are chosen – hence the function choice in Eq. (13).

The symmetry between the two spatial directions suggests that the trigonometric/hyperbolic function choice should be symmetric between the two directions. Moreover, since two trigonometric functions are needed, symmetry and the fact that the interface conditions are integral, leads to the fact that the entire expansion be represented by the first function, and the other three be complementary through odd/even and trigonometric/hyperbolic flips. Hence, there are only two possibilities for the first eigenfunction, namely $\sin\sinh$ and $\sin\cosh$. It so happens that the first choice leads to a coupled scheme, which rigorous scheme stability analysis as well as numeric examination show to be algorithmically unstable. The second choice, on the other hand, produces a formally stable scheme, which also decouples, enabling the analytical solution of the entire magnetic problem, as shown next.

In a discretized scheme the integral continuity conditions are the physically correct ones, since only integral continuity conditions allow the fulfillment of the Ampère and magnetic Gauss laws over arbitrary loops along subcell lines within the domain. Consequently, the proposed scheme is the only possible first order consistent and stable discrete integration scheme for the magnetic equation with a given transverse current in a composite.

The subcell interface continuity equations corresponding to the potential given in Eq. (13) are obtained by integrating ${\bm{H}}^{(\beta\gamma)}=-\nabla\phi^{h(\beta\gamma)}+{\bm{H}}^{p(\beta\gamma)}$, in conjunction with Eq. (13), over the subcell edges (the integrated values needed for interface continuity resulting in satisfaction of the integral Ampère and Gauss laws for a loop in the discretized specimen). Consequently, the coefficients $\hat{h}^{(\beta\gamma)}_{1},...,\hat{h}^{(\beta\gamma)}_{4}$ can be established as follows:

[TABLE]

where the superscript and subscript on $\bar{H}$ denote the normal to the surface and component, respectively, i.e.,

[TABLE]

and where

[TABLE]

The solution for the coefficients is given by:

[TABLE]

It can be easily verified that the regularity condition requires that: $\pi/4<\sqrt{c^{(\beta\gamma)}_{1}}h_{\beta}/2<\pi/2$ and $\pi/4<\sqrt{c^{(\beta\gamma)}_{2}}l_{\gamma}/2<\pi/2$.

This can be satisfied by the choice: $c^{(\beta\gamma)}_{1}=9\pi^{2}/[16(h^{(\beta)})^{2}],c^{(\beta\gamma)}_{2}=9\pi^{2}/[16(l^{(\gamma)})^{2}]$.

The relations between the integrated opposite subcell edges provide:

[TABLE]

Face-integrated continuity of interfaces, along with Eq. (54) yields:

[TABLE]

Solving the recursion equation (67) analytically gives:

[TABLE]

The last relation in (80) implies that components 3 and 4 can be also taken at the opposite global boundary, since the magnetic flux densities should be uniform in the direction in which they point. Thus symmetric global boundary conditions are assumed for the normal magnetic flux density components, i.e., a symmetric specimen and hence skew-symmetric tangential components’ global boundary conditions.

The assumption is reasonable, since normal magnetic flux components penetrate the composite remaining unaltered passing through different constituents, due to the fact that normal magnetic flux continuity across interfaces is permeability-independent (unlike for tangential magnetic flux components).

Next, in order to obtain the unknown single-index functions in the first two components in the right-hand-side of Eq. (80), we write them explicitly for the two opposite global boundaries (for each of the first two components there), and then solve for the unknowns, using skew-symmetry of the tangential field for a symmetric-cross-section conducting wire.

This yields:

[TABLE]

Substituting these relations in Eq. (80), we get the face-integrated field components explicitly, depending only on the global boundary conditions:

[TABLE]

Hence, the final expressions of the coefficients $\hat{h}^{(\beta\gamma)}_{1},...,\hat{h}^{(\beta\gamma)}_{4}$ take the form:

[TABLE]

where the normal components of the global boundary conditions for magnetic field intensity are obtained by integrating the Biot-Savart law both over each conducting subcell and each global boundary subcell edge:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Just like in Librescu et al. (2003), the global boundary conditions are presently obtained by summing up the contributions of all the current-bearing fibers, along with the Biot-Savart law. The difference is that here, two additional integrations were performed, one over each conducting fiber area, to account for the integral contribution without approximation, and the other is integration along a global boundary subcell edge, which is needed as the boundary conditions are enforced in the integral sense, for reasons explained earlier.

Finally, the magnetic field intensity distribution components in subcell $(\beta\gamma)$ are given by:

[TABLE]

[TABLE]

With the established electric and magnetic field intensities, $\bm{E}$ and ${\bm{H}}^{(\beta\gamma)}$, in each subcell, the ponderomotive force variation in the subcell, as given in Eq. (9), can be expressed as follows:

[TABLE]

5 The Thermal Field

In the following, an approximate solution for the temperature distribution in the composite is established using a hybrid HOT – finite difference approach. From Eqs. (3), (6a), (7) and (10), the temperature field within each subcell should satisfy Poisson’s equation:

[TABLE]

Let the temperature $\theta^{(\beta\gamma)}$ within subcell $(\beta\gamma)$ be decomposed into the homogeneous and particular parts:

[TABLE]

such that

[TABLE]

A solution to Eq. (101b) is:

[TABLE]

The homogeneous solution can be chosen in the form:

[TABLE]

This expression is constructed in a manner similar to the one used in the magnetic case – four linearly independent harmonic eigenfunctions, corresponding to to the four local boundary conditions for each subcell, two on each edge, as expected for a Poisson equation. The only difference is that this time there are two ’even-even’ functions, in addition to two ’odd-odd’ functions. This is due to the fact that in this case one condition is for heat flux continuity, which is an integral condition, but the other condition is for temperature continuity, which has no integral version and is a macroscopically local requirement. Moreover, no differentiation is required, just value setting.

These two features lead to the choice of two ’even-even’ eigenfunctions, as a means to correctly represent bidirectional thermal self-coupling in the resulting integration scheme. Now, evaluating the temperature at the centers of the edges $y^{(\beta)}_{2}=h_{\beta}/2,y^{(\gamma)}_{3}=0$ and $y^{(\beta)}_{2}=0,y^{(\gamma)}_{3}=l_{\gamma}/2$, and, in addition, integrating the normal heat flux densities along these edges, yields, in conjunction with Eq. (100), the following relation:

[TABLE]

with these auxiliary definitions employed:

[TABLE]

A regularity condition, corresponding to a non-vanishing determinant of the matrix in Eq. (120), is satisfied by the symmetric choice: $\sqrt{\lambda^{(\beta\gamma)}_{1}}h_{\beta}/2=\sqrt{\lambda^{(\beta\gamma)}_{2}}l_{\gamma}/2=\pi/4$, for which the determinant, $\Delta$, of the matrix in Eq. (5), satisfies the inequality: $\Delta\leq-3.7205$, guaranteeing regularity.

Next, evaluating the temperature at the centers of the opposite pair of adjacent edges, $y^{(\beta)}_{2}=-h_{\beta}/2,y^{(\gamma)}_{3}=0$ and $y^{(\beta)}_{2}=0,y^{(\gamma)}_{3}=-l_{\gamma}/2$, and, in addition, integrating the normal heat flux densities along these edges, yields, in conjunction with Eq. (103), the following, second set of relations:

[TABLE]

Consequently, the temperature and normal heat flux ’jumps’ between the opposite edges of subcell $(\beta\gamma)$, can be expressed as follows:

[TABLE]

with the following auxiliary matrix used:

[TABLE]

where it can be shown that: det$\left({\bm{M}}^{(\beta\gamma)}\right)$ = 1.

Rearranging Eq. (147), the final form of a set of $4N_{\beta}N_{\gamma}$ subcell equations relating opposite edges’ values for a given subcell is obtained:

[TABLE]

In addition to the above ’subcell’ relations, interfacial conditions between neighboring subcells must be satisfied, as follows:

[TABLE]

Finally, the boundary conditions at the external boundaries must be imposed. We follow here Librescu et al. (2003) by assuming convection boundary conditions, which take the form:

[TABLE]

where $\Lambda^{(\gamma)}_{\beta}$ and $\Lambda^{(\beta)}_{\gamma}$ are the convection coefficients prescribed at the boundaries $X_{2}=\pm H/2$ and $X_{3}=\pm L/2$, respectively. In general, the convection coefficient is defined by $\Lambda=\kappa_{0}\text{Nu}/L_{max}$ where Nu, $\kappa_{0}$ and $L_{max}$ are the Nusselt number, thermal conductivity of ambient air ($\kappa_{0}=0.024Wm^{-1}K^{-1}$) and the length of the largest side of the specimen, respectively.

In every subcell there are 8 unknowns given by the 4 temperatures at the centers of the edges and the 4 integrals along the subcell edges of the normal heat flux densities. Thus, there are $8N_{\beta}N_{\gamma}$ unknown variables in the considered composite. In the same time, Eq. (173) provides $4N_{\beta}N_{\gamma}$ relations. In addition, the interfacial conditions (5) form $2(N_{\beta}-1)N_{\gamma}+2N_{\beta}(N_{\gamma}-1)$ relations. Finally, the external boundary conditions, (5), provide another $2(N_{\beta}+N_{\gamma})$ relations. Consequently, there are $8N_{\beta}N_{\gamma}$ relations, which can be solved for the mid-edges’ temperatures and total normal heat fluxes. These enable the determination of the coefficients $\hat{\theta}_{1},...,\hat{\theta}_{4}$, by employing Eq. (138) or (156). Last, the temperature change $\theta^{(\beta\gamma)}$ within any subcell $(\beta\gamma)$, as produced by the application of an external electric field, $E_{0}$, can be determined from Eqs. (102) and (103). The derived numerical integration scheme appears to be stable (in terms of the recursion mapping encompassed in the outlined overall matrix equation comprised of Eqs. (173-5)).

6 The Mechanical Field

The method of solution of the mechanical problem is also similar to the HOT approach that has been presented in chapter 11 of Aboudi et al. (2013), which deals with functionally graded materials. Presently, however, the incorporation of the ponderomotive force and the temperature field need to be dealt with. Thus, the solution of the mechanical problem is based on satisfying the governing equations, continuity and boundary conditions in the average (integral) sense. This implies that the equilibrium equation is solved in strong form (i.e., not by numerical variational minimization) yet in the volume-integral sense. Furthermore, the four three-component force vectors acting on the edges of a rectangular two-dimensional subcell, are balanced by the three-component body force vector acting on the volume of the subcell. This results in three equations in each subcell. In addition, the continuity conditions are enforced in the integral sense, that is forces (normalized by edge areas) acting on the edges of adjacent subcells are set equal, which gives six more equations for each subcell. Finally, the displacements continuity is enforced. In principle, displacements could be made continuous at the centers of the edges, that is locally, just as was done in the case of the temperatures in the thermal problem. However, unlike the temperature, which is a purely local, intensive state variable in the material, the displacement is in fact a spatially integrable quantity. Consequently, displacement continuity is enforced here for the center-of-edge displacements, that is for the integrals of the displacement field along the edge, normalized by the edge length, rather than for the local displacement of the centers of the edges. This gives another six equations, leading to a total of fifteen for a subcell. This calls for fifteen subcell variables used to interpolate the displacement vector field within the subcell, as a quadratic polynomial in two variables (suitable for the plane deformation assumption), without the bilinear term.

The displacement components within subcell $(\beta\gamma)$ are, therefore, expanded into second-order polynomials, as follows:

[TABLE]

[TABLE]

[TABLE]

The resulting strain components in subcell $(\beta\gamma)$ are given by:

[TABLE]

By averaging the equilibrium equations (4) over the subcell volume (area in the two-dimensional case considered here), one obtains, in conjunction with Eqs. (6) and (8), the following relations:

[TABLE]

where the volume-average (area-average in the two-dimensional case) of the ponderomotive force-to-unit-volume, used in Eq. (6), can be expressed as:

[TABLE]

and the volume (area) average of the principle components of the thermal stress involving the temperature gradient (resulting from the divergence of the thermal stress tensor, which is spherical in the thermally isotropic case considered here) is given by:

[TABLE]

By utilizing the expressions given for the ponderomotive force in Eqs. (4-98), the volume average (area-average here) of $\bar{F}^{(\beta\gamma)}_{i}$ can be readily determined. As for the components of the volume (area) average of the temperature gradient, these can be determined from the expression for $\theta^{(\beta\gamma)}$ given by Eqs. (100), (102) and (103), assuming the expansion coefficients are solved for.

It should be noted that instead of writing the equilibrium equations in terms of edge-averaged tractions expressed through stress-strain relations via the coefficients expanding the displacement field, the equivalent direct volume (area) averaging on the equilibrium equations was employed, making use of the divergence theorem applied to a subcell.

Now, with the stress tensor given be Eq. (8), the line-average stress components along the subcell edges can be expressed as follows:

[TABLE]

By substituting Eqs. (8) and (6) in Eq. (6), the following expressions are obtained:

[TABLE]

where the definition of edge-averaged thermal stresses is introduced, as follows:

[TABLE]

By recalling Eq. (103), edge-averaged thermal stresses can be obtained explicitly, as:

[TABLE]

The line-averaged displacement components along the subcell edges are similarly defined as:

[TABLE]

Substituting the displacement expansions given in Eqs.(176)-(178) into Eq. (6) gives:

[TABLE]

Next, static condensation of the equilibrium equations with three of the coefficients of the expansion of the displacement field vector is performed. The remaining twelve expansion coefficients are consequently related to the twelve edge-averaged displacements, which are then tied to the twelve displacement/traction continuity conditions of each subcell.

First, manipulation of pairs in Eqs. (6), results in:

[TABLE]

and

[TABLE]

Next, substituting $W^{(\beta\gamma)}_{i(20)}$ and $W^{(\beta\gamma)}_{i(02)}$ into Eq. (6) yields the following expressions for $W^{(\beta\gamma)}_{i(00)}$:

[TABLE]

By substituting the expressions for $W^{(\beta\gamma)}_{i(mn)}$ given in Eq. (6)-(6) into the edge-averaged traction components in Eq. (6), the following matrix equation is obtained for subcell $(\beta\gamma)$:

[TABLE]

where ${\bm{K}}^{(\beta\gamma)}$ is a $12\times 12$ matrix, the elements of which depend on the dimensions of subcell $(\beta\gamma)$ and the material properties in this subcell, and where the ponderomotive force and temperature distributions’ contributions are given by:

[TABLE]

Evidently, the total number of unknowns $U^{\pm(\beta\gamma)}_{2i},U^{\pm(\beta\gamma)}_{3i}$ in all subcells is $12N_{\beta}N_{\gamma}$.

Now, continuity of tractions (forces) between subcells is expressed in the following equations:

[TABLE]

whereas continuity of edge-averaged displacements reads:

[TABLE]

These conditions form $12N_{\beta}N_{\gamma}-6N_{\beta}-6N_{\gamma}$ equations.

Finally, assuming that the external boundaries of the composite specimen are traction-free, the global boundary conditions add the following set of equations:

[TABLE]

where $\beta=1,...,N_{\beta}$ and $\gamma=1,...,N_{\gamma}$. These conditions form additional $6(N_{\beta}+N_{\gamma})$ equations, which along with the aforementioned interface continuity conditions provide the total of required $12N_{\beta}N_{\gamma}$ equations. The solution of this sparse system of algebraic equations provides the variables $U^{\pm(\beta\gamma)}_{2i}$ and $U^{\pm(\beta\gamma)}_{3i}$, from which the displacements, strains and stresses can be determined at any point in the composite.

The final step in this procedure is proper placement of the composite specimen in space, the so-called displacement pinning, used in problems with all-free global boundary conditions to prevent matrix singularity in the final equation system, related to rigid body motion. In the considered case, in order to prevent rigid-body displacement, the entire edge-averaged displacement vector of the left and bottom edges of the bottom left subcell is set to zero:

[TABLE]

Accordingly, static condensation of the overall system of equations is performed, eliminating these six displacement-related variables, along with six equations they are involved in.

7 Applications

The presented theory applies to composites with any number of fibers of arbitrary cross-section shape that should be distributed in a symmetric fashion. The theory is next illustrated for the case of a composite of a rectangular cross-section consisting of an epoxy matrix reinforced by two symmetrically allocated iron fibers of round cross-section, see Fig. 1(c). The iron volume ratio is: $v_{f}$ = 0.7. The composite is divided into $N_{\beta}$ = 100, $N_{\gamma}$ = 50 identical square subcells. The cross-sectional dimensions are macroscopic, with: H = 0.1 m and L = 0.05 m. The material properties of the iron fibers and the epoxy matrix are given in Table 1. In regard with the thermal problem, the convective boundary conditions were applied once with a Nusselt number of Nu = 1, which corresponds to nearly pure heat conduction from the composite to the air through a boundary layer of the order of magnitude of the cross-section, and once for a Nusselt number of Nu = 1000, corresponding to forced convection into ambient air, with the maximum typical convection coefficient, as given in Holman (1976). Clearly, in practical applications, more precise knowledge of the thermal problem in the air surrounding the specimen would be beneficial. Here, however, only the solid-mechanical problem is treated, and the convection coefficient is assumed to be known in advance.

The volume-averages of the $x_{2}$ components of the ponderomotive body force-per-unit-volume in the upper ($f$ = 1) and lower ($f$ = 2) fiber are defined by:

[TABLE]

where $\bar{F}^{(\beta\gamma)}_{2}$ is as defined in Eq. (181). It should be interesting to consider the pure particle-in-vacuum Lorentz force defined in a subcell $(\beta\gamma)$ by:

[TABLE]

This forms a counterpart of the ponderomotive force defined in Eq. (181). Furthermore, the average of the second components of the Lorentz force in the upper ($f=1$) and lower ($f=2$) fiber can be defined by:

[TABLE]

Following one of the perspectives of the present work, the forces defined above are used to illustrate the contribution of the magnetization-induced-force to the total ponderomotive force.

The results of the calculations for the specific case outlined above, in terms of magnetic field intensities, magnetic flux densities, ponderomotive force components, temperature, geometric strain tensor components, and equivalent stress distributions, are given in color-maps in Figs. 2,3,5,6 and 7 for applied electric field of: $E_{0}$ = 0.003 V/m. Fig. 4 presents the comparison of the average ponderomotive force-per volume in one fiber using once only the Lorentz part and once the total Maxwell-Minkowski model for the considered regime, for different values of the applied electric field. Figs. 8 and 9 exhibit the strain response of the composite for increasing values of the applied electric field in the range rendering the constitutive laws employed linear as assumed.

Table 1. Material constants for the iron fibers and epoxy matrix:

[TABLE]

The constant scalars $E$, $\nu$, $\alpha$, $\sigma$, $\mu_{r}$ and $\kappa$ denote Young’s modulus, Poisson’s ratio, coefficient of thermal expansion, electric conductivity, relative permeability and thermal conductivity, respectively.

8 Analysis of Results

Fig. 2(a) shows a color-map of the distribution of $H_{3}$ for $E_{0}$ = 0.003 V/m. One notes the continuity in the $x_{2}$ direction across material interface between iron and epoxy. This is expected since $x_{2}$ is the direction of variation of the tangential component of the magnetic field intensity, which is just $H_{3}$. Fig. 2(b) shows a color-map of the distribution of $B_{3}$. This time the nearly identical values of $H_{3}$ near a material interface pointing in the direction of $x_{2}$ are multiplied by permeabilities differing by orders of magnitude. Accordingly, the appearance of tangential continuity vanishes. In addition, one may clearly see the magnetic flux density variation with $x_{2}$, which shows negative values in the lower fiber, zero in the middle of the cross-section and positive values in the upper fiber, not too different from the distribution for the case of a single fiber, where the magnetic flux density varies linearly with distance from the center.

Fig. 3 reveals quasi-linear variation of the ponderomotive force with the coordinates, corresponding to $E_{0}$ = 0.003 V/m, going through zero in the symmetry axes, implying self-contraction (material on the right experiences a negative force acting rightwards, material on the left bears positive force rightwards, and the same goes for the vertical direction) of the widths of the fibers and a vertical attraction force between the two fibers, as expected from experimental observation of two equal-sign-current-bearing wires, with the additional result that the force density vanishes at the center, which is obligatory for an attraction force that has to change sign in the middle. Also, as expected, the force vanishes in the matrix, since there is no current or magnetization there.

Fig. 4 shows the variation of the vertical component of the ponderomotive force with applied electric field, where the dashed line corresponds to calculation with artificially omitted magnetization (in the expression for the ponderomotive force). One finds the expected quadratic dependence. The interesting effect here is how magnetization builds up to decrease the force experienced by the material, thus diminishing internal magnetic energy increase caused by external work, minimizing free energy, which is expected. The interesting part is the extent of this backlash, which reaches almost half the size of the primary phenomenon. This is of course due to the large magnetic susceptibility of iron.

Fig. 5 shows a color-map of the distribution of the temperature deviation from ambient temperature for $E_{0}$ = 0.003 V/m. One can acknowledge the almost uniform distribution within the fibers, resulting from the relatively high thermal conductivity of iron, and the relatively pronounced non-uniform temperature profile in the epoxy, corresponding to its poor thermal conduction, the high temperature in the interface with the fiber, resulting from Joule heating, and the convection at the global boundary. Fig. 5(a) shows the case of convection by conduction, with Nu = 1, where the maximum temperature increase is about four degrees kelvin, which is small enough to still assume constant conductivities.

Fig. 5(b) shows the distributions for Nu = 1000, representing the upper limit of the typical values in forced convection by air. The corresponding temperature deviation there is around one-fiftieth of degree kelvin, which is reasonable for a case where generated Joule heating is very effectively removed by convection at the boundary. Finally, the smoothness of the distribution, although not surprising for a Poisson equation with a constant inhomogeneous part, implies that the discretization was sufficient, at least for the thermal problem.

Figs. 6-9 are dedicated to the mechanical response of the partially conductive composite to applied electric field $E_{0}$ = 0.003 V/m, which produces two competing effects:

the ponderomotive force, which causes transverse fiber-contraction and intra-fiber attraction that lead to compressive strain, and, in the same time heating, causing thermal expansion.

In Fig. 6(a), concerning $\epsilon_{22}$ for Nu = 1, although one might have expected maximum compressive strain in the narrow piece of epoxy locked between the two fibers, the strain is actually mostly expansive there. The reason is the fact that the temperature increase there is as high as it gets within a fiber, whereas the ponderomotive force vanishes there, since it is in the axis of symmetry where it changes sign. Overall compressive strain does exists. It is vertically around the centers of the fibers, closer to the global boundary, where the ponderomotive force is maximal, there is little epoxy to be thermally expanded, but enough to be influenced by compressive stresses from the iron. Fig. 6(b) shows the distribution of $\epsilon_{22}$, which does not appear counterintuitive. The maximum expansion is at the boundaries around the level of the centers of the fibers. The ponderomotive force there causes only displacement toward the center of the specimen, without strain. Thermal expansion therefor dominates. Fig. 6(c) shows $\epsilon_{23}$, which appears to have maximum values at directions positioned forty-five degrees to the principle axes of the specimen, as expected in linear isotropic thermo-elasticity. Fig. 6(d) shows the second-invariant (von Mises) equivalent of the deviatoric stress in each subcell. Due to multiplication of strain by elastic moduli, it is, as expected, a simple map, showing high stress at the stiff iron fibers and lower stress in the epoxy. The stresses ratio is an order of magnitude, though the moduli ratio is two orders of magnitude. This owes to the fact that the strains in the epoxy are an order of magnitude higher, which in turn results from thermal domination in the strains, along with an order of magnitude higher thermal expansion coefficient in the epoxy.

Fig. 7 shows the mechanical response for $E_{0}$ = 0.003 V/m for the case of forced convection with Nu = 1000. This case is representative of the purely magnetic effect on strain. The temperature increase is so small, due to efficient heat transfer, that nearly the entire specimen is in compression in the direction of $x_{2}$, as shown in Fig. 7(a). Perhaps the most interesting feature of the entire example is the short stripe of epoxy just between the two fibers, where the compressive strain $\epsilon_{22}$ is maximal. The advantage of Fig. 7(a) is that it is representative of the macroscopic physics in the Nu$\gg$1 case: two current-bearing wires embedded in insulation, attracted to each other, thus compressing the epoxy between them. Fig. 7(b) shows $\epsilon_{33}$, in which the prominent feature is the expansive strain in the epoxy circumventing the fibers in the region between them, where the ponderomotive force is small in absolute value and the coefficient of thermal expansion is large. Fig. 7(c) shows the typical picture of extremal shear stresses in forty-five degrees to the symmetry axes, only unlike in Fig. 6(c), the larger strains are in the middle of the specimen, since the strains in the edges are more cylindrical, given that the ponderomotive force is compressive in both directions and the thermal strain is an order of magnitude smaller and cannot change the sign of the ratio of the 22 and 33 strain components near the specimen edges. Fig. 7(d) is very interesting although intuitive. It shows distinctive stress concentration in the fibers around the area of their contact with the layer of epoxy that separates them. This is what one would expect in the case of purely attractive forces, which is effectively the situation when there is efficient forced convection removing all current-generated heat. One notes the nice round isostress contourlines of atmospheric magnitude in the proximity of the near-contact region.

Figs. 8 and 9 show overall specimen horizontal and vertical strain-response curves for varying applied electric field, varied in a range retaining the constitutive assumptions reliable. Fig. 8 shows the normal stresses in both directions for increasing $E_{0}$ values for Nu = 1, after averaging the strain over the specimen volume (which is justified, strain being volume integrable in the small-strain limit). The quadratic nature of the response is preserved, the ponderomotive force being quadratic in current and the equilibrium equation being linear in strain. Both volume-averaged strain components are positive, signaling that the principle effect for Nu = 1 is thermal expansion. One notes that $\epsilon_{33}$ is somewhat larger than $\epsilon_{22}$, which should be attributed to the fact that in the direction of $x_{2}$ there is more iron to produce overall larger ponderomotive force. As can be seen from Fig. 3, $\underset{x_{2},x_{3}}{\text{max}}\ \bar{F}_{2}(x_{2},x_{3})\approx 2\ \underset{x_{2},x_{3}}{\text{max}}\ \bar{F}_{3}(x_{2},x_{3})$. Consequently, the negative part of $\epsilon_{22}$ is larger than in $\epsilon_{33}$, and thus it is smaller by a finite observable amount.

Fig. 8 is interesting from the second perspective of the present work, which is computational micromechanical analysis of the steady-state response of a conductive-fiber–insulative-matrix composite to electric field, in the numerical scale of electrical-engineering applications (that is, where an insulated large-scale cable with two wires is considered – even though usually not both wires in a cable bear current), using the strong-form approach to the solution of the multi-physics boundary value problem.

Fig. 9 serves the two additional perspectives of this research. One is the construction and illustration of a computational apparatus supporting a sensing/actuating device. Due to the large Nusselt number, thermal expansion is negligible and thus the average strain response of the specimen for Nu = 1000 shown in Fig. 9 is negative, that is, there is macroscopic compression, as one would expect from a two-wire insulated cable. The reason why this result is more illustrative of the idea of an actuator is also the fact that both average compression and expansion can be obtained if the currents in the wires are of equal or opposite signs. Although it is assumed throughout the paper that the electric field is given, the same apparatus works if instead the current in each wire is given. Then, by applying opposite sign currents, expansion $\epsilon_{22}$ values can be obtained, as well as compression for equal sign currents. This, together with the fact that air cooling can be applied easier in the case of a device than just in a general electric cable application, makes the result in Fig. 9 be representative of the third perspective of the paper as stated in the introduction, i.e, that of device engineering (sensor/actuator).

Finally, the dashed curve in Fig. 9 represents the first perspective of the paper as stated in the introduction, which is the examination of the quantitative effect of accounting for the magnetization part of the ponderomotive force in the mechanical response in the case of a ferromagnetic conductor (in the linear regime).

Clearly, neglecting magnetization results in overestimation of the average compression strain of the specimen. This is rather intuitive. The surprising part is that although, as shown in Fig. 4, the total ponderomotive force is about half the Lorentz force part, Fig. 9 shows that the average strain when neglecting the magnetization force part of the ponderomotive force, is not twice larger than the actual strain, but about five times larger (in absolute value). The reason for this apparent discrepancy is in the interplay with thermal expansion, which is still present in the iron fibers even for Nu = 1000. In the case of calculation with the total ponderomotive force, the magnetic part of $\epsilon_{22}$ in the fibers nearly cancels-out the thermal strain there. The result is that $\epsilon_{22}$ in the fibers can add nearly nothing to the volume-averaged value. If the magnetization part of the force is neglected, the magnetic strain in the fibers overcomes the thermal expansion, and the total negative-strain contribution from the fibers adds to the average value. This effect is stronger than doubling, since the nearly zero strain exists in the realistic calculation not only in the fibers but almost in the entire specimen except for in a small layer between the fibers. Finally, the strain in the aforementioned layer is also at least doubled when only the ’particle-in-vacuum’ Lorentz force is considered.

9 Conclusions

In this article, a theory was presented, in which a computational procedure for micromechanical steady state response analysis of a thermally conducting linear elastic composite with electrically conductive ferromagnetic fibers embedded in an insulating matrix, was derived. The crux of the work is the construction of a stable computational scheme for the derivation of the ponderomotive force induced by an applied electric field. The expression for the force, along with thermal analysis is consequently used in conjunction with mechanical analysis. The mechanical analysis discretizes the domain into subcells, every one of which is filled with homogeneous material. The equations of elasticity are solved for the array of subcells in strong form – that is not by minimizing a variation of a functional, as in the finite element approach, but by directly solving the boundary value problem. Quadratic expansion of the displacement vector within each subcell is employed. This approach, although producing a stable solution scheme for the elastic problem, does not allow solving the magnetic problem, where quadratic or any other standard interpolation produces a singular matrix relation between the interpolation coefficients and the edge-averaged values. Due to this fact, a different approach was taken here, in which linear combination of point-wise analytical solutions of the differential operator representing the magnetic problem was employed, consequently establishing the coefficients by enforcing local boundary conditions in the integral sense. It should be noted that in the case of the magnetic problem, additional benefit was gained from constructing the solution from point-wise analytical solutions of the differential equations. This gain is embodied in the fact that the introduction of an auxiliary gauge condition, such as the Coulomb gauge, became unnecessary, since the magnetic field was not assumed to be the curl of an auxiliary vector field, but was rather constructed as the superposition of the Ampère solution in a conductor and the gradient of a harmonic function. This is, of course, only possible in the linear case.

In the thermal problem, an approach similar to the one taken in the magnetic case was chosen, although quadratic interpolation and fulfillment of the differential equation only in the integral sense could have been an alternative possible route.

With the ponderomotive force obtained analytically for each subcell, and the temperature field obtained analytically up to four coefficients for each subcell determined numerically, the mechanical problem was solved consistently, enforcing mechanical equilibrium on the subcell level and requiring continuity of the integral quantities of surface forces on subcell boundaries, as well as spatially integrated displacements.

The numerical example solved with the derived theory showed reasonably intuitive results, which in itself is a certain validation of the method. It was demonstrated that the approach is suitable for the multi-physics analysis of a current-bearing two-wire cable, where ponderomotive and thermal effects interplay is captured quite reasonably. Furthermore, it is the assertion of this paper that the derived theory can be used as the computational apparatus for, and in fact suggests a sensing/actuating device with quadratic mechanical response to either uniform electric field, if a straight long wire is considered, or the large-scale gradient of electric field, if a two-wire thin ringed cable is considered. The present analysis is limited to steady-state response of such a device. Future work may extend the analysis to include transient response.

Regarding the response to applied electric field, it may also be noted that the ability to respond to the large-scale gradient of electric field, rather than to its local value, which may be obtained in a thin ringed device (owing to the fact that current will be induced in the ring only in a nonuniform electric field – the local smoothness is required, though, for the electric problem to be solved by an electric field uniform throughout the composite’s cross-section), may be useful in sensing/applications where high sensing resolution is required.