Given:
![{\displaystyle {\begin{matrix}{\text{(1)}}\qquad u_{,xx}&=&\alpha u_{,t}\\{\text{(2)}}\qquad u_{,xxxx}&=&\alpha u_{,tt}\end{matrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c04a01290741f741f4e4b65bd3e9cbfa689a58cb)
What physical processes do the two PDEs model?
Equation (1) models the one dimensional heat problem where
represents temperature at
at time
and
represents the thermal diffusivity.
Equation (2) models the transverse vibrations of a beam. The variable
represents the displacement of the point on the beam corresponding to position
at time
.
is
where
is Young's modulus,
is area moment of inertia,
is mass density, and
is area of cross section.
Write out the PDEs in elementary partial differential notation
![{\displaystyle {\begin{aligned}u_{,xx}&=\alpha u_{,t}\qquad \equiv \qquad {\frac {\partial ^{2}u}{\partial x^{2}}}=\alpha {\frac {\partial u}{\partial t}}\\u_{,xxxx}&=\alpha u_{,tt}\qquad \equiv \qquad {\frac {\partial ^{4}u}{\partial x^{4}}}=\alpha {\frac {\partial ^{2}u}{\partial t^{2}}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b4ae18e5f5d8e1a770f095cbc160d67136aaec25)
Determine whether the PDEs are elliptic, hyperbolic, or parabolic. Show how you arrived at your classification.
Equation (1) is parabolic for all values of
.
We can see why by writing it out in the form given in the notes on Partial differential equations.
![{\displaystyle (1){\frac {\partial ^{2}u}{\partial x^{2}}}+(0){\frac {\partial ^{2}u}{\partial x\partial t}}+(0){\frac {\partial ^{2}u}{\partial t^{2}}}+(0){\frac {\partial u}{\partial x}}+(-\alpha ){\frac {\partial u}{\partial t}}+(0)u=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f1b32a8b00bca1d0409763d21ea13db7d916204c)
Hence, we have
,
,
,
,
,
, and
. Therefore,
![{\displaystyle {b^{2}-4ac=0-(1)(0)=0\Rightarrow {\text{parabolic}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1f543025fffc43054cee1d5d7258d6f06ccae22c)
Equation (2) is hyperbolic for positive values of
, parabolic when
is zero, and elliptic when
is less than zero.
Given:
![{\displaystyle {\frac {d}{dx}}\left(a(x){\frac {du(x)}{dx}}\right)+c(x)u(x)=f(x)\quad {\mbox{for }}\quad x\in (O,L)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1fc3d78c74fdbaad37687e7c55d5cfdf573f4488)
- Heat transfer: :
![{\displaystyle kA{\frac {d^{2}T}{dx^{2}}}+ap\beta T=f}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4723b8bba45dbe52b77ac6547f83063957ae6e9b)
- Flow through porous medium: :
![{\displaystyle \mu {\frac {d^{2}\phi }{dx^{2}}}=f}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fccbb80c3aef63af2a1c253231662cd5cdd8b0ff)
- Flow through pipe: :
![{\displaystyle {\frac {1}{R}}{\frac {d^{2}P}{dx^{2}}}=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ccddef1bb65bf7916ead52903b4eb0f15ebf3a10)
- Flow of viscous fluids: :
![{\displaystyle \mu {\frac {d^{2}v_{z}}{dx^{2}}}=-{\frac {dP}{dx}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/44a5d69771959e6ac1330a656fd3843903a48f5f)
- Elastic cables: :
![{\displaystyle T{\frac {d^{2}u}{dx^{2}}}=f}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1104c8a0d77ea8b5408aa1ae2b9988e7fecd6fcd)
- Elastic bars: :
![{\displaystyle EA{\frac {d^{2}u}{dx^{2}}}=f}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a92b2bdeaad145db8b53cbf3291586db794c243b)
- Torsion of bars: :
![{\displaystyle GJ{\frac {d^{2}\theta }{dx^{2}}}=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/315750514c6df8b59f7d730e9ff3c21609df9aa9)
- Electrostatics: :
![{\displaystyle \varepsilon {\frac {d^{2}\phi }{dx^{2}}}=\rho }](https://wikimedia.org/api/rest_v1/media/math/render/svg/f4f77e022744351fee048c1b1fd018edccc619ac)
Given:
The residual over an element is
![{\displaystyle R^{e}(x,c_{1}^{e},c_{2}^{e},\dots ,c_{n}^{e}):={\cfrac {d}{dx}}\left(a~{\cfrac {du_{h}^{e}}{dx}}\right)+c~u_{h}^{e}-f(x)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9f8d41cfc74abcf6c6dcaad3640df15003316bc2)
where
![{\displaystyle u_{h}^{e}(x)=\sum _{j=1}^{n}c_{j}^{e}\varphi _{j}^{e}(x)~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fc4b8e7c04091ff8b531051549ed054b54bdefaa)
In the least squares approach, the residual is minimized in the
following sense
![{\displaystyle {\frac {\partial }{\partial c_{i}}}\int _{x_{a}}^{x_{b}}[R^{e}(x,c_{1}^{e},c_{2}^{e},\dots ,c_{n}^{e})]^{2}~dx=0,\qquad i=1,2,\dots ,n~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6150c420def861ff65898d5d81664147d5171800)
What is the weighting function used in the least squares approach?
We have,
![{\displaystyle {\frac {\partial }{\partial c_{i}}}\int _{x_{a}}^{x_{b}}[R^{e}]^{2}~dx=\int _{x_{a}}^{x_{b}}2R^{e}{\frac {\partial R^{e}}{\partial c_{i}}}~dx=0~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/47739e6a69829de3f627a720f93d267835345600)
Hence the weighting functions are of the form
![{\displaystyle {w_{i}(x)={\frac {\partial R^{e}}{\partial c_{i}}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2025fe683e76247c6486d3c55228d3e75cef2e12)
Develop a least-squares finite element model for the problem.
To make the typing of the following easier, I have gotten rid of the subscript
and the superscript
. Please note that these are implicit in what follows.
We start with the approximate solution
![{\displaystyle u(x)=\sum _{j=1}^{n}c_{j}\varphi _{j}(x)~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8b7a89c934ef603493986e3be7289cc28f673709)
The residual is
![{\displaystyle R={\cfrac {d}{dx}}\left[a~\left(\sum _{j}c_{j}{\cfrac {d\varphi _{j}}{dx}}\right)\right]+c~\left(\sum _{j}c_{j}\varphi _{j}\right)-f(x)~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/72416b5d53bd417de587652db451a2b57d90d87e)
The derivative of the residual is
![{\displaystyle {\frac {\partial R}{\partial c_{i}}}={\cfrac {d}{dx}}\left[a~\left(\sum _{j}{\frac {\partial c_{j}}{\partial c_{i}}}{\cfrac {d\varphi _{j}}{dx}}\right)\right]+c~\left(\sum _{j}{\frac {\partial c_{j}}{\partial c_{i}}}\varphi _{j}\right)-f(x)~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cade660aa945709af2fef8b2b9b1c1ca5bddff3d)
For simplicity, let us assume that
is constant, and
.
Then, we have
![{\displaystyle R=a~\left(\sum _{j}c_{j}{\cfrac {d^{2}\varphi _{j}}{dx^{2}}}\right)-f~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/26137aef6f965a13267dde6f0cd1529897bddfd3)
and
![{\displaystyle {\frac {\partial R}{\partial c_{i}}}=a~\left(\sum _{j}{\frac {\partial c_{j}}{\partial c_{i}}}{\cfrac {d^{2}\varphi _{j}}{dx^{2}}}\right)~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/62eef81825f49613c7f58efe83bcba11c598994b)
Now, the coefficients
are independent. Hence, their derivatives
are
![{\displaystyle {\frac {\partial c_{j}}{\partial c_{i}}}={\begin{cases}0&{\text{for}}~~i\neq j\\1&{\text{for}}~~i=j\end{cases}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3dc246a4f875de9900ae0cce02dcc7b5bacc7c4c)
Hence, the weighting functions are
![{\displaystyle {\frac {\partial R}{\partial c_{i}}}=a~{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/81b32b59d72756e674fd5ee5fae170759226be26)
The product of these two quantities is
![{\displaystyle R{\frac {\partial R}{\partial c_{i}}}=\left(a\sum _{j}c_{j}{\cfrac {d^{2}\varphi _{j}}{dx^{2}}}-f\right)\left(a{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}\right)=\sum _{j}\left[\left(a^{2}{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}{\cfrac {d^{2}\varphi _{j}}{dx^{2}}}\right)c_{j}\right]-af{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/23f955c9827d82da02741f7ee6315210a83ed3c2)
Therefore,
![{\displaystyle \int _{x_{a}}^{x_{b}}R{\frac {\partial R}{\partial c_{i}}}~dx=0\implies \sum _{j}\left[\left(\int _{x_{a}}^{x_{b}}a^{2}{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}{\cfrac {d^{2}\varphi _{j}}{dx^{2}}}~dx\right)c_{j}\right]-\int _{x_{a}}^{x_{b}}af{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}~dx=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/143b07ace71a4d1a4c83e03903389306f3cc0a4e)
These equations can be written as
![{\displaystyle K_{ij}c_{j}=f_{i}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/de12efe2c21c4dd80eb2eeb69d2e7739214fc308)
where
![{\displaystyle {K_{ij}=\int _{x_{a}}^{x_{b}}a^{2}{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}{\cfrac {d^{2}\varphi _{j}}{dx^{2}}}~dx}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/580319294d1612a728858cf2f9f0ab5b59082b68)
and
![{\displaystyle {f_{i}=\int _{x_{a}}^{x_{b}}af{\cfrac {d^{2}\varphi _{i}}{dx^{2}}}~dx~.}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3afb7545f81de27d86619bb3ceb277c388e62201)
- Discuss some functions
that can be used to approximate ![{\displaystyle u}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c3e6bb763d22c20916ed4f0bb6bd49d7470cffd8)
For the least squares method to be a true "finite element" type of method we have to use finite element shape functions. For instance, if each element has two nodes we could choose
![{\displaystyle \varphi _{1}={\cfrac {x_{b}-x}{h}}\qquad {\text{and}}\qquad \varphi _{2}={\cfrac {x-x_{a}}{h}}~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/adf14c651e1d91405859207c120b4887835edae7)
Another possiblity is set set of polynomials such as
![{\displaystyle \varphi =\{1,x,x^{2},x^{3}\}~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d858c815600ec6c2ba6492b6647a416f8a30e4a5)
Given:
![{\displaystyle -{\frac {d}{dx}}\left(k{\frac {dT}{dx}}\right)=q\quad {\mbox{ for }}\quad x\in (0,L)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1c7a814e0fdea52c859d9c8273e18b999eb6ddd4)
where
is the temperature,
is the thermal conductivity, and
is the heat generation. The Boundary conditions are
![{\displaystyle {\begin{matrix}T(0)&=&T_{0}\\k{\frac {dT}{dx}}+\beta (T-T_{\infty })+{\hat {q}}{\bigg |}_{x=L}&=&0\end{matrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/221bcb0209b2dac70b86f102b91beca23960b2f7)
Formulate the finite element model for this problem over one element
First step is to write the weighted-residual statement
![{\displaystyle 0=\int _{x_{a}}^{x_{b}}w_{i}^{e}\left[-{\frac {d}{dx}}\left(k{\frac {dT_{h}^{e}}{dx}}\right)-q\right]dx}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7f9d28aad695a58be69873277c9297c837673947)
Second step is to trade differentiation from
to
, using
integration by parts. We obtain
![{\displaystyle 0=\int _{x_{a}}^{x_{b}}\left(k{\frac {w_{i}^{e}}{dx}}{\frac {dT_{h}^{e}}{dx}}-w_{i}^{e}q\right)dx-\left[w_{i}^{e}k{\frac {dT_{h}^{e}}{dx}}\right]_{x_{a}}^{x_{b}}{\text{(3)}}\qquad }](https://wikimedia.org/api/rest_v1/media/math/render/svg/eaf53a02a1e1fe6dd0c49c65e55da713ea4a6a13)
Rewriting (3) in using the primary variable
and
secondary variable
as described in the
text book to obtain the weak form
![{\displaystyle 0=\int _{x_{a}}^{x_{b}}\left(k{\frac {w_{i}^{e}}{dx}}{\frac {dT_{h}^{e}}{dx}}-w_{i}^{e}q\right)dx-w_{i}^{e}(x_{a})Q_{a}^{e}-w_{i}^{e}(x_{b})Q_{b}^{e}{\text{(4)}}\qquad }](https://wikimedia.org/api/rest_v1/media/math/render/svg/7dfa56e29964cb1fcef58bed4be11d551e69b1b0)
From (4) we have the following
![{\displaystyle {\begin{matrix}K_{ij}^{e}&=&\int _{x_{a}}^{x_{b}}k{\frac {N_{i}^{e}}{dx}}{\frac {dN_{j}^{e}}{dx}}dx\\f_{i}^{e}&=&\int _{x_{a}}^{x_{b}}qN_{i}^{e}dx\\K_{ij}^{e}T_{j}^{e}&=&f_{i}^{e}+N_{i}^{e}(x_{a})Q_{a}^{e}+N_{i}^{e}(x_{b})Q_{b}^{e}\end{matrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/76873eaf54d086c1d560e654803e66d27c55d7d8)
Use ANSYS to solve the problem of heat conduction in an insulated rod. Compare the solution with the exact solution. Try at least three refinements of the mesh to confirm that your solution converges.
The following inputs are used to solve for the solution of this problem:
/prep7
pi = 4*atan(1) !compute value for pi
ET,1,LINK34!convection element
ET,2,LINK32!conduction element
R,1,1!AREA = 1
MP,KXX,1,0.01!conductivity
MP,HF,1,25 !convectivity
emax = 2
length = 0.1
*do,nnumber,1,emax+1 !begin loop to creat nodes
n,nnumber,length/emax*(nnumber-1),0
n,emax+2,length,0!extra node for the convetion element
*enddo
type,2 !switch to conduction element
*do,enumber,1,emax
e,enumber,enumber+1
*enddo
TYPE,1 !switch to convection element
e,emax+1,emax+2!creat zero length element
D,1,TEMP,50
D,emax+2,TEMP,5
FINISH
/solu
solve
fini
/post1
/output,p5,txt
prnsol,temp
prnld,heat
/output,,
fini
ANSYS gives the following results
Node
|
Temperature
|
1
|
50.000
|
2
|
27.590
|
3
|
5.1793
|
Node
|
Heat Flux
|
1
|
-4.4821
|
4
|
4.4821
|
Write you own code to solve the problem in part 2.
The following Maple code can be used to solve the problem.
> restart;
> with(linalg):
> L:=0.1:kxx:=0.01:beta:=25:T0:=50:Tinf:=5:
> # Number of elements
> element:=10:
> # Division length
> ndiv:=L/element:
> # Define node location
> for i from 1 to element+1 do
> x[i]:=ndiv*(i-1);
> od:
> # Define shape function
> for j from 1 to element do
> N[j][1]:=(x[j+1]-x)/ndiv;
> N[j][2]:=(x-x[j])/ndiv;
> od:
> # Create element stiffness matrix
> for k from 1 to element do
> for i from 1 to 2 do
> for j from 1 to 2 do
> Kde[k][i,j]:=int(kxx*diff(N[k][i],x)*diff(N[k][j],x),x=x[k]..x[k+1]);
> od;
> od;
> od;
> Kde[0][2,2]:=0:Kde[element+1][1,1]:=0:
> # Create global matrix for conduction Kdg and convection Khg
> Kdg:=Matrix(1..element+1,1..element+1,shape=symmetric):
> Khg:=Matrix(1..element+1,1..element+1,shape=symmetric):
> for k from 1 to element+1 do
> Kdg[k,k]:=Kde[k-1][2,2]+Kde[k][1,1];
> if (k <= element) then
> Kdg[k,k+1]:=Kde[k][1,2];
> end if;
> od:
> Khg[element+1,element+1]:=beta:
> # Create global matrix Kg = Kdg + Khg
> Kg:=matadd(Kdg,Khg):
> # Create T and F vectors
> Tvect:=vector(element+1):Fvect:=vector(element+1):
> for i from 1 to element do
> Tvect[i+1]:=T[i+1]:
> Fvect[i]:=f[i]:
> od:
> Tvect[1]:=T0:
> Fvect[element+1]:=Tinf*beta:
> # Solve the system Ku=f
> Ku:=multiply(Kg,Tvect):
> # Create new K matrix without first row and column from Kg
> newK:=Matrix(1..element,1..element):
> for i from 1 to element do
> for j from 1 to element do
> newK[i,j]:=coeff(Ku[i+1],Tvect[j+1]);
> od;
> od;
> # Create new f matrix without f1
> newf:=vector(element,0):
> newf[1]:=-Ku[2]+coeff(Ku[2],T[2])*T[2]+coeff(Ku[2],T[3])*T[3]:
> newf[element]:=Tinf*beta:
> # Solve the system newK*T=newf
> solution:=multiply(inverse(newK),newf);
> exact := 50-448.2071713*x;
> with(plots):
> # Warning, the name changecoords has been redefined
> data := [[ x[n+1], solution[n]] <math>n=1..element]:
> FE:=plot(data, x=0..L, style=point,symbol=circle,legend="FE"):
> ELAS:=plot(exact,x=0..L,legend="exact",color=black):
> display({FE,ELAS},title="T(x) vs x",labels=["x","T(x)"]);