doc/user/heatedblock.tex

%
% $Id$
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%           Copyright 2003-2007 by ACceSS MNRF
%       Copyright 2007 by University of Queensland
%
%                http://esscc.uq.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 examine the deformation of a linear elastic body caused by expansion through a heat distribution. We want 
a 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 trough 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 a 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 trying 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 we may use $Y\hackscore{i}=\lambda+\frac{2}{3} \mu) \; \alpha \; T\hackscore{i}$. 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}
This Lame equation is in fact symmetric, given the difference in $D$ and $d$ as compared to the scalar case.
The \LinearPDE class is notified of 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 damage. 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	ksteube	1316	%
2	jgs	110	% $Id$
3	gross	625	%
4	ksteube	1316	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5	gross	625	%
6	ksteube	1316	% Copyright 2003-2007 by ACceSS MNRF
7			% Copyright 2007 by University of Queensland
8			%
9			% http://esscc.uq.edu.au
10			% Primary Business: Queensland, Australia
11			% Licensed under the Open Software License version 3.0
12			% http://www.opensource.org/licenses/osl-3.0.php
13			%
14			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15			%
16
17	gross	578	\section{Elastic Deformation}
18			\label{ELASTIC CHAP}
19	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
20			a displacement field $u\hackscore{i}$ which solves the momentum equation
21	gross	578	\index{momentum equation}:
22	gross	579	\begin{eqnarray}\label{HEATEDBLOCK general problem}
23	gross	578	- \sigma\hackscore{ij,j}=0
24			\end{eqnarray}
25			where the stress $\sigma$ is given by
26	gross	579	\begin{eqnarray}\label{HEATEDBLOCK linear elastic}
27	gross	578	\sigma\hackscore{ij}= \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
28			- (\lambda+\frac{2}{3} \mu) \; \alpha \; (T-T\hackscore{ref})\delta\hackscore{ij} \;.
29			\end{eqnarray}
30			In this formula $\lambda$ and $\mu$ are the Lame coefficients, $\alpha$ is the
31			temperature expansion coefficient, $T$ is the temperature distribution and $T_{ref}$ a reference temperature. Note that
32	gross	579	\eqn{HEATEDBLOCK general problem} is similar to eqn{WAVE general problem} introduced in section~\Sec{WAVE CHAP} but the
33	gross	578	inertia term $\rho u\hackscore{i,tt}$ has been dropped as we assume a static scenario here. Moreover, in
34			comparison to the \eqn{WAVE stress}
35	gross	579	definition of stress $\sigma$ in \eqn{HEATEDBLOCK linear elastic} an extra term is introduced
36	woo409	757	to bring in stress due to volume changes trough temperature dependent expansion.
37	gross	568
38	gross	578	Our domain is the unit cube
39	gross	593	\begin{eqnarray} \label{HEATEDBLOCK natural location}
40	gross	578	\Omega=\{(x\hackscore{i} \| 0 \le x\hackscore{i} \le 1 \}
41			\end{eqnarray}
42			On the boundary the normal stress component is set to zero
43	gross	579	\begin{eqnarray} \label{HEATEDBLOCK natural}
44	gross	578	\sigma\hackscore{ij}n\hackscore{j}=0
45			\end{eqnarray}
46			and on the face with $x\hackscore{i}=0$ we set the $i$-th component of the displacement to $0$
47	gross	579	\begin{eqnarray} \label{HEATEDBLOCK constraint}
48	gross	578	u\hackscore{i}(x)=0 & \mbox{ where } & x\hackscore{i}=0 \;
49			\end{eqnarray}
50			For the temperature distribution we use
51	gross	579	\begin{eqnarray} \label{HEATEDBLOCK temperature}
52			T(x)= T\hackscore{0} e^{-\beta \\|x-x^{c}\\|};
53	gross	578	\end{eqnarray}
54	gross	579	with a given positive constant $\beta$ and location $x^{c}$ in the domain. Later in \Sec{MODELFRAME} we will use
55	ksteube	1316	$T$ from a time-dependent temperature diffusion problem as discussed in \Sec{DIFFUSION CHAP}.
56
57	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
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			-A\hackscore{ijkl} u\hackscore{k,l}){,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			+\delta\hackscore{ik} \delta\hackscore{jl}
79			\delta\hackscore{il} \delta\hackscore{jk}) \\
80			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	ksteube	1316	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			\frac{1}{6} ((\sigma\hackscore{00}-\sigma\hackscore{11})^2+
108			(\sigma\hackscore{11}-\sigma\hackscore{22})^2
109			(\sigma\hackscore{22}-\sigma\hackscore{00})^2)
110			+ \sigma\hackscore{01}^2+\sigma\hackscore{12}^2+\sigma\hackscore{20}^2}
111			\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_0exp(-betalength(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))+lamtrace(g)*kronecker(mydomain)-sigma0
154			sigma_mises=sqrt(((sigma[0,0]-sigma[1,1])2+(sigma[1,1]-sigma[2,2])2+ \
155			(sigma[2,2]-sigma[0,0])**2)/6. \
156			+sigma[0,1]2 + sigma[1,2]2 + sigma[2,0]**2)
157			#... output ...
158			saveVTK("deform.xml",disp=u,stress=sigma_mises)
159			\end{python}
160
161			\begin{figure}
162	gross	599	\centerline{\includegraphics[width=\figwidth]{figures/HeatedBlock.eps}}
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
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision