Test Case of Thermal Dilatation in Axisymmetrical Coordinates Static in Vectorial Form

1. Introduction

This page presents the test case of Thermal Dilatation implement with CFPDEs with constant temperature without volumic force in static and axisymmetrical case, solve in vectorial form. This simulation have goal to verify the implementation of thermal dilatation.

2. Run the calculation

The command line to run this case is :

    mpirun -np 16 feelpp_toolbox_coefficientformpdes --config-file=thermal-dilatation.cfg --cfpdes.gmsh.hsize=5e-3

3. Data Files

The case data files are available in Github here :

CFG file - Edit the file
JSON file - Edit the file
GEO file - Edit the file

4. Equation

The domain of resolution is \(\Omega_c^{axis}\) with bounds \(\Gamma_c^{axis}\), \(\Gamma_{D \hspace{0.05cm} elas}^{axis}\) the bound of elastic Dirichlet conditions, \(\Gamma_{N \hspace{0.05cm} elas}^{axis}\) the bound of elastic Neumann conditions such that \(\Gamma_c^{axis} = \Gamma_{D \hspace{0.05cm} elas}^{axis} \cup \Gamma_{N \hspace{0.05cm} elas}^{axis}\).

Differential Formulation of Coupled Elasticity Equation and AV-Formulation in Stationary and Axisymmetrical case

\[\text{(Elasticity-AV Axis Static)} \\ \left\{ \begin{matrix} - \nabla \cdot \bar{\bar{\sigma}} = 0 & \text{on } \Omega_c^{axis} & \text{(Elas-1)} \\ \mathbf{u} = 0 & \text{on } \Gamma_{D \hspace{0.05cm} elas}^{axis} & \text{(D elas)} \\ \bar{\bar{\sigma}} \cdot \mathbf{n} = 0 & \text{on } \Gamma_{N \hspace{0.05cm} elas}^{axis} & \text{(N elas)} \end{matrix} \right.\]

With :

Lamé’s coefficients : \(\lambda = \frac{E \, v}{(1-2 v) \, (1+v)}\), \(\mu = \frac{E}{2 \, (1+v)}\) with
- Young modulus \(E\)
- Poisson’s coefficient \(v\)

Displacement \(\mathbf{u} = \begin{pmatrix} u_x \\ u_y \\ u_z \end{pmatrix}\)

Tensor of small deformation \(\bar{\bar{\epsilon}} = \frac{1}{2} \left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right)\)

Stress Tensor with Hooke Law : \(\bar{\bar{\sigma}} = \bar{\bar{\sigma}}^E = \lambda \, Tr(\bar{\bar{\epsilon}})+2\mu \bar{\bar{\epsilon}} - (3\lambda+2\mu) \alpha_T \left( T - T_0 \right) \) with coefficient of linear dilatation \(\alpha_T\), the temperature \(T\) and the rest temperature \(T_0\)

And notations :

Divergence of tensor : \(\nabla \cdot \bar{\bar{\sigma}} = \begin{pmatrix} \nabla \cdot \bar{\bar{\sigma}}_{r,:} + \frac{\bar{\bar{\sigma}}_{rr} - \bar{\bar{\sigma}}_{\theta \theta}}{r} \\ \nabla \cdot \bar{\bar{\sigma}}_{z,:} + \frac{\bar{\bar{\sigma}}_{rz}}{r} \end{pmatrix}\)
Scalar produce of tensor : \(\bar{\bar{\sigma}} \cdot \mathbf{n} = \begin{pmatrix} \bar{\bar{\sigma}}_{r,:} \cdot \mathbf{n}^{axis} \\ \bar{\bar{\sigma}}_{z,:} \cdot \mathbf{n}^{axis} \end{pmatrix} = \begin{pmatrix} \bar{\bar{\sigma}}_{rr} \, n_r + \bar{\bar{\sigma}}_{rz} \, n_z \\ \bar{\bar{\sigma}}_{zr} \, n_r + \bar{\bar{\sigma}}_{zz} \, n_z\end{pmatrix}\)

5. Geometry

The geometry is a torus of the conductor in cartesian coordinates \((x,y,z)\) or rectangle in axisymmetric coordinates \((r,z)\), surrounded by air.

Axisymmetrical cut of Geometry

Geometry in Three Dimensions

The geometrical domains are :

Conductor : the torus is composed by conductor material
- Interior : interior surface of conductor ring
- Exterior : exterior surface of conductor ring
- Upper : upper of conductor ring
- Bottom : bottom of conductor ring
Air : the air surround Conductor
- zAxis : a bound of Air, correspond to \(Oz\) axis (\(\{(z,r), \, z=0 \}\))
- infty : the rest of bound of Air

Symbol

Description

value

unit

\(r_{int}\)

interior radius of torus

\(75e-3\)

\(r_{ext}\)

exterior radius of torus

\(100.2e-3\)

\(z_1\)

half-height of torus

\(25e-3\)

6. Boundary Conditions

We impose the boundary conditions :

Strong Dirichlet : \(\mathbf{u} = 0\) on Upper and Bottom. The Dirichlet condition represents the embedding of mechanical part.
Neumann : \(\bar{\bar{\sigma}} \cdot \mathbf{n} = 0\) on Interior and Exterior. The Neumann condition represents the freedom of displacement.

On JSON file, the boundary conditions are writed :

Boundary conditions on JSON file

   "BoundaryConditions":
    {
        "elastic":
        {
            "Dirichlet":
            {
                "Bottom":
                {
                    "expr":"{0,0}"
                },
                "Upper":
                {
                    "expr":"{0,0}"
                }
            }
        }
    }

7. Weak Formulation

We obtain :

Elastic equations in Axisymmetric and Stationnary Case

\[\tiny{ \text{(Weak Static Elasticity Axis Vectoria)} \\ \left\{ \begin{align*} \int_{\Omega_c^{axis}} \left( r \left( \lambda \left( \frac{\partial u_r}{\partial r} + \frac{1}{r} u_r + \frac{\partial u_z}{\partial z} \right) - (3\lambda + 2 \mu) \, \alpha_T \, \left( T-T_0 \right) \right) \left( \frac{\partial \xi_r}{\partial r} + \frac{1}{r} \xi_r \right) + r \mu \left[ 2 \left( \frac{\partial u_r}{\partial r} \frac{\partial \xi_r}{\partial r} + \frac{1}{r^2} u_r \xi_r \right) + \left( \frac{\partial u_r}{\partial z} \frac{\partial \xi_r}{\partial z} + \frac{\partial u_z}{\partial r} \frac{\partial \xi_r}{\partial z} \right) \right] \\ + r \mu \left( \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r} \right) \frac{\partial \xi_z}{\partial r} + r \left[ \lambda \left( \frac{\partial u_r}{\partial r} + \frac{1}{r} u_r + \frac{\partial u_z}{\partial z} \right) + 2\mu \frac{\partial u_z}{\partial z} - (3\lambda + 2 \mu) \, \alpha_T \, \left( T-T_0 \right) \right] \frac{\partial \xi_z}{\partial z} \right) dr dz \\ = \int_{\Omega_c^{axis}}{r \, \mathbf{F} \cdot \mathbf{ξ} \ dr dz} + \int_{\Gamma_{N \hspace{0.05cm} elas}^{axis}}{ r \, \mathbf{g} \cdot \mathbf{ξ} \ dr dz } + \int_{\Gamma_{R \hspace{0.05cm} elas}^{axis}}{ r \left( \mathbf{g} - k \mathbf{u} \right) \cdot \mathbf{ξ} \ dr dz } & \hspace{0.5cm} \text{(Weak Static Elasticity Axis)} \\ \text{for} \ \mathbf{ξ} = \begin{pmatrix} \xi_r \\ \xi_z \end{pmatrix}_{axis} \in H^1_{u_0}(\Omega_c) \end{align*} \right. }\]

8. Parameters

The parameters of problem are :

Symbol

Description

Value

Unit

\(E\)

Young Modulus

\(2.1e6\)

\(Pa\)

\(v\)

Poisson’s coefficient

\(0.33\)

\(dimensionless\)

\(\lambda\)

Lame’s coefficient

\(\frac{E \, v}{(1-2v)(1+v)}\)

\(Pa\)

\(\mu\)

Lame’s coefficient

\(\frac{E}{2 (1+v)}\)

\(Pa\)

\(T\)

temperature

\(293\)

\(K\)

\(T_0\)

rest temperature

\(303\)

\(K\)

On JSON file, the parmeters are writed :

Parameters on JSON file

    "Parameters": {
        "E": 2.1e6,
        "nu": 0.33,
        "mu": "E/(2*(1+nu)):E:nu",
        "lambda":"E*nu/((1+nu)*(1-2*nu)):E:nu",
        "T0":293,
        "T":303,
        "alpha_T":"17e-6",
        "sigma_T":"-E/(1-2*nu)*alpha_T*(T-T0):E:nu:alpha_T:T:T0",
        "F":"{0,0}"
    }

9. Coefficient Form PDEs

We use the application Coefficient Form PDEs. The coefficient associate to Weak Formulation are :

Coefficient

Description

Expression

\(c\)

diffusion coefficient

\(r \mu\)

\(\gamma\)

conservative flux source term

\(\begin{pmatrix} - \lambda \left( r \frac{\partial u_r}{\partial r} + u_r + r \frac{\partial u_z}{\partial z} \right) - \mu \, r \, \frac{\partial u_r}{\partial r} - r \sigma_T & - \mu \, r \, \frac{\partial u_z}{\partial r} \\ - \mu \, r \, \frac{\partial u_r}{\partial z} & - \lambda \left( r \, \frac{\partial u_r}{\partial r} + u_r + r \, \frac{\partial u_z}{\partial z} \right) - \mu \, r \, \frac{\partial u_z}{\partial z} - r \sigma_T \end{pmatrix}\)

\(f\)

source term

\(\begin{pmatrix} r F_r - \lambda \frac{\partial u_r}{\partial r} - (\lambda+2\mu) \frac{u_r}{r} - \lambda \frac{\partial u_z}{\partial z} - \sigma_T \\ r F_z \end{pmatrix}\)

On JSON file, the coefficients are writed :

CFPDEs coefficients on JSON file

    "Materials":
    {
        "Conductor":
        {
            "elastic_c":"x*mu:x:mu",
            "elastic_gamma":"{-lambda*(x*elastic_grad_u_00+elastic_u_0+x*elastic_grad_u_11)-x*mu*elastic_grad_u_00-x*sigma_T, -x*mu*elastic_grad_u_10, -x*mu*elastic_grad_u_01, -lambda*(x*elastic_grad_u_00+elastic_u_0+x*elastic_grad_u_11)-x*mu*elastic_grad_u_11-x*sigma_T}:x:lambda:mu:elastic_u_0:elastic_grad_u_00:elastic_grad_u_01:elastic_grad_u_10:elastic_grad_u_11:sigma_T",
            "elastic_f":"{x*F_0-lambda*elastic_grad_u_00-(lambda+2*mu)*elastic_u_0/x-lambda*elastic_grad_u_11-sigma_T, x*F_1}:x:F_0:F_1:mu:lambda:elastic_u_0:elastic_grad_u_00:elastic_grad_u_01:elastic_grad_u_11:sigma_T",

            "sigma_rr":"(lambda+2*mu)*elastic_grad_u_00+lambda*elastic_u_0/x+lambda*elastic_grad_u_11+sigma_T:lambda:mu:x:elastic_u_0:elastic_grad_u_00:elastic_grad_u_11:elastic_grad_u_22:sigma_T",
            "sigma_thth":"lambda*elastic_grad_u_00+(lambda+2*mu)*lambda*elastic_u_0/x+lambda*elastic_grad_u_11+sigma_T:lambda:mu:x:elastic_u_0:elastic_grad_u_00:elastic_grad_u_11:elastic_grad_u_22:sigma_T",
            "sigma_zz":"lambda*elastic_grad_u_00+lambda*lambda*elastic_u_0/x+(lambda+2*mu)*elastic_grad_u_11+sigma_T:lambda:mu:x:elastic_u_0:elastic_grad_u_00:elastic_grad_u_11:elastic_grad_u_22:sigma_T",
            "sigma_rz":"mu*(elastic_grad_u_01+elastic_grad_u_10):mu:elastic_grad_u_01:elastic_grad_u_10",

            "sigma_1":"lambda/x*elastic_u_0+(lambda+mu)*(elastic_grad_u_00+elastic_grad_u_11)+sigma_T+mu*sqrt((elastic_grad_u_00-elastic_grad_u_11)*(elastic_grad_u_00-elastic_grad_u_11)+4*(elastic_grad_u_01+elastic_grad_u_10)*(elastic_grad_u_01+elastic_grad_u_10)):x:lambda:mu:elastic_u_0:elastic_grad_u_00:elastic_grad_u_11:elastic_grad_u_01:elastic_grad_u_10:sigma_T",
            "sigma_2":"lambda*elastic_grad_u_00+(lambda+2*mu)/x*elastic_u_0+lambda*elastic_grad_u_11+sigma_T:x:lambda:mu:elastic_u_0:elastic_grad_u_00:elastic_grad_u_11:sigma_T",
            "sigma_3":"lambda/x*elastic_u_0+(lambda+mu)*(elastic_grad_u_00+elastic_grad_u_11)+sigma_T-mu*sqrt((elastic_grad_u_00-elastic_grad_u_11)*(elastic_grad_u_00-elastic_grad_u_11)+4*(elastic_grad_u_01+elastic_grad_u_10)*(elastic_grad_u_01+elastic_grad_u_10)):x:lambda:mu:elastic_u_0:elastic_grad_u_00:elastic_grad_u_11:elastic_grad_u_01:elastic_grad_u_10:sigma_T"
        }
    }

10. Numeric Parameters

Mesh size \(1 \, mm\)

Mesh of Geometry

11. Run the calculation

11.1. cfpdes

11.2. CATIA

12. Comparison

I run with CATIA (software of conception) the same case. We can compare the result :

Deformation result of dilatation compute with cfpdes : the opaque part is the initial geometry and the transparent part is the geometry after dilatation (with scale=1000)

Deformation result of dilatation compute with catia

Deformation result of dilatation compute with cfpdes (with scale=1000) : the colours correspond to displacement

Deformation result of dilatation compute with catia : the colours correspond to displacement

The comparison isn’t reliable (visual comparison), the result of CFPDEs and CATIA have the form and the same order (near of center-exterior of torus, the displacement is around \(7.2 \, mm\) for the two results).

With comparison with CATIA, we have some confidence in this method.