User:Eml5526.s11.team2/hwk4

=Problem 4.1: Solve the PDE using a half Fourier Basis with the Weighted Residual Method =

Given: PDE, Boundary Conditions, and Basis Function
PDE:
 * {| style="width:100%" border="0"

$$  \displaystyle 2\frac{d^2u}{dx^2}+3=0 $$     (4.1.1)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Boundary Conditions:
 * {| style="width:100%" border="0"

$$  \displaystyle u(x=1)=0 $$     (4.1.2)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle -\frac{du}{dx}(x=0)=4 $$     (4.1.3)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Half Fourier Basis:
 * {| style="width:100%" border="0"

$$  \displaystyle b_j(x)={1, sin(x), sin(2x),...,sin(jx)} $$     (4.1.4)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Find: Unknown Coefficients in $$ u^h(x) $$ expansion to approximate $$ u(x) $$ of the PDE.
We want to approximate the solution of Eq.4.1.1 by a series expansion.
 * {| style="width:100%" border="0"

$$  \displaystyle u^h(x)=\sum_{i=0}^n{d_ib_i(x)} $$     (4.1.5) To solve for the unknown coefficients $$ d_i, $$ we must develop the same number of linearly independent equations as we have unknowns. The number of unknowns we generate with the expansion of $$ u^h(x) $$ is entirely up to us, i.e. is it dependent on the number of degrees of freedom we decide to use for the problem formulation. The following solutions are examples of using 3, 5 and 7 degrees of freedom.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Develop 2 equations which satisfy the given Boundary Conditions.
First let's expand Eq.4.1.1 out to n=2 terms.
 * {| style="width:100%" border="0"

$$  \displaystyle u^h(x)=d_0b_0+d_1b_1+d_2b_2 $$     (4.1.6) Now consider Eq.4.1.4, and substitute the appropriate terms into Eq.4.1.6.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle u^h(x)=d_0(1)+d_1\Big(sin(x)\Big)+d_2\Big(sin(2x)\Big) $$     (4.1.7) Let's rewrite this equation in terms of a product between 2 vectors.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle u^h(x)=\left[1\quad sin(x) \quad sin(2x)\right ]^T\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix} $$     (4.1.8)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

In order for us to solve for the unknown $$d_i's$$ in Eq.4.1.8, we must have 3 linearly independent equations describing linear combinations of the $$d_i's.$$Two of the equations are easily obtained directly from the given boundary conditions in Eq.4.1.3 and Eq.4.1.4. We want our approximation to closely resemble our exact solution, therefore we must require identical behavior between Eq.4.1.6 and Eq.4.1.1 at the boundaries.
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{align} u^h(x=1)=d_0(1)+d_1sin(1)+d_2sin(2)&=0\\ \end{align} $$     (4.1.9)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \frac{du}{dx}(x=0)=d_1cos(0)+2d_2cos(0)=-4 $$     (4.1.11)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Use the Weighted Residual Method to develop a $$ 3^{rd} $$ linearly independent equation.
Now we just need one more linearly independent equation and we can then readily solve for the $$d_i's.$$ To develop another set of equations which are linearly independent of Eq.4.1.10 and Eq.4.1.12, we take the inner product between, Eq.4.1.15 and our basis set of functions Eq.4.1.4. Before getting to the inner product first, let's define our linear differential operator as the following
 * {| style="width:100%" border="0"

$$  \displaystyle P(\cdot)=\frac{d^2}{dx^2}+\frac{3}{2}=0 $$     (4.1.13) Now, using our operator in Eq.4.1.13, we can operate on our approximation in Eq.4.1.8.
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{align} P(u^h)=\frac{d^2}{dx^2}\left[1\quad sin(x)\quad sin(2x)\right ]^T\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}+\frac{3}{2}&=0\\ \left[0\quad -sin(x)\quad -4sin(2x)\right ]^T\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}+\frac{3}{2}&=0 \end{align} $$     (4.1.14) Which we can then rewrite Eq.4.1.14 as
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \left[0\quad -sin(x)\quad -4sin(2x)\right ]^T\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}=-\frac{3}{2} $$     (4.1.15) Now, we are in a position to take the inner product of Eq.4.1.15(our approximation) with the basis set in Eq.4.1.4 with $$ n=2.$$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \left \langle b_j, P(u^h) \right \rangle =\int_0^1b_jP(u^h)=0 $$     (4.1.16)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \int_0^1\begin{bmatrix} 1\\sin(x)\\sin(2x)\end{bmatrix} \left[0\quad -sin(x)\quad -4sin(2x)\right ]^T\,dx\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}=-\frac{3}{2}\int_0^1\begin{bmatrix} 1\\sin(x)\\sin(2x)\end{bmatrix}\,dx $$     (4.1.17)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \int_0^1\begin{bmatrix} 0&-sin(x)&-4sin(2x)\\0&-sin^2(x)&-4sin(x)sin(2x)\\0&-sin(x)sin(2x)&-4sin^2(2x)\end{bmatrix}\,dx\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}=-\frac{3}{2}\int_0^1\begin{bmatrix} 1\\sin(x)\\sin(2x)\end{bmatrix}\,dx $$     (4.1.18)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{bmatrix} 0&-1&2cos(2)-2\\0&sin(2)/4-1/2&2sin(3)/3-2sin(1)\\0&sin(3)/6-sin(1)/2&sin(4)/2-2\end{bmatrix}\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}=\begin{bmatrix} -3/2\\-3/2\\-3sin^2(1)/2\end{bmatrix} $$     (4.1.19) We can use any one of the above equations as our $$ 3^{rd} $$ equation. For simplicity, let's choose the $$ 1^{st} $$ equation in the above matrix expression. This would be
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0" cellpadding="5" cellspacing="0"

