/[escript]/trunk/doc/user/heatedblock.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3989 - (show annotations)
Tue Sep 25 02:21:54 2012 UTC (6 years, 10 months ago) by jfenwick
File MIME type: application/x-tex
File size: 7689 byte(s)
More copyright fixes.
pyvisi traces removed.
Some install doco
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2012 by University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Open Software License version 3.0
8 % http://www.opensource.org/licenses/osl-3.0.php
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[0])*[1.,0.,0.] \
150 +whereZero(x[1])*[0.,1.,0.] \
151 +whereZero(x[2])*[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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26