Elastic and Thermoelastic Responses of Orthotropic Half-Planes

Abstract

A unified approach is presented for constructing explicit solutions to the plane elasticity and thermoelasticity problems for orthotropic half-planes. The solutions are constructed in forms which decrease the distance from the loaded segments of the boundary for any feasible relationship between the elastic moduli of orthotropic materials. For the construction, the direct integration method was employed to reduce the formulated problems to a governing equation for a key function. In turn, the governing equation was reduced to an integral equation of the second kind, whose explicit analytical solution was constructed by using the resolvent-kernel algorithm.

1. Introduction

The design of engineering applications and advancement of material science are driven projection and implementation of new composite materials with averaged properties meeting specific requirements [1]. These requirements are usually formulated through both empirical evidence and a priori analysis implemented by inverse and optimization algorithms [2,3,4]. Such algorithms, in turn, are based on the complex analysis of primary responses of multi-phase materials to the collective action of physical impacts, accounting for a wide range of operational and constructional features [5,6] (e.g., anisotropy of materials, fiber packing, and stratification). Moreover, the efficiency in developing optimization algorithms and procedures for the inverse identification of effective material parameters for advanced composite materials are directly related to the solutions of problems in solid mechanics [7,8,9].

A homogenization-based approach is widely used to predict the effective mechanical or thermal properties in solid mechanics problems. Due to the large number of layers forming, for instance, a stratified composite, and the need to meet the continuity conditions for the temperature and the components of displacement, stress, and heat-flux vectors, the implementation of approaches based on the classical theory of thermoelasticity becomes problematic regarding the related boundary value problems for composite materials. Therefore, it is necessary to use specific approximate models that meet the original interface conditions on the surfaces connecting the individual components. These models are created by applying various averaging procedures, e.g., the thick plate theory approach [10], mixture theory approach [11], asymptotic homogenization [12], matrix methods [13], and tolerance theory [14]. Another approach is homogenization based on non-asymptotic techniques derived from non-standard analysis, which was proposed in [15] and adapted for composites of layered structure [16]. Within the framework of this approach, the homogenized model with microlocal parameters is now widely used in composite materials research [17,18]. The new approach to this model has also found application in modeling layered composites with functionally graded properties for thermoelastic problems [19].

An interesting approach for handling composite materials consisting of a base layer covered with thin multilayer coatings of multifunctional purpose is based on the method of generalized boundary condition [20,21]. This approach captures the effect of thin multilayer coating and avoids computational complications due to the mismatch of the geometrical and thermophysical properties on the layers’ interfaces.

Efficient analysis of the elastic and thermoelastic response of composite materials exhibiting anisotropic features significantly depends on the accuracy of accounting for the dissimilarity of the material moduli in different spatial directions [22,23,24]. Hence, the methods for such analysis of isotropic solids may largely fail when attempted directly on anisotropic solids. If, for example, the eigenfunction method was used (separation of variables) to studying local effects in semi-infinite elastic composites, the desired solution is to decay when moving away from the loaded area according to the Saint-Venant principle. In contrast to the isotropic case, however, the eigenvalues of such problems for anisotropic solids will no longer be defined uniquely but depend on the interrelation between the elastic moduli involved into the coefficients of characteristic equations [6,25]. This affects the type and character of corresponding eigenfunctions and complicates them, ensuring their vanishing behavior at distant points. The latter drawback limits the applicability of such solutions for the entire spectra of feasible anisotropic moduli.

In this paper, we extend a technique that avoids the latter complication to represent a solution for an elastic orthotropic half-plane in a unified form that explicitly depends on the locally applied loadings and ensures decaying at a distance from the loaded zone. The technique is based on applying the direct integration method [6], which implies the reduction of a corresponding boundary value problem to the governing equations for individual stress tensor components by using the relationships derived via the integration of the equilibrium equations. This method has been efficiently used to analyze numerous direct and inverse boundary value problems in different coordinate systems [26,27]. In [28], this method was used to develop a technique for solving a three-dimensional elasticity problem for a transversely isotropic half-space. This technique reduces governing equations of the problem to integral equations of the second kind. With the application of the resolvent-kernel algorithm for solving the derived integral equations, solutions are constructed in forms that are irrespective of the interdependences between the elastic material moduli and decays on the distance from the loaded zones. Moreover, the solution is given in the form of explicit expressions on the loadings. These features make it quite attractive for the identification, inverse, and optimization algorithms [29].

2. Formulation of the Problem

