3D thermo-mechanical multifield problems

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. xxx (xxxx) xxx www.elsevier.com/locate/cma

Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems Ehsan Rabizadehb , Amir Saboor Bagherzadehc , Cosmin Anitescub , Naif Alajlana , Timon Rabczuka ,∗ a Department of Computer Engineering, King Saud University, Riyadh, Saudi Arabia Institute of Structural Mechanics (ISM), Bauhaus-Universität Weimar, Marienstraße 15, D-99423 Weimar, Germany c Department of Applied Mathematics, Faculty of Mathematical Sciences, Ferdowsi University of Mashhad, Mashhad, Iran b

Received 31 January 2019; received in revised form 21 September 2019; accepted 29 September 2019 Available online xxxx

Abstract In this article, we present a 3D adaptive method for thermoelastic problems based on goal-oriented error estimation where the error is measured with respect to a pointwise quantity of interest. We developed a method for a posteriori error estimation and mesh adaptation based on dual weighted residual (DWR) method relying on the duality principles and consisting of an adjoint problem solution. Here, we consider the application of the derived estimator and mesh refinement to two-/three dimensional (2D/3D) thermo-mechanical multifield problems. In this study, the goal is considered to be given by singular pointwise functions, such as the point value or point value derivative at a specific point of interest (PoI). An adaptive algorithm has been adopted to refine the mesh to minimize the goal in the quantity of interest. The mesh adaptivity procedure based on the DWR method is performed by adaptive local h-refinement/coarsening with allowed hanging nodes. According to the proposed DWR method, the error contribution of each element is evaluated. In the refinement process, the contribution of each element to the goal error is considered as the mesh refinement criterion. In this study, we substantiate the accuracy and performance of this method by several numerical examples with available analytical solutions. Here, 2D and 3D problems under thermo-mechanical loadings are considered as benchmark problems. To show how accurately the derived estimator captures the exact error in the evaluation of the pointwise quantity of interest, in all examples, considering the analytical solutions, the goal error effectivity index as a standard measure of the quality of an estimator is calculated. Moreover, in order to demonstrate the efficiency of the proposed method and show the optimal behavior of the employed refinement method, the results of different conventional error estimators and refinement techniques (e.g., global uniform refinement, Kelly, and weighted Kelly techniques) are used for comparison. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Thermoelastic multifield problems; Finite element method (FEM); Dual weighted residual (DWR) based error estimation; Goal-oriented error estimation (GOEE); Pointwise quantity of interest (QoI); Adaptive mesh refinement (AMR)

∗ Corresponding author.

E-mail address: [email protected] (T. Rabczuk). https://doi.org/10.1016/j.cma.2019.112666 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

2

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Nomenclature As long as no other meaning is explicitly given to a certain quantity within the text, its meaning corresponds to the following notation list. Acronyms/Abbreviations AISI BC DIN DOF PoI TOL UNS W-Kelly

American iron and steel institute Boundary condition Deutsches Institut f¨ur Normung (German institute for standardization) Degree of freedom Point of interest Tolerance Unified numbering system Weighted Kelly

General notations a, A a, A a, A a, A

Scalars (with regular lightface italic font) Vectors or matrices (with boldface upright font) Second-order tensors (with boldface italic font) Fourth-order tensors (with boldface upright sans-serif font)

Roman letters A a(•)(•) Bu Bθ C0 Ck C d du dθ dim D D D0 (e) eu eθ eex χ eχ e∗~ e∗~ ex F Fb {F}(e)

System of governing equations Bilinear form of A Displacement-strain interpolation matrix; Displacement derivative matrix Temperature-gradient interpolation matrix; Temperature derivative matrix Space of continuous functions. Space of functions with continuous derivatives up to order k Stiffness or elasticity tensor Number of spatial dimensions Degrees of freedom of displacement field at each node Degrees of freedom of temperature field at each node Problem spatial dimension; Number of spatial dimensions Elasticity matrix [Vector-valued] Displacement trial function space [Vector-valued] Displacement test function space Element of “e” Exact error of the displacement field in the primal problem Exact error of the temperature change field in the primal problem Exact discretization error of the primal problem Exact discretization error of the primal problem Error of the dual problem; Exact error of the dual problem Exact error of the dual problem Scalar field Body force vector per unit volume Global mechanical-thermal force vector of element (e)

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

{F}h G(•, •) Hk Ih I I Isym J J (•) Jε (•) [K](e) [K]h l(•) L 2 (Ω ) L 2 (∂Ω ) n (e) n Ne Ni Nn Nθ Nu NV N∗ Pin Pout Ph q q¯n Qv rm ru, rθ Ru , Rθ R t t tn ¯tn T T0 T¯ T T0 u u¯ U

3

Global mechanical-thermal force vector associated with FEM Green’s function Hilbert space of order k Point interpolant of dual solution from dual FE space into primal FE space Second order identity tensor Fourth order identity tensor Symmetric fourth order identity tensor Discrete differential operator Goal functional Regularized form of the functional J (•) Stiffness matrix of element (e) Global stiffness matrix associated with FEM Linear form (functional) of A Lebesgue spaces of square-integrable functions on the domain Ω Lebesgue spaces of square-integrable functions on the boundary ∂Ω Number of nodes of element (e) Outward-pointing boundary normal unit vector Number of the elements Element basis (shape) function Lagrangian shape function associated with node n Shape function vector of the temperature change scalar field Shape function matrix of the displacement vector field Dimension of the function space V h = dim(V h ) Set of natural numbers without zero = Z+ = {1, 2, . . .} Internal pressure External pressure Finite element partition Heat flux vector Prescribed heat flux boundary condition Volumetric heat source Mean radius Edge residuals Cell residuals Set of real numbers Wall thickness Surface traction vector Normal traction vector on boundary with normal unit vector of n Prescribed traction boundary condition vector Absolute temperature Reference temperature Prescribed temperature boundary condition [Scalar-valued] Temperature change trial function space [Scalar-valued] Temperature change test function space Displacement vector Prescribed displacement boundary condition vector Nodal values of the displacement vector field u

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

4

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

V V V0 V′ W k, p x Z≥0

Mixed (product) function space for displacement and temperature change Trial (ansatz) functional space Test functional space Dual vector space of V Sobolev functional space of order k and integrability p Coordinate vector; Point coordinates Set of non-negative integers = {0, 1, 2, . . .}

Greek letters α αT β χ χ

δi j {∆}h {∆}(e) ϵ ϵM ϵθ εJ ε ex J ε es J ε ηJ η J% Γ γ Γe ΓD ΓN ~ˆ h κ κ λ, µ Ω Ω ∂Ω Ωe ϕ φj ρ(χ h )(•) ψ ψj σ θ

Thermal expansion tensor Linear thermal expansion coefficient Stress temperature tensor Thermo-mechanical pair; displacement-temperature change pair := (u, θ) Trial function of the primal problem Kronecker delta Global vector of unknowns associated with FEM Nodal values vector of element (e); Vector of unknowns of element (e) Strain tensor Elastic or mechanical strain tensor Thermal strain tensor Goal error; Exact goal error Exact goal error Estimated goal error Error tolerance Relative error in the quantity of interest J Percent relative error in the quantity of interest J Boundary Stress temperature modulus Boundary of element (e) Dirichlet boundary Neumann boundary Numerical dual solution calculated with higher accuracy than primal solution Thermal conduction coefficient; Coefficient of thermal conductivity Second-order thermal conductivity tensor Lam´e constants Problem domain; Lipschitz domain; Simulation domain without boundary Closure of the domain Ω ; Simulation domain Smooth boundary of domain Ω ; Lipschitz boundary of domain Ω Area of element (e) Temperature change test function Global basis function for the temperature finite-dimensional space T h Residual functional; Primal residual Displacement test function Global basis function for the displacement finite-dimensional space D h Cauchy stress tensor Temperature change

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

5

Nodal values of the temperature change field θ Test function of the primal problem Thermo-mechanical test function pair := (ψ, ϕ)

Θ ϑ ϑ Operators : % \ \ J• K ⌊•⌋ ∂i

Frobenius inner product or double contraction operator “Modulo”/remainder operator Complement or set difference operator Integer division operator Jump operator; Jump across the inter-element edges/face Floor function Partial derivative with respect to xi := ∂∂xi

∂x

Partial derivative with respect to xi := ∁ T

,i

mod tr(•) ∇ T (•) ∇ s (•) ∇ θs ∇ us ∇ ∇(•) ∇· (•) ∇×(•) A ∆ dim(•) ⟨•, •⟩ Hom(•, •) ∩ (•, •)Ω (•, •)∂Ω (•, •) (•, •)Γi ∑ span {. . . } ⊗ ∪ ⨄

∂ ∂ xi

[Set] Complement operator Transpose operator Partial derivative with respect to x := ∂∂xi Mod operator; “Modulo”/remainder operator Trace of (•) Transpose gradient operator := (∇(•))T Symmetric gradient operator := 21 (∇ ⊗ (•) + (•) ⊗ ∇) Temperature symmetric gradient operator Displacement symmetric gradient operator Del or nabla operator, denoting spatial derivation Gradient of (•) Divergence of (•) Curl of (•) Assembly operator Laplacian operator := ∇ · ∇ = ∇ 2 Dimension of the function space • Duality pairing Internal Hom functor Set intersection operator Standard L 2 (Ω ) scalar product Standard L 2 (∂Ω ) scalar product Inner product defined on Ω := (•, •)Ω Inner product defined on Γi ⊂ ∂Ω Discrete summation operator Linear span Tensorial product Set union operator Disjoint union operator

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

6

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Sub- and Superscripts h(d)

D N

d indicates the approximation degree of the numerical solution

h

Dirichlet subscript Neumann subscript Continuity order of FEM solution; Element continuity order

C (e) h h

i ij i jkl q t θ u

x m

The

associated with the element (e)

Finite-dimensional subspace of the infinite-dimensional function space Finite element approximation of the exact solution Vector component in xi direction i jth component of the second-order tensor i jklth component of the fourth-order tensor Heat flux superscript Traction superscript Temperature superscript Displacement superscript Vector component in x direction Tensor order superscript; tensor of order m

Miscellaneous symbols := ≈ ≡ O ̸∈ ̸⊆ ⊥a ∀ ∅ •| Γ ∈ ⊂

Definition sign Approximation sign; Approximately equal Equivalence sign Order of approximation sign; Order of accuracy No set membership sign; Not member of Not subset sign; Not subset of Orthogonality sign (with respect to bilinear form a) “For all” sign Empty set symbol Evaluation sign denoting • evaluated at Γ Set membership sign meaning being an element of a set; Member of Subset sign; Subset of

1. Introduction In recent years, substantial attention has been devoted to multifield problems and their numerical analysis. Many realistic engineering problems are categorized as multifield problems which involve the interaction among diverse fields. One important category of multifield problems is thermoelasticity, where in addition to the mechanical field, thermal field and their effects on each other are considered. Thermoelasticity theory which relates two separate and independent developed theories of elasticity and heat conduction, investigates the interaction between displacement and temperature fields [1,2]. In the classical theory of thermoelasticity, the Fourier heat diffusion equation, which describes the temperature distribution, is a parabolic partial differential equation and does not contain any displacement or strain terms. It is well known that the Fourier heat equation predicts an infinite propagation speed for heat waves [3]. The infinite speed assumption for the Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

7

thermal disturbances propagation by conventional heat conduction theory is appropriate for most of the engineering applications (some exceptions are cases with extremely short time responses like thermal shock problems or very low, near zero, temperatures [4]). In literature, considerable effort has been focused on the study of thermoelastic problems. In 1950, Mindlin and Cheng provided an investigation on steady-state thermal stresses in a half-space [5], then in 1957, a detailed analysis of steady-state problem for half-space was published by Sternberg and McDowell [6]. Publishing the first thermoelasticity book by Melan and Parkus in 1953, which is about the steady-state thermoelasticity [7], was the important development in this field. In 1957, Nowacki [8] investigated a half-space domain with discontinuous temperature boundary conditions on its surface and a few years later, in 1960, a half-space and a layer under steadystate thermoelastic condition has been considered by Sneddon and Lockett [9]. Further references to the previous studies in thermoelasticity can be found in [10,11] and in the books of Boley and Weiner [12], Nowacki [13,14] and Parkus [15]. More recent researches in the field of thermoelasticity can be found in [16–21]. Thermoelasticity has been one of the important branches of solid mechanics and numerous developments have been aimed towards analyzing them numerically. In general, computational errors are intrinsic in numerical analysis and it is not possible to avoid numerical errors in the results. Although, many investigations have been devoted to the numerical simulation of thermo-mechanical problems, fewer attempts have been directed toward analyzing the error and error estimation in thermoelastic applications. The accuracy and reliability of the solution is one of the main concerns of each numerical analysis, and different error estimators have been developed to evaluate the errors. Generally, error estimation techniques can be classified as a priori or a posteriori methods [22]. In this article, the a posteriori error estimation in which the results of analysis are used to assess the accuracy of the solution, is considered. Babuˇska and Rheinboldt, who published a paper on finite element error estimation for two point elliptic boundary value problem in 1978 [23], were the first researchers who investigated a posteriori error estimation method in finite element applications. They also had a substantial contributions to the adaptive mesh refinement by relating the error to each element. For a review of the different a posteriori error estimation techniques we refer to Gr¨atsch and Bathe [24], Ainsworth and Oden [22], or Verf¨urth [25]. Traditionally, according to the procedure used to obtain the estimates, the a posteriori error estimation techniques can be categorized into two major branches of recovery based and residual based error estimators. The recovery type a posteriori error estimators assess the error based on the difference between the (discontinuous) gradient field of the numerical results and the relevant recovered gradient field obtained by application of a recovery operator [26]. In 1987, Zienkiewicz and Zhu introduced a simple recovery type error estimation (known as Z -Z or Z 2 error estimator), which was based on the comparison of smoothed solution gradients with the gradients of the solution themselves [27]. Later, in 1992, Zienkiewicz and Zhu modified and developed their error estimator by introducing so called superconvergent patch recovery (SPR) method [28,29]. In order to improve the error estimation and adapt the method for different applications, several modifications have been proposed to the conventional SPR method (e.g., SPR-EB by Wiberg et al. [30] and also by Blacker and Belytschko [31] in 1994; SPR-C [32], SPRXFEM [33], and SPR-CX [34] by R´odenas et al. in 2007, 2008, and 2010, respectively; WSPR by Moslemi and Khoei [35] in 2009). The other category of a posteriori error estimation, which is based on a residual quantity (e.g., equation residual, (discontinuous) derivative residual, material constitutive law residual) [26], was originally introduced by Babuˇska and Rheinboldt [23]. Further development and application of the element residual method in error estimation were performed by Demkowicz et al. [36,37], Bank et al. [38,39], Ainsworth and Oden [40,41]. Later, in 1989, an extensive study of error residual methods is reported by Oden et al. [42]. For further detail, the reader is referred to [43–49] and references therein. In conventional global error estimates, usually the errors are estimated in energy norm, which gives a general presentation on the overall accuracy of an approximate solution. However, in many real-life numerical simulations, such an information is not sufficient and from the viewpoint of engineering purposes local errors or error in particular quantities of interest are of main interest. The quantities of interest or goal is a functional which could be presented in terms of primal solution over certain subdomains, lines or at specific points of the solution domain. In order to achieve this objective, recently, a new branch in a posteriori error estimation theory has emerged as goal-oriented error estimation (GOEE). The GOEE approach, which is usually based on the use of duality techniques, estimates the error in quantities of interest (QoI) directly by relating the goal error to the solutions of primal and adjoint Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

