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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26