Consider a plane thermoelasticity problem for an elastic orthotropic half-plane $(x,y)\in (-\infty ,\infty )\times [0,\infty )$ in the dimensionless Cartesian coordinate system $(x,y)$. The problem is governed [22,30] by the equations of equilibrium in terms of stresses

$$\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xy}}{\partial y}=0,\hspace{1em}\frac{\partial {\sigma }_{xy}}{\partial x}+\frac{\partial {\sigma }_{yy}}{\partial y}=0$$

(1)

and the strain-compatibility equation

$$\frac{{\partial }^2{\epsilon }_{xx}}{\partial y^2}+\frac{{\partial }^2{\epsilon }_{yy}}{\partial x^2}=\frac{{\partial }^2{\epsilon }_{xy}}{\partial x\partial y}$$

(2)

Here, ${\sigma }_{xx}$, ${\sigma }_{yy}$, ${\sigma }_{xy}$ and ${\epsilon }_{xx}$, ${\epsilon }_{yy}$, ${\epsilon }_{xy}$ are the in-plane normal and tangential components of the stress- and strain-tensor, respectively. These components are related via the constitutive equations of the Duhamel–Neumann law for orthotropic materials [31], which takes the following form

$${\epsilon }_{xx}=a_{11}{\sigma }_{xx}+a_{12}{\sigma }_{yy}+{\alpha }_{1}T,\hspace{1em}{\epsilon }_{yy}=a_{12}{\sigma }_{xx}+a_{22}{\sigma }_{yy}+{\alpha }_{2}T,\hspace{1em}{G}_{xy}{\epsilon }_{xy}={\sigma }_{xy}$$

(3)

under both the plane strain and plane stress hypotheses. Here,