8

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

problems [50–54]. In recent years, recovery based GOEE and adaptivity techniques have been implemented in thermoelasticity problems by Rabizadeh et al. [55,56]. One of the most well-known techniques to derive GOEE is given by dual weighted residual (DWR) method which relies on duality argument. So-called DWR method, which is based on a variational formulation of the problem, derives the weighted a posteriori error estimates with respect to goal error functional by employing a suitable adaptation of global “duality arguments” (known as “Aubin–Nitsche-Trick”). DWR error estimators consist of residual evaluations weighted by adjoint sensitivity measures. The solution of the global dual (auxiliary) problem which provides the weight functions, is a measure for the sensitivity of primal solution with respect to the goal functional [57–60]. Some of the other adaptivity techniques can be referred as [61–63]. In numerical calculations, after estimating the error, in order to reduce the error, a mesh refinement procedure is performed. In engineering applications, three dimensional (3D) simulation of thermoelastic materials plays a significant role. Generally, in 3D problems, the computational cost increases enormously by global uniform mesh refinement. In order to overcome this problem, mesh adaptivity techniques have been introduced and developed. There are different adaptive mesh refinement (AMR) techniques, but in the majority of them the method tries to adapt the accuracy of the solution within certain regions with higher error contribution according to an error indicator. The classical adaptive refinement algorithms control the error measured in energy norm. Although, controlling the global error also reduces the goal error, traditional mesh adaptivity techniques do not necessarily optimally control the local errors to enhance the accuracy in quantity of interest approximation. In order to optimize the refinement methods for controlling the errors measured in quantity of interest, goal-oriented AMR techniques have been developed. Unlike the classical adaptive finite element formulation which consider the global energy norm as error control criteria, the goal-oriented associated adaptivity techniques are based on the concept of controlling error in terms of special goal-oriented criteria. The AMR strategies can be categorized into three main classes of h-, p-, and r -adaptivity methods. h-adaptive mesh refinement/coarsening refers to local refinement and/or coarsening of the mesh. In p-refinement method which results in p-nonconformity in mesh, the solution is enhanced by locally varying the polynomial degree of the shape functions. Although in h and p adaptive methods the initial structure of the mesh is unchanged, r -adaptivity increases the accuracy by relocating or moving the elements which can be interpreted as a re-meshing technique. In engineering applications, the h-refinement is the most popular. In this contribution, an extension of the dual weighted residual (DWR) method to the analysis of three dimensional thermoelastic problems is presented. The underlying theoretical framework is that of the DWR method goal-oriented error estimation for finite element discretization of general variational equations, developed by Becker and Rannacher in [57,59]. Here we consider its application to thermoelastic problems. In order to reduce the error in the quantity of interest efficiently and optimally in 2D and 3D problems, we perform goal-oriented adaptive refinement schemes to target the error in the quantity of interest directly. The 3D thermoelastic FEM simulation and the corresponding dual weighted residual error estimation with the concomitant goal-oriented mesh adaptivity are accomplished in C++ programming language within the framework of the deal.II FE library [64,65]. In the numerical results, it is substantiated that application of dual weighted residual based h-adaptivity procedure results in the considerable reduction in the computational costs. Overall, we conclude that the employed DWR approach is particularly designed to achieve high solution accuracy at a minimal computational cost. In this study, industrial cases of thick-walled cylindrical and spherical pressure vessels are considered as the benchmark problems and implementation of the developed method investigated. The analyses of thick-walled pressure vessels subjected to internal pressure and temperature change are important in solid mechanics and engineering applications. Spherical and cylindrical pressure vessels are widely used in petrochemical industries to store the cold pressurized gas or liquidized gas like propane, butane, etc. The spherical storage vessels contain more volume per unit surface area and the stress resistance is uniform. The smaller surface area per unit volume of the spherical vessels in comparison with other vessel shapes, which minimizes the heat transfer between surroundings and the inside medium, makes spherical vessels the optimal shape for storing cold pressurized gas in petrochemical industries. The remainder of this paper is organized as follows. In Section 2, the system of governing equations of the thermoelastic problem is set up. Section 3 focuses on the variational formulation of the problem and accordingly extracting the relevant bilinear and linear forms. The finite element discretization is discussed in Section 4. Section 5 is devoted to the main objective of this article. This section, which presents the goal-oriented error estimation, Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

9

consists of two subsections. The first part is dedicated to formulating the dual weighted residual method for the thermoelastic problems. Subsequently, in Section 5.2, the pointwise quantities of interests are discussed. The numerical results provided in Section 6 demonstrate the efficiency, reliability, and optimality of the proposed method for four model problems with analytical solutions. Finally, concluding remarks are outlined in Section 7. The analytical solutions of the studied examples are presented in Appendices A.1, A.2. 2. Governing equations In the thermoelastic problems, the governing equations consist of equation of motion (Cauchy momentum equation) and energy balance equation. In this article, the studied thermoelastic problem is assumed to be in steady state condition. Considering this assumption, the classical (strong) formulation of the general governing equations for the linear thermo-mechanical problem ignoring the inertia and time derivative terms, on the domain Ω , reads as follows [3,66–69]: ∇· σ + Fb = 0,

in Ω

(1)

∇· q − Q v = 0,

in Ω

(2)

with boundary conditions (BCs) ¯ u|Γ u = u, := σ · n|Γ t = t¯n , T |Γ θ = T¯ ,

(3)

qn |Γ q = q · n|Γ q = q¯n ,

(6)

tn |Γ t

(4) (5)

where σ , Fb , q, Q v , u, t, and T are Cauchy stress tensor, body force vector per unit volume, heat flux vector, volumetric heat source, displacement vector, surface traction vector, and absolute temperature, respectively. Let consider domain Ω ⊂ Rd with d ∈ {2, 3} the number of spatial dimension, as the problem domain. We assume the domain Ω to be a simply connected, open, bounded set with sufficiently (piecewise) smooth boundary ∂Ω = Ω \Ω . Accordingly, Γ u , Γ t , Γ θ and Γ q are boundaries of the domain Ω with prescribed boundary conditions of “displacement”, “traction”, “temperature” and “heat flux”, referred by “u”, “t”, “θ”, and “q” superscripts, respectively, and n ∈ Rd is the outward-pointing unit vector normal to Γ [3,55]. The boundary conditions in terms of the known tractions on the boundary directly follow from Cauchy’s formula. The stress tensor σ , can be calculated from the stress–strain constitutive equation or the generalized Hook’s law. Therefore, for an elastic material under thermal effects we have σ = C : ϵ M = C : (ϵ − ϵ θ ) ( ) 1 (7) = C : (∇ s u − αθ ) = C : ∇u + ∇ Tu + βθ, 2 where C is the fourth-rank tensor of material properties, called the stiffness or elasticity tensor; α is the second order thermal expansion tensor, β = −C : α is the stress temperature tensor, and θ = T − T0 is the temperature change above the reference temperature (T0 ) at which strains and stresses are zero. ϵ denotes the strain tensor which is the symmetric gradient of the displacement vector, and consists of the elastic strain (ϵ M ) and thermal ( or Tmechanical ) strain (ϵ θ ) tensors. The operators “:”, “∇ s ”, “∇ T ” and the superscript “ ” are the Frobenius inner product 1 2 3 (double contraction), symmetric gradient, transpose gradient, and transpose operators, respectively. Eq. (7) is the general stress–strain constitutive relation which is valid for a heterogeneous anisotropic linear elastic solid continuum. In this study, the material is assumed to be homogeneous isotropic linear thermoelastic material.4 Consequently, by substituting the stress tensor into the equilibrium equation (1) yields the Navier equation (also known as the Navier–Cauchy equation) for a homogeneous isotropic linear thermoelastic material µ∆u + (λ + µ)∇∇ · u − α T (3λ + 2µ)∇T = −Fb , 1 2 3 4

and

(8)

( ) : 2 B ≡ tr A · B T = Ai j Bi j and 4 C : 2 D ≡ Ci jkl Dkl . 1 ∇ (•) := 2 (∇ ⊗ (•) + (•) ⊗ ∇). ∇ T(•) := (∇(•))T . For homogeneous isotropic linear elastic material C = λI ⊗ I + 2µIsym , α = α T I, and κ = κ I, where I and I are second order identity symmetric fourth order identity tensors, respectively. 2A s

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

10

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

where λ and µ are the Lam´e constants, α T is the linear thermal expansion coefficient, “∆” is the Laplacian (= ∇ · ∇ = ∇ 2 ) operator and γ = α T (3λ + 2µ) is known as the stress temperature modulus. The other governing equation of the thermoelastic problem is the heat balance equation (2). For heterogeneous anisotropic linear material, the classical Duhamel’s generalization of the heat transfer [70] leads to: q = −κ · ∇T,

(9)

where κ is the second order thermal conductivity tensor and the minus sign shows that the heat flow occurs from a higher to a lower temperature. By substituting (9) in the energy conservation Eq. (2) and considering a homogeneous isotropic linear thermoelastic solid continuum (see footnote 4), the heat conduction and balance equations reduce to q = −κ∇T,

(10)

κ∆T = −Q v ,

(11)

respectively; where Eq. (10) is known as the Fourier’s heat conduction law and κ is the thermal conduction coefficient. In the static thermoelastic problem, the equilibrium equation (8) is coupled but the heat transfer equation (11) depends only on the temperature and therefore can be solved separately (staggered approach). However, we opt to solve the equations simultaneously (monolithic algorithms). Therefore, the system of governing equations A for pairs χ := {u, T }, which models a stationary coupled thermoelasticity, is defined as { } µ∆u + (λ + µ)∇∇ · u − α T (3λ + 2µ)∇T + Fb A(χ ) := = 0, (+BCs). (12) κ∆T + Q v A incorporated with BCs determines the pair χ := {u, θ} of displacement vector u and scalar temperature change θ of a coupled linear thermoelastic problem with a homogeneous isotropic material. 3. Monolithic variational formulation Let us consider so-called Lipschitz domain Ω ⊂ Rd (d = 2, 3) with Lipschitz boundary ∂Ω = Ω \Ω , as the computational domain. The strong form of the governing equations in the compact monolithic form written as a linear system reads J χ = f, (+BCs),

in Ω ,

(13)

where the linear mapping J : V → V ′ is a discrete differential operator and f ∈ V ′ := Hom(V, F), with V ′ being the corresponding dual vector space of V, and F being a scalar field. The mixed (product) function space V for the displacement u and temperature change θ is defined as V := D × T , where d 1 1 H1u,Γ u (Ω ) = H u (Ω ) × · · · × H u (Ω ) ¯ u,Γ ¯ u,Γ ¯   

D(Ω ) :=

d times

= T (Ω ) :=

[

H1¯

θ ,Γ

]d

{ ⏐ } ⏐ = u ⏐ u ∈ H1 (Ω )d ; u|Γ u = u¯ , { ⏐ } ⏐ 1 ¯ θ (Ω ) = θ ⏐ θ ∈ H (Ω ); θ |Γ θ = θ ,

1 Hu,Γ u (Ω ) ¯

(14) (15)

with Hk (Ω )d being the Hilbert space of order k defined on Ω d . Fig. 1 represents schematic of an elastic body under thermo-mechanical loadings. In this figure, different Dirichlet and Neumann boundaries of both thermal and mechanical fields are depicted. Calculus of variations shows that the strong form (13) with mentioned boundary conditions, is equivalent to finding χ ∈ V such that for all ϑ ∈ V 0 := H10,Γ u (Ω )d × H1 θ (Ω ) holds: 0,Γ

[J χ ] (ϑ) = f (ϑ) ,

(16)

where ϑ is called test function and vanishes on the Dirichlet boundaries where the solution is already known. By defining the bilinear form, Eq. (16) leads to a more generic form of a weak formulation. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

11

Fig. 1. Schematic of an elastic body in thermo-mechanical multifield problems.

Accordingly, for the convenience, we introduce the following canonical notation for the natural variational (weak) formulation of A(χ ) based on the mixed (product) function space V for the displacement u and temperature change θ: Find χ ∈ V(Ω ) := D(Ω ) × T (Ω ), such that a (χ ) (ϑ) = l (ϑ) ∀ϑ := (ψ, ϕ) ∈ V 0 (Ω ), (17) ⏐ } ⏐ ¯ θ|Γ θ = θ¯ , (u, θ) ⏐ u ∈ W 1,2 (Ω )d , θ ∈ W 1,2 (Ω ); u|Γ u = u, ⏐ { } ⏐ V 0 (Ω ) :=D 0 (Ω ) × T0 (Ω ) = (ψ, ϕ) ⏐ ψ ∈ W 1,2 (Ω )d , ϕ ∈ W 1,2 (Ω ); ψ|Γ u = 0, ϕ|Γ θ = 0 , V(Ω ) := D(Ω ) × T (Ω ) =

