 Contents of /trunk/doc/user/heatedblock.tex

Revision 3989 - (show annotations)
Tue Sep 25 02:21:54 2012 UTC (7 years ago) by jfenwick
File MIME type: application/x-tex
File size: 7689 byte(s)
 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2012 by University of Queensland 4 5 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Open Software License version 3.0 8 9 % 10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 % Development since 2012 by School of Earth Sciences 12 % 13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 14 15 \section{Elastic Deformation} 16 \label{ELASTIC CHAP} 17 In this section we want to examine the deformation of a linear elastic body caused by expansion through a heat distribution. 18 We want a displacement field $u_{i}$ which solves the momentum 19 equation\index{momentum equation}: 20 \begin{eqnarray}\label{HEATEDBLOCK general problem} 21 - \sigma_{ij,j}=0 22 \end{eqnarray} 23 where the stress $\sigma$ is given by 24 \begin{eqnarray}\label{HEATEDBLOCK linear elastic} 25 \sigma_{ij}= \lambda u_{k,k} \delta_{ij} + \mu ( u_{i,j} + u_{j,i}) 26 - (\lambda+\frac{2}{3} \mu) \; \alpha \; (T-T_{ref})\delta_{ij} \;. 27 \end{eqnarray} 28 In this formula $\lambda$ and $\mu$ are the Lame coefficients, $\alpha$ is the 29 temperature expansion coefficient, $T$ is the temperature distribution and $T_{ref}$ a reference temperature. 30 Note that \eqn{HEATEDBLOCK general problem} is similar to \eqn{WAVE general problem} 31 introduced in \Sec{WAVE CHAP} but the inertia term $\rho u_{i,tt}$ 32 has been dropped as we assume a static scenario here. 33 Moreover, in comparison to the \eqn{WAVE stress} definition of stress $\sigma$ 34 in \eqn{HEATEDBLOCK linear elastic} an extra term is introduced to bring in 35 stress due to volume changes through temperature dependent expansion. 36 37 Our domain is the unit cube 38 \begin{eqnarray} \label{HEATEDBLOCK natural location} 39 \Omega=\{(x_{i}) | 0 \le x_{i} \le 1 \} 40 \end{eqnarray} 41 On the boundary the normal stress component is set to zero 42 \begin{eqnarray} \label{HEATEDBLOCK natural} 43 \sigma_{ij}n_{j}=0 44 \end{eqnarray} 45 and on the face with $x_{i}=0$ we set the $i$-th component of the displacement to $0$: 46 \begin{eqnarray} \label{HEATEDBLOCK constraint} 47 u_{i}(x)=0 & \mbox{ where } & x_{i}=0 \; 48 \end{eqnarray} 49 For the temperature distribution we use 50 \begin{eqnarray} \label{HEATEDBLOCK temperature} 51 T(x)= T_{0} e^{-\beta \|x-x^{c}\|} 52 \end{eqnarray} 53 with a given positive constant $\beta$ and location $x^{c}$ in the domain. 54 55 %Later in \Sec{MODELFRAME} we will use 56 % $T$ from a time-dependent temperature diffusion problem as discussed in \Sec{DIFFUSION CHAP}. 57 When we insert \eqn{HEATEDBLOCK linear elastic} we get a second order system 58 of linear PDEs for the displacements $u$ which is called the Lame equation\index{Lame equation}. 59 We want to solve this using the \LinearPDE class. 60 For a system of PDEs and a solution with several components the \LinearPDE class takes PDEs of the form 61 \begin{equation}\label{LINEARPDE.SYSTEM.1 TUTORIAL} 62 -(A_{ijkl} u_{k,l})_{,j}=-X_{ij,j} \; . 63 \end{equation} 64 $A$ is a \RankFour and $X$ is a \RankTwo. 65 We show here the coefficients relevant for the problem we are trying to solve. 66 The full form is given in \eqn{LINEARPDE.SYSTEM.1}. 67 The natural boundary conditions\index{boundary condition!natural} take the form 68 \begin{equation}\label{LINEARPDE.SYSTEM.2 TUTORIAL} 69 n_{j} A_{ijkl} u_{k,l}=n_{j}X_{ij} 70 \end{equation} 71 while constraints\index{constraint} take the form 72 \begin{equation}\label{LINEARPDE.SYSTEM.3 TUTORIAL} 73 u_{i}=r_{i} \mbox{ where } q_{i}>0 74 \end{equation} 75 $r$ and $q$ are each a \RankOne. 76 We can easily identify the coefficients in \eqn{LINEARPDE.SYSTEM.1 TUTORIAL}: 77 \begin{eqnarray}\label{LINEARPDE ELASTIC COEFFICIENTS} 78 A_{ijkl}=\lambda \delta_{ij} \delta_{kl} + \mu ( 79 \delta_{ik} \delta_{jl} 80 + \delta_{il} \delta_{jk}) \\ 81 X_{ij}=(\lambda+\frac{2}{3} \mu) \; \alpha \; (T-T_{ref})\delta_{ij} \\ 82 \end{eqnarray} 83 The characteristic function $q$ defining the locations and components where constraints are set is given by: 84 \begin{equation}\label{HEATEDBLOCK MASK} 85 q_{i}(x)=\left\{ 86 \begin{array}{cl} 87 1 & x_{i}=0\\ 88 0 & \mbox{otherwise}\\ 89 \end{array} 90 \right. 91 \end{equation} 92 Under the assumption that $\lambda$, $\mu$, $\beta$ and $T_{ref}$ 93 are constant we may use $Y_{i}=(\lambda+\frac{2}{3} \mu) \; \alpha \; T_{i}$. 94 However, this choice would lead to a different natural boundary condition 95 which does not set the normal stress component as defined in \eqn{HEATEDBLOCK linear elastic} to zero. 96 97 Analogously to the concept of symmetry for a single PDE, we call the PDE 98 defined by \eqn{LINEARPDE.SYSTEM.1 TUTORIAL} symmetric\index{symmetric PDE} if 99 \begin{eqnarray}\label{LINEARPDE.SYSTEM.SYMMETRY TUTORIAL} 100 A_{ijkl} =A_{klij} \\ 101 \end{eqnarray} 102 This Lame equation is in fact symmetric, given the difference in $D$ and $d$ as compared to the scalar case. 103 The \LinearPDE class is notified of this fact by calling its \method{setSymmetryOn} method. 104 105 After we have solved the Lame equation we want to analyse the actual stress distribution. 106 Typically the von-Mises stress\index{von-Mises stress} defined by 107 \begin{equation} 108 \sigma_{mises} = \sqrt{ 109 \frac{1}{2} ((\sigma_{00}-\sigma_{11})^2 110 + (\sigma_{11}-\sigma_{22})^2 111 + (\sigma_{22}-\sigma_{00})^2) 112 + 3( \sigma_{01}^2+\sigma_{12}^2+\sigma_{20}^2) } 113 \end{equation} 114 is used to detect material damage. 115 Here we want to calculate the von-Mises stress and write it to a file for visualization. 116 117 The following script, which is available in \file{heatedblock.py} in the 118 \ExampleDirectory, solves the Lame equation and writes the displacements and 119 the von-Mises stress\index{von-Mises stress} into a file \file{deform.vtu} in 120 the \VTK file format\index{scripts!\file{diffusion.py}}: 121 \begin{python} 122 from esys.escript import * 123 from esys.escript.linearPDEs import LinearPDE 124 from esys.finley import Brick 125 from esys.weipa import saveVTK 126 #... set some parameters ... 127 lam=1. 128 mu=0.1 129 alpha=1.e-6 130 xc=[0.3, 0.3, 1.] 131 beta=8. 132 T_ref=0. 133 T_0=1. 134 #... generate domain ... 135 mydomain = Brick(l0=1., l1=1., l2=1., n0=10, n1=10, n2=10) 136 x=mydomain.getX() 137 #... set temperature ... 138 T=T_0*exp(-beta*length(x-xc)) 139 #... open symmetric PDE ... 140 mypde=LinearPDE(mydomain) 141 mypde.setSymmetryOn() 142 #... set coefficients ... 143 C=Tensor4(0., Function(mydomain)) 144 for i in range(mydomain.getDim()): 145 for j in range(mydomain.getDim()): 146 C[i,i,j,j]+=lam 147 C[i,j,i,j]+=mu 148 C[i,j,j,i]+=mu 149 msk=whereZero(x)*[1.,0.,0.] \ 150 +whereZero(x)*[0.,1.,0.] \ 151 +whereZero(x)*[0.,0.,1.] 152 sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain) 153 mypde.setValue(A=C, X=sigma0, q=msk) 154 #... solve pde ... 155 u=mypde.getSolution() 156 #... calculate von-Mises stress 157 g=grad(u) 158 sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0 159 sigma_mises=sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \ 160 (sigma[2,2]-sigma[0,0])**2)/2. \ 161 +3*(sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)) 162 #... output ... 163 saveVTK("deform.vtu", disp=u, stress=sigma_mises) 164 \end{python} 165 166 \begin{figure} 167 \centerline{\includegraphics[width=\figwidth]{HeatedBlock}} 168 \caption{von-Mises Stress and Displacement Vectors} 169 \label{HEATEDBLOCK FIG 2} 170 \end{figure} 171 172 \noindent Finally, the results can be visualized by calling 173 \begin{verbatim} 174 mayavi2 -d deform.vtu -f CellToPointData -m VelocityVector -m SurfaceMap 175 \end{verbatim} 176 Note that the filter \text{CellToPointData} is applied to create a smoother 177 representation of the von-Mises stress. 178 \fig{HEATEDBLOCK FIG 2} shows the results where the colour of the vertical 179 planes represent the von-Mises stress and a horizontal plane of arrows shows 180 the displacements vectors. 181