A new convolution variational boundary element technique for design sensitivity analysis and topology optimization of anisotropic thermo-poroelastic structures

Abstract The main objective of this paper is to propose a new efficient time domain boundary element technique for design sensitivity analysis and topology optimization of anisotropic thermo-poroelastic structures. The spatial regularization scheme based on integration by parts is performed in Laplace domain in conjunction with time-domain convolution variational boundary element method (CVBEM). The generating optimization problem is solved using the Moving Asymptotes Method (MAM) with the adjoint variable method to obtain optimal material distribution, influence of anisotropy on the thermo-poroelastic stresses sensitivities. The results from the numerical model demonstrate the validity, efficiency and accuracy of the proposed technique.

Understanding the effects and behavior of the elastic stress wave propagation is very important in the design of many engineering optimization problems.Analytical solutions are very difficult or even impossible to achieve.Therefore, the problem can be solved approximately by numerical methods which has been developed and successfully implemented the boundary element method to solve a wide range of general elastic problems.The boundary element method (BEM) ( Dom ınguez, 1993;Fahmy, 2011Fahmy, , 2012aFahmy, , 2012bFahmy, , 2012cFahmy, , 2012dFahmy, , 2012e, 2013aFahmy, , 2013bFahmy, , 2013cFahmy, , 2017aFahmy, , 2017bFahmy, , 2018aFahmy, , 2018bFahmy, , 2018cFahmy, , 2019aFahmy, , 2019bFahmy, , 2019cFahmy, , 2019d, 2019e ) , 2019e ) has been developed and successfully implemented to solve our current complex problem.Generally, Laplacedomain fundamental solutions are easier to obtain than time-domain fundamental solutions for engineering and scientific problems (Wang & Achenbach, 1994, 1995).consequently, the CVBEM is very useful for problems where time-domain fundamental solutions are not known, because it utilizes the Laplacedomain fundamental solutions of the governing equations of the problem.So, CVBEM expands the range of engineering problems that can be solved with the classical time-domain BEM.
The aim of this article is to develop a novel methodology for solving thermo-poroelastic topology optimization problems.The spatial regularization procedure based on integration by parts is performed in Laplace domain in conjunction with the time-domain convolution variational boundary element method (CVBEM).The spatial discretization is performed using the collocation scheme.An implicit differentiation method has been used to compute the temperature, displacements and pore pressure sensitivities with respect to the control variables.The adjoint variable method has been used to calculate sensitivities for objective function and constraints to allow for Moving Asymptotes Method (MAM) to solve the resulting topology optimization problem and find the optimal material distribution and influence of anisotropy on the optimal design.
The effects of anisotropy on the thermo-poroelastic stresses sensitivities are very pronounced.The optimum design of an elliptical sandwich structure has been considered to demonstrate the validity, efficiency and accuracy of the proposed technique.

Formulation of the problem
First, we consider the region g with boundary C as shown in Figure 1.
According to Biot (1956aBiot ( , 1956b) ) and Schanz (2001), the governing equations of anisotropic thermo-poroelastic structures can be described as According to Fourier's law, the heat flux ɋ p is defined by where r is the total stress tensor, e is the linear strain tensor, C ajlg is the constant elastic moduli, q ¼ q s 1 À / ð Þþ /q f is the bulk density, q s is the solid density, q f is the fluid density, u is the displacement field in the solid, u f is the displacement field in the fluid, h is the temperature, B are stress-temperature coefficients, q is the fluid-specific flux, F are the total body forces, tr denotes the trace, A ¼ / 1 þ Q=R ð Þis Biot's effective stress coefficient, Q and R are the solid-fluid coupling parameters, p is the pore pressure, f is the fluid volume variation measured in unit reference volume, , 1985).
The equations of motion which describe our problem can be written in matrix form as follows (Bonnet, 1987;Messner & Schanz, 2011) where the operator Bx and the tractions tg are defined as where B e x is the linear elastostatics operator and T e x is the traction derivative with respect to x: The thermo-poroelastic medium occupies the region X that bounded by a closed curve C that is subdivided into two non-intersective parts C D and C N : The

