doc/user/heatedblock.tex

% $Id$
%
%           Copyright © 2006 by ACcESS MNRF
%               \url{http://www.access.edu.au
%         Primary Business: Queensland, Australia.
%   Licensed under the Open Software License version 3.0
%      http://www.opensource.org/licenses/osl-3.0.php
%
\section{Elastic Deformation}
\label{ELASTIC CHAP}
In this section we want to discuss the deformation of a linear elastic body caused by expansion through a heat distribution. We want 
to displacement field $u\hackscore{i}$ which solves the momentum equation
\index{momentum equation}:
\begin{eqnarray}\label{HEATEDBLOCK general problem}
 - \sigma\hackscore{ij,j}=0
\end{eqnarray}
where the stress $\sigma$ is given by
\begin{eqnarray}\label{HEATEDBLOCK linear elastic}
 \sigma\hackscore{ij}= \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
 - (\lambda+\frac{2}{3} \mu)  \; \alpha  \;  (T-T\hackscore{ref})\delta\hackscore{ij} \;.
\end{eqnarray}
In this formula $\lambda$ and $\mu$ are the Lame coefficients, $\alpha$ is the 
temperature expansion coefficient, $T$ is the temperature distribution and $T_{ref}$ a reference temperature. Note that 
\eqn{HEATEDBLOCK general problem} is similar to eqn{WAVE general problem} introduced in section~\Sec{WAVE CHAP} but the
inertia term $\rho u\hackscore{i,tt}$ has been dropped as we assume a static scenario here. Moreover, in 
comparison to the \eqn{WAVE stress}
definition of stress $\sigma$ in \eqn{HEATEDBLOCK linear elastic} an extra term is introduced 
to bring in stress due to volume changes through temperature dependent expansion.   

Our domain is the unit cube 
\begin{eqnarray} \label{HEATEDBLOCK natural location}
\Omega=\{(x\hackscore{i} | 0 \le x\hackscore{i} \le 1 \}
\end{eqnarray}
On the boundary the normal stress component is set to zero
\begin{eqnarray} \label{HEATEDBLOCK natural}
\sigma\hackscore{ij}n\hackscore{j}=0
\end{eqnarray}
and on the face with $x\hackscore{i}=0$ we set the $i$-th component of the displacement to $0$
\begin{eqnarray} \label{HEATEDBLOCK constraint}
u\hackscore{i}(x)=0 & \mbox{ where } & x\hackscore{i}=0 \; 
\end{eqnarray}
For the temperature distribution we use 
\begin{eqnarray} \label{HEATEDBLOCK temperature}
T(x)= T\hackscore{0} e^{-\beta \|x-x^{c}\|}; 
\end{eqnarray}
with a given positive constant $\beta$ and location $x^{c}$ in the domain. Later in \Sec{MODELFRAME} we will use
$T$ from an time-dependent temperature diffusion problem as discussed in \Sec{DIFFUSION CHAP}.
   
When we insert~\eqn{HEATEDBLOCK linear elastic} we get a second oder system of linear PDEs for the displacements $u$ which is called
the Lame equation\index{Lame equation}. We want to solve
this using the \LinearPDE class to this. For a system of PDEs and a solution with several components the \LinearPDE class 
takes PDEs of the form
\begin{equation}\label{LINEARPDE.SYSTEM.1 TUTORIAL}
-A\hackscore{ijkl} u\hackscore{k,l}){,j}=-X\hackscore{ij,j} \; .
\end{equation}
$A$ is of \RankFour and $X$ is of \RankTwo. We show here the coefficients relevant 
for the we traying to solve. The full form is given in~\eqn{LINEARPDE.SYSTEM.1}. 
The natural boundary conditions \index{boundary condition!natural} take the form:
\begin{equation}\label{LINEARPDE.SYSTEM.2 TUTORIAL}
n\hackscore{j} A\hackscore{ijkl} u\hackscore{k,l}=n\hackscore{j}X\hackscore{ij} \;.
\end{equation}
Constraints \index{constraint} take the form
\begin{equation}\label{LINEARPDE.SYSTEM.3 TUTORIAL}
u\hackscore{i}=r\hackscore{i} \mbox{ where } q\hackscore{i}>0
\end{equation}
$r$ and $q$ are each \RankOne. 
We can easily identify the coefficients in~\eqn{LINEARPDE.SYSTEM.1 TUTORIAL}:
\begin{eqnarray}\label{LINEARPDE ELASTIC COEFFICIENTS}
A\hackscore{ijkl}=\lambda \delta\hackscore{ij} \delta\hackscore{kl} + \mu ( 
+\delta\hackscore{ik} \delta\hackscore{jl}
\delta\hackscore{il} \delta\hackscore{jk}) \\
X\hackscore{ij}=(\lambda+\frac{2}{3} \mu) \;  \alpha \; (T-T\hackscore{ref})\delta\hackscore{ij} \\
\end{eqnarray}
The characteristic function $q$ defining the locations and components where constraints are set is given by:
\begin{equation}\label{HEATEDBLOCK MASK}
q\hackscore{i}(x)=\left\{ 
\begin{array}{cl}
1  & x\hackscore{i}=0  \\ 
0  & \mbox{otherwise}   \\
\end{array}
\right. 
\end{equation}
Under the assumption that $\lambda$, $\mu$, $\beta$ and $T\hackscore{ref}$ 
are constant setting $Y\hackscore{i}=\lambda+\frac{2}{3} \mu) \; \alpha \; T\hackscore{i}$ seems to be also possible. However,
this choice would lead to a different natural boundary condition which does not set the normal stress component as defined
in~\eqn{HEATEDBLOCK linear elastic} to zero.

Analogously to concept of symmetry for a single PDE, we call the PDE defined by~\eqn{LINEARPDE.SYSTEM.1 TUTORIAL} symmetric if
\index{symmetric PDE}
\begin{eqnarray}\label{LINEARPDE.SYSTEM.SYMMETRY TUTORIAL}
A\hackscore{ijkl} =A\hackscore{klij} \\
\end{eqnarray}
Note that different from the scalar case now the coefficients $D$ and $d$ have to be inspected. It is easy to see that 
the Lame equation in fact is symmetric. The \LinearPDE class is notified by this fact by calling its \method{setSymmetryOn} method.

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
\begin{equation}
\sigma\hackscore{mises} = \sqrt{
\frac{1}{6} ((\sigma\hackscore{00}-\sigma\hackscore{11})^2+
               (\sigma\hackscore{11}-\sigma\hackscore{22})^2
               (\sigma\hackscore{22}-\sigma\hackscore{00})^2)
+  \sigma\hackscore{01}^2+\sigma\hackscore{12}^2+\sigma\hackscore{20}^2}
\end{equation}
is used to detect material damages. Here we want to calculate the von--Mises and write the stress to a file for visualization.

\index{scripts!\file{diffusion.py}}:
The following script which is available in \file{heatedbox.py} in the \ExampleDirectory solves the Lame equation
and writes the displacements and the von--Mises stress\index{von--Mises stress} into a file \file{deform.xml} in the \VTK file format: 
\begin{python}
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Brick
#... set some parameters ...
lam=1.
mu=0.1
alpha=1.e-6
xc=[0.3,0.3,1.]
beta=8.
T_ref=0.
T_0=1.
#... generate domain ...
mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
x=mydomain.getX()
#... set temperature ...
T=T_0*exp(-beta*length(x-xc))
#... open symmetric PDE ...
mypde=LinearPDE(mydomain)
mypde.setSymmetryOn()
#... set coefficients ...
C=Tensor4(0.,Function(mydomain))
for i in range(mydomain.getDim()):
  for j in range(mydomain.getDim()):
     C[i,i,j,j]+=lam
     C[j,i,j,i]+=mu
     C[j,i,i,j]+=mu
msk=whereZero(x[0])*[1.,0.,0.] \
   +whereZero(x[1])*[0.,1.,0.] \
   +whereZero(x[2])*[0.,0.,1.]
sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)
mypde.setValue(A=C,X=sigma0,q=msk)
#... solve pde ...
u=mypde.getSolution()
#... calculate von-Misses stress
g=grad(u)
sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0
sigma_mises=sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \
                  (sigma[2,2]-sigma[0,0])**2)/6. \
                 +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)
#... output ...
saveVTK("deform.xml",disp=u,stress=sigma_mises)
\end{python}

\begin{figure}
\centerline{\includegraphics[width=\figwidth]{figures/HeatedBlock.eps}}
\caption{von--Mises Stress and Displacement Vectors.}
\label{HEATEDBLOCK FIG 2}
\end{figure}

Finally the the results can be visualize by calling
\begin{python}
mayavi -d deform.xml -f CellToPointData -m VelocityVector -m SurfaceMap &
\end{python}
Note that the filter \text{CellToPointData} is applied to create smooth representation of the 
von--Mises stress. \fig{HEATEDBLOCK FIG 2} shows the results where the vertical planes showing the 
von--Mises stress and the horizontal plane shows the vector representing displacements.

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