Numerical analysis of a problem in micropolar thermoviscoelasticity

1.   Introduction

It seems that Voigt ^[1] was the first researcher who noted the polar structure of several materials as crystalline solids. Later, Cosserat brothers ^[2] proposed a theory of micropolar materials. In such theory, each material point has six degrees of freedom: three for the displacements and three for the orientation usually called directions. After some time without relevant contributions in this field, in 1964 Eringen and Suhubi ^[3,4] studied the nonlinear theory of micropolar elasticity which took into account the microstructural motions, and later Eringen ^[5] developed its linear version. A large number of engineering materials may be modeled more realistically using these micropolar continua, where the classical theory of elasticity is inadequate. Some particular examples can be, for instance, the large-scale development and utilization of composite, reinforced and coarse-grained materials. The basic idea of micropolar elasticity is the fact that the granular character of the medium is taken into account, describing deformations by a microrotation. Using this theory, solids support couple stresses in addition to force stresses. Nowacki ^[6], Tauchert et al.^[7], Tauchert ^[8] and other researchers extended the theory by introducing the thermal effects (see also ^[9,10,11] for a good review on micropolar themoelasticity). Since then, a large number of papers have been published dealing with mathematical issues of this type of materials (see ^{[12,13,14,15,16,17,18,19]}, among others ^{[20,21,22,23,24,25]}).

In this work, we continue the research started in ^[26], where isotropic and homogeneous micropolar thermoviscoelastic materials were considered. The existence and uniqueness of solutions were shown by using the theory of linear semigroups. The analyticity of the solutions was also obtained supposing that the supply terms were absent. Here, our aim is to study this problem from the numerical point of view, introducing a fully discrete approximation, proving a discrete stability property and a priori error estimates, and performing some two-dimensional numerical results.

The paper is structured in the following way. The thermomechanical model and its weak form are presented in Section 2. We recall an existence and uniqueness result obtained by Magaña and Quintanilla ^[26]. Then, in Section 3 fully discrete approximations are introduced, based on the finite element method for the spatial approximation and the backward Euler scheme to discretize the time derivatives. A discrete stability property is shown and a priori error estimates are provided, from which, under some regularity conditions, the linear convergence is derived. Finally, two-dimensional numerical examples are shown in Section 4.

2.   The micropolar thermoviscoelastic model

Let us denote by $\Omega\subset \mathbb{R}^d$, $d = 2, 3$, a domain with boundary $\Gamma = \partial \Omega$ assumed to be Lipschitz, and occupied by our thermoviscoelastic material with micropolar properties. As usual, the spatial variable is represented by ${\boldsymbol{x}}\in \Omega$, and for the time we use the notation $t\in [0, T_f],$ $T_f > 0$ being the final time. To save cumbersome expressions, we will remove the dependence of our functions with respect to these two variables. Moreover, a subscript after a comma denotes the partial derivative with respect to the corresponding variable (i.e., $f_{i, j} = \frac{\partial f_i}{\partial x_j}$), and the time derivative is denoted by a superposed dot for the first order or two superposed dots for the second order. Repeated indices mean summation over those indices.

Let ${\boldsymbol{u}} = (u_i)_{i = 1}^d$, ${\boldsymbol{\phi}} = (\phi_i)_{i = 1}^d$ and $T$ be the displacement field, the microrotation vector and the temperature, respectively.

Since we assume that the material is homogeneous and isotropic, the micropolar thermoviscoelastic problem is written as follows (see ^[26] for details).

Problem P. Find the displacement field ${\boldsymbol{u}}:\Omega\times [0, T_f]\rightarrow \mathbb{R}^d$, the microrotation vector ${\boldsymbol{\phi}}:\Omega\times [0, T_f]\rightarrow \mathbb{R}^d$ and the temperature $T:\Omega\times [0, T_f]\rightarrow \mathbb{R}$ such that, for $i = 1, \ldots, d$,

In Eqs (2.1)–(2.3), $\Delta$ denotes the classical Laplace operator, $\epsilon_{ijk}$ represents the classical alternating symbol, $\rho$, $J$, $a$, $\mu$, $\sigma$, $\lambda$, $\mu_v$, $\sigma_v$, $\lambda_v$, $b$, $\gamma$, $\alpha$, $\beta$, $b^*$, $\gamma_v$, $\alpha_v$, $\beta_v$, $\kappa$ and $\kappa^*$ are constitutive coefficients and ${\boldsymbol{F}}^1 = (F_i^1)_{i = 1}^d$, ${\boldsymbol{F}}^2 = (F_i^2)_{i = 1}^d$ and $R$ denote the supply terms. Moreover, in order to simplify the writing, from now on we assume that $T_0 = 1$.

We note that in the case that $\sigma,$ $\sigma_{\nu}$, $\gamma$, $\alpha$, $\beta$, $\gamma_{\nu}$, $\alpha_{\nu}$, $\beta_{\nu}$, $b^*$ an $\kappa^*$ vanish we recover the classical theory of thermoviscoelasticity (see, for instance, ^[27,28,29]).

As it is pointed out in ^[26], the mechanical internal energy and the dissipation are given by

where tensors ${\boldsymbol{e}} = (e_{ij})_{i, j = 1}^d$ and $\dot{ {\boldsymbol{e}}} = (\dot{e}_{ij})_{i, j = 1}^d$ are defined as

Since we must impose that both functions are positive definite, we will assume the following conditions on the constitutive coefficients:

We note that conditions (2.7) are slightly stronger than those required in ^[26]; however, we have modified them for the sake of simplicity in the calculations of the next section. See Appendix A for the derivation of the relations of the viscous coefficients.

The aim of this paper is to study numerically Problem P and so, we first introduce its variational formulation. Therefore, let us define the spaces $Y = L^2(\Omega)$, $H = [L^2(\Omega)]^d$ and $Q = [L^2(\Omega)]^{d\times d}$, with their respective scalar product $(\cdot, \cdot)_Y$, $(\cdot, \cdot)_H$ and $(\cdot, \cdot)_Q$, being the corresponding norms $\|\cdot\|_Y$, $\|\cdot\|_H$ and $\|\cdot\|_Q$. We also consider the variational spaces $E$ and $V$ as follows,

The corresponding scalar products and norms are denoted as $(\cdot, \cdot)_E$ and $(\cdot, \cdot)_V$, and $\|\cdot\|_E$ and $\|\cdot\|_V$, respectively.

Using the classical procedure and the prescribed boundary conditions, we arrive to the variational form of Problem P written in terms of the temperature $T$, the velocity field ${\boldsymbol{v}} = (v_i)_{i = 1}^d = \dot{ {\boldsymbol{u}}} = (\dot{u}_i)_{i = 1}^d$ and the microrotation speed ${\boldsymbol{\psi}} = (\psi_i)_{i = 1}^d = \dot{ {\boldsymbol{\phi}}} = (\dot{\phi}_i)_{i = 1}^d$.

Problem VP. Find the temperature $T: [0, T_f]\rightarrow E$, the velocity field ${\boldsymbol{v}}: [0, T_f]\rightarrow V$ and the microrotation speed ${\boldsymbol{\psi}}: [0, T_f]\rightarrow V$ such that $T(0) = T^0$, ${\boldsymbol{v}}(0) = {\boldsymbol{v}}^0$, ${\boldsymbol{\psi}}(0) = {\boldsymbol{\psi}}^0$, and, for a.e., $t\in (0, T_f)$ and for all ${\boldsymbol{w}}, {\boldsymbol{\xi}}\in V, \, \zeta\in E$,

where $\nabla$ and ${\rm{div}}\;$ represent the usual gradient and divergence operators, and the displacement field ${\boldsymbol{u}}$ and the microrotation ${\boldsymbol{\phi}}$ are obtained from the respective equations:

In Problem VP we have used the following bilinear operators to simplify the writing:

We note that, using conditions (2.7), the dissipation $D$ is positive definite and so, it leads to the following property that we will use in the next section:

where $M^*$ is a given positive constant.

The following well posedness result was proved in ^[26].

Theorem 1. Let conditions (2.7) hold. Then, there exists a unique solution to Problem VP with the following regularity:

3.   Finite element approximations: An a priori error analysis

In this section we will study, from the numerical point of view, a fully discrete approximation of Problem VP. First, in order to provide the spatial approximation, we construct the finite element spaces $E^h$ and $V^h$ given by

where, as usual, we assume that $\overline{\Omega}$ is a polyhedral domain and we denote by ${ \mathcal{T}}^h$ a regular triangulation of $\overline{\Omega}$ in the sense of ^[30]. $P_1(Tr)$ represents the space of linear polynomials in $Tr$ and parameter $h > 0$ is the spatial discretization size.

The discretization of the time derivatives is obtained with a uniform partition of the time interval $[0, T_f]$, denoted by $0 = t_0 < t_1 < \ldots < t_N = T_f$, being $k = T_f/N$ the time step size. As a usual notation, if $f$ is a continuous function we denote $f_n = f(t_n)$, and for the sequence $\{z_n\}_{n = 0}^N$, let $\delta z_n = (z_n-z_{n-1})/k$ be its divided differences.

Applying the well-known implicit Euler scheme, we have the fully discrete approximation of Problem VP as follows.

Problem VP$^{hk}$. Find the discrete temperature $T^{hk} = \{T_n^{hk}\}_{n = 0}^N\subset E^h$, the discrete velocity field ${\boldsymbol{v}}^{hk} = \{ {\boldsymbol{v}}_n^{hk}\}_{n = 0}^N\subset V^h$ and the discrete microrotation speed ${\boldsymbol{\psi}}^{hk} = \{ {\boldsymbol{\psi}}_n^{hk}\}_{n = 0}^N\subset V^h$ such that $T_0^{hk} = T^h_0$, ${\boldsymbol{v}}_0^{hk} = {\boldsymbol{v}}^h_0$, ${\boldsymbol{\psi}}_0^{hk} = {\boldsymbol{\psi}}^h_0$, and, for $n = 1, \ldots, N$ and for all ${\boldsymbol{w}}^h, {\boldsymbol{\xi}}^h\in V^h, \, \zeta^h\in E^h$,

where the discrete displacement field ${\boldsymbol{u}}^{hk}$ and the discrete microrotation vector ${\boldsymbol{\phi}}^{hk}$ are updated from the respective equations:

In the previous fully discrete problem, the discrete initial conditions $T_0^{h}$, ${\boldsymbol{u}}_0^{h}$, ${\boldsymbol{v}}_0^{h}$, ${\boldsymbol{\phi}}_0^{h}$ and ${\boldsymbol{\psi}}_0^{h}$ are approximations of the initial conditions $T^{0}$, ${\boldsymbol{u}}^{0}$, ${\boldsymbol{v}}^{0}$, ${\boldsymbol{\phi}}^{0}$ and ${\boldsymbol{\psi}}^{0}$ defined as

Here, we denote by $\mathcal{P}_1^{h}$ and $\mathcal{P}_2^{h}$ the interpolation operators over the finite element spaces $V^h$ and $E^h$, respectively (see, e.g., ^[30]).

Using Lax-Milgram lemma and assumptions (2.7), it is straightforward to show the existence of a unique discrete solution to Problem VP$^{hk}$.

We have the following discrete stability property.

Theorem 2. Under the assumptions of Theorem 1, then it follows that the sequences $\{ {\boldsymbol{u}}^{hk}, {\boldsymbol{v}}^{hk}, {\boldsymbol{\phi}}^{hk}, {\boldsymbol{\psi}}^{hk}, T^{hk}\}$, generated by Problem $VP^{hk}$, satisfy the stability estimate:

where $M$ is a positive constant which is independent of the discretization parameters $h$ and $k$.

Proof. Taking ${\boldsymbol{w}}^h = {\boldsymbol{v}}_n^{hk}$ as a test function in discrete variational Eq (3.3) we have

Thus, keeping in mind that

where, from now on, $\epsilon > 0$ is assumed small enough and we have used Cauchy-Schwarz inequality and the arithmetic-geometric mean inequality

we find that

Proceeding in a similar form with the terms of the discrete microrotation speed and also keeping in mind that

we have

Finally, we obtain the estimates for the discrete temperature $T_n^{hk}$. Taking $\zeta^h = T_n^{hk}$ as a test function in discrete variational Eq (3.5) it follows that

and taking into account that

we obtain

Combining the previous estimates and taking into account property (2.12), we have

Multiplying them by $k$ and summing up the resulting equation it follows that

Finally, the stability estimates are a direct consequence of the application of a discrete version of Gronwall's inequality (see, e.g., ^[31]), the properties of the interpolation operators $P_1^h$ and $P_2^h$ (see ^[30]) and the regularities on the initial conditions.

In the rest of the section, we will obtain some a priori error estimates on the numerical errors ${\boldsymbol{v}}_n- {\boldsymbol{v}}_n^{hk}$, ${\boldsymbol{u}}_n- {\boldsymbol{u}}_n^{hk}$, ${\boldsymbol{\psi}}_n- {\boldsymbol{\psi}}_n^{hk}$, ${\boldsymbol{\phi}}_n- {\boldsymbol{\phi}}_n^{hk}$ and $T_n-T_n^{hk}$. We have the following.

Theorem 3. Let the assumptions of Theorem 1 still hold. If we denote by $({\boldsymbol{v}}, {\boldsymbol{\psi}}, T)$ the solution to problem $VP$ and by $({\boldsymbol{v}}^{hk}, {\boldsymbol{\psi}}^{hk}, T^{hk})$ the solution to problem $VP^{hk}$, then we have the following a priori error estimates, for all ${\boldsymbol{w}}^h = \{ {\boldsymbol{w}}_j^{h}\}_{j = 0}^N, \, {\boldsymbol{\xi}}^h = \{ {\boldsymbol{\xi}}_j^{h}\}_{j = 0}^N\subset V^h$ and $\zeta^h = \{\zeta_j^{h}\}_{j = 0}^N\subset E^h$,

where $M$ is again a positive constant which does not depend on parameters $h$ and $k$.

Proof. First, we obtain the error estimates on the velocity field. Thus, subtracting variational Eq (2.8) at time $t = t_n$ for a test function ${\boldsymbol{w}} = {\boldsymbol{w}}^h\in V^h\subset V$ and discrete variational Eq (3.3) we have, for all ${\boldsymbol{w}}^h\in V^h$,

and therefore, we have, for all ${\boldsymbol{w}}^h\in V^h$,

Using the estimates

it follows that, for all ${\boldsymbol{w}}^h\in V^h$,

where $\epsilon > 0$ is assumed small enough and we have used several times inequality Eq (3.7).

Proceeding in an analogous form, taking into account that

we obtain the following error estimates for the microrotation vector, for all ${\boldsymbol{\xi}}^h\in V^h$,

where $\epsilon > 0$ is assumed again small enough.

Finally, we derive the error estimates on the temperature field. Then, subtracting variational Eq (2.10) at time $t = t_n$ for a test function $\zeta = \zeta^h\in E^h\subset E$ and discrete variational Eq (3.5), we find that, for all $\zeta^h\in E^h$,

Therefore, it follows that, for all $\zeta^h\in E^h$,

Keeping in mind that

using inequality Eq (3.7) several times we find that, for all $\zeta^h\in E^h$,

where $\epsilon > 0$ is small enough.

Combining estimates (3.9)–(3.11) we have, for all ${\boldsymbol{w}}^h, \, {\boldsymbol{\xi}}^h\in V^h$ and $\zeta^h\in E^h$,

Multiplying the above estimates by $k$ and summing up to $n$ we have

Keeping in mind the estimates:

we use a discrete version of Gronwall's inequality ^[31] and we conclude a priori error estimates (3.8).

From estimates (3.8), if we assume suitable regularity conditions on the continuous solution, we can obtain the convergence order of the approximations. So, we have the following result.

Corollary 1. Under the assumptions of Theorem 3 and the additional regularity conditions:

we conclude the linear convergence of the approximations given by Problem VP$^{hk}$; that is, there exists a constant $M > 0$, independent of parameters $h$ and $k$, such that

4.   Numerical results

In this last section, we describe some two-dimensional numerical results obtained in the simulation of an academical example to show the accuracy of the finite element approximation, and the generation of the microrotations and temperature due to a surface mechanical force.

Since the problem is two-dimensional, we present below the corresponding version of Problem P. Defining the solutions

after some calculations we find the following problem:

where $\alpha = 1, 2$ and the two-dimensional Laplacian operator is defined as $\Delta^* = \frac{\partial}{\partial x_1} +\frac{\partial}{\partial x_2}.$

The numerical scheme provided by the fully discrete problem $VP^{hk}$ was implemented on a 3.2 Ghz PC using MATLAB, and a typical run (using parameters $h = \sqrt{2}/32$ and $k = 0.001$) took about $145$ seconds of CPU time.

4.1.   First example: numerical convergence and discrete energy decay

As an academical example, in order to show the numerical behavior of the fully discrete approximations we solve a two-dimensional problem with the following data:

the homogeneous Dirichlet boundary conditions, for all ${\boldsymbol{x}} = (x_1, x_2)\in \Gamma$,