$$\begin{array}{c}{a}_{11}=\frac{{\tilde{a}}_{11}}{E_x},\hspace{1em}{a}_{22}=\frac{{\tilde{a}}_{22}}{E_y},\hspace{1em}{a}_{12}=-\frac{{\tilde{a}}_{12}}{E_y}=-\frac{{\tilde{a}}_{21}}{E_x},\\ {\tilde{a}}_{11}=\left\{\begin{array}{l}1,\\ 1-{\nu }_{zx}{\nu }_{xz},\end{array}\right.\hspace{1em}{\tilde{a}}_{22}=\left\{\begin{array}{l}1,\\ 1-{\nu }_{zy}{\nu }_{yz},\end{array}\right.\\ {\tilde{a}}_{12}=\left\{\begin{array}{l}{\nu }_{yx},\\ {\nu }_{yx}+{\nu }_{zx}{\nu }_{yz},\end{array}\right.\hspace{1em}{\tilde{a}}_{21}=\left\{\begin{array}{l}{\nu }_{xy},\\ {\nu }_{xy}+{\nu }_{xz}{\nu }_{zy},\end{array}\right.\\ {\alpha }_1=\left\{\begin{array}{l}{\alpha }_x,\\ {\alpha }_x+{\alpha }_z{\nu }_{zx},\end{array}\right.\hspace{1em}{\alpha }_2=\left\{\begin{array}{l}{\alpha }_y,\\ {\alpha }_y+{\alpha }_z{\nu }_{zy},\end{array}\right.\end{array}$$

(4)

$E_x$ and $E_y$ are the Young moduli in the directions $x$ and $y$, respectively; ${\nu }_{j{l}}$ are the Poisson ratios describing the contraction in the ${l}$-direction at tension in the j-direction, ${l},j=\{x,y,z\},\hspace{0.33em}{l}\ne j$; $G_{xy}$ is the shear modulus for the coordinate plane $(x,y)$; ${\alpha }_x,$ ${\alpha }_y$, and ${\alpha }_z$ are the linear thermal expansion coefficients in the directions $x,y,$ and $z$, respectively; $T(x,y)=T^{\ast }(x,y)-T_0^{\ast }(x,y)$, $T^{\ast }(x,y)$ is the temperature field, which can be found from a corresponding heat conduction problem [31], and $T_0^{\ast }(x,y)$ is the reference temperature of the non-stressed state. In the foregoing expressions with braces, the upper and lower lines correspond to the plane stress and plane strain cases, respectively. The orthotropic elastic modules in (4) comply with the following symmetry condition [22]

$$E_j{\nu }_{{l}j}=E_{{l}}{\nu }_{j{l}},\hspace{1em}{l},j=\{x,y,z\},\hspace{1em}\hspace{1em}\hspace{0.33em}{l}\ne j.$$

Assume boundary $y=0$ of the half-plane to be acted upon by locally distributed normal and tangential force loadings imposed by the following boundary conditions:

$${\sigma }_{yy}=-p(x),\hspace{1em}{\sigma }_{xy}=q(x),\hspace{1em}\hspace{1em}y=0,\hspace{1em}|x|\hspace{0.33em}<\infty ,$$

(5)

where $p(x)$ and $q(x)$ are the given functions of some local distribution profiles, i.e., $\underset{|x|\hspace{0.33em}\to \hspace{0.33em}\infty }{\lim }p(x)=$$\underset{|x|\hspace{0.33em}\to \hspace{0.33em}\infty }{\lim }q(x)=0$. Note that for the correct formulation of the boundary value problem, the force loadings (5) and temperature $T(x,y)$ are to meet the necessary conditions derived in [32].

3. Reduction to Governing Equations

By following the strategy of the direct integration method [6,26], consider the normal stress ${\sigma }_{yy}$ and the in-plane total stress $\sigma ={\sigma }_{xx}+{\sigma }_{yy}$ to be the key functions. By eliminating the tangential stress from equilibrium Equation (1), we derive the following governing equation for the key functions:

$$\Delta {\sigma }_{yy}=\frac{{\partial }^2\sigma }{\partial x^2}$$

(6)

where $\Delta ={\partial }^2/\partial x^2+{\partial }^2/\partial y^2$ is the in-plane Laplace differential operator. This equation implies general equilibrium in terms of the stress-tensor components and thus is irrespective of the material properties.

Making use of the constitutive Equation (3) along with the equilibrium ones (1), we can represent the compatibility Equation (2) in terms of stresses, as follows:

$$\Delta 〈a_{11}\sigma +{\alpha }_{1}T〉=({\alpha }_1-{\alpha }_2)\frac{{\partial }^{2}T}{\partial x^2}+\left(2{\beta }_1-\frac{1}{G_{xy}}\right)\frac{{\partial }^2{\sigma }_{yy}}{\partial y^2}-{\beta }_2\frac{{\partial }^2{\sigma }_{yy}}{\partial x^2}$$

(7)

Here, ${\beta }_1=a_{11}-a_{12},$ ${\beta }_2=a_{22}-a_{11}$. Within the context of expressions (4), the coefficients in governing Equation (7) involve the orthotropic moduli of the material. Note that this equation can also be obtained by assuming the material moduli constant in an analogous equation [33] derived for inhomogeneous orthotropic material. In the case of isotropic material, Equation (7) can naturally be reduced to the corresponding governing equation derived in [32].

When acting upon Equation (7) by the differential operator ${\partial }^2/\partial x^2$ and making use of Equation (6), we can derive the following fourth-order differential equation:

$$\frac{{\partial }^4{\sigma }_{yy}}{\partial y^4}+2a_1\frac{{\partial }^4{\sigma }_{yy}}{\partial x^2\partial y^2}+a_2\frac{{\partial }^4{\sigma }_{yy}}{\partial x^4}=-\frac{{\partial }^2}{\partial x^2}\left({\alpha }_1^0\frac{{\partial }^{2}T}{\partial y^2}+{\alpha }_2^0\frac{{\partial }^{2}T}{\partial x^2}\right)$$

(8)

for the stress ${\sigma }_{yy}$, solely, where the coefficients are given as

$$a_1=\frac{1}{a_{11}}\left(a_{12}+\frac{1}{2G_{xy}}\right),\hspace{1em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{a}_2=\frac{a_{22}}{a_{11}},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{\alpha }_j^0=\frac{{\alpha }_j}{a_{11}},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}j=1,2$$

(9)

By using the second equilibrium Equation (1) together with the second boundary condition (5), we can derive a condition for the derivative of the key function in the form as follows:

$$\frac{\partial {\sigma }_{yy}}{\partial y}=-\frac{dq(x)}{dx},\hspace{1em}\hspace{1em}y=0,\hspace{1em}|x|\hspace{0.33em}<\infty .$$

(10)

In such a manner, the original thermoelasticity problem is reduced to a governing Equation (8) with coefficients (9) expressed through the orthotropic elastic moduli (4) under the boundary conditions presented by the first condition (5) and condition (10) for the key function ${\sigma }_{yy}$. After this function is found by solving the mentioned boundary value problem, the in-plane total stress $\sigma$ can be determined from Equation (6). The normal stress ${\sigma }_{xx}$ can then be found via the formula ${\sigma }_{xx}=\sigma -{\sigma }_{yy}$. The tangential stress ${\sigma }_{xy}$ can be resulted as follows

$${\sigma }_{xy}=\frac{q(x)}{2}-\frac{1}{2}{\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{\partial {\sigma }_{xx}(x,\eta )}{\partial x}\mathrm{sgn}(y-\eta )\mathrm{d}\eta }=-\frac{1}{2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{\partial {\sigma }_{yy}(\xi ,y)}{\partial y}\mathrm{sgn}(x-\xi )\mathrm{d}\xi }$$

based on the integration of the equilibrium Equation (1) within the context of boundary condition in (5).

4. Solution of the Governing Equations

To solve the formulated problem, we employ the integral Fourier transform with respect to variable $x$:

$$\overline{f}(y)=\overline{f}(y;s)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f(x,y)\exp (-isx)\mathrm{d}x}$$

(11)

where $s$ is the transform parameter, $i^2=-1$, and $f(x,y)$ is an arbitrary function whose transform (11) exists in at least the generalized sense [34].

Having applied transform (11) to Equations (8) and (10) along with the first condition in (5), we arrive at the following ordinary differential equation:

$$\frac{{\mathrm{d}}^4{\overline{\sigma }}_{yy}}{\mathrm{d}{y}^4}-2s^2{a}_1\frac{{\mathrm{d}}^2{\overline{\sigma }}_{yy}}{\mathrm{d}{y}^2}+s^4{a}_2{\overline{\sigma }}_{yy}=s^2\left({\alpha }_1^0\frac{{\mathrm{d}}^2\overline{T}}{\mathrm{d}{y}^2}-s^2{\alpha }_2^0\overline{T}\right)$$

(12)

accompanied with the conditions

$${\overline{\sigma }}_{yy}=-\overline{p},\hspace{1em}\frac{\mathrm{d}{\overline{\sigma }}_{yy}}{\mathrm{d}y}=-is\overline{q},\hspace{1em}\hspace{1em}y=0.$$

(13)

The characteristic equation of (12) takes the following form

$${\mu }^4-2a_1{(s\mu )}^2+s^4{a}_2=0$$

(14)

with $\mu$ being the eigenvalues, which can be given as

$$\begin{array}{c}{\mu }_j={(-1)}^j{\lambda }_{1}s,\hspace{1em}\hspace{1em}{\mu }_{2+j}={(-1)}^j{\lambda }_{2}s,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}j=1,2,\hspace{0.17em}\hspace{0.17em}\\ {\lambda }_1=\sqrt{a_1-\sqrt{a_1^2-a_2}}={\lambda }_1^{\mathrm{Re}}+i{\lambda }_1^{\mathrm{Im}},\hspace{1em}\hspace{1em}{\lambda }_2=\sqrt{a_1+\sqrt{a_1^2-a_2}}={\lambda }_2^{\mathrm{Re}}+i{\lambda }_2^{\mathrm{Im}}.\end{array}$$

(15)

As it can be seen from formulae in (15), the type of the eigenvalues critically depends on the interrelation between the coefficients $a_1$ and $a_2$, which is verbalized by expressions for ${\lambda }_1$ and ${\lambda }_2$. Thus, in order to be able of capturing all possible features of the solution to Equation (12) for various orthotropic materials and ensuring its vanishing with $y\to \infty$, it is necessary to exhaustively analyze all possible cases of interrelation between coefficients of Equation (14), which within the context of formules (4) and (9), implies the analysis of interdependences of the orthotropic material moduli. Based on the theoretical and experimental studies [22,35], the following constraints can be established for the elastic moduli involved into expressions for the coefficients (9):

$$\begin{array}{c}{E}_{{l}}>0,\hspace{1em}{G}_{12}>0,\hspace{1em}{l}=1,2,3,\hspace{1em}|{\nu }_{13}+{\nu }_{12}{\nu }_{23}|\hspace{0.33em}<\frac{1}{\sqrt{E_2{E}_3}}\sqrt{E_1-{\nu }_{12}{E}_2}\sqrt{E_2-{\nu }_{23}{E}_3},\\ |{\nu }_{12}|\hspace{0.33em}<\sqrt{\frac{E_1}{E_2}},\hspace{1em}|{\nu }_{21}|\hspace{0.33em}<\sqrt{\frac{E_2}{E_1}},\hspace{1em}|{\nu }_{23}|\hspace{0.33em}<\sqrt{\frac{E_2}{E_3}},\hspace{1em}|{\nu }_{32}|\hspace{0.33em}<\sqrt{\frac{E_3}{E_2}}.\end{array}$$

(16)

In view of the equations in (4), (9), and (16), it can be concluded that

$$a_2>0$$

The feasible dependences (i.e., the ones satisfying conditions (16)) of the real and imaginary parts ${\lambda }_j^{\mathrm{Re}}/\sqrt[4]{a_2}$ and ${\lambda }_j^{\mathrm{Im}}/\sqrt[4]{a_2}$, $j=1,2,$ of eigenvalues (15) on parameter $t=a_1/\sqrt{a_2}$ are shown in Figure 1. The corresponding expressions for ${\lambda }_j^{\mathrm{Re}}/\sqrt[4]{a_2}$ and ${\lambda }_j^{\mathrm{Im}}/\sqrt[4]{a_2}$ are given in

Table A1. Expressions for the real and imaginary parts ${\lambda }_j^{\mathrm{Re}}/\sqrt[4]{a_2}$ and ${\lambda }_j^{\mathrm{Im}}/\sqrt[4]{a_2}$, $j=1,2,$ of eigenvalues (15) on parameter $t=a_1/\sqrt{a_2}$ for the material moduli satisfying physical constraints (16).

$\mathit{t}={\mathit{a}}_1/\sqrt{{\mathit{a}}_2}$	${\mathit{\lambda }}_1^{\mathbf{Re}}/\sqrt[4]{{\mathit{a}}_2}$	${\mathit{\lambda }}_1^{\mathbf{Im}}/\sqrt[4]{{\mathit{a}}_2}$	${\mathit{\lambda }}_2^{\mathbf{Re}}/\sqrt[4]{{\mathit{a}}_2}$	${\mathit{\lambda }}_2^{\mathbf{Im}}/\sqrt[4]{{\mathit{a}}_2}$
$t<-1$	0	0	$\sqrt{|t|+\sqrt{t^2-1}}$	$\sqrt{|t|-\sqrt{t^2-1}}$
$t=-1$	0	0	1	1
$-1<t<0$	$\sin \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{|t|}\right)$	$\sin \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{|t|}\right)$	$-\cos \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{|t|}\right)$	$\cos \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{|t|}\right)$
$t=0$	$1/\sqrt{2}$	$1/\sqrt{2}$	$-1/\sqrt{2}$	$1/\sqrt{2}$
$0<t<1$	$\cos \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{t}\right)$	$\cos \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{t}\right)$	$-\sin \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{t}\right)$	$\sin \left(\frac{1}{2}\arctan \frac{\sqrt{1-t^2}}{t}\right)$
$t=1$	1	1	0	0
$1<t$	$\sqrt{t-\sqrt{t^2-1}}$	$\sqrt{t+\sqrt{t^2-1}}$	0	0

