/[escript]/branches/stage3.0/doc/user/heatedblock.tex
ViewVC logotype

Annotation of /branches/stage3.0/doc/user/heatedblock.tex

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26