Boundary element implementation for the temperature field
Using the boundary element method as described in Fahmy (2012a), we obtain The source term in (11) has been approximated using a series f q and coefficients a q as ð Thus, the representation formula can be expressed as The temperature T and the heat flux q and their corresponding particular solutions T q and q q can be approximated as where h, ɋ, h ɋ , ɋ ɋ and U are matrices.By implementing point collocation procedure to (13) and using ( 14), we obtain the following system Using ( 16) into (15) we have By implementing a point collocation procedure to Eq. ( 16) as in Fahmy (2012a), we get From Eq. ( 18) the following expression can be derived Which can be substituted into (17) producing where The nodal vectors are subdivided into known k and unknown u superscripts Thus, from ( 20) we obtain the following matrix equation Making use of ( 23), the unknown fluxes q u ðtÞ can be written as 24), the second row of (23) can be expressed as follows where By using the finite difference method, we can write (25) as where 4. BEM implementation for traction and displacement fields According to problem (8), the representation formula of the unknown field ûg can be described as where Based on the anisotropic fundamental solutions of Wang andAchenbach (1994, 1995), we assumed that the fundamental solution Û r ð Þ and the corresponding traction T y can be written in Laplace domain as follows (Messner & Schanz, 2011).
where the solid displacement fundamental solution Ûs r ð Þ can be expressed as with Eq. ( 31) can be written as (Schanz, 2001) According to the regularization procedure in this article, the solid displacement fundamental solution Ûs r ð Þ can be divided into a singular part where The remaining parts of the fundamental solution Û r ð Þ can be calculated as (Schanz, 2001) Now, using Eqs.( 28) and ( 29), we obtain (Steinbach, 2008) lim where and the double layer operator can be written as By using Eqs.( 39)-( 42), the boundary integral equation can be expressed as By implementing inverse Laplace transformation, the boundary integral equation can be obtained as According to Messner and Schanz (2011), the poroelastodynamic fundamental solution can be expressed as follows (45) By using Stokes theorem, the differentiable vector aðyÞ can be written as which can be expressed as follows ð By using (47) we can derive the following formula (Messner & Schanz, 2011) By implementing Eq. ( 48) to a vector a ¼ vu we obtain (Kielhorn, 2009) By using ( 34) and ( 45), we can derive T s in the following form According to Martins and Hwang (2013), we can write T e y as follows which can be expressed using the Eq. ( 34) as follows By implementing (29) to (53), we have By applying ( 49) and ( 50) to (54), and using the Gaussian quadrature rule with the Duffy transformation, we get (Wang & Achenbach, 1994, 1995) K which is an equivalent form of (29) under closed boundary assumption.
The second term of the right side of (55) can therefore be manipulated to obtain where According to Fahmy (2019a), we can write the following limit lim with the regularized double layer kernel function (46).By augmenting Ûs s to Ûs and using (50) we can write (55) as follows According to the time discretization with discrete times t n ¼ nDt Dt > 0 ð Þ, the convolution integral can be written as Thus, we can write By using Cauchy's integral formula, we can calculate the integration weights x n which introduced by Lubich (1988aLubich ( , 1988b) ) as follows By using polar coordinate transformation z ¼ Re Àiu , the integral Eq. ( 62) can be approximated as By using ( 63), we can write (61) in the following form According to the procedure of Kielhorn (2009), we get which can be expressed in Laplace domain as follows Let the boundary C ¼ oX is discretized into N e surface triangles boundary elements s e as (Figure 2(a)) We can now describe the two subspaces on C h as where the unknown Neumann datum and unknown Dirichlet datum are approximated as where u a i k ½ are continuous polynomial shape functions, w b j k ½ are piecewise discontinuous polynomial shape functions and k ¼ 1, 2, 3, 4 are the poroelastic degrees of freedom.Substituting ( 71) and ( 72) into (67), we obtain where According to Maier, Diligenti, and Carini (1991) and considering matrices instead of integral operators, we can write (73) in compact form as follows where According to Maier, Diligenti, and Carini (1990) with consideration of f ¼ Ly and taking into account (61), thus, it is easily to prove that Ly, y 0 ¼ Ly 0 , y for any y and y 0 (76) Due to the symmetry of L in space and time, we can write the following variational formula (see Maier  , 1990, 1991) where x is a field point, n is a load point and Now, we evaluate the consequent variation of the functional (77) due to variation dy of the variables By using ( 75) and ( 76), we can write for any dy Thus, y is the boundary solution