of Appendix A.

The observed variation in the type and multiplicity of the eigenvalues for all feasible dispersions of the orthotropic material moduli affects significantly the form of a general solution to Equation (12), which is critical for assuring its vanishing when moving away from the loaded boundary $y=0$ of the half-plane towards its depth. To the best of our knowledge, the exhaustive coverage of all possible cases of vanishing analytical solutions remained beyond the circle of interest in the vast majority of papers dedicated to the mechanical analysis of semi-infinite orthotropic solids.

The cumbersomeness of the exhaustive analyzes of all possible cases of interdependences between the anisotropic material moduli was demonstrated in [25] when solving a three-dimensional problem in a half-space made of a relatively simpler type of material, i.e., a transversely isotropic one, possessing a smaller number of independent elastic moduli compared to an orthotropic material. In order to avoid this cumbersomeness, it was suggested in [28] to use an alternative way for constructing analytical solutions to such problems by reducing them to integral equations of the second kind with subsequent application of the resolvent-kernel method. Herein, we extend the latter technique for solving the boundary value problem (12) and (13) for the case of orthotropic material.

By following the strategy of paper [28], let us represent Equation (12) in the form as follows:

$$\frac{{\mathrm{d}}^4{\overline{\sigma }}_{yy}}{\mathrm{d}{y}^4}-2s^2{a}^2\frac{{\mathrm{d}}^2{\overline{\sigma }}_{yy}}{\mathrm{d}{y}^2}+s^4{a}^4{\overline{\sigma }}_{yy}=2s^2(a_1-a^2)\frac{{\mathrm{d}}^2{\overline{\sigma }}_{yy}}{\mathrm{d}{y}^2}+s^2\left({\alpha }_1^0\frac{{\mathrm{d}}^2\overline{T}}{\mathrm{d}{y}^2}-s^2{\alpha }_2^0\overline{T}\right),$$