(4.1.20)
 * style="width:10%" |
 * style="width:80%" |
 * style="width:10%" | <p style="text-align:right">
 * style="width:10%" | <p style="text-align:right">
 * }

Now we are set to expression our formulated problem in Matrix Form.

Display in Matrix Form

 * {| style="width:100%" border="0"

$$  \displaystyle \begin{bmatrix} 1&sin(1)&sin(2)\\0&1&2\\0&-1&2\Big(cos(2)-1 \Big)\end{bmatrix}\begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}=\begin{bmatrix} 0\\-4\\-3/2\end{bmatrix} $$     (4.1.21)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Solve for unknown d.
Our stiffness matrix is
 * {| style="width:100%" border="0"

$$  \displaystyle \underline{\mathbf{K}}=\begin{bmatrix} 1&sin(1)&sin(2)\\0&1&2\\0&-1&2\Big(cos(2)-1 \Big)\end{bmatrix} $$     (4.1.22) Using MATLAB to aid in our solution we find the inverse of our stiffness matrix to be
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \underline{\mathbf{K}}^{-1}=\begin{bmatrix} 1&0&0\\1327/1577&1&-1\\401/441&2&-2644/1443\end{bmatrix} $$     (4.1.23) Therefore, we find $$ {\mathbf{d}} $$ to be
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle {\mathbf{d}} = {\mathbf{K}}^{-1}{\mathbf{f}} $$     (4.1.24)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{bmatrix} d_0\\d_1\\d_2\end{bmatrix}=\begin{bmatrix} 1&0&0\\1327/1577&1&-1\\401/441&2&-2644/1443\end{bmatrix} \begin{bmatrix}0\\-4\\-3/2\end{bmatrix} $$     (4.1.25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Construct,Plot and find the error between the approximate and exact solutions for u(x).
A plot of both the approximate and exact solutions can be found below.



We can see that our approximate solution is though very close to exact, does not converge as quickly as we saw the polynomial basis function converge. This to be due to the slower converging half Fourier Basis when approximating a polynomial, versus using a polynomial basis. When using a polynomial basis of degree n, it will provide us with the exact solution of any polynomial of degree n or less. This can be easily proven through the use of the Integral Mean Value Theorem. Let's move on to using 5 degrees of freedom (n=4) and see our or approximate solution can improve.
 * {| style="width:100%" border="0"

$$  \displaystyle \Delta \mathbf{E_{abs}}=u^h(x=0.5)-u(x=0.5) $$     (4.1.27)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Or in other words approximately 3%. Not to bad but lets see if we can get better.

Construct,Plot and find the error between the approximate and exact solutions for u(x)..
Being that we have extensively explained the mathematics for case 2 above, we will just show the MATLAB code that I have written to solve this problem while focus more on the results of increasing the degrees of freedom, i.e. the number of coefficients that we expand to in our approximate solution. This not only makes it easier on the writer, but on the reader as well. The following is the Plot for the case n=4.



We see that we have a little further away with the 2 extra degrees of freedom. Let's calculate our error below.
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{align} \Delta \mathbf{E_{abs}}=u^h(x=0.5)-u(x=0.5) \end{align} $$     (4.1.29) Or in other words we have increased our error with an increase in expansion terms of our approximate solution.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Computational Difficulties
This case proved to be too computational demanding for my computer. I let the code, which worked for both cases 2 and 4 above run for 45 min and still MATLAB said it was busy. Therefore, I force quit my program and can conclude our best result is obtained with the n=2 case.

=Problem 4.2: Heat Conduction Fish & Belytschko Page 73 Problem 3.5 =

Given: The Following Boundary Conditions
Where (4.2.1) is the essential boundary condition defined at $$ \displaystyle x = 0 $$ and (4.2.2) is the natural boundary condition defined at $$ \displaystyle x = 10 $$ on the domain $$ \Omega = ]0,10[ \quad $$

Find:The Weak Form for 1-D Heat Conduction From the Strong Form
Obtain the weak form for the equation of the heat conduction with the boundary conditions stated above. The condition on the natural boundary condition is a convection condition.

Solution: Obtaining the Weak Form
The strong form has the following equation for a 1-D heat conduction body.

Assuming a steady-state from for the heat conduction and a constant area ( 4.2.3) reduces to

Multiplying by an arbitrary weighting function $$ \displaystyle w(x) $$ and integrating with respect to x across the domain leads to the following form

Making use of Integration by Parts on the first term of (4.2.5) leads to the following equation.

Using the requirement of the essential boundary condition on the weighting function satisfying the essential boundary location implies that $$ \displaystyle w(0) = 0 $$ and also using the natural boundary condition (4.2.2) $$ k \frac{\partial T}{\partial x}\bigg|_{x=0} = hT $$ leads to the following

Therefore the weak form can be represented finding a $$ \displaystyle T(x) $$ such that for all smooth $$ \displaystyle T(x) $$ with $$ \displaystyle T(0) = 100 $$ (4.2.8) is satisfied.

= Problem 4.3: Construct weak form of heat problem and find solution =

Given: Strong Form and Boundary Conditions
Statement:- Given the strong form for the heat conduction problem in a circular plate:
 * {| style="width:100%" border="0"

$$  \displaystyle k\frac{d}\left( {r\frac} \right) + rs = 0,0 < r \leqslant R $$ (4.3.1) natural boundary condition :$$ \frac\left( {r = 0} \right) = 0$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

essential boundary condition :$$T\left( {r = R} \right) = 0$$

where R is the total radius of the plate, s is the heat source per unit length along the plate radius, T is the temperature and k is the conductivity. Assume that k, s and R are given:

a. Weak Form
Construct the weak form for the above strong form.

b. Solution Using Quadratic Polynomial Trial Function
Use quadratic trial (candidate) solutions of the form $$ \displaystyle T = {\alpha _0} + {\alpha _1}r + {\alpha _2}{r^2}$$ and weight functions of the same form to obtain a solution of the weak form.

c. Solve the Differential Equation
Solve the differential equation with the boundary conditions and showthat the temperature distribution along the radius is given by:


 * {| style="width:100%" border="0"

$$  \displaystyle T = \frac{s}\left( {{R^2} - {r^2}} \right) $$     (4.3.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

a. Weak Form Derivation
Strong Form Can be written as:


 * {| style="width:100%" border="0"

$$  \displaystyle p\left( T \right) = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle p\left( T \right) = k\frac{d}\left( {r\frac} \right) + rs $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }

$$w\left( r \right)$$ is weight function Weighted Residual Form can be written as


 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^R {w\left( r \right)p\left( T \right)dr = 0} $$     (4.3.3)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^R {w\left( r \right)\left[ {k\frac{d}\left( {r\frac} \right) + rs} \right]dr} = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^R {w\left( r \right)k\frac{d}\left( {r\frac} \right)dr + } \int\limits_0^R {w\left( r \right)rsdr} = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }

Write $$w\left( r \right) = w$$ Integrating the first term by parts.


 * {| style="width:100%" border="0"

$$  \displaystyle \left[ {wkr\frac} \right]_0^R - \int\limits_0^R {\left( {\frackr\frac} \right)} dr + \int\limits_0^R {wrsdr} = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }

Assuming $$w\left( R \right) = 0$$ and using natural boundary condition $$\frac\left( {r = 0} \right) = 0$$


 * {| style="width:100%" border="0"

$$  \displaystyle 0 - 0 - \int\limits_0^R {\left( {\frackr\frac} \right)} dr + \int\limits_0^R {wrsdr} = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^R {\left( {\frackr\frac} \right)} dr = \int\limits_0^R {wrsdr} $$     (4.3.4)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Writing in standard form:


 * }

b. Solution using quadratic polynomial trial function
let $$T$$ and $$w$$ are of the following
 * {| style="width:100%" border="0"

$$  \displaystyle T = {\alpha _0} + {\alpha _1}r + {\alpha _2}{r^2} $$ and $$ \displaystyle w = {\beta _0} + {\beta _1}r + {\beta _2}{r^2} $$     (4.3.6)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Now using essential Boundary Condition i.e $$T\left( R \right) = 0$$


 * {| style="width:100%" border="0"

$$  \displaystyle {\alpha _0} + {\alpha _1}R + {\alpha _2}{R^2} = 0 $$     (4.3.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle {\alpha _0} = - \left( {{\alpha _1}R + {\alpha _2}{R^2}} \right) $$     (4.3.8)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Now using Natural boundary condition i.e. $$\frac\left( {r = 0} \right) = 0$$


 * {| style="width:100%" border="0"

$$  \displaystyle {\alpha _1} + 2{\alpha _2}\left( 0 \right) = 0 $$     (4.3.9)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle {\alpha _1} = 0 $$     (4.3.10)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle {\alpha _0} = - {\alpha _2}{R^2} $$     (4.3.11)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Now using homogeneous boundary condition $$w\left( R \right) = 0$$
 * {| style="width:100%" border="0"

$$  \displaystyle {\beta _0} + {\beta _1}R + {\beta _2}{R^2} = 0 $$     (4.3.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle {\beta _0} = - \left( {{\beta _1}R + {\beta _2}{R^2}} \right) $$     (4.3.13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Also


 * {| style="width:100%" border="0"

$$  \displaystyle \frac = {\alpha _1} + 2{\alpha _2}r $$     (4.3.14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

As $${\alpha _1} = 0$$
 * {| style="width:100%" border="0"

$$  \displaystyle \frac = 2{\alpha _2}r $$     (4.3.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Now differentiating $$w$$ with respect to $$r$$
 * {| style="width:100%" border="0"

$$  \displaystyle \frac = {\beta _1} + 2{\beta _2}r $$     (4.3.16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Putting values of $$w,\frac,\frac$$ in the standard weak form.


 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^R {kr\left( {{\beta _1} + 2{\beta _2}r} \right)} \left( {2{\alpha _2}r} \right)dr = \int\limits_0^R {sr\left( {{\beta _0} + {\beta _1}r + {\beta _2}{r^2}} \right)} dr $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle k\int\limits_0^R {\left( {2{\beta _1}{\alpha _2}{r^2} + 4{\beta _2}{\alpha _2}{r^3}} \right)dr} = s\int\limits_0^R {\left( {{\beta _0}r + {\beta _1}{r^2} + {\beta _2}{r^3}} \right)dr} $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \left[ {\frac{3} + k{\beta _2}{\alpha _2}{r^4}} \right]_0^R = \left[ {\frac{2} + \frac{3} + \frac{4}} \right]_0^R $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }

Putting value of $$$$ from equation


 * {| style="width:100%" border="0"

$$  \displaystyle \frac{3} + k{\beta _2}{\alpha _2}{R^4} = \frac{2} + \frac{3} + \frac{4} $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }

Rearranging terms


 * {| style="width:100%" border="0"

$$  \displaystyle {\beta _1}\left[ {\frac{3} + \frac{2} - \frac{3}} \right] + {\beta _2}\left[ {k{\alpha _2}{R^4} + \frac{2} - \frac{4}} \right] = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }

As $${\beta _1},{\beta _2}$$ are arbitrary they can't be zero, putting bracketed terms equal to zero both equation give same value of $${\alpha _2}$$.


 * {| style="width:100%" border="0"

$$  \displaystyle {\alpha _2} = \frac $$     (4.3.17)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using equation (4.3.10) $${\alpha _0}$$
 * {| style="width:100%" border="0"

$$  \displaystyle {\alpha _0} = - {\alpha _2}{R^2} = \frac $$     (4.3.18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Putting these values in $$T = {\alpha _0} + {\alpha _1}r + {\alpha _2}{r^2}$$


 * {| style="width:100%" border="0"

$$  \displaystyle T = \frac - \frac $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }

Now checking this solution by putting in given strong form given in equation.....


 * {| style="width:100%" border="0"

$$  \displaystyle k\frac{d}\left( {r\frac{d}\left( {\frac{s}\left( {{R^2} - {r^2}} \right)} \right)} \right) + rs = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle k\frac{d}\left( { - \frac} \right) + rs = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle k\left( { - \frac{k}} \right) + rs = 0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle - \cancel{sr} + \cancel{sr} = 0 $$     (4.3.20)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Comments
As we are not getting any particular interval of 'r' by putting the weak form solution in the given differential equation, so this solution is valid on given domain of 'r'.

c. Exact solution for given differential equation
Starting from 4.3.1 and rearranging we get


 * {| style="width:100%" border="0"

$$  \displaystyle d\left( {r\frac} \right) = - \frac{s}{k}rdr $$     (4.3.21)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Integrating this we get
 * {| style="width:100%" border="0"

$$  \displaystyle r\frac = - \frac{s}{r^2} + {c_1} $$     (4.3.22)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using the boundary condition and solving
 * {| style="width:100%" border="0"

$$  \displaystyle \frac(r = 0) = 0 $$     (4.3.23)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \Rightarrow {c_1} = 0 $$     (4.3.24)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Next we rearrange again and integrate
 * {| style="width:100%" border="0"

$$  \displaystyle dT = - \frac{s}rdr $$     (4.3.25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle T = - \frac{s}{r^2} + {c_2} $$     (4.3.26)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Again using the boundary condition to solve for the constant of integration
 * {| style="width:100%" border="0"

$$  \displaystyle T(R) = 0 $$     (4.3.27)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \Rightarrow {c_2} = \frac{s}{R^2} $$     (4.3.28)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

We arrive at the exact solution to the strong form which was found from the weak form
 * {| style="width:100%" border="0"

$$  \displaystyle T(r) = \frac{s}\left( {{R^2} - {r^2}} \right) $$     (4.3.29)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Solution from weak form and exact solution from the differential equation are same for given heat conduction problem.

=Problem 4.4: Deriving the Weak Form from the Strong Form for the Torsion of a Bar Fish & Belytschko Page 73 Problem 3.7 =

Given: The Strong Form
Given the strong form for the circular bar in torsion (Fig. 1):

Where the natural boundary condition is

Where the essential boundary condition is



Where the essential boundary condition $$ \displaystyle m(x) $$ is a distributed moment per unit length, M is the torsion moment, $$ \displaystyle \phi $$ is the angle of rotation, $$ \displaystyle G $$ is the shear modulus and $$ \displaystyle J $$ is the polar moment of intertia given by $$ \displaystyle J = \frac{\pi C^4}{2} $$, where $$ \displaystyle C $$ is the radius of the circular shaft.

b. Determine Constants of Integration
Assume that $$ \displaystyle m(x) = 0 $$ and integrate the differential equation given above. Find the integration constants using the boundary conditions.

(a) Obtaining the Weak Form
Multiplying (4.4.1) by an arbitrary weighting function $$ \displaystyle w(x) $$ and integrating with respect to x across the domain leads to the following form

Making use of Integration by Parts on the first term of (4.4.4) leads to the following equation.

Using the requirement of the essential boundary condition on the weighting function satisfying the essential boundary location implies that $$ \displaystyle w(0) = 0 $$ and also using the natural boundary condition (4.4.2) $$ JG \frac{\partial \phi}{\partial x}_{x=l} = M $$ leads to the following

Therefore the weak form can be represented finding a $$ \displaystyle \phi(x) $$ such that for all smooth $$ \displaystyle \phi(x) $$ with $$ \displaystyle \phi(0) = 0 $$ (4.4.7) is satisfied.

(b)Obtaining the Constants of Integration
Applying equation (4.4.1) and with $$ \displaystyle m(x) = 0 $$ leads to the following form.

Integrating (4.4.8) with respect to x generates a constant of integration and provides the following form.

Applying the natural boundary condition from (4.4.2) leads to the following

Solving (4.4.9) for $$ \displaystyle \frac{\partial \phi}{\partial x} $$ and substituting in $$ \displaystyle C_1 = -M $$ leads to the following

Integrating (4.4.11) with respect to x generates a constant of integration and provides the following form, which is the well known equation for angle of twist of a bar.

Applying the essential boundary condition from (4.4.3) to (4.4.12) leads to the following

Therefore the constants and the angle of twist are as follows.

=Problem 4.5: Constraint Breaking Solution (CBS) for Multiple Families of Basis Functions=

Given: Multiple Families of Basis Functions on a Domain
For

The possible basis function families:

e) Exponential Family:
Where the above families satisfy the the CBS with a domain of

1) Show the Following and Plot $$ \displaystyle b_j $$ and $$ \displaystyle \bar{b_j} $$ for $$ \displaystyle j=1,2,3$$
i) $$ \displaystyle \overline  (x) = const $$ ii) And plot $$\displaystyle {b_j} $$ and $$\displaystyle \overline   $$ for $$\displaystyle j=1,2,3 $$

a) Polynomial Family
Satisfying the CBS for the zero term:

Satisfying the CBS for the remaining terms:

The CBS family solution:

Plot of basis functions and CBS basis functions vs. the domain of -2 to 4:



b) Cosine Family
Satisfying the CBS for the zero term:

Satisfying the CBS for the remaining terms:

The CBS family solution:

Plot of basis functions and CBS basis functions vs. the domain of -2 to 4:



c) Sine Family
Satisfying the CBS for the zero term:

Satisfying the CBS for the remaining terms:

The CBS family solution:

Plot of basis functions and CBS basis functions vs. the domain of -2 to 4:



d) Fourier Family
Satisfying the CBS for the zero term:

Satisfying the CBS for the remaining terms:

The CBS family solution:

Plot of basis functions and CBS basis functions vs. the domain of -2 to 4:



e) Exponential Family
Satisfying the CBS for the zero term:

Satisfying the CBS for the remaining terms:

The CBS family solution:

Plot of basis functions and CBS basis functions vs. the domain of -2 to 4:



Conclusion
It seems that the strength or effective representation of the basis family would be expected to follow the trend: (exponential family) << (polynomial family)<< (sine family) ~ (cosine family) < (fourier family). Although this trend may not follow for every case(small j's and linear or simple curvature), it demonstrates that as j becomes significantly large, the fourier family of basis functions tends to become more accurate.

Part 2: Show that the exponential family is linearly independent on the domain
To be linearly independent, the gram matrix of the basis functions must not equal zero:

The gram matrix consists of the dot products of basis functions as follows:

However, with a given domain, the Gram Matrix can also be written as:

where,

A few examples:

etc...

The gram matrix of the exponential family becomes:

The determinate of the exponential gram matrix does not zero:

Therefore, the exponential family of basis functions is linearly independent on the domain

=Problem 4.6: Elastic bar with a variable distributed spring force p(x)=

Given: bar of length l, cross-sectional area A(x), Young’s modulus E(x) with body force b(x)
Consider an elastic bar with a variable distributed spring p(x) along its length as shown in Figure 1. The distributed spring imposes an axial force on the bar in proportion to the displacement. Consider a bar of length l, cross-sectional area A(x), Young’s modulus E(x) with body force b(x) and boundary conditions as shown in Figure 1

a. Derivation of Strong Form
Let $$\displaystyle \vec{n}(x) $$ be defined as the unit vector having positive orientation in the direction of $$ \sigma $$ on the element we are considering. This means that the sign convention of our unit vectors will be

Which corresponds to the convention depicted in the free-body diagram shown above. Let us also define $$ \displaystyle \sigma(x) $$ as the force for a given area as,

Upon the rearrangement of terms we arrive at

Classical mechanics allows us to write the sum of the forces in static system is zero

Using the free-body diagram, let us now consider all of the forces acting on the elastic bar $$\displaystyle b(x),p(x), N(x), $$ and $$\displaystyle  N(x+dx) $$. Balancing the forces acting on our finite element results in the following relation,

Next we will substitute (1.5) and (1.6) into (1.10) to get

Recall that a Taylor's Series expansion (See:Taylor Series) about $$ a $$ near $$ x $$ has the following form:

Next we will need to determine a relationship for the displaced area and stress, namely $$\displaystyle \sigma(x+dx)$$ and $$ \displaystyle   A(x+dx)$$. To accomplish this we will performing a Taylor Series expansion choosing $$\displaystyle a=x+dx $$ and neglect all higher order terms.

Taking the product of (1.13) and (1.14) gives us

Substitute (1.15) into (1.11) results in

This can be simplified by canceling the like terms and dividing through by $$ \displaystyle dx $$ as shown below,

to get the following simplified relation,

Next we must apply the product rule in the reverse directions. Meaning we will consider (1.19) the expansion obtained by applying the product rule and recombine terms to get the unexpanded form shown here,

tensile stress according to Hooke's Law

Substituting (1.22) into our simplified and unexpanded form of our force balance shown in (1.20) we arrive at the solution,
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{d}[E(x)\cdot A(x)\frac] + b(x)-p(x) = 0 $$     (4.6.19) Standard form of the strong form can be written as
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

b. Derivation of weak form
Now strong form can be written as
 * {| style="width:100%" border="0"

$$  \displaystyle P\left( u \right) = 0 $$     (4.6.21)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Where
 * {| style="width:100%" border="0"

$$  \displaystyle p\left( u \right) = \frac{d}\left[ {E\left( x \right)A\left( x \right)\frac} \right] + \left( {b\left( x \right) - p\left( x \right)} \right) $$     (4.6.22)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Let $${w\left( x \right)}$$ is a weight function, then writing Weighted Residual Form
 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^l {w\left( x \right)\left[ {P\left( u \right)} \right]dx = 0} $$     (4.6.23)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^l {w\left( x \right)\left[ {\frac{d}\left[ {E\left( x \right)A\left( x \right)\frac} \right] + \left( {b\left( x \right) - p\left( x \right)} \right)} \right]dx = 0} $$     (4.6.24)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Writing $$w\left( x \right) = w,E\left( x \right) = E,A\left( x \right) = A,b\left( x \right) = b,p\left( x \right) = p$$
 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^l {w\left[ {\frac{d}\left[ {EA\frac} \right] + \left( {b - p} \right)} \right]dx = 0} $$     (4.6.25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Expanding integral
 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^l {w\frac{d}\left[ {EA\frac} \right]dx} + \int\limits_0^l {w\left( {b - p} \right)dx}  = 0 $$     (4.6.26)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$  \displaystyle \left[ {wEA\frac} \right]_0^l - \int\limits_0^l {EA\frac\fracdx} + \int\limits_0^l {w\left( {b - p} \right)dx}  = 0 $$     (4.6.27)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Now assuming $$w\left( 0 \right) = 0$$ and also from natural boundary condition $$E(l)\frac = \bar t$$ above integral will take form as
 * {| style="width:100%" border="0"

$$  \displaystyle w\left( l \right)A\bar t - 0 - \int\limits_0^l {EA\frac\fracdx} + \int\limits_0^l {w\left( {b - p} \right)dx}  = 0 $$     (4.6.28)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Rearranging terms


 * {| style="width:100%" border="0"

$$  \displaystyle \int\limits_0^l {EA\frac\fracdx} = \int\limits_0^l {w\left( {b - p} \right)dx}  + w\left( l \right)A\bar t $$ (4.6.29)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Now writing in standard format

Reference
=Problem 4.7: Calculix User's Guide Containing Illustrative Examples For Building FEA FBDs=

First Step: Download
Calculix is a open source finite element suite. This software contains its own pre- and post processing algorithms which enables the user to prepare and solve a model after defining the appropriate boundary conditions. The solver, CCX, was written and is maintained by Guido Dhondt. The pre- and post processor, CGX, was written and is maintained by Klaus Wittig. Because CalculiX is open source users can download the latest version from Convergent Mechanical Solution's website. Convergent has done an excellent job integrating an installation wizard with the distribution which makes the installation process a breeze for experts and novice alike.


 * Users intending on downloading CalculiX with a Windows operating system should visit: http://www.bconverged.com/products.php
 * The built-in GUI makes it simple to follow the on screen prompts for a quick pain-free installation.
 * Users will be prompted to select a working directory which will be used as the starting file path whenever opening CalculiX editor.
 * The source code, Linux builds, documentation, tests, and other resources for Calculix : http://www.calculix.de/

Second Step: Test
Before starting Calculix we need to accomplish one more step.


 * {| style="width:100%" border="0"

$$  \displaystyle [Start]--[All Programs]--[bConverged]--[Calculix]--[Test Calculix *.*] $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Using Calculix

 * To open Calculix in the user friendly windows version go to


 * {| style="width:100%" border="0"

$$  \displaystyle [Start]--[All Programs]--[bConverged]--[Calculix]--[Calculix Command] $$
 * style="width:95%" |
 * style="width:95%" |
 * }


 * Doing this will open the Calculix command window.
 * The command window will automatically open in the working direction that we have defined during the installation process.
 * Next we introduce the program command line parameters which serve as options when creating or edit a file.


 * {|style="width:100%" border="0"

$$  \displaystyle Command: cgx [-a|-b|-c|-duns2d|-duns3d|-isaac2d|-isaac3d|-foam|-step|-stl]  filename $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Now we proceed to show an example of a pre-built FEA problem. Examples can be download from Calculix website or you can access them once bConverged has been installed and the pre-build samples have been tested. Because of the amount of sample files CalculiX must test, something initiated by the user after installation, the time will very for the completion of the initial tests. To find the examples simply go to the program directory containing bConverged. The location will very slightly depending on if you are running a 32 or 64-bit version of windows.

Since we are in the working directory when opening CalculiX we will need to copy the appropriate examples to our working directory to open them. Or you can simply change your current directory to a directory of your choosing based on where your want the examples to be placed. In the command line type ' cgx -b sphere.fbd '. Important !!!: Always make sure that the mouse is positioned on the graph window before writing commands. This point is extremely important because failing to do so will result in commands not being entered in the command window. Initially this is somewhat cumbersome because we have trained ourselves to point the mouse to the location where we will be inserting text. After a few minutes you will get the feel for doing this. Position your command window such that you can simultaneously view your entered commands and geometric figure.

Command line $$ \displaystyle \Rightarrow $$ C:\..........\examples\basic>cgx -b sphere.fbd



By using mouse you can rotate, zoom in out and out, or reposition the figure. Rotating the object is controlled by the left mouse button, zooming in and out by the middle mouse button and translating the object is controlled by the right mouse button. Before moving on try these operations with the mouse in order to become well acquainted with how these movements orient the figure.

Learning By Example
Lets learn more by doing an example on Calculix. We can broaden our knowledge on commands by doing examples. We can take modelling disc as an example.

Building Your First Model

 * Run Calculix Command
 * Type cgx -b discmodel.fbd (You probably get an error message. Because there is no exists file name: discmodel)
 * To start with modelling we have to remember the construction procedure.
 * First we will construct nodes, lines, surfaces, and the bodies.
 * Lets look at point command.


 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} Command:\quad 'pnt' <name(char<9)>& |'!' [<x> <y> <z>]| \\ & [   ]| \\                                  &  [<P1> <P2>  ]| \\ & [<setname(containing nodes)>]
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$
 * }

Creating Points

 * In order to type commands in to the command window please remember to press and release to the graph window with left button mouse button.
 * Type pnt y  -0.00000        1.00000        0.00000  (recognizing that py is a dummy with three coordinates on xyz plane)
 * Type pnt P0  -0.00000       -0.00000        0.00000
 * Type pnt P001 0.70711       -0.00000       -0.70711
 * Type pnt P003 -0.00000      -0.00000       -1.00000
 * Type pnt P005 -0.70711      -0.00000       -0.70711
 * Type pnt P006 -1.00000      -0.00000        0.00000
 * Type pnt P009 -0.70711       -0.00000       0.70711
 * Type pnt P00A 0.00000       -0.00000        1.00000
 * Type pnt P00G 0.70711       -0.00000        0.70711
 * Type pnt P00I 1.00000       -0.00000       -0.00000




 * Even though you have entered the points into the console you will not be able to see the point attributes until you enter a command to display them.
 * To see the points with attributes(labels) enter the following command.
 * Type plot pa all (get the same view with graph)
 * The general command used here is plot and has all of the available options shown here.
 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} Command: 'plot' & ['n'|'e'|'f'|'p'|'l'|'s'|'b'|'S'|'L']&['a'|'d'|'p'|'q'|'v'] -> \\ & 'w'|'k'|'r'|'g'|'b'|'y'|'m'
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$
 * }

Creating Lines

 * Now we will look at how to generate lines.
 * Type line L001 P00I P001 p0 4
 * This will create the line connecting p001 to p00I.
 * To show this line
 * Type plus la all
 * The plus command is used when one wants to retain the previous plot and add addition
 * {| style="width:100%" border="0"

$$  \displaystyle Command: 'line' <name(char<9)>|'!' <p1> <p2> <cp|seq>
 * style="width:95%" |
 * style="width:95%" |

$$
 * }




 * As you see above we just created arc for a portion of disc.
 * Even if the command's name is line we also create arcs with the same command.
 * L001 is a dummy name.
 * P00I and P001 are limits.
 * p0 is the center.
 * And 4 is the number of divisions.
 * To see the line attributes use the following:


 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} Command: 'plus'  & ['n'|'e'|'f'|'p'|'l'|'s'|'b'|'S'|'L']&['a'|'d'|'p'|'q'|'v'] -> \\ & 'w'|'k'|'r'|'g'|'b'|'y'|'m'|'i'
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$
 * }


 * Plus performs the same job as plot.
 * Plus retains all previous visuals along with the initial plot filters.
 * Type line L002 P001 P003 p0 4
 * Type line L003 P003 p0 8
 * Type line L004 p0 P00I 8
 * Type line L005 P003 P005 p0 4
 * Type line L006 P005 P006 p0 4
 * Type line L007 P006 p0 8
 * Type line L009 P006 P009 p0 4
 * Type line L00A P009 P00A p0 4
 * Type line L00C P00A p0 8
 * Type line L00G P00A P00G p0 4
 * Type line L00I P00G P00I p0 4

Creating Surfaces
To finish the model we need to add surfaces. To do so we will use the following generalized syntax.


 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} Command: 'gsur' &<name(char<9)>|'!' '+|-' 'BLEND| ' '+|-' <line|lcmb> '+|-' -> \\ &<line|lcmb>
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$
 * }


 * This keyword is used to define or redefine a surface.
 * Surfaces are defined with boundary lines.
 * Type gsur A001 + BLEND - L003 - L002 - L001 - L004
 * Type gsur A002 + BLEND - L007 - L006 - L005 + L003
 * Type gsur A003 + BLEND - L00C - L00A - L009 + L007
 * Type gsur A004 + BLEND + L004 - L00I - L00G + L00C



Creating Meshes

 * In the final part of this example we will create the meshes.
 * The options tr3, qu4, and qu8 are for meshing surfaces.
 * The options he8 and he20 are used to mesh bodies.


 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} Command: 'elty' &[ ] ['be2'|'be3'|'tr3'|'tr3u'|'tr6'|'qu4'|-> \\ &'qu8'|'he8'|'he8f'||'he8i'|'he8r'|'he20'|'he20r']
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$
 * }


 * Type elty all QU4
 * Type mesh all
 * We just have meshed the surface with an elements of qu4.


 * Type plot m all



Basic Models
Shown in the figures below are the example models that come installed with the CalculiX suite. We were able to regenerate these models by applying similar commands to the one present in the above tutorial. After playing around with this software we were finally able to obtain the right combination of commands the display the colorful models with attributes shown below.

Turbine: Illustrating the Use of Colors and Surface Element Attribute Specifications


=Problem 4.8: The Mass Matrix in Weak Form Using F&B problem 3.4=

Given: PDE and Corresponding Parameters
where the weak form is A trial solution for u(x) is suggested as and w(x) as

Find: The Mass Matrix in Weak Form
with n=3 and for dynamics:

Find $$\underline {F}(t)$$

State the Weak Form of the dynamic PDE for elasticity
Our given basis is The discretized forms of our trial solution and our weighting functions respectively are

Substitute back the operator definitions
=Contributing Members & Referenced Lecture=