Design sensitivity analysis and topology optimization
The considered optimization problem can be defined as follows (Matsumoto, Tanaka, & Ogawa, 2003) min where r is the residual of state equation, N is the design variables number, % c is a fraction of solid volume.Also, we assumed that q 1 , q 2 ð Þ ¼ 0:01, 0:99 ð Þ : Eq. ( 67) can be written as The objective function may be defined as The augmented objective function can be described as (Igumnov & Markov, 2016;Schanz, 2001) where A and B are nonsymmetric matrices and k is the adjoint variable vector According to Schanz (2001), the sensitivities can be calculated as By implementing the Moving Asymptotes Method (MAM) with the adjoint variable method (Itz a, Viveros, & Parra, 2016) a standard filter densities q e and the sensitivities for the BEM analysis can be computed as follows (Bourdin, 2001;Bruns & Tortorelli, 2001) Where w is the weighting function, x e are the element coordinates, N e denotes the elements set whose distance from e is lower than R and t i is the volume of element.

Numerical results and discussion
The proposed technique in the current work should be applicable to any anisotropic thermo-poroelastic problem.The models of Braun, Ghabezloo, Delage, Sulem, and Conil (2018) and Giot, Granet, Faivre, Massoussi, and Huang (2018) may be considered as special cases of our current complex and general study.To illustrate the numerical results calculated by the proposed technique presented in the current paper, the material considered as an anisotropic poroelastic material is a prismatic solid, and the physical parameters for this material are (Braun et al., 2018): The elasticity tensor of anisotropic prismatic structure where / ¼ 0:15, k ¼ 0:987 Â 10 À15 m 2 , q s ¼ 1600 kg=m 3 , q F ¼ 1113 kg=m 3 , p ¼ 25 MPa, Q=R ¼ 0:65, C ¼ 0:66, l ¼ 0:001, and Dt ¼ 0:0007 s: The conditions have been considered in the calculations are as follows Figure 2(b) shows the considered structure geometry.Table 1 shows the structure optimal design using CVBEM considered in this paper.that the anisotropy has a signicant effect on the Thermo-poroelastic stresses sensitivities.Also, the results are highly sensitive to assumed initial conditions.The displacement and pore pressure variations with the time are plotted in Figures 9 and 10 for the analytical, FEM and CVBEM methods to demonstrate the validity of our implemented technique.It is noticed from these figures that the results of the CVBEM are in very excellent agreement with those obtained using the analytical method of Braun et al. (2018) and FEM of Giot et al. (2018).Our results thus confirm that our technique is efficient and precise.

Conclusion
The CVBEM technique is proposed for transient studies of thermo-poroelastic optimization problems, where the regularization can be performed using integration by parts to remove the strong singularity and improve the efficiency of our proposed method.
The CVBEM is very useful for problems where timedomain fundamental solutions are not known, because it utilizes the Laplace-domain fundamental solutions of the governing equations of the problem.Also, our proposed technique reduces the computational time and memory.For these mentioned reasons, our proposed CVBEM is very effective, accurate and powerful technique.The accuracy, efficiency and validity of our proposed CVBEM technique were confirmed by comparing the obtained results of our considered problem with the corresponding analytical results of Braun et al. (2018) and with the FEM results of Giot et al. (2018), it can be noticed that the CVBEM results show excellent agreement with the analytical and numerical FEM results, our results thus confirm that the proposed CVBEM technique is very effective, accurate and powerful technique for anisotropic thermo-poroelastic optimization problems.Understanding the design sensitivity and optimization of anisotropic thermo-poroelastic structures are very important in many fields such as geothermal engineering, petroleum engineering, geomechanics engineering, geotechnical engineering, biomechanics engineering, material science and material engineering etc.

Figure 1 .
Figure 1.Geometry of the fundamental solution integration.

Figure 3 .
Figure 3. Variation of the r 11 sensitivity with time t:

Figure 10 .
Figure 10.Variation of the pore pressure p with time t:

Table 1 .
Investigation of the influence of anisotropy in an elliptical sandwich structure.

Table . 2
. Comparison of computer resources needed for FEM and CVBEM modelling of elliptical sandwich structure.