(17)

where $a^4=a_2>0$. A general solution to Equation (17) that is vanishing at $y\to \infty$ can be given as

$$\begin{gathered} {\overline{\sigma }}_{yy}(y)=(A_1+y{A}_2)\exp (-|s|ay) \\ -\frac{|s|}{2}\frac{a_1-a^2}{a}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\overline{\sigma }}_{yy}(\eta )\left(1-|s|a|y-\eta |\right)\exp (-|s|a|y-\eta |)\mathrm{d}\eta } \\ -\frac{|s|}{4a^3}{\displaystyle \underset{0}{\overset{\infty }{\int }}\overline{T}(\eta )(\left(1+|s|a|y-\eta |\right){\alpha }_2^0+a^2\left(1-|s|a|y-\eta |\right){\alpha }_1^0)\exp (-|s|a|y-\eta |)\mathrm{d}\eta }, \end{gathered}$$

(18)

where $A_j$, $j=1,2$, are constants of integration. These constants can be determined by putting (18) into conditions (13). As a result, we arrive at the following equation for the key stress ${\overline{\sigma }}_{yy}$:

$${\overline{\sigma }}_{yy}(y)=-\hspace{0.33em}(\left(1+|s|ay\right)\overline{p}+isy\overline{q})\exp (-|s|ay)+\Theta (y)+{\displaystyle \underset{0}{\overset{\infty }{\int }}{\overline{\sigma }}_{yy}(\eta )\mathcal{K}(y,\eta )\mathrm{d}\eta }$$

(19)

Here,

$$\mathcal{K}(y,\eta )=\frac{|s|}{2}\frac{a_1-a^2}{a}(\left(1-|s|a\eta +|s|ay(3-2|s|a\eta )\right)\exp (-|s|a(y+\eta ))$$

$$-\hspace{0.33em}\left(1-|s|a|y-\eta |\right)\exp (-|s|a|y-\eta |))$$

$$\Theta (y)=\frac{|s|}{4a^3}{\displaystyle \underset{0}{\overset{\infty }{\int }}\overline{T}(\eta )(}\left\{\left(1\hspace{0.33em}+\hspace{0.33em}|s|a\eta \hspace{0.33em}+|s|ay(1+2|s|a\eta )\right){\alpha }_2^0\right.$$

$$\left.+\hspace{0.33em}{a}^2\left(1-|s|a\eta +|s|ay(3-2|s|a\eta )\right){\alpha }_1^0\right\}\exp (-|s|a(y+\eta ))$$

$$-\left\{\left(1+|s|a|y-\eta |\right){\alpha }_2^0+a^2\left(1-|s|a|y-\eta |\right){\alpha }_1^0\right\}\exp (-|s|a|y-\eta |))\mathrm{d}\eta$$

Equation (19) is an integral equation of the second kind. An analytical solution to this equation expressing the key function ${\overline{\sigma }}_{yy}(y)$ through the mapping functions of the force loadings $\overline{p}$ and $\overline{q}$ on the boundary and the temperature $\overline{T}(y)$, can be constructed by means of the resolvent-kernel algorithm [36] to obtain the following:

$${\overline{\sigma }}_{yy}(y)=-\overline{p}{\phi }_p(y)-is\overline{q}{\phi }_q(y)+{\phi }_T(y)$$

(20)

where

$${\phi }_p(y)=\left(1+|s|ay\right)\exp (-|s|ay)+{\displaystyle \underset{0}{\overset{\infty }{\int }}\left(1+|s|a\eta \right)\mathcal{R}(y,\eta )\exp (-|s|a\eta )\mathrm{d}\eta }$$

$${\phi }_q(y)=y\exp (-|s|ay)+{\displaystyle \underset{0}{\overset{\infty }{\int }}\eta \mathcal{R}(y,\eta )\exp (-|s|a\eta )\mathrm{d}\eta }$$

$${\phi }_T(y)=\Theta (y)+{\displaystyle \underset{0}{\overset{\infty }{\int }}\Theta (\eta )\mathcal{R}(y,\eta )\mathrm{d}\eta }$$

and the resolvent kernel is constructed as

$$\mathcal{R}(y,\eta )={\displaystyle \sum _{n=0}^{\infty }{\mathcal{K}}_{n+1}(y,\eta )}$$

(21)

Here,

$${\mathcal{K}}_1(y,\eta )=\mathcal{K}(y,\eta ),\hspace{1em}\hspace{1em}{\mathcal{K}}_{n+1}(y,\eta )={\displaystyle \underset{0}{\overset{\infty }{\int }}{\mathcal{K}}_1(y,\tau ){\mathcal{K}}_n(\tau ,\eta )\mathrm{d}\tau },\hspace{1em}n=1,2,\dots$$