{

(18) (19)

where a(•)(•) is the bilinear form and l(•) is the linear form (functional) of A; the spaces V and V 0 are trial (ansatz) and test functional spaces, respectively; W k, p (Ω )d is the Sobolev functional space of order k and integrability p defined on Ω d . The trial space V which satisfies the prescribed essential (Dirichlet) boundary conditions, contains members of V 0 shifted by the value of Dirichlet condition. Using the abstract framework, the bilinear and linear forms associated with this monolithic (mixed) variational formulation are a (χ ) (ϑ) :=[J χ ] (ϑ) = ⟨J χ , ϑ⟩, (20) l (ϑ) :=

f (ϑ) = ⟨f, ϑ⟩,

(21)

where the product ⟨•, •⟩ denotes the duality pairing (product) between linear functionals from the dual space V ′ of V and the space V 0 . In the following, (•, •)Ω and (•, •)∂Ω denote the standard L 2 (Ω ) and L 2 (∂Ω ) scalar products6 of functions, vectors, and tensors, respectively. Now we apply this abstract approach to the thermoelastic problems. Considering the Riesz representation theorem on Hilbert spaces and definitions (20) and (21), the bilinear form a and linear form l associated with the thermoelastic problem, are obtained by solving variational form (16) based on the system of governing equations (12) as 5

(J χ , ϑ)Ω = (f, ϑ)Ω , ∫ T χ ϑ J dΩ = ϑ T f dΩ ,

∫ Ω

∫ Ω

(22) (23)

Ω

(

ψ

ϕ

[ ) µ∆ + (λ + µ)∇∇· 0

−α T (3λ + 2µ)∇ κ∆

]{ } ∫ ( u ψ dΩ = − θ Ω

ϕ

{ } ) Fb dΩ . Qv

We may rewrite Eq. (24) as ∫ ∫ ( ) ( ) µψ · ∆u + (λ + µ)ψ · ∇∇·u − γ ψ · ∇θ + κϕ∇·∇θ dΩ = − ψ · Fb + ϕ Q v dΩ . Ω 5 6

(24)

(25)

Ω

∫ For simplicity ⟨•, •⟩ is used instead of ⟨•, •⟩V ′ ,V 0 , where ⟨ω, υ⟩V ′ ,V 0 ≡ (ω, υ) L 2 (Ω) := Ω ω · υ dΩ , ω ∈ V ′ , υ ∈ V 0 . 2 2 L (Ω ) and L (∂Ω ) are the Lebesgue spaces of square-integrable functions on the domain Ω and boundary ∂Ω , respectively.

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

12

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Performing integration by parts for the terms with second-order derivatives and making use of the divergence theorem (Gauss–Ostrogradsky theorem) we arrive at: ∫ ( ) µ∇ψ : ∇u + λ∇·ψ∇·u + µ∇ψ : ∇uT − γ ∇·ψθ + κ∇ϕ · ∇θ dΩ Ω ∫ ∫ ∫ ( ) = ψ · Fb + ϕ Q v dΩ + (ψ · σ ) · n dΓ − (ϕq) · n dΓ . (26) Ω

Γt

Γq

By introducing the Cauchy’s formula into the traction surface integral of Eq. (26) and substituting the Neumann boundary conditions (4) and (6), we may rewrite the right side as ∫ ( ) µ∇ψ : ∇u + λ∇·ψ∇·u + µ∇ψ : ∇uT − γ ∇·ψθ + κ∇ϕ · ∇θ dΩ Ω ∫ ∫ ∫ ( ) = ψ · Fb + ϕ Q v dΩ + ψ · ¯tn dΓ − ϕ q¯n dΓ . (27) Ω

Γt

Γq

Accordingly, the bilinear form a : V × V 0 → R which maps a pair of trial and test functions to a real number, and the functional l : V 0 → R a mapping from the test functions space to the real number space becomes a (χ ) (ϑ) = a ((u, θ)) ((ψ, ϕ)) := ( ) ( ) ( ) ( ) ( ) µ ∇u, ∇ψ Ω + λ ∇·u, ∇·ψ Ω + µ ∇uT , ∇ψ Ω − γ θ, ∇·ψ Ω + κ ∇θ, ∇ϕ Ω ,

(28)

l (ϑ) = l ((ψ, ϕ)) := ( ) ( ) ( ) ( ) Fb , ψ Ω + Q v , ϕ Ω + ¯tn , ψ Γ t − q¯n , ϕ Γ q , (29) ∫ ∫ where (•, •)Ωi := Ωi• · • dΩ and (•, •)Γi := Γi• · • dΓ denote the inner products defined on Ωi ⊂ Ω and Γi ⊂ ∂Ω , respectively, with “·” representing the scalar product between dual quantities (simple or double index saturation operation between vectors or tensors). Hereafter, we omit the subscript Ωi whenever the set Ωi is identical with the domain Ω , on which the system of differential equations is posed. 4. Finite element discretization As mentioned in Section 3, the governing system of equations (12) can be formulated as mixed variational form (17). As the function spaces V and V 0 are infinite dimensional, it is only possible to find the solution in special cases or in specially constructed problems. In order to overcome this problem, the infinite-dimensional (continuous) variational problem should be restricted to a finite-dimensional (discrete) variational problem. Therefore, in general cases, numerical methods are used to discretize the governing equations and obtain the approximate solutions. For discretizing the current thermoelastic problem, we use a standard Galerkin finite element method based on the trial space V and test space V 0 . The monolithic variational form (17) of the thermoelastic problem is the basis of the finite element Galerkin discretization. By applying Galerkin finite element method for seeking the solution within a finite dimensional subspace V h ⊂ V we may numerically calculate the best approximation χ h ∈ V h to the solution. By restricting the solution and variations to a finite dimensional subspace V h , the discretization of the problem can be easily obtained. Accordingly, the problem of finding the approximate solution reads: Find χ h ∈ V h := D h × T h , such that a(χ h )(ϑ h ) = l(ϑ h ) ∀ϑ h ∈ V 0h ,

(30)

{ } D h = span ψ1 , ψ2 , . . . , ψ NV ,

(31)

{ } T h = span φ1 , φ2 , . . . , φ NV ,

(32) N

V where D h and T h construct the same size7 finite dimensional spaces with bases denoted by Ψ := {ψ j (x)} j=1 and NV h Φ := {φ j (x)} j=1 , respectively so that NV = dim(V ).

7

NV = dim(V h ) = dim(D h ) = dim(T h ).

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

13

In finite element method, in order to define the finite element space V h , the domain Ω is partitioned into a finite Ne set of non-overlapping subdomains P h = {Ωi }i=1 with disjoint interiors,8 where Ne denotes the total umber of the h subdomains (elements) of the partition P . Once the domain Ω is partitioned into cells, it is possible to define local function space V eh on each element (e) to constitute the discretized global function space V h . In the so-called ansatz method, to approximate the solution pair χ := {u, θ} in element (e), we employ the ansatzs: ⏐ ⏐ ⏐ ⏐ (e) u(x)⏐ ≈ uh (x)⏐ = N(e) ∈ D eh , x ∈ Ωe , N(e) (33) u (x)U u : Ωe → R, (e) (e) ⏐ ⏐ ⏐ ⏐ (e) θ (x)⏐ ≈ θ h (x)⏐ = N(e) ∈ Teh , x ∈ Ωe , N(e) (34) θ (x)Θ θ : Ωe → R, (e)

(e)

(e)

(e)

where U and Θ are the nodal values vector of the displacement vector field u, and the temperature change scalar 9 field θ , respectively. Similarly, N(e) u is the element shape function matrix for the vector field of displacement and (e) Nθ is the element shape function (see footnote 9) vector for the scalar field of temperature change. Ωe indicates the area of the element (e). After partitioning the domain, one may build the global system of equations by substituting ansatzs (33) and (34) and relevant test functions in the discrete variational problem (30) in each element and assembling the element matrices in to the global matrices. Considering the fact that a(•)(•) and l(•) are linear, the Ritz–Galerkin10 discretization leads to a linear system of equations [K]h {∆}h = {F}h ,

(35)

where [K]h , {∆}h , and {F}h are the global stiffness matrix, vector of the unknowns, and mechanical-thermal force vector, respectively [71,72]. By assembly process, the global matrices can be obtained as Ne

Ne

[K]h = A [K](e) ,

{∆}h = A {∆}(e) ,

e=1

e=1

Ne

{F}h = A {F}(e) , e=1

(36)

where A is the assembly operator, [K](e) , {∆}(e) and {F}(e) are the stiffness matrix, the nodal values vector and the mechanical-thermal force vector of the element (e), respectively and are defined as ] [ (e) Kuu −K(e) (e) uθ , (37) [K] = 0 K(e) θθ { (e) } U {∆}(e) = , (38) Θ (e) { (e) } Fu {F}(e) = , (39) F(e) θ with [Kuu ](e) = [Kuθ ](e) =

{Fθ }(e)

T

B(e) D(e) B(e) u u dΩe ,

∫Ωe ∫Ωe

T

B(e) γ (e) N(e) u θ dΩe , T

B(e) κ (e) B(e) θ θ dΩe , Ωe ∫ ∫ T (e) T (e) ¯tn dΓe , = N(e) F dΩ + N(e) e u u b Γe Γe ∫ ∫ T (e) (e) T (e) = Nθ Q v dΓe − N(e) q¯n dΓe , θ

[Kθθ ](e) = {Fu }(e)

∫

Γe

N 8 ⨄e Ω

i

= Ω,

and Ωm

(40) (41) (42) (43) (44)

Γe

⋂

Ωn = ∅ for m ̸= n, where Ω and Ω are the simulation domain, and the simulation domain without the

i=1

boundary, respectively. 9 Element basis function. 10 In the standard Ritz–Galerkin method, basis functions for both ansatzs and test functions are assumed to be the same. In this case, the uh and θ h approximations are called the Ritz projections of u and θ, respectively. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

14

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

where U(e) , Θ (e) , D(e) , Ωe and Γe are the nodal values of the displacement vector field u, the nodal values of the temperature change field θ , the elasticity matrix, the area and the boundary of the element (e), respectively. B(e) u (e) (e) and B(e) θ contain the derivatives of the element shape functions Nu and Nθ , respectively. The displacement-strain (Bu ) and temperature-gradient (Bθ ) interpolation matrices and the relevant shape functions are defined as s (e) B(e) u = ∇ u Nu ,

(45)

B(e) θ

(46)

=

∇ θs

N(e) θ ,

s

with ∇ being the symmetric gradient and (Nu(e) )i j := Nn δi j˜

i = 1, . . . , du , j = 1, . . . , . . . , du × n (e) ,

(47)

(Nθ(e) )i j

i = 1, . . . , dθ , j = 1, . . . , . . . , dθ × n ,

(48)

n = ( j − 1)\dim + 1, j˜ = ( j − 1)%dim + 1,

(49)

:= Nn δi j˜

(e)

(50)

where Nn , n (e) , and du/θ are the Lagrangian shape function associated with node n of element (e), number of element nodes, and number of the degrees of freedom (DOFs) of the corresponding field at each node, respectively. Operator “\” is the integer division11 and operator “%” is the “modulo”/remainder12 operator. 5. Goal-oriented error estimation (GOEE) Due to the nature of numerical analysis methods, presence of error in results is inevitable, therefore in any numerical simulation, the main concern is the accuracy of the approximation. In many cases, the objectives of the numerical analysis is not only estimating the preliminary variables, but also calculating a functional of the primal solution. This functional which is presented based on preliminary variables, is called goal or quantity of interest (QoI). In such problems, the estimation of goal error is of interest rather than the solution error itself. The goal-oriented error estimation (GOEE) has been developed in order to quantify and control the local error in a quantity of interest (QoI). The goal might be expressed as a local or integral, continuous or noncontinuous functional based on the primal solution. Goal-oriented error estimation and mesh adaptation are one of the challenging topics in numerical analysis especially when dealing with pointwise quantities of interest due to the discontinuity of the goal functional. The main goal of this section is to employ the dual weighted residual (DWR) method for the adaptive solution of static coupled thermoelastic problems with pointwise goal of the primal solution value or its derivatives at a particular point. The considered QoI is discontinuous and only defined at specific points in the domain. 5.1. Dual weighted residual (DWR) method The DWR method estimates the error in quantities of interest based on a duality technique. In this technique, in addition to the primal problem (the original problem), a corresponding global dual or adjoint problem also needs to be solved. Solving the primal and dual problems will lead to estimating the error in the QoI. The auxiliary (dual) problem, which provides (local) sensitivity measures with respect to the goal functional, is defined based on the primal problem and the quantity of interest (QoI). Although the bilinear forms of the primal and dual problems are the same, the goal functional is considered as the functional of the dual problem [60,73]. In the following, we will apply and extend the abstract theory of the DWR method introduced by Becker and Rannacher [57,59], to error control in the approximation of thermoelastic problems. Let us assume the primal and relevant discretized variational formulations of the current problem to be (17) and (30), respectively. We rewrite these equations as: a (χ ) (ϑ) = l (ϑ) ∀ϑ ∈ V 0 , a(χ h )(ϑ h ) = l(ϑ h ) ∀ϑ h ∈ V 0h , 11 12

(51) (52)

a\b ≡ ⌊a/b⌋. a%b := a mod b ≡ a − b × ⌊a/b⌋.

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

15

where the solution spaces are χ ∈ V and χ h ∈ V h . Since the χ h is the finite element approximation of the exact solution χ , the exact discretization error of the primal problem in the primal variable field is defined as: ex eex − χh. χ =χ

(53)

Accordingly the Galerkin orthogonality13 property in the primal problem reads: ∀ϑ h := (ψ h , ϕ h ) ∈ V 0h .

h a(eex χ )(ϑ ) = 0

(54)

Now, consider the linear output functional J (•) : V → R defined on the space V of admissible functions for the primal problem. The aim of the goal-oriented error estimation is to calculate the error of J . The exact error in the QoI can be written as: ex ε ex − J h = J (χ ex ) − J (χ h ). J = J

(55)

Considering the fact that the goal functional J is assumed to be linear, relations (53) and (55) lead to: χ ex ) − J (χ h ) = J (χ ex − χ h ) = J (eex ε ex χ ). J = J(

(56)

Since there is not any analytical solution available for general coupled thermoelastic problems, the exact solution χ and the subsequent exact goal value cannot be computed, therefore we approximate the exact error by an estimated error calculated based on the DWR method. The basic part of the DWR method for ‘goal-oriented’ a posteriori error estimation is defining the dual problem. For any given (linear) goal functional J (•) defined on V 0 , we introduce the continuous and subsequently the discrete dual problems as: [59,73] a (υ) (~) = J (υ) h

h

h

∀υ ∈ V 0 , h

a(υ )(~ ) = J (υ ) ∀υ ∈

(57)

V 0h ,

(58)

where a(•)(•) is the bilinear associated with the primal problem and the test functions are chosen from the corresponding test function space. The ~ := {~ u , ~θ } ∈ V denotes the solution of the adjoint problem, which are the sensitivities that measure the influence of the error functional J ; hence, the exact error of the dual problem is defined as: e∗~ ex = ~ ex − ~ h ,

(59)

Accordingly, the Galerkin orthogonality property in the dual problem reads: a(υ h )(e∗~ ex ) = 0

∀υ h ∈ V 0h ,

(60)

By definition, the primal and dual solutions are mutually adjoint to each other on V in the sense that: J (χ ) = a (χ ) (~) = l (~) . Inserting as special test function υ := eχ primal error

(61) 14

into (57) yields:

a (eχ ) (~) = J (eχ ) .

(62)

Considering the Galerkin orthogonality property (54) of the primal problem, we can rewrite (62) as: a(eχ )(~ − ϑ h ) = a(χ ex − χ h )(~ − ϑ h ) = J (eχ ) ∀ϑ h ∈ V 0h ,

(63)

Similarly, substituting the dual error e∗~ as the test function in the primal problem (51) and considering the Galerkin orthogonality property (60) of the dual problem along with Eq. (63) carries over the mutual duality of the primal and dual solution (relation (61)) to the corresponding errors as J (eχ ) = a(eχ )(~) = a(eχ )(e∗~ ) = a(χ )(e∗~ ) = l(e∗~ ).

(64)

In multifield problems, the quantity of interest can depend on either or all of the fields. Therefore, in thermoelasticity, which we deal with a scalar-vector valued problem, the goal functional error can be expressed in terms of displacement and temperature change errors as: εJ = J (eχ ) = J ((eu , eθ )). 13 14

(65)

h In the sense of the bilinear form a, the error is orthogonal to the test space V 0h (eex χ ⊥a V 0 ). Henceforth, for simplicity we drop the “ex” superscript form errors (e = eex , e∗ = e∗ ex , and ε = ε ex ).

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

16

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Considering Eq. (63), the error in the quantity of the interest can be calculated as: εJ = J ((eu , eθ )) = a(eχ )(~ − ϑ h ) = a((eu , eθ ))(~ − ϑ h ) = a((eu , eθ ))((~ u − ψ h , ~θ − ϕ h )),

(66)

where eu = u − uh and eθ = θ − θ h are the exact errors of the displacement and temperature change fields of the primal problem. By substituting the bilinear form (28) we obtain: ( ) ( ) εJ =µ ∇eu , ∇(~ u − ψ h ) + λ ∇ · eu , ∇·(~ u − ψ h ) ( ) ( ) ( ) +µ ∇euT , ∇(~ u − ψ h ) − γ eθ , ∇·(~ u − ψ h ) + κ ∇eθ , ∇(~θ − ϕ h ) . (67) Replacing the primal errors eu and eθ based on the displacement and temperature change, we obtain: [ ( ) ( ) εJ = µ ∇u, ∇(~ u − ψ h ) + λ ∇ · u, ∇·(~ u − ψ h ) )] ) ( ) ( ( +µ ∇uT , ∇(~ u − ψ h ) − γ θ, ∇·(~ u − ψ h ) + κ ∇θ, ∇(~θ − ϕ h ) − [ ( ) ( ) µ ∇uh , ∇(~ u − ψ h ) + λ ∇ · uh , ∇·(~ u − ψ h ) ( ) ( ) ( )] T +µ ∇uh , ∇(~ u − ψ h ) − γ θ h , ∇·(~ u − ψ h ) + κ ∇θ h , ∇(~θ − ϕ h ) .

(68)

Splitting the global integrations over Ω into the contributions of the cells e ∈ P h and implementing cell-wise integration by parts, considering the divergence theorem and the constitutive equation (7), leads to: εJ =

Ne { [ ∑ ( ) ( ) − ∇·σ (u) , ~ u − ψ h Ωe − κ ∇·∇θ, ~θ − ϕ h Ωe e=1

( ) ( ) ] + σ (u) ·n, ~ u − ψ h ∂Ωe + κ∇θ ·n, ~θ − ϕ h ∂Ωe − [ ( ) ( ) − ∇·σ (uh ) , ~ u − ψ h Ωe − κ ∇·∇θ h , ~θ − ϕ h Ωe ( ) ( ) ]} + σ (uh ) ·n, ~ u − ψ h ∂Ωe + κ∇θ h ·n, ~θ − ϕ h ∂Ωe .

(69)

According to the governing equations (1) and (11), the first two integrals of Eq. (69) can be modified by ∇·σ (u) = −Fb and κ∇·∇θ = −Q v , respectively. Moreover, because the displacement and temperature fields have continuous first order derivatives inside the domain Ω , their normal derivative on one cell’s boundary cancels with that on its neighbor, where the normal unit vector has the opposite sign (n| Γ ∩∂Ω ′ = −n| Γ ∩∂Ωe ). In other words, e because σ (u) and κ∇θ are continuous, the third and fourth terms vanish on interior boundaries by summation on all cells. But, at the Dirichlet boundary (Γ D ) of the domain, where there is no neighbor cell with which these terms could cancel, the weights ~ u −ψ h and ~θ −ϕ h can be chosen as zero, since ~ has zero boundary values, and ϑ h can be chosen to have the same. On the other hand, at the Neumann boundary15 (Γ N ) of the domain, according to the Cauchy’s formula and Fourier’s law (10) these term can be replaced by the relevant Neumann boundary condition. Considering all these facts, we obtain the error representation εJ =

Ne {[ ∑ (

Fb , ~ u − ψ h

) Ωe

] ( ) ( ) ( ) + Q v , ~θ − ϕ h Ωe + ¯tn , ~ u − ψ h ∂Ωe ∩Γ t + q¯n , ~θ − ϕ h ∂Ωe ∩Γ q −

[ ( e=1 ) ( ) ]} ) ( ) ( − ∇·σ (uh ) , ~ u − ψ h Ωe − κ ∇·∇θ h , ~θ − ϕ h Ωe + σ (uh ) ·n, ~ u − ψ h ∂Ωe + κ∇θ h ·n, ~θ − ϕ h ∂Ωe . (70) Comparing Eq. (70) with bilinear form (28) and functional (29) shows that the summation of the terms inside the first bracket and inside the second bracket on all cells result in the functional and the bilinear form, respectively. By introducing the notation ρ(χ h )(•) := l(•) − a(χ h )(•) for the residual of the Galerkin approximation χ h (primal residual) we have ( ) ( )( ) εJ = l (~ u − ψ h , ~θ − ϕ h ) − a (uh , θ h ) (~ u − ψ h , ~θ − ϕ h ) (71) = l(~ − ϑ h ) − a(χ h )(~ − ϑ h ) =: ρ(χ h )(~ − ϑ h ). 15

Γ N ≡ Γ D∁ := ∂Ω \Γ D .

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

17

The residual ρ(χ h )(•) := l(•) − a(χ h )(•) of the Galerkin approximation χ h (primal residual) can be interpreted as a functional on the solution space V 0 . Finally, by combining the element integral terms of Eq. (70) and considering the primal residual16 (71), the QoI error reads as follows: εJ = J (eχ ) = ρ(χ h )(~ − ϑ h ) =

Ne { ∑ (

Fb + ∇·σ (uh ) , ~ u − ψ h

) Ωe

( ) + Q v + κ∇·∇θ h , ~θ − ϕ h Ωe

e=1

+ ¯tn , ~ u − ψ h (

)

∂Ωe ∩Γ t

( ) ( ) ( ) } + q¯n , ~θ − ϕ h ∂Ωe ∩Γ q − σ (uh ) ·n, ~ u − ψ h ∂Ωe − κ∇θ h ·n, ~θ − ϕ h ∂Ωe .

(72)

Depending on where the element boundary is, the last two integrals of Eq. (72) can have three different cases. In case the element boundary is located at the relevant Neumann boundary of the domain, these integrals can be merged with the relevant Neumann condition terms. If the element boundary coincides with the relevant Dirichlet boundary, mentioned integrals vanish due to the zero value of the corresponding weights. On the interior element boundaries (Γ ̸⊆ ∂Ω ), unlike the continuous terms which vanish by summation on all cells, the last two terms of Eq. (72) cannot be eliminated by summation and should be approximated by jump function. Since the solution derivatives (subsequently stress and heat flux) in the standard Lagrangian conformal finite element formulation based on the assumed displacement/heat approach and elements with C 0 continuity, are not required to be continuous at the element boundaries, the last two integrals of Eq. (72) do not have the same value between two adjacent element to cancel each other. These terms can be approximated by the averaging between two adjacent elements. Hence, for the two neighboring elements e, e′ ∈ P h with common edge Γ , the stress on the element edge is approximated by: ] ( ) ) ( ) 1 [( σ (uh ) ·n Γ ⊂∂Ω ′ + σ (uh ) ·n Γ ⊂∂Ωe σ (uh ) ·n Γ ⊂∂Ωe ≈ e 2 ] ⏐ ⏐ ] ⏐ ⏐ 1[ 1[ σ (uh ) ⏐Γ ′ ·n′ + σ (uh ) ⏐Γe ·n = − σ (uh ) ⏐Γ ′ − σ (uh ) ⏐Γe ·n = e e 2 2 1 = − Jσ (uh ) K Γe · n, (73) 2 where, Γe′ := Γ ∩ ∂Ωe′ and Γe := Γ ∩ ∂Ωe refer to the common edge in elements (e′ ) and (e), respectively and n′ is the normal unit vector of edge Γe′ pointing from (e′ ) to (e). In above averaging, considering the equal weights 12 for two elements, is due to the equal distribution assumption of the boundary contributions onto the two neighboring elements sharing the common edge. The operator J•K denotes the jump across the inter-element edges (in 2-D) or faces (in 3-D) as: J•K Γ ∩∂Ωe := •| Γ ∩∂Ω ′ − •| e

Γ ∩∂Ωe

.

(74)

Similarly,17 the heat flux on the common edge of the two adjacent elements is approximated by: ( ) 1 κ∇θ h ·n Γ ⊂∂Ωe ≈ − Jκ∇θ h K Γe · n. (75) 2 Eventually, by substituting (73) and (75) in (72) and defining the cell and edge residuals, we obtain the goal error representation εJ = J (eχ ) = ρ(χ h )(~ − ϑ h ) =

Ne { ∑ (

R u (χ h ), ~ u − ψ h

) Ωe

( ) + R θ (χ h ), ~θ − ϕ h Ωe

e=1

( ) ( ) } + r u (χ h ), ~ u − ψ h ∂Ωe + r θ (χ h ), ~θ − ϕ h ∂Ωe ,

(76)

with cell residuals R u and R θ and edge residuals r u and r θ , as smoothness indicators, defined by R u (χ h )|Ωe := Fb + ∇·σ (uh ) , R θ (χ h )|Ωe := Q v + κ∇·∇θ h ,

(77) (78)

16 In this article, all residuals are evaluated with respect to the primal problem. Henceforth, for simplicity we use “residual” instead of “primal residual”. ( ) 17 In general • · n ≈ − 12 J•K Γe · n. Γ ⊂∂Ωe Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

18

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

u

r (

χh

)|Γ ⊂∂Ωe

r (χ h )|Γ ⊂∂Ωe θ

⎧ 1 ⎪ ⎪ 2 Jσ (uh ) K·n ⎨ := 0 ⎪ ⎪ ⎩¯ tn − σ (uh ) ·n ⎧1 h ⎪ ⎨ 2 Jκ∇θ K·n := 0 ⎪ ⎩ q¯n − κ∇θ h ·n

if Γ ⊂ ∂Ωe \∂Ω , if Γ ⊂ Γ u ,

(79)

if Γ ⊂ Γ t , if Γ ⊂ ∂Ωe \∂Ω , if Γ ⊂ Γ θ , if Γ ⊂ Γ q .

(80)

These residuals are weighted by ~ u − ψ h and ~θ − ϕ h , which are the influence (sensitivity) factors and indicate how important the residuals on each cell is for the evaluation of the error in the goal functional. In order to evaluate the goal error, the weights ~ − ϑ h should be calculated. Indeed, since ϑ h is an arbitrary discrete test function, for instance, it can be set to be the point interpolation of the dual solution, ϑ h = Ih ~, with the interpolant Ih : V 0 → V 0h . Hence, the dual solution yields local sensitivity information with respect to the given quantity of interest. Considering these weights, the exact representation of the error in the QoI would be εJ = J (eχ ) = ρ(χ h )(~ − Ih ~) =

Ne { ∑ (

R u (χ h ), ~ u − Ih ~ u

) Ωe

( ) + R θ (χ h ), ~θ − Ih ~θ Ωe

e=1

( ) ( ) } + r u (χ h ), ~ u − Ih ~ u ∂Ωe + r θ (χ h ), ~θ − Ih ~θ ∂Ωe ,

(81)

The dual solution ~ is only analytically known for some very special cases and in general cases, which the dual solution is not available, the exact evaluation of the goal error is not possible. Consequently, we can compute the dual solution numerically and approximate ~ through a discrete function ~ˆ h , obtained from solving the discretized dual problem (58). Accordingly, we obtain an approximate error representation εJ = J (eχ ) ≈ a(eχ )(~ˆ h − Ih ~ˆ h ) = ρ(χ h )(~ˆ h − Ih ~ˆ h ).

(82)

For evaluation of the error from Eq. (82), computation of the dual problem with the same accuracy as the primal problem (~ˆ h = ~ h ∈ V 0h ) would not be applicable. The difficulty is that, since ~ h and subsequently ~ h − Ih ~ h lie in V 0h and according to the Galerkin orthogonality, approximating the dual solution with the same accuracy (the same mesh and element type) as the primal one yields J (eχ ) ≈ a(eχ )(~ h − Ih ~ h ) = ρ(χ h )(~ h − Ih ~ h ) ≡ 0.

(83)

Hence, the dual solution must be approximated with higher accuracy and ~ˆ has to be from a larger space than the primal finite element space. In order to overcome this problem and get an approximation ~ˆ h ̸∈ V 0h , methods such as global-higher order approximation (using a higher order finite element), higher resolution approximation (using finer mesh), or local higher-order approximations (using a patch-wise higher-order interpolation), can be adopted [59,73]. Although the last approach is the cheapest, in this study, for simplicity we opt the first method and compute the dual problem with higher polynomial degree than the primal one. To solve the adjoint problem, we employed a global-higher order finite element of degree d which d is greater than p, the approximation degree of the primal solution. By considering d > p, the weights ~ h(d) − Ih ~ h(d) would be finer that the trial space of the primal solution and the residual would not be orthogonal, therefore, the overall error estimate would not be zero. Ultimately, we obtain a computable error representation h

εJ = J (eχ ) ≈ ρ(χ h( p) )(~ h(d) − Ih ~ h(d) )

= a(eχ )(~ h(d) − Ih ~ h(d) ) = a(χ − χ h( p) )(~ h(d) − Ih ~ h(d) ) Ne { ∑ ( u h( p) ) ( ) = R (χ ), ~ uh(d) − Ih ~ uh(d) Ωe + R θ (χ h( p) ), ~θh(d) − Ih ~θh(d) Ωe e=1

( ) ( ) } + r u (χ h( p) ), ~ uh(d) − Ih ~ uh(d) ∂Ωe + r θ (χ h( p) ), ~θh(d) − Ih ~θh(d) ∂Ωe =εJes ,

d > p,

(84)

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

19

where interpolant Ih projects the numerically approximated dual solution with degree d into the finite element space of the primal solution with degree p. Eq. (84) shows that the error in the QoI can be expressed in terms of cell and edge residual weighted by factors provided by dual solution. Thus, the adjoint solution describes the dependencies of error in the QoI, which is concentrated at some cells K , to the local properties of the data (i.e., the residuals of cells K ′ ). Considering these features, the dual solution can be interpreted as a “generalized” Green’s function G(K , K ′ ) [73,74]. Usually, in error estimation instead of the absolute error, the relative error (η) or the percent relative error (η% ) is used [75]. In goal-oriented error estimation, the exact relative error in the quantity of interest can be defined as: |ε J | |J (eχ )| (85) = ; ηexJ % = ηexJ × 100. ηexJ := |J (χ ex )| |J (χ ex )| For the problems without analytical solution, the estimated relative error is calculated instead of the exact relative error. An estimation of the relative error is calculated by |εJes | |εJes | ηexJ ≈ ηesJ := , (86) = es |J (χ )| |J (χ h ) + ε es | J

χ es

χh

ees χ

is the estimated primal field. where = + In conclusion, for evaluation of the goal error, after solving the primal problem numerically, the dual problem should be solved by FEM with higher degree than the primal one. Then, considering Eq. (84), the error contribution of each element is evaluated by calculating the cell and edge residuals weighted by the (local) adjoint sensitivity factors. Finally, all element errors are summed to approximate the total error in the QoI. Furthermore, Eq. (84) is a cell-wise quantity, therefore in the employed refinement step it is used as a mesh refinement criterion where based on each element contribution in the goal error, mesh is refined/coarsened adaptively with allowing hanging nodes. 5.2. Pointwise error estimation In many real cases, the quantity of interest is “singular”; in the other words, the functional J (•) is not defined on W 1,2 (Ω ) properly but only on some subspaces containing the solution and on the finite element space [57]. One of the typical examples of such cases is the pointwise functional. In this study our main objective of implementing DWR a posteriori error estimation is to find an estimate for the pointwise error. Pointwise error estimation aims at assessing the accuracy of the solution or function of the solution at a given point x0 ∈ Ω [50,51]. This article is devoted to the study of a particular case of pointwise error estimation to evaluate the error between the exact and approximated solutions or their derivative at a given point of interest (PoI). Point value error: To evaluate the accuracy of the primal solution, the quantity of interest would be QoI := χ i (x0 )

i ∈ {1, . . . , dim + 1}, x0 ∈ Ω .

(87)

According to the Sobolev Embedding Theorem,18 since χ ∈ V(Ω ) := T (Ω ) × T (Ω ) and Ω ⊂ Rd , the solution may not be continuous for geometrical dimensions d ⩾ 2, therefore χ i may not be defined/well-defined at x0 [51]; subsequently, the linear goal functional J (•) corresponding to the solution value at x0 is not necessarily bounded, J (χ i ) := χ i (x0 )

i ∈ {1, . . . , dim + 1}, x0 ∈ Ω .

(88)

1,2

In the other words, while the functions of W are continuous in 1D, for the continuity (i.e., having meaningful point values) in 2D/3D, the functions of W 2,2 or higher orders are required [22,50]. These statements implies that when the solution convergence is obtained in H1 , only the point value continuity in 1D is guaranteed and in higher geometrical dimensional domains (2D/3D) the H1 convergence does not ensure the point value continuity [76]. In order to circumvent the continuity issue of the quantity of the interests defined in regions of smaller dimension than the domain dimension (i.e., points/lines in 2D and points/lines/surfaces in 3D), two methods have been 18 From the Sobolev Embedding Theorem we know that if 2 m > d, then Hm ⊂ C 0 , therefore in 1D, H1 ⊂ C 0 (so the solution is continuous) and in 2D/3D, H2 ⊂ C 0 .

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

20

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

suggested: using regularization by Bangerth and Rannacher (see [73, chap. 3.3]); and using so-called mollification by Prudhomme and Oden (see [77, chap. 2] and [22, chap. 8.5]). In the current study we prefer to employ the regularization approach for pointwise quantities of interest. By considering the regularized functional Jε , the corresponding dual solution behaves like a regularized Green function. In spite of the dependency of the type of regularization on the particular situation, for prescribed error tolerance T O L, the regularization should be such that [57] |J (χ i ) − Jε (χ i )| ≈ ε,

ε := T O L ,

(89)

where Jε (•) is the regularized form of the functional J (•). Accordingly, to estimate the error at a given point x0 ∈ Ω , instead of the original functional J we consider the regularized functional Jε as the quantity of the interest, ∫ χ i (x) dx = χ i (x0 ) + O(ε 2 ) Jε (χ i ) := |Bε (x0 )|−1 Bε (x0 )

= J (χ i ) + O(ε 2 ),

(90)

where Bε := {x ∈ Ω : |x − x0 | < ε} is the ε-ball centered around the point x0 and ε := T O L denotes the accuracy we want to achieve. Considering the regularized goal functional instead of the pointwise one in GOEE, avoids over-refinement near the singularities which means less computational work. Point value derivative error In order to estimate the pointwise error in the first derivative of the solution, the quantity of interest and subsequently the goal functional would be J (χ i ) ≡ QoI := ∂ j χ i (x0 )

i, j ∈ {1, . . . , dim + 1}, x0 ∈ Ω .

(91)

Since the pointwise derivatives are generally not defined at the element interfaces or not necessarily well-defined and bounded at x0 ∈ Ω , the corresponding dual solution may not exist in the sense of H1 (Ω ). Therefore, an approach similar to the point value error estimation, should be taken to estimate error in the point value derivative. Thus, the original functional J should be regularized according to ∫ Jε (χ i ) := |Bε (x0 )|−1 ∂ j χ i (x) dx = ∂ j χ i (x0 ) + O(ε 2 ) Bε (x0 )

= J (χ i ) + O(ε 2 ),

(92)

where again Bε is the ε-ball centered around the point x0 and ε := T O L. By considering the current regularized functional Jε , the resultant dual solution behaves like a regularized derivative Green function. In this case, |J (χ i ) − Jε (χ i )| ≈ ε and the corresponding dual solution would be well-defined in H1 (Ω ). Considering this regularization avoids over-refinement near the singularities and reduces computational cost. 6. Numerical results and discussion In this section, in order to substantiate the developed DWR goal-oriented error estimation method, 2D and 3D numerical examples with available analytical solutions are studied. Comparing the results with analytical solution and other error estimators and refinement techniques demonstrates the efficiency and reliability of the proposed adaptive algorithm. In this section, four numerical examples are investigated and in each example, in addition to the proposed goal-oriented mesh adaptivity process, for comparison, global uniform, Kelly [78] and weighted Kelly (WKelly) refinement techniques are also implemented. In the presented method, the models are discretized by h-nonconforming unstructured meshes19 and goal-oriented based adaptive h-refinement using hanging nodes is implemented to enhance the accuracy of analyses. In this study only 1-irregular mesh, allowing maximum one disconnected (hanging) node over an edge, is considered. In the refinement process, in order to obtain a smoother mesh than just the bare minimum, Bank’s [79] “one-irregular” and “k-neighbor” rules are considered [64,65]. 19

Initially conforming structured but h-nonconforming unstructured after adaptation.

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

21

Fig. 2. Model geometry and boundary conditions of the thick-walled cylinder (example 6.1).

Employed adaptive local h-refinement/coarsening with allowing hanging nodes saves CPU time very much and reduces computational cost considerably. Unlike, in the global uniform refinement process, at each refinement step, all elements are uniformly subdivided into 4 and 8 elements in 2D and 3D meshes, respectively. All problems reported here are analyzed considering homogeneous isotropic material and in order to assume constant (temperature-independent) material properties, it is assumed that the temperature does not depart from the initial values very much. 6.1. Thick-walled cylinder subjected to internal pressure and temperature change with the goal of point value In this example, in order to verify the accuracy of the presented error estimation method, a thick-walled cylindrical pressure vessel under thermo-mechanical loadings is studied. We employed DWR GOEE and adaptive mesh refinement methods to estimate a pointwise error and verified the accuracy of the implemented methods by analytical solution. The model geometry and boundary conditions are illustrated in Fig. 2. Generally, the cylinders with wall thickness to radius ratio greater than 0.1 is classified as “thick-walled” cylinders. Here, we selected the thickness to mean radius ratio of 0.66. The geometry parameters of the problem are 2t t = ≈ 0.66, (93) rin = 5 m, rout = 10 m, rm rin + rout where, t is the wall thickness and rm is the mean radius. The vessel is designed to store a cold pressurized gas with operating inside pressure Pin = 25 bar and temperature Tin = 230 K (≈ −43 ◦ C). The outside pressure and temperature assumed to be Pout = 1 bar (≈ 1 atm) and Tout = 290 K, respectively. All other geometrical parameters and boundary conditions are presented in Fig. 2. The analytical solution of this example is stated in Appendix A.1. We suppose the vessel material to be stainless Steel Grade 304 (AISI 304, DIN 1.4301, UNS S30400) with properties [80,81] E =200 GPa,

ν =0.27,

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

22

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 3. Initial mesh used to analyze the thick-walled cylinder of problem 6.1.

λ =92.4341 GPa, K =15 W/(m K),

µ =78.7402 GPa, α T =15.0 × 10−6 1/K.

(94)

In this example, our main interest is the displacement in x direction (u x ) at a specific point of interest (PoI). Therefore, we employed DWR GOEE method to estimate the error of u x(PoI ) and refined adaptively based on the presented DWR method. As depicted in Fig. 2, the point of interest is at rin + rout = 7.5 m, θ = 45◦ . (95) r PoI = 2 In each error estimation step, the primal model is discretized by an unstructured mesh with isoparametric bilinear quadrilateral (Q4) Lagrangian elements and solved numerically. The dual problem is solved on the primal mesh but with higher order elements. Based on the dual and primal solutions and according to Eq. (84), the goal error in each element is calculated. Then, based on each element contribution in the goal error, mesh is refined/coarsened adaptively with allowing hanging nodes. Afterwards, the new refined mesh is considered as the primal mesh for the next error estimation step. This process is repeated until we reach the desired accuracy. In this example in order to reach goal error order 10−9 , there were 21 refinement steps. Fig. 3 presents the initial mesh while Figs. 4(a)–4(f) show the refined meshes after 5th, 7th, 10th, 13th, 15th, and 17th refinement steps, respectively. The refined meshes of Fig. 4 show that, in order to minimize the error of the displacement in x direction at a specific point, the refinement process mainly tends to refine the elements near PoI and the boundary fixed in x direction (please note the subtle interplay between resolving around the point of interest and the fixed boundary). Figs. 5(a)–5(d) show contour plots of the temperature change and displacement in r, x and y directions for the primal problem in the last refined mesh, respectively. Contours of Fig. 5, show a very good agreement with exact solutions of Appendix A.1. As mentioned in Section 5.2, the dual solution can be interpreted as a Green function. Figs. 6(a)–6(d) show the regularized Green function for evaluating u x(PoI ) in 5th, 7th, 13th, and 15th refinement steps, respectively. In these figures, which show the dual solution, the refined mesh is elevated by the x-component of the dual solution. According to the analytical solution (see Appendix A.1), the exact value of the goal in this example is 0.02635804. Considering the goal analytical value, the exact errors in different refinement steps of the current example are calculated and compared with the DWR estimated goal error in Fig. 7. This comparison shows that even for such Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

23

Fig. 4. Refined meshes of problem 6.1 based on the pointwise quantity of interest u x(PoI ) .

a complicated quantity of interest as the point value, the true and estimated errors are actually very close to each other. Moreover, the plotted errors show a monotonic decrease by increasing the total number of the DOFs (except for the initial very coarse meshes), which indicates the stable convergence of the refinement process. As depicted in Fig. 7, in the full logarithmic scale, the estimated and exact errors fit straight lines with slopes of −1.048 and −1.031, respectively. In order to show the optimality of the employed mesh refinement, the complexity O(N −1 ) is also indicated in this graph. In Fig. 8, the achieved error levels for different error estimation and mesh refinement criteria are compared. The selected methods for the comparison are DWR, global uniform refinement, Kelly, and weighted Kelly techniques. In this graph, the estimated errors in different refinement steps of each method are plotted in log–log scale. In this example, the weights in the weighted Kelly indicator are chosen to be 1/(r 2 +0.12 ), where r is the distance to the PoI. From these results it can be concluded that the estimated goal √ error in DWR method is proportional to O(1/N ), while in global uniform refinement is proportional to O(1/( N × log(N ))). This comparison shows that among the selected error estimation and refinement techniques, the employed DWR method with convergence rate 1.048 is the optimum and most efficient one. In error estimation, in order to assess the accuracy of the employed error estimator, error effectivity index is used. The error effectivity index, which represents the degree of over or underestimation on the current mesh, is defined as the ratio between the estimated and exact errors and should be ideally close to 1.0 [82]. For goal-oriented error estimation, the error effectivity index is defined by θε :=

εesJ εJ

,

(96)

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

24

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 5. Solution contour plots of the primal problem 6.1 in the last refined mesh.

Fig. 6. Elevated representation of the regularized Green function for evaluating u x(PoI ) in problem 6.1.

where θε is the goal error effectivity index. Since the QoI in this problem is linear, the goal error effectivity index can be calculated as ε es ε esJ ε esJ θε := J = = . (97) εJ J (eχ ) J (χ ex ) − J (χ h ) Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

25

Fig. 7. Estimated and exact errors of the problem 6.1 for DWR method.

Fig. 8. Errors of the problem 6.1 estimated by different refinement techniques.

It is also possible to represent the effectivity in the QoI defined as θ QoI :=

J (χ h ) + ε esJ , J (χ ex )

(98)

which also should desirably be close to 1. χ h ) can be calculated from the FEM results. However, In Eqs. (97) and (98), ε es J can be obtained from (84) and J ( J (χ ex ) can only be computed for problems with an available analytical solution. The exact value of the QoI in this example is J (χ ex ) = 0.02635804 and the exact error in the quantity of interest is obtained by ε J = J (χ ex − χ h ) = J (χ ex ) − J (χ h ) h = u ex x (PoI ) − u x

(PoI )

= 0.02635804 − u hx

(PoI ) .

(99)

Accordingly, the goal error effectivity index of the DWR method is calculated and plotted in Fig. 9. As depicted in this plot, in spite of the fact that the implemented method overestimates the error in the very coarse meshes of this example (the initial and first two refined meshes), by further refinements, the θε lies in the range of 0.9–1.0 and approaches to the ideal value of 1. Evolution of the goal error effectivity index shows the very good accuracy and the efficiency of the implemented DWR and refinement methods. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

26

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 9. The goal error effectivity index of DWR method in the problem 6.1, considering u x(PoI ) as the quantity of interest. Table 1 Goal error estimation results of the problem 6.1 using DWR method. Refinement step

DOF (Primal)

DOF (Dual)

J (u hx )

εJ

εJes

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

45 114 189 336 594 1017 1716 3024 4908 7656 13080 20766 32463 52548 80457 119727 195228 289923 411621 679521 997137 1415901

135 387 672 1239 2220 3888 6588 11781 19185 30033 51636 82131 128697 208863 319980 476658 778170 1156281 1642452 2713281 3982848 5656251

0.026090624 0.026313937 0.026291983 0.026319282 0.026341035 0.026349482 0.026353596 0.026356242 0.026356684 0.026356941 0.026357519 0.026357667 0.026357795 0.026357903 0.026357953 0.026357972 0.026358006 0.026358024 0.026358023 0.026358031 0.026358038 0.026358039

2.6742×10−4 4.4107×10−5 6.6061×10−5 3.8762×10−5 1.7009×10−5 8.5620×10−6 4.4481×10−6 1.8022×10−6 1.3599×10−6 1.1028×10−6 5.2541×10−7 3.7719×10−7 2.4921×10−7 1.4078×10−7 9.1348×10−8 7.2445×10−8 3.7922×10−8 1.9659×10−8 2.0978×10−8 1.3458×10−8 6.3420×10−9 4.6342×10−9

3.7193×10−4 5.9139×10−5 7.9554×10−5 3.6045×10−5 1.4703×10−5 7.4484×10−6 4.0967×10−6 1.7019×10−6 1.2772×10−6 1.0466×10−6 4.9680×10−7 3.6352×10−7 2.3838×10−7 1.3701×10−7 8.8622×10−8 7.0281×10−8 3.6946×10−8 1.9167×10−8 2.0211×10−8 1.3034×10−8 6.0964×10−9 4.5025×10−9

The goal error estimation results obtained for the current problem using dual weighted residual estimator and the corresponding adapted meshes are shown in Tables 1 and 2 (the point value and estimated error results in these tables are rounded and the other presented data are calculated based on the original results with higher precision). The presented results indicate clearly that by implementing the current DWR estimator and refinement methods, the goal error effectivity index (θε ) and effectivity in the QoI (θ QoI ) values approach to the desired value of 1. 6.2. Thick-walled cylinder subjected to internal pressure and temperature change with the goal of point value derivative As the second example, the problem investigated in the previous section has been studied with a different pointwise quantity of interest. In this problem, we considered singular point value derivative of ∂u x /∂ x at the same Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

27

Table 2 Relative error and effectivity indices of DWR method in the problem 6.1. Refinement step

DOF (Primal)

ηex J

ηes J

θε

θ QoI

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

45 114 189 336 594 1017 1716 3024 4908 7656 13080 20766 32463 52548 80457 119727 195228 289923 411621 679521 997137 1415901

1.0146×10−2 1.6734×10−3 2.5063×10−3 1.4706×10−3 6.4532×10−4 3.2484×10−4 1.6876×10−4 6.8375×10−5 5.1592×10−5 4.1838×10−5 1.9934×10−5 1.4310×10−5 9.4549×10−6 5.3409×10−6 3.4657×10−6 2.7485×10−6 1.4387×10−6 7.4585×10−7 7.9587×10−7 5.1060×10−7 2.4061×10−7 1.7582×10−7

1.4055×10−2 2.2424×10−3 3.0167×10−3 1.3676×10−3 5.5786×10−4 2.8260×10−4 1.5543×10−4 6.4569×10−5 4.8457×10−5 3.9709×10−5 1.8848×10−5 1.3792×10−5 9.0440×10−6 5.1982×10−6 3.3622×10−6 2.6664×10−6 1.4017×10−6 7.2719×10−7 7.6678×10−7 4.9450×10−7 2.3129×10−7 1.7082×10−7

1.390809194 1.340798869 1.204261472 0.929894020 0.864388923 0.869936510 0.921017105 0.944335350 0.939247641 0.949105772 0.945558733 0.963761386 0.956546350 0.973273406 0.970160855 0.970131146 0.974273941 0.974992293 0.963444234 0.968472439 0.961274720 0.971590440

1.003965019 1.000570290 1.000511936 0.999896903 0.999912487 0.999957751 0.999986671 0.999996194 0.999996866 0.999997871 0.999998915 0.999999481 0.999999589 0.999999857 0.999999897 0.999999918 0.999999963 0.999999981 0.999999971 0.999999984 0.999999991 0.999999995

point of interest (PoI) as the quantity of interest. In order to optimally minimize the goal error, we implemented the presented DWR goal error estimator and accordingly the adaptive mesh refinement processes. The efficiency of the method is demonstrated by comparing the results with analytical solution. The exact value of the goal in this example is J (∂x u x ) = 0.00506189 (see Appendix A.1). In this example, 23 refinement steps are required to reach the goal error level of 10−9 . While the initial mesh is the same as the problem 6.1 (Fig. 3), the refined meshes after the 5th, 10th, 13th, 17th, 19th, and 21st refinement steps are presented in Figs. 10(a)–10(f), respectively. In this problem, the exact solution is same as the example 6.1 and only the meshes are different, therefore the primal solutions are roughly the same but the dual solutions are quite different. The regularized derivative Green function for evaluating ∂u x /∂ x at the point of interest corresponding to the 8th, 9th, 10th, and 11th refinement steps are presented in Figs. 11(a)–11(d), respectively. According to the analytical solution, the exact error in the quantity of interest is calculated by ε J = J (χ ex − χ h ) = J (χ ex ) − J (χ h ) h = ∂x u ex x (PoI ) − ∂x u x

(PoI )

= 0.00506189 − ∂x u hx

(PoI ) .

(100)

Subsequently, the goal error estimation results of the presented DWR method in each refinement step of the problem 6.2 are presented in Table 3. Checking the estimated error, we see that for the DWR method, the estimated error has the convergence rate of 0.941 and is proportional to O(1/N ). Considering the goal analytical value, the goal error effectivity index which is a standard measure of the quality of an estimator, is calculated and plotted in Fig. 12. This plot shows that the presented DWR error estimation method, calculates the goal error with very good accuracy even for a complicated singular pointwise quantity of interest like point value derivative at a single point. In this example, although DWR underestimates the actual error, after the first three very coarse meshes, the θε lies in the range of 0.9–1.0, which proves the accuracy and efficiency of the implemented error estimator and mesh refinement methods. Notwithstanding the irregular behavior of the error effectivity index in this problem, by refining the mesh, fluctuations are decreased and the mean value tends to the ideal value of 1. In addition to presented goal error effectivity index, calculating the θ QoI shows that the proposed technique accurately captures the true error in the evaluation of the point value derivative and a single point. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

28

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 10. Refined meshes of the problem 6.2 based on the pointwise quantity of interest

Fig. 11. Dual solution, regularized derivative Green function for evaluating

⏐

∂u x ⏐ ∂ x ⏐ PoI ,

⏐

∂u x ⏐ ∂ x ⏐ PoI .

in problem 6.1.

Similar to the previous example, comparing the estimated error among DWR, Kelly, weighted Kelly (with (r + 0.12 )−1 weights), and global uniform refinement methods shows that the DWR method has the optimal convergence rate and better efficiency among them. As depicted in Fig. 13, the convergence rate of the DWR method is 0.941. 2

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

29

Table 3 Goal error estimation results of the problem 6.2 using DWR method. Refinement step

No. of active cells

DOF (Primal)

DOF (Dual)

J (u hx )

εJes

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

8 14 41 56 89 143 209 371 560 911 1418 2381 3710 5981 9674 15473 24302 38978 61193 94163 148955 229571 347537 536279

45 72 174 228 336 549 753 1317 1980 3099 4740 7818 12033 19188 30807 48879 76410 122094 190863 292566 462159 711156 1074432 1653183

135 231 612 819 1233 2073 2883 5094 7737 12168 18663 30900 47715 76164 122562 194667 304527 487047 761796 1168302 1846083 2841594 4294062 6607830

0.0050051216 0.0050069059 0.0050402112 0.0050507257 0.0050575234 0.0050569003 0.0050589620 0.0050583819 0.0050593138 0.0050603200 0.0050611874 0.0050615044 0.0050616025 0.0050617190 0.0050618248 0.0050618030 0.0050618341 0.0050618681 0.0050618646 0.0050618724 0.0050618809 0.0050618831 0.0050618849 0.0050618854

3.4514×10−5 4.1773×10−5 1.9317×10−5 1.0346×10−5 4.0426×10−6 4.5137×10−6 2.6801×10−6 3.3633×10−6 2.3348×10−6 1.4727×10−6 6.4821×10−7 3.5504×10−7 2.7068×10−7 1.5859×10−7 5.8888×10−8 8.2488×10−8 5.1078×10−8 1.9436×10−8 2.3159×10−8 1.5239×10−8 7.2957×10−9 5.3219×10−9 3.5168×10−9 2.9817×10−9

Fig. 12. Evolution of the goal error effectivity index of the problem 6.2 considering

⏐

∂u x ⏐ ∂ x ⏐ PoI

as the quantity of interest.

6.3. Thick-walled sphere subjected to internal pressure and temperature change with the goal of point value In this example we implement the presented DWR method to estimate the goal error in a thick-walled sphere under thermo-mechanical loading. In order to prove the accuracy of the implemented error estimation method, the results are evaluated by analytical solution. The analytical solution of this problem is stated in Appendix A.2. Let us consider a sphere under pressure/temperature loads as illustrated in Fig. 14. The geometry parameters of the problem are t 2t rin = 5 m, rout = 10 m, = ≈ 0.66, (101) rm rin + rout Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

30

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 13. Errors of the problem 6.2 estimated by different refinement techniques.

Fig. 14. Model geometry and boundary conditions of the thick-walled sphere (example 6.3).

Similar to the example 6.1, the operation conditions are Tin = 230 K, Tout = 290 K,

Pin = 25 bar, Pout = 1 bar,

(102)

and the material considered to be the same as problem 6.1 with properties (94). In this example, we considered the displacement in x direction (u x ) at point of interest (PoI) as the quantity of interest. As depicted in Fig. 14, the point of interest is assumed to be at rin + rout r PoI = = 7.5 m, ϕ = 45◦ , θ = 60◦ . (103) 2 According to the analytical solution (Appendix A.2), the exact value of the quantity of the interest is 0.01831743. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

31

Fig. 15. Initial mesh used to analyze the thick-walled sphere of problem 6.3.

Fig. 16. Refined meshes of the problem 6.3 based on the pointwise quantity of interest

⏐

∂u x ⏐ ∂ x ⏐ PoI .

In order to reach the error in the QoI of order 10−6 , we require 9 refinement steps. In each step, the primal model is discretized by an unstructured mesh with isoparametric trilinear hexahedral (H8) Lagrangian elements and solved numerically. Similar to the 2D examples, the dual problem is solved on the primal mesh but with higher order elements. The initial mesh and the refined meshes of 7th, 8th, and 9st refinement steps are presented in Figs. 15 and 16(a)–16(c), respectively. The refined meshes of Fig. 16, show that in this example, the presented refinement method refines the areas near PoI and corners of the fixed boundary in the x direction more than the other areas. This point is also clearly depicted in the spherical slice of the refined meshes shown in Fig. 17. In these figures, the point at the middle of the slice is the point of the interest. The counter plots of the temperature change, displacement in x and r directions for the primal problem in the last refined mesh are plotted in Figs. 18(a)–18(c), respectively. Due to the symmetry of the solution, the couture plots of the displacement in y and z directions are similar to the x direction. Contours of Fig. 18, show a very good agreement with exact solutions of Appendix A.2. The comparison between the estimated error by DWR method and the exact error calculated based on the exact goal value of 0.01831743 is shown in Fig. 19. This comparison shows the good accuracy of the employed method even for a singular QoI in a 3D problem. In this plot, the estimated error shows a monotonic decrease by increasing the total number of the DOFs with convergence rate of 0.807. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

32

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 17. Spherical slice of the refined meshes of the problem 6.3 at radius r PoI .

Fig. 18. Solution contour plots of the primal problem 6.3 in the last refined mesh.

Fig. 19. DWR estimated and exact errors in the QoI of the problem 6.3.

For assessing the agreement between exact and estimated error curves, the so-called effectivity index parameter is calculated and presented in Fig. 20. Figs. 19 and 20 show that for the coarse meshes (initial and first three refined meshes), the presented method overestimates the true error (ε esJ > ε J , θε > 1), but with further refinements the exact error is underestimated. In Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

33

Fig. 20. Evolution of the goal error effectivity index of the problem 6.3. Table 4 Goal error estimation results of the problem 6.3 using DWR method. Refinement step

DOF (Primal)

DOF (Dual)

J (u hx )

εJ

εJes

0 1 2 3 4 5 6 7 8 9

228 456 1300 2224 4588 8036 16380 34792 67696 180744

1220 2704 8404 14788 31356 56456 116004 253752 504792 1355592

0.018086918 0.018077492 0.018198115 0.018266274 0.018267977 0.018285846 0.018305114 0.018307450 0.018313350 0.018314369

2.3051×10−4 2.3994×10−4 1.1932×10−4 5.1158×10−5 4.9455×10−5 3.1586×10−5 1.2317×10−5 9.9815×10−6 4.0815×10−6 3.0621×10−6

3.8655×10−4 3.7329×10−4 1.5822×10−4 5.8697×10−5 4.8563×10−5 2.9590×10−5 1.1807×10−5 9.2985×10−6 3.7934×10−6 2.8906×10−6

this example, although in the first three refinement steps, DWR method underestimates the exact error, with more refinements the θε lies in the range of 0.9–1.0, which proves the accuracy and efficiency of the presented error estimator and mesh refinement methods. These goal error estimation results considering the DWR method, are also presented in Tables 4 and 5 (It should be noted that the point value and estimated error results in these table are rounded and the other presented data are calculated based on the original results with higher precision). In these tables, the first row corresponds to the initial mesh. These results show that, after the fourth refinement step, the average value of the goal error effectivity index (θε ) and effectivity in the QoI (θ QoI ) are 0.9470574 and 0.9999586, respectively. In order to evaluate the efficiency of the implemented adaptivity method, the estimated errors are compared with the results of other mesh refinement strategies. This comparison is presented in Fig. 21. In the weighted Kelly method, the weights are considered to be (r 2 + 0.12 )−1 . The results show that the employed DWR method with convergence rate of 0.9072 has the best efficiency among selected methods. 6.4. Thick-walled sphere subjected to internal pressure and temperature change with the goal of point value derivative In this part, as the last example, the problem of Section 6.3 is studied with goal of the point value derivative ∂u x /∂ x at the same point of interest (PoI). The DWR method has been implemented to estimate the error of the current pointwise singular goal and refine additively to optimally minimize the goal error. In order to show the accuracy and efficiency of the proposed method, the results have been compared with the exact solution. According to the analytical solution Appendix A.2, the exact value of the goal in this example is J (∂x u x ) = 0.00407859. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

34

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx Table 5 Relative error and effectivity indices of DWR method in the problem 6.3. Refinement step

DOF (Primal)

ηex J

ηes J

θε

θ QoI

0 1 2 3 4 5 6 7 8 9

228 456 1300 2224 4588 8036 16380 34792 67696 180744

1.2584×10−2 1.3099×10−2 6.5138×10−3 2.7928×10−3 2.6999×10−3 1.7243×10−3 6.7243×10−4 5.4492×10−4 2.2282×10−4 1.6717×10−4

2.0925×10−2 2.0231×10−2 8.6193×10−3 3.2031×10−3 2.6513×10−3 1.6156×10−3 6.4461×10−4 5.0765×10−4 2.0709×10−4 1.5780×10−4

1.676923077 1.555748206 1.326035290 1.147366608 0.981962737 0.936818751 0.958595492 0.931580080 0.929411765 0.943975642

1.008518668 1.007279727 1.002123740 1.000411572 0.999951302 0.999891054 0.999972158 0.999962717 0.999984272 0.999990634

Fig. 21. Errors of the problem 6.3 estimated by different refinement techniques.

In this example, the initial mesh is considered to be the same as the problem 6.3. In order to reach the error level of 10−7 , 9 refinement steps on the initial mesh are required. Figs. 22(a)–22(d) show the refined meshes after the 5th, 6th, 7th, and 8st refinement steps, respectively. In order to depict the refinement process around the point of the interest (PoI), the slice of the refined meshes of Fig. 22 at azimuthal angle ϕ = 45◦ are presented in Fig. 23 Since the exact solution of this problem is the same as the example 6.3 (see Appendix A.2) and only the meshes are different, the primal solution of this example is roughly the same as the previous one. Despite the equality in the primal solution, the dual solutions are quite different. Figs. 24(a) and 24(b) show the elevated representation of ⏐ x⏐ the dual solutions for evaluating u x(PoI ) (problem 6.3) and ∂u (problem 6.4), respectively. These plots relate ∂ x PoI to the dual solutions on the ϕ = 45◦ slice plane in the last refined meshes of the corresponding problems. In Table 6, the goal error estimation results of the presented DWR method are presented in each refinement step. Studying the estimated error, shows a convergence rate of 0.752 for the implemented DWR and refinement methods. Considering to the analytical solution, the true error in the quantity of interest is calculated by ε J = J (χ ex − χ h ) = J (χ ex ) − J (χ h ) h = ∂x u ex x (PoI ) − ∂x u x

(PoI )

= 0.00407859 − ∂x u hx

(PoI ) .

(104)

Based on the formula (104) and the results of Table 6, for DWR method, the true error in the quantity of the interest is calculated and compared with the estimated error in Fig. 25. This comparison indicates that the exact and estimated errors show very good accordance even for a singular goal of the point value derivative in 3D. Moreover, Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 22. Refined meshes of the problem 6.4 to optimize the error of the pointwise goal of

35

⏐

∂u x ⏐ ∂ x ⏐ PoI .

Table 6 Goal error estimation results of the problem 6.4 using DWR method. Refinement step

No. of active cells

DOF (Primal)

DOF (Dual)

J (u hx )

εJes

0 1 2 3 4 5 6 7 8 9

24 66 206 570 1207 2243 3916 7808 16565 34317

228 604 1352 3652 7000 12296 21256 40448 82084 164412

1220 3516 8652 24884 48780 87708 153256 297952 613376 1242428

0.0040080330 0.0040145317 0.0040611361 0.0040643595 0.0040729112 0.0040755768 0.0040758990 0.0040772540 0.0040776165 0.0040780227

5.1673×10−5 5.6670×10−5 1.5183×10−5 1.3162×10−5 5.1874×10−6 2.7976×10−6 2.4849×10−6 1.2573×10−6 9.1346×10−7 5.3869×10−7

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

36

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 23. Slice of the refined meshes of the problem 6.4 at azimuthal angle ϕ = 45◦ .

Fig. 24. Elevated representation of the dual solutions on the ϕ = 45◦ slice plane of the 9th refinement steps (scaled differently).

as depicted in Fig. 25, the estimated errors show a monotonic decrease by increasing the total number of the DOFs, which indicates the stable convergence of the refinement process. In order to quantify this comparison and evaluate the accuracy of error estimate, the goal error effectivity index, the ratio between the two curves of Fig. 25, is calculated and plotted in Fig. 26(a). In this example, similar to Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

37

Fig. 25. Estimated and exact errors of the problem 6.4 for DWR method.

0

Fig. 26. The standard measures of the DWR goal error estimator for the problem 6.4: (a) Goal error effectivity index (θε ); (b) Effectivity in the quantity of interest (θ QoI ).

the 2D case, the DWR method underestimates the true goal error and after the third refinement step θε lies in the acceptable range of 0.9–1.0. Fig. 26(a) denotes that further mesh refinements lead the error effectivity index to the ideal value of 1. Furthermore, θε fluctuations are reduced while approaching to the ideal value. Moreover, as depicted in Fig. 26(b), θ QoI also keeps heading to the desired value of 1 by further refinements which also implies the accuracy of the proposed method. Similar to the previous examples, comparing the DWR results with Kelly, weighted Kelly (with (r 2 + 0.12 )−1 weights), and global uniform refinement methods indicates the optimal convergence rate and better efficiency of the presented DWR method. As shown in Fig. 27, the convergence rate of the DWR method is 0.752. 7. Conclusion In this paper we presented an adaptive thermo-mechanical finite element formulation based on dual weighted residual (DWR) based goal-oriented error estimation (GOEE). The DWR a posteriori error estimation is employed to estimate the singular pointwise quantity of interest in 2D and 3D problems and the results have been evaluated by comparing with analytical solution. In investigated examples, two singular goals of point value and point value derivative are analyzed and meshes are refined adaptively to minimize the error in the quantity of interest. In the refinement process, local h-refinement allowing hanging node in 2D and 3D is adopted to enhance the accuracy of analysis and quantity of interest. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

38

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 27. Errors of the problem 6.4 estimated by different refinement techniques.

The method has been applied to several numerical problems with analytical solutions and the estimated errors are compared with the exact errors. These comparisons show the very good accordance of the estimated and true goal errors even for singular pointwise quantities of interest in 2D and 3D problems. In order to quantify these comparisons, the goal error effectivity index is calculated. In all examples, after few initial coarse meshes, the goal error effectivity index lies in the acceptable range of 0.9–1.0 and by further refinements approaches to the ideal value of 1. In estimating the point value goal error, the employed DWR method in the initial steps of the refinement process, which the meshes are still coarse, overestimates the error but with further refinements the exact error is underestimated. In the examples with point value derivative goal, the DWR method underestimated the actual goal error. The evolution of error effectivity index in problems with point value goal is smoother than the problems with point value derivation goals. Nevertheless, in examples with point value derivative as the quantity of interest, by refining the mesh, the irregular behavior and fluctuations of the goal error effectivity index are decreased and the mean value tends to the desired value of 1. Moreover, in each example, apart from the analytical solution, the proposed goal-oriented mesh adaptivity process have been compared with global uniform, Kelly and weighted Kelly (W-Kelly) refinement techniques. This comparison shows that among the selected error estimation and refinement techniques, the employed DWR method is the most optimum and efficient one and the global uniform refinement has the least convergence rate. On average, in the studied examples, the DWR, Kelly, and W-Kelly methods show, respectively, 0.68%, 0.42%, and 0.26% higher convergence rates compared to the global uniform refinement. This efficiency enhancement leads to the considerable computational time saving; for example in the problem 6.1, by implementing the proposed DWR method it is possible to reach the error tolerance of 10−6 by a mesh with only about 2% DOFs of a globally uniform refined mesh and this difference grows exponentially by demanding higher estimation accuracies. Therefore, we can sum up that implementing the developed DWR method into adaptive finite element algorithm leads to economical computational meshes. In conclusion, in all examples: the DWR errors show a monotonic decrease by increasing the total number of the DOFs, which indicates the stable convergence of the refinement process; the true and estimated errors are actually very close to each other; the evolution of the goal error effectivity index shows the very good accuracy of the implemented DWR methods; and the higher convergence rate of the employed method compared to the other methods prove the efficiency of the employed refinement technique. Acknowledgment The authors extend their appreciation to the Distinguished Scientist Fellowship Program (DSFP) at King Saud University, Saudi Arabia for funding this work. Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

39

Appendix A.1. Analytical solution of a thick-walled cylinder subjected to internal pressure and temperature change In order to extract the exact solution of a thick-walled cylinder under thermo-mechanical loading in 2D, a general case of a very long cylinder with internal pressure and temperature of Pi and Ti and external pressure and temperature of Po and To is considered. Due to the axisymmetric conditions, the strain component ϵr θ is zero. For this case, the strain–displacement and stress–strain relations considering plane strain conditions are [3] du r (r ) 1 = [σrr − ν(σθ θ + σzz )] + α T θ T , dr E u r (r ) 1 ϵθθ := = [σθ θ − ν(σrr + σzz )] + α T θ T , r E σzz = ν (σrr + σθ θ ) − Eα T θ T , ϵrr :=

(105) (106) (107)

where ϵi j and σi j are strain and stress tensor components in cylindrical coordinates, and i, j ∈ {r, √ θ, z} which are radial, azimuthal and axial directions in cylindrical coordinate system, respectively, with r = x 2 + y 2 and θ = arctan(y/x). In above relations θ T (r ) = T (r ) − T0 refers to the temperature change, which for the model studied here is [3] ln( rro ) ∆θ T ( ro ) ∆θ = θ + ln , (108) θ T (r ) = θo + o T ln( rroi ) ln(c) r where c = ro /ri , ∆θ T = θi − θo , θi = Ti − T0 , and θo = To − T0 . Following the approach stated in [3], the radial displacement and stress tensor components of exact solution of the thick-walled cylinder under thermo-mechanical loadings, we obtain [( ( )) ( 2 ) ] 1 − ν + ln rro 1 − 2ν θo ro 1 1+ν + 2(1 − ν) r+ α T ∆θ T u r (r ) = + 2(1 − ν) 1 − c2 ∆θ T ln(c) 1 − c2 r [( ) ( ) ] 1 + ν ro2 1 r2 r2 + 1 + (1 − 2ν) 2 Po − 1 + (1 − 2ν) 2 Pi , (109) E 1 − c2 r ro ri ] [( ) ( ) ( 2 ) −ν + ln rro θo 1 1+ν 1 − 2ν ro + 2(1 − ν) α T ∆θ T ϵrr (r ) = + − 2(1 − ν) 1 − c2 ∆θ T ln(c) 1 − c2 r 2 [( ) ( ) ] 1 1 + ν ro2 1 1 1 − + (1 − 2ν) P − (1 − 2ν) − Pi , (110) o E 1 − c2 ro2 r 2 ri2 r 2 [ ( )] ( ) ln rro 1 E ro2 σrr (r ) = 1− 2 − α T ∆θ T 2(1 − ν) 1 − c2 r ln(c) [( ) ( ) ] ro2 1 1 1 1 + − Po − 2 − 2 Pi , (111) 1 − c2 ro r r2 r2 [i ] ( ) ( ) ln rro 1 ro2 1 E 1+ 2 − + α T ∆θ T σθθ (r ) = 2(1 − ν) 1 − c2 r ln(c) ln(c) [( ) ( ) ] ro2 1 1 1 1 + + P − + Pi . (112) o 1 − c2 ro2 r 2 ri2 r 2 Considering the vector and tensor transformation from cylindrical coordinate system to the Cartesian one, the displacement in x direction and its derivative with respect to x, respectively, are u x (r, θ) = u r (r ) cos(θ ), ∂u x (r, θ) ϵx x (r, θ) := = ϵrr (r ) cos2 (θ ) + ϵθ θ (r ) sin2 (θ ). ∂x

(113) (114)

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

40

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

A.2. Analytical solution of a thick-walled sphere subjected to internal pressure and temperature change Consider a general case of a thick-walled sphere subjected to thermo-mechanical loading with constant pressure and temperature of Pi and Ti applied on the inner surface of r = ri and pressure and temperature of Po and To applied on the outer surface of r = ro . Due to the symmetrical loading and geometry in both polar θ and azimuthal φ directions the strain components of ϵr θ , ϵθ φ , and ϵr φ are zero and ϵθ θ = ϵφφ . For this case, the strain–displacement and stress–strain relations are [3,66] du r (r ) 1 ϵrr := = [σrr − 2νσθ θ ] + α T θ T , (115) dr E 1 u r (r ) = [(1 − ν)σθ θ − νσrr ] + α T θ T , ϵθθ = ϵφφ := (116) r E where ϵi j and σi j are strain and stress tensor components in spherical coordinates, and i, j ∈ √ {r, θ, φ} which are radial, polar, and azimuthal directions in spherical coordinate system, respectively, with r = x 2 + y 2 + z 2 , θ = arccos(z/r ), and φ = arctan(y/x). In above relations θ T (r ) = T (r ) − T0 refers to the temperature change, which for the model studied here is [3] 1 − ( rro ) ∆θ T ( ro ) θ T (r ) = θo − 1− , (117) ro ∆θ T = θo − 1 − ( ri ) 1−c r where c = ro /ri , ∆θ T = θi − θo , θi = Ti − T0 , and θo = To − T0 . By solving the Navier equation (8) in the spherical coordinate system, considering the temperature distribution (117) and thermal/mechanical boundary conditions, we obtain [3,83] ] [ ri2 ro2 ri r o 2ν(ri + ro ) 1+ν 2 2 u r (r ) = r α T ∆θ T r + ri r o + ro − 2 − 2(1 − ν) ro3 − ri3 i r 1+ν + ϵrr (r ) = + σrr (r ) = − σθθ (r ) = −

α T (ro3 θo − ri3 θi ) 1 + ν ri3 ro3 (Po − Pi ) 1 − 2ν ro3 Po − ri3 Pi − − r, 2E E ro3 − ri3 r 2 (ro3 − ri3 ) ro3 − ri3 ] [ α T (ro3 θo − ri3 θi ) 1 + ν ri ro ν(ri + ro ) ri2 ro2 + α ∆θ + − T T 3 3 1 − ν ro3 − ri 1+ν r ro3 − ri3 1 + ν ri3 ro3 (Po − Pi ) 1 − 2ν ro3 (Po − ri3 Pi ) − , E E r 3 (ro3 − ri3 ) ro3 − ri3 ] [ ri ro E 2 2 1 2 2 1 (ri + ro ) − (ri + ri ro + ro ) + (ri + ro ) 3 α T ∆θ T 1 − ν ro3 − ri3 r r ( ( 3) 3 3 3) r r ro r 1 − i3 Po + 3 i 3 1 − o3 Pi , 3 3 r r r o − ri r o − ri [ ] E ri r o 2 2 1 2 2 1 2(ri + ro ) − (ri + ri ro + ro ) − (ri + ro ) 3 α T ∆θ T 2(1 − ν) ro3 − ri3 r r ( ( ) ) ri3 ri3 ro3 ro3 1 + 3 Po + 3 1 + 3 Pi . 2r 2r ro3 − ri3 ro − ri3

(118)

(119)

(120)

(121)

Considering the vector and tensor transformation from spherical coordinate system to the Cartesian one, the displacement in x direction and its derivative with respect to x, respectively, are u x (r, θ, φ) = u r (r ) cos(φ) sin(θ), ∂u x (r, θ, φ) ϵx x (r, θ, φ) := = (ϵrr (r ) − ϵθ θ (r )) cos2 (φ) sin2 (θ ) + ϵθ θ (r ). ∂x

(122) (123)

References [1] Witold Nowacki, Dynamic Problems of Thermoelasticity (in Polish), PWN-Polish Scientific Publishers, Warszawa, 1966. [2] Witold Nowacki, Dynamic Problems of Thermoelasticity (English Edition), Noordhoff International Publishing, Leyden, 1975. [3] Richard B. Hetnarski, M. Reza Eslami, Thermal stresses - Advanced theory and applications, Solid Mech. Appl. 158 (2009) 1–591.

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

41

[4] P. Hosseini-Tehrani, M.R. Eslami, Boundary element analysis of finite domains under thermal and mechanical shock with the lord-shulman theory, J. Strain Anal. Eng. Des. 38 (1) (2003) 53–64. [5] Raymond D. Mindlin, David H. Cheng, Thermoelastic stress in the semi-infinite solid, J. Appl. Phys. 21 (1950) 931–933. [6] Eli Sternberg, E.L. McDowell, On the steady-state thermoelastic problem for the half-space, Quart. Appl. Math. 14 (1955) 381–391, A technical report to the Office of Naval Research, Department of the Navy. [7] E. Melan, H. Parkus, Wärmespannungen infolge stationärer Temperaturfelder, Springer, 1953. [8] W. Nowacki, A three-dimensional thermoelastic problem with discontinuous boundary conditions, Arch. Mech. Stos. 9 (1957) 319–324. [9] I.N. Sneddon, F.J. Lockett, On the steady-state thermoelastic problem for the half-space and the thick plate, Quart. J. Appl. Math. 18 (1960) 145–153. [10] P. Chadwick, Thermoelasticity, the dynamical theory, Prog. Solid Mech. 1 (1960) 263–328. [11] A. Boley Bruno, Thermal stresses, in: Structural Mechanics. Proceedings of 1st Symposium on Naval Structural Mechanics, Pergamon Press, New York, 1960, pp. 378–406. [12] Bruno A. Boley, Jerome H Weiner, Theory of Thermal Stresses, John Wiley and Sons, New York, 1960. [13] Witold Nowacki, Thermoelasticity, Addison-Wesley Publishing Co., Reading, Massachusetts, 1962. [14] Witold Nowacki, Thermoelasticity, second ed., Pergamon Press, PWN-Polish Scientific Publishers, Oxford, Warszawa, 1962. Pergamon Press, PWN-Polish Scientific Publishers, Oxford, Warsaw, 1986. [15] H. Parkus, Thermoelasticity, Blaisdell Publishing Co., Waltham, Mass, 1968. [16] P. Phung-Van, Cuong-Le Thanh, H. Nguyen-Xuan, M. Abdel-Wahab, Nonlinear transient isogeometric analysis of FG-CNTRC nanoplates in thermal environments, Compos. Struct. 201 (2018) 882–892. [17] Junyan Ni, Xiaowei Wang, Jianming Gong, Magd Abdel Wahab, Thermal, metallurgical and mechanical analysis of circumferentially multi-pass welded P92 steel pipes, Int. J. Press. Vessels Pip. 162 (2018) 164–175. [18] Loc V. Tran, M. Abdel Wahab, Seung-Eock Kim, An isogeometric finite element approach for thermal bending and buckling analyses of laminated composite plates, Compos. Struct. 179 (2018) 35–49. [19] Junyan Ni, Xiaowei Wang, Jianming Gong, Magd Abdel Wahab, A multi-phase model for transformation plasticity using thermodynamics-based metallurgical algorithm, Int. J. Mech. Sci. 148 (2018) 135–148. [20] Loc V. Tran, P. Phung-Van, Jaehong Lee, M. Abdel Wahab, H. Nguyen-Xuan, Isogeometric analysis for nonlinear thermomechanical stability of functionally graded plates, Compos. Struct. 140 (2016) 655–667. [21] P. Phung-Van, Loc V. Tran, A.J.M. Ferreira, H. Nguyen-Xuan, M. Abdel-Wahab, Nonlinear transient isogeometric analysis of smart piezoelectric functionally graded material plates based on generalized shear deformation theory under thermo-electro-mechanical loads, Nonlinear Dynam. 87 (2) (2017) 879–894. [22] Mark Ainsworth, J. Tinsley Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley, New York, 2000, ID: 761321857. [23] Ivo Babuška, W.C. Rheinboldt, A-posteriori error estimates for the finite element method, Internat. J. Numer. Methods Engrg. 12 (10) (1978) 1597–1615. [24] T. Grätsch, K.-J. Bathe, A posteriori error estimation techniques in practical finite element analysis, Comput. Struct. 83 (4) (2005) 235–265. [25] R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Wiley-Teubner, Chichester-Stuttgart, 1996. [26] Weimin Han, A Posteriori Error Analysis Via Duality Theory : with Applications in Modeling and Numerical Approximations, Springer, New York, 2005. [27] O.C. Zienkiewicz, J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineerng analysis, Internat. J. Numer. Methods Engrg. 24 (1987) 337–357. [28] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique, Internat. J. Numer. Methods Engrg. 33 (1992) 1331–1364. [29] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity, Internat. J. Numer. Methods Engrg. 33 (1992) 1365–1382. [30] Nils-Erik Wiberg, Fethi Abdulwahab, Saulius Ziukas, Enhanced superconvergent patch recovery incorporating equilibrium and boundary conditions, Internat. J. Numer. Methods Engrg. 37 (20) (1994) 3417–3440. [31] T. Blacker, T. Belytschko, Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements*, Internat. J. Numer. Methods Engrg. 37 (1994) 517–536. [32] J.J. Ródenas, M. Tur, F.J. Fuenmayor, A. Vercher, Improvement of the superconvergent patch recovery technique by the use of constraint equations: The SPR-C technique, Internat. J. Numer. Methods Engrg. 70 (2007) 705–727. [33] J.J. Ródenas, O.A. González-Estrada, J.E. Tarancón, F.J. Fuenmayor, A recovery-type error estimator for the extended finite element method based on singular+smooth stress field splitting, Internat. J. Numer. Methods Engrg. 76 (4) (2008) 545–571. [34] J.J. Ródenas, O.A. González-Estrada, P. Díez, F.J. Fuenmayor, Accurate recovery-based upper error bounds for the extended finite element framework, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2607–2621. [35] H. Moslemi, A.R. Khoei, 3D adaptive finite element modeling of non-planar curved crack growth using the weighted superconvergent patch recovery method, Eng. Fract. Mech. 76 (11) (2009) 1703–1728. [36] L. Demkowicz, J.T. Oden, T. Strouboulis, Adaptive finite elements for flow problems with moving boundaries. part i: Variational principles and a posteriori estimates, Comput. Methods Appl. Mech. Engrg. 46 (2) (1984) 217–251. [37] L. Demkowicz, J.T. Oden, T. Strouboulis, An adaptive p-version finite element method for transient flow problems with moving boundaries, Finite Elem. Fluids 6 (1985) 291–305, Chichester, England and New York, Wiley-Interscience. [38] R.E. Bank, Analysis of a local a posteriori error estimate for elliptic equations, in: Accuracy Estimates and Adaptive Refinements in Finite Element Computations, John Wiley & Sons Ltd., 1986, pp. 119–128.

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

42

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

[39] R.E. Bank, A. Weiser, Some a posteriori error estimators for elliptic partial differential equations, in: Mathematics of Computation, Vol. 44, 1985, p. 283. [40] Mark Ainsworth, J. Tinsley Oden, A procedure for a posteriori error estimation for h-p finite element methods, Comput. Methods Appl. Mech. Engrg. 101 (1992) 73–96. [41] Mark Ainsworth, J. Tinsley Oden, A unified approach to a posteriori error estimation using element residual methods, Numer. Math. 65 (1993) 23–50. [42] J.T. Oden, L. Demkowicz, W. Rachowicz, T.A. Westerman, Toward a universal hp-adaptive finite element strategy. Part 2: A posteriori error estimation, Comput. Methods Appl. Mech. Engrg. 77 (1989) 113–180. [43] Pedro Díez, Núria Parés Mariné, Antonio Huerta, Recovering lower bounds of the error by postprocessing implicit residual a posteriori error estimates, Internat. J. Numer. Methods Engrg. 56 (10) (2003) 1465–1488. [44] Pedro Díez, Núria Parés Mariné, Antonio Huerta, Accurate upper and lower error bounds by solving flux-free local problems in stars, Rev. Eur. Elem. Finis 13 (5-6-7) (2004) 497–507. [45] T. Gerasimov, M. Rüter, E. Stein, An explicit residual-type error estimator for Q1-quadrilateral extended finite element method in two-dimensional linear elastic fracture mechanics, Internat. J. Numer. Methods Engrg. 90 (April) (2012) 1118–1155. [46] C. Carstensen, S.A. Funken, Constants in Clément-interpolation error and residual-based a posteriori estimates in finite element methods, East-West J. Math. 8 (3) (2000) 153–175. [47] C. Johnson, P. Hansbo, Adaptive finite element methods in computational mechanics, Comput. Methods Appl. Mech. Engrg. 101 (1992) 143–181. [48] E. Stein, T. Gerasimov, M. RüÌ´Lter, Error-controlled adaptive multiscale analysis for crack initiation and propagation in brittle materials, Adapt. Model. Simul. 26 (2013). [49] T. Gerasimov, E. Stein, P. Wriggers, Constant-free explicit error estimator with sharp upper error bound property for adaptive FE analysis in elasticity and fracture, Internat. J. Numer. Methods Engrg. (2014) pages on–line. [50] S. Prudhomme, J.T. Oden, On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors, Comput. Methods Appl. Mech. Engrg. 176 (1–4) (1999) 313–331. [51] J.T. Oden, S. Prudhomme, Goal-oriented error estimation and adaptivity for the finite element method, Comput. Math. Appl. 41 (2001) 735–756. [52] S. Korotov, P. NeittaanmaÌ´Lki, S. Repin, A posteriori error estimation of goal-oriented quantities by the superconvergence patch recovery, J. Numer. Math. 11 (1) (2003) 33–59. [53] S. Korotov, A posteriori error estimation of goal-oriented quantities for elliptic type BVPs, J. Comput. Appl. Math. 191 (2) (2006) 216–227. [54] O.A. González-Estrada, J.J. Ródenas, Stéphane P.A. Bordas, E. Nadal, F.J. Fuenmayor, Locally equilibrated superconvergent patch recovery for efficient error estimation in quantities of interest, in: 6th Eur. Congr. Comput. Methods Appl. Sci. Eng. ECCOMAS 2012, 2012, pp. 1–28. [55] Ehsan Rabizadeh, Amir Saboor Bagherzadeh, Timon Rabczuk, Adaptive thermo-mechanical finite element formulation based on goal-oriented error estimation, Comput. Mater. Sci. 102 (2015) 27–44. [56] Ehsan Rabizadeh, Amir Saboor Bagherzadeh, Timon Rabczuk, Goal-oriented error estimation and adaptive mesh refinement in dynamic coupled thermoelasticity, Comput. Struct. 173 (2016) 187–211. [57] Roland Becker, Rolf Rannacher, A feed-back approach to error control in finite element methods: basic analysis and examples, East-West J. Numer. Math. 4 (4) (1996) 237–264. [58] Roland Becker, Rolf Rannacher, Weighted a posteriori error control in fe methods, in: ENUMATH 97, 1998, pp. 621–637. [59] Roland Becker, Rolf Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numer. 10 (2001) 1–102. [60] R. Rannacher, Adaptive finite element methods for partial differential equations, 2003, ArXiv Mathematics e-prints. [61] P.R. Budarapu, R. Gracie, S.W. Yang, X. Zhaung, T. Rabczuk, Efficient coarse graining in multiscale modeling of fracture, Theor. Appl. Fract. Mech. 69 (2014) 126–143. [62] S.W. Yang, P.R. Budarapu, Mahapatra D.R., Bordas S.P.A., G. Zi, T. Rabczuk, A meshless adaptive multiscale method for fracture, Comput. Mater. Sci. 96 (Part B) (2015) 382–395. [63] P.R. Budarapu, J. Reinoso, M. Paggi, Concurrently coupled solid shell-based adaptive multiscale method for fracture, Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365, https://www.sciencedirect.com/science/article/pii/S0045782516313883. [64] W. Bangerth, R. Hartmann, G. Kanschat, deal. II – a general purpose object oriented finite element library, ACM Trans. Math. Softw. (TOMS) 33 (4) (2007) 24/1–24/27. [65] D. Arndt, W. Bangerth, D. Davydov, T. Heister, L. Heltai, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin, D. Wells, The deal. II library, version 8.5, J. Numer. Math. 25 (3) (2017) 137–146. [66] Allan F. Bower, Applied Mechanics of Solids, first ed., CRC Press, 2009. [67] Sanjay Bagade, N.W. Khobragade, Thermoelastic Problems in Physics and Engineering: Mathematical Approach to Thermoelasticity, LAP LAMBERT Academic Publishing, 2012. [68] William Alan Day, Heat Conduction Within Linear Thermoelasticity, Springer, Cop., Berlin, 1985. [69] Dorin IeSan, ¸ Thermoelastic Models of Continua, Springer, London, 2011. [70] Jean-Marie Constant Duhamel, Sur les équations générales de la propagation de la chaleur dans les corps solides dont la conductibilité n’est pas la même dans tous les sens, J. Ec. Polytech. Paris 13 (21) (1832) 356–399. [71] M.R. Eslami, A First Course in Finite Element Analysis, first ed., Amirkabir University Press, Tehran, 2003. [72] K.J. Bathe, Finite Element Procedures in Engineering Analysis, in: Prentice-Hall civil engineering and engineering mechanics series, Prentice-Hall, 1982.

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.

E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

43

[73] Wolfgang Bangerth, Rolf Rannacher, Adaptive Finite Element Methods for Differential Equations, Birkhäuser Verlag, Basel Boston, 2003. [74] Donald Estep, Michael Holst, Mats Larson, Generalized green’s functions and the effective domain of influence, SIAM J. Sci. Comput. 26 (4) (2005) 1314–1339. [75] F.J. Fuenmayor, J.L. Oliver, Criteria to achieve nearly optimal meshes in the h-adaptive finite element mehod, Internat. J. Numer. Methods Engrg. 39 (1996) 4039–4061. [76] Dominique Chapelle, Klaus-Jürgen Bathe, The Finite Element Analysis of Shells : Fundamentals, in: Computational Fluid and Solid Mechanics, Springer, Berlin, 2011. [77] John Tinsley Oden, Junuthula Narasimha Reddy, An Introduction to the Mathematical Theory of Finite Elements, John Wiley & Sons, New York, 1976. [78] D.W. Kelly, J.P. De S.R. Gago, O.C. Zienkiewicz, I. Babuska, A posteriori error analysis and adaptive processes in the finite element method: Part I—error analysis, Internat. J. Numer. Methods Engrg. 19 (11) (1983) 1593–1619. [79] Randolph E. Bank, The efficient implementation of local mesh refinement algorithms, Adapt. Comput. Methods Partial Differ. Equ. (1983) 74–81. [80] Philip D. Harvey (Ed.), American society for metals, in: Engineering Properties of Steel, American Society for Metals, Metals Park, Ohio, 1982. [81] Donald Peckner, I.M. Bernstein, Handbook of Stainless Steels, McGraw-Hill, New York, 1977. [82] Mark Ainsworth, J. Tinsley Oden, A posteriori error estimation in finite element analysis, Comput. Methods Appl. Mech. Eng. 142 (1–2) (1997) 1–88. [83] William S. Slaughter, The Linearized Theory of Elasticity, Springer Science + Business Media, LLC, New York, 2013.

Please cite this article as: E. Rabizadeh, A.S. Bagherzadeh, C. Anitescu et al., Pointwise dual weighted residual based goal-oriented a posteriori error estimation and adaptive mesh refinement in 2D/3D thermo-mechanical multifield problems, Computer Methods in Applied Mechanics and Engineering (2019) 112666, https://doi.org/10.1016/j.cma.2019.112666.