the initial conditions, for all $(x_1, x_2)\in (0, 1)\times (0, 1)$ and $\alpha = 1, 2$,

and the supply terms, for all $(x_1, x_2, t)\in (0, 1)\times (0, 1)\times (0, 1)$,

In this case, the exact solution to the above two-dimensional problem can be easily calculated and it has the form, for $(x_1, x_2, t)\in [0, 1]\times [0, 1]\times [0, 1]$:

In order to analyze the numerical convergence, several uniform partitions of the domain have been performed dividing $\Omega = (0, 1)\times(0, 1)$ into $2(nd)^2$ triangles (that is, the spatial discretization parameter $h$ equals $\frac{\sqrt{2}}{nd}$). Thus, the approximation errors estimated by

are presented in Table 1 for several values of the discretization parameters $h$ and $k$. Moreover, the evolution of the error depen\-ding on the parameter $h+k$ is plotted in Figure 1. The convergence of the approximations is clearly found, and the linear convergence, stated in Corollary 1, seems to be achieved, providing a good accuracy.

If we assume now that there are not supply terms, and we use the final time $T_f = 5$, the following data:

and the initial conditions, for all $(x_1, x_2, t)\in (0, 1)\times (0, 1)\times (0, 1)$,

The discrete energy $E_n^{hk}$ is given by

where the discrete version of tensor ${\boldsymbol{e}}$ is defined as

Taking the time discretization parameter $k = 0.01$, the evolution in time of the discrete energy given above is plotted in Figure 2 (in both natural and semi-log scales) for different values of $\lambda$, $\mu$ and $\sigma$. As can be seen, it converges to zero and an exponential decay seems to be achieved in both cases.

4.2.   Second example: application of a surface mechanical force

In this second example, our aim is to show the coupling effect in a thermo-viscoelastic problem with a prescribed mechanical force. Therefore, we consider the same square domain clamped on its left vertical part $\{0\}\times [0, 1]$. We also assume that $\phi_3$ and $T$ vanish on the whole boundary, and we use the following expression for the surface mechanical force ${\boldsymbol{f}}$ applied on the horizontal upper boundary $[0, 1]\times\{1\}$:

The rest of the boundary is assumed traction-free. Moreover, we employ the following data:

and null initial conditions for all the variables. Taking the time discretization parameter $k = 0.01$, in Figure 3 we plot the norm of the displacement field over the deformed mesh at final time. As expected, the bending of the body has been produced due to the applied force and the clamping boundary condition.

Moreover, in Figure 4 the microrotations and the temperature are shown at final time. They have been generated by the deformation of the body, reaching the lowest (microrotations) and highest (temperature) values in its interior with a quadratic form. We also see the influence of the null boundary condition.

A.   Derivation of the relations for the viscous coefficients

In conditions (2.7) we have proposed some assumptions on the constitutive constants to guarantee that the internal energy and the dissipation are positive definite. The assumptions on the coefficients $\lambda$, $\mu$, $\sigma$, $\gamma$, $\alpha+\beta$, $\mu_v+\lambda_v$, $\mu_v+\sigma_v$ and $\sigma_v$ are standard and we can find them in the book of Eringen ^[10]; however, we do not know any reference where the assumptions on $\beta_v$, $\gamma_v$, $\alpha_v$, $\kappa$, $b^*$ and $\kappa^*$ are obtained.

To clarify this aspect we need to impose that the matrices

are positive definite.

The eigenvalues of the first matrix are $\beta_v+\gamma_v$ and $3\alpha_v+\beta_v+\gamma_v$. Therefore, they should be positive.

The eigenvalues of the second matrix are $\beta_v+\gamma_v$ and

Last two eigenvalues are positive if and only if $\gamma_v-\beta_v+\kappa > 0$ and $(\gamma_v-\beta_v+\kappa)^2 > 2(b^*+ \kappa^*)^2+(\beta_v+\kappa-\gamma_v)^2.$ This last condition is equivalent to

Thus, our conditions become $\kappa > 0$ and $2\kappa(\gamma_v-\beta_v) > (b^*+\kappa^*)^2.$

Acknowledgments

This paper is part of the projects PGC2018-096696-B-I00 and PID2019-105118GB-I00, funded by the Spanish Ministry of Science, Innovation and Universities and FEDER "A way to make Europe".

Conflict of interest

The authors declare there is no conflicts of interest.