Note that for many practical cases, the resolvent kernel (21) can be found in a closed analytical form [36]. Otherwise, the practical computation of the resolvent can be implemented using the approximate formula

$$\mathcal{R}(y,\eta )\approx {\mathcal{R}}_N(y,\eta )={\displaystyle \sum _{n=0}^N{\mathcal{K}}_{n+1}(y,\eta )}$$

(22)

where $N$ is an integer, which can be evaluated either by a numerical experiment or through an analytical procedure [37].

After constructing the stress ${\overline{\sigma }}_{yy}(y)$ in the form (20), we can use the following formulae:

$${\overline{\sigma }}_{xx}(y)=-\frac{1}{s^2}\frac{{\mathrm{d}}^2{\overline{\sigma }}_{yy}(y)}{\mathrm{d}{y}^2},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{\overline{\sigma }}_{xy}(y)=\frac{i}{s}\frac{\mathrm{d}{\overline{\sigma }}_{yy}(y)}{\mathrm{d}y},$$

(23)

derived from Equation (1) in the mapping domain of the integral transform (11). By substituting (20) into (23), we can represent the stresses ${\overline{\sigma }}_{xx}(y)$ and ${\overline{\sigma }}_{xy}(y)$ in a similar to (20) form, as follows:

$${\overline{\sigma }}_{xx}(y)=-\overline{p}{\Psi }_p(y)-is\overline{q}{\Psi }_q(y)+{\Psi }_T(y)$$

(24)

$${\overline{\sigma }}_{xy}(y)=-\overline{p}{\vartheta }_p(y)-is\overline{q}{\vartheta }_q(y)+{\vartheta }_T(y)$$

(25)

The expressions for functions ${\Psi }_p(y)$, ${\Psi }_q(y)$, ${\Psi }_T(y)$, and ${\vartheta }_p(y)$, ${\vartheta }_q(y)$, and ${\vartheta }_T(y)$ are given in Appendix B.

After constructing all the stress-tensor components in the mapping domain of the transform (11) in the form (20), (24), and (25), we can restore them in the physical domain by implementing the inverse transform

$$f(x,y)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\overline{f}(y;s)\exp (isx)\mathrm{d}s}$$

The constructed expressions (20), (24), and (25) allow for expressing the stress tensor components in an orthotropic half-plane through the applied force loadings and the temperature field explicitly. They appear to be uniform for any feasible case of material properties obeying the constraints (16) and ensuring the stress tensor components exhibit decaying behavior at the infinitely distant points of the half-plane.

5. Numerical Example and Discussion

Consider numerical implementation of the proposed solution algorithm for several case studies of the orthotropic materials with properties presented in

Table 1. Material properties of some orthotropic composites [38].

	Material	${\mathit{E}}_{\mathit{x}}, \mathbf{GPa}$	${\mathit{E}}_{\mathit{y}}, \mathbf{GPa}$	${\mathit{\nu }}_{\mathit{x}\mathit{y}}$	${\mathit{G}}_{\mathit{x}\mathit{y}}, \mathbf{GPa}$
1	Carbon Fiber Composite UP-1	15.70	14.80	0.67	42.00
2	Carbon Fiber Textile Reinforced Plastic UOL-300-1A	150.70	8.00	0.29	4.00

for the case of plane stress adopted in expressions for elastic coefficients in (4). The considered orthotropic composite materials are based on carbon fiber plastics which are of practical implementation for aircraft engineering. The computations are performed for the smooth normal impact of the half-plane boundary (5), where $p(x)=p_0\exp (-\kappa x^2)$ and $q(x)=0$, $p_0=\mathrm{const}$, $\kappa =\mathrm{const}$. It can be shown that the eigenvalues ${\lambda }_1$ and ${\lambda }_2$ in formulae (15) are complex conjugated for material 1 and real and dissimilar for material 2 from Table 1.

The full-field distributions of the dimensionless stress tensor components ${\sigma }_{yy}/p_0$, ${\sigma }_{xx}/p_0$, and ${\sigma }_{xy}/p_0$ are depicted in Figure 2 for material 1 (a, c, e) and 2 (b, d, f) demonstrating the effect of anisotropy in the stress patterns for the considered type of loading. The computations have been performed by making use of the exact solutions to Equation (12).

Figure 3 demonstrates the dimensionless stress ${\sigma }_{yy}/p_0$ computed at $x=0$ and $\kappa =1$ by unified formula (20) for different values of $N$ in the expression for approximate resolvent kernel (22) used instead of the exact one (21) for the orthotropic materials with properties given in Table 1. As we can see, keeping three terms in the series approximating the resolvent kernel is enough for achieving engineering accuracy. The approximation towards the exact solution depends on the character of the eigenvalues. For the case of complex-conjugated eigenvalues (material 1), the succession is single sided, while it is flipping from side to side for the case of real eigenvalues (material 2).

Consider the application of the developed solution technique for the computation of thermal stresses in orthotropic half-planes in the case of plane strain (bottom lines under the braces in expressions (4)) for the materials with properties given in

Table 2. Elastic and thermoelastic material moduli of some orthotropic composites.

	Material	${\mathit{E}}_{\mathit{x}}, \mathbf{GPa}$	${\mathit{E}}_{\mathit{y}}, \mathbf{GPa}$	${\mathit{\nu }}_{\mathit{y}\mathit{x}}$	${\mathit{\nu }}_{\mathit{z}\mathit{x}}$	${\mathit{\nu }}_{\mathit{z}\mathit{y}}$	${\mathit{G}}_{\mathit{x}\mathit{y}}, \mathbf{GPa}$	${\mathit{\alpha }}_{\mathit{x}}\text{×} {10}^5, 1/\mathbf{K}$	${\mathit{\alpha }}_{\mathit{y}}\text{×} {10}^5, 1/\mathbf{K}$	${\mathit{\alpha }}_{\mathit{z}}\text{×} {10}^5, 1/\mathbf{K}$
3	Glass-Epoxy Composite [39]	13.7	55.9	0.277	0.4	0.68	5.6	4.0	1.0	4.0
4	Boron-Epoxy Fiber-Reinforced Composite [40]	2.068	20.68	0.25	0.25	0.25	1.034	5.1	0.84	5.1

. The analysis of the case of plane strain involves a greater number of material moduli compared to the plane-stress case considered above.

Assuming the heat conduction moduli to be the same in all the spatial directions, let the thermal stresses be induced by the steady-state temperature field, with the Fourier mapping (11) given in the form

$$\overline{T}(y)=-s^2\frac{{\tau }_0\sqrt{\pi }}{4}\exp \left(-|s|y-\frac{s^2}{4}\right),\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{\tau }_0=\mathrm{const}$$

(26)

This temperature field was derived and discussed in detail in [32], emphasizing it fits the correctness conditions for the solutions to stationary thermoelasticity problems in half-planes [6,32]. For this case study, the boundary of the half-plane is assumed to be free of external force loadings, i.e., $p(x)=q(x)=0$ in (5).

The effect of material anisotropy in the thermal stresses is illustrated in Figure 4 depicting the full-field distributions of the stress ${\sigma }_{yy}$ for materials 3 and 4 (see Table 2). Note that eigenvalues (15) for both materials are real and dissimilar. We can observe the similarity in the character of stress distribution for the considered orthotropic materials. However, the magnitude of thermal stresses is significantly different.

The convergence of solutions when implementing the approximate formula for the resolvent kernel (22) instead of the exact one for the computation of thermal stress (20) is demonstrated in Figure 5. Similar to the preceding case studies for the half-plane subject to normal force loadings, we can observe good convergence of the approximate solution to the exact one, i.e., keeping three terms in the series presenting the resolvent kernel is sufficient for achieving the desired accuracy. For material 4 (Figure 5b), the convergence is slightly slower.

6. Conclusions

Advanced design and processing of machinery, equipment, and structures elements are associated with the development of technologies for creating composite materials with properties optimized by specific target criteria. This encourages complex studies of physicochemical processes and thermal stresses affecting the strength and durability of materials and structures under the complex action of the dominant types of loading, such as the force and temperature fields, taking into account a wide range of operational and structural features.

A technique for solving the plane elasticity and thermoelasticity problems for an orthotropic homogeneous half-plane is developed based on the direct integration method. The key features of the presented technique are the following.

• This technique allows for constructing an analytical solution in the form that is unified for all possible interdependencies between the elastic moduli of an orthotropic material.
• For all possible cases of orthotropic materials, the solution fits the Saint-Venant principle, i.e., decreases when receding the loaded zone of the half-plane boundary.
• The solution is constructed in the form of explicit expressions of the stress tensor components on the applied force and thermal loadings.

These features make the constructed solutions to be a perfect candidate for implementation with the engineering analysis within the inverse problems, optimization of the stress patterns, identification of effective material parameters, etc.

As a certain disadvantage of the proposed solution scheme, it can be noted that the practical computation of the solution involves approximated formula (22) for the resolvent kernel. Nevertheless, the considered case studies have shown rapid convergence of the approximate solution to the exact one.