User:Eml5526.s11.team01/Hwk4

=Problem 4.1: Solving using WRF with a change in Basis Function= Solving a differential equation using Weight Residual Form.

Given
The partial differential equation is:

The boundary conditions are:

$$ u\left(1\right) = 0$$

$$\frac{du}{dx}\left(0\right) = -4 $$

As weighting function should be used:

Find

 * 1. Find an approximate solution of the form $$U^{h}=\sum_{i=0}^{n}d_{i}b_{i}$$ for $$\displaystyle n=2$$
 * 2. Find two equations that enforce the boundary conditions
 * 3. Project the weight residues
 * 4. Display the equations in matrix form
 * 5. Solve for $$\displaystyle d $$
 * 6. Construct $$\displaystyle U^{h}$$ and plot $$\displaystyle U^{h}$$ and $$\displaystyle u$$
 * 7. Repeat 1 to 6 for $$\displaystyle n=4$$ and $$\displaystyle n=6$$

Find an approximate solution with n=2
When n = 2 we have:
 * {| style="width:100%" border="0"

$$d_{i}=\left[ \begin{matrix} d_{0} & d_{1} & d_{2} \\ \end{matrix} \right]$$
 * $$\displaystyle (4.1.3) $$
 * }
 * }

and


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

$$b_{i}=\left[ \begin{matrix} 1 & \sin \left( x \right) & \sin \left( 2x \right) \\ \end{matrix} \right]$$
 * $$\displaystyle (4.1.4) $$
 * }
 * }

Therefore:


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


 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |

$$U^{h}=d_{0}+d_{1}\sin \left( x \right)+d_{2}\sin \left( 2x \right)$$
 * $$\displaystyle (4.1.5) $$
 * style= |
 * }

Find two equations that enforce the boundary conditions
The boundary conditions are:
 * {| style="width:100%" border="0"


 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |

$$u\left( 1 \right)=0=d_{0}+d_{1}\sin \left( 1 \right)+d_{2}\sin \left( 2 \right)$$ Differentiating
 * $$\displaystyle (4.1.6) $$
 * style= |
 * }
 * {| style="width:100%" border="0"



$$\frac{du}{dx}=0+d_{1}\cos \left( x \right)+2d_{2}\cos \left( 2x \right)$$
 * $$$$
 * }
 * }

And evaluating


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

$$\frac{du}{dx}(0)=-4=d_{1}+2d_{2}$$
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * $$\displaystyle (4.1.7) $$
 * style= |
 * }

Project the weight residues
In this case changing the function $$u\left( x \right)$$ by the approximated function $$\displaystyle U^{h}$$ in the PDE $$\left(4.1.1\right)$$, we found that $$P\left( u \right)$$ is:


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

$$P\left( U^{h} \right)=2\frac{\partial ^{2}U^{h}}{\partial x^{2}}+3=-2d_{1}\sin \left( x \right)-8d_{2}\sin \left( 2x \right)+3=0$$
 * $$\displaystyle (4.1.8) $$
 * }
 * }

$$\left(4.1.8\right)$$ can be written as:
 * {| style="width:100%" border="0"



$$P\left( U^{h} \right)=\left[ \begin{matrix} 0 & -2d_{1}\sin \left( x \right) & -8d_{2}\sin \left( 2x \right) \\ \end{matrix} \right]\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]+3$$
 * $$\displaystyle (4.1.9) $$
 * }
 * }

Projecting the residue we have:


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

$$ \int_{a}^{b}b(x)P\left ( U^{h} \right )= 0 $$ After substitution of $$\left(4.1.4\right) and \left(4.1.9\right) into \left(4.1.10\right)$$, the following equation is obtained
 * $$\displaystyle (4.1.10) $$
 * }
 * }
 * {| style="width:100%" border="0"


 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |

$$\int_{0}^{1}{\left[ \begin{matrix} 1 \\   \sin \left( x \right)  \\ \sin \left( 2x \right) \\ \end{matrix} \right]}\left[ \begin{matrix} 0 & 2d_{1}\sin \left( x \right) & 8d_{2}\sin \left( 2x \right) \\ \end{matrix} \right]dx\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=3\int_{0}^{1}{\left[ \begin{matrix} 1 \\   \sin \left( x \right)  \\ \sin \left( 2x \right) \\ \end{matrix} \right]}dx$$
 * $$\left(4.1.11\right)$$
 * style= |
 * }

Display the equations in matrix form
Performing the product and then the integration indicated in $$\left(4.1.1\right)$$, we have
 * {| style="width:100%" border="0"


 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |

$$\left[ \begin{matrix} 0 & 2\sin \left( x \right) & 8\sin \left( 2x \right) \\ 0 & 2\sin ^{2}\left( x \right) & 8\sin \left( x \right)\sin \left( 2x \right) \\ 0 & 2\sin \left( x \right)\sin \left( 2x \right) & 8\sin ^{2}\left( x \right) \\ \end{matrix} \right]\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=\left[ \begin{matrix} 3 \\   1.3791  \\   2.1242  \\ \end{matrix} \right]$$
 * $$\left(4.1.12\right)$$
 * style= |
 * }

The latter equation is clearly of the form $$\displaystyle Kd=F$$. We can observe that matrix K is not symmetric.

Solve for $$\displaystyle d $$
In equation $$\left(4.1.12\right)$$ we have three equations; but, as we need to enforce the boundary conditions we take the only one of them and solve it together with equations $$\left(4.1.6\right)$$ and $$\left(4.1.7\right)$$ From these equations we now have the system:


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



$$\left[ \begin{matrix} 1 & .8415 & .9093 \\   0 & 1 & 2  \\   0 & .7944 & 4.7568  \\ \end{matrix} \right]\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=\left[ \begin{matrix} 1 \\   -4  \\   2.1242  \\ \end{matrix} \right]$$
 * $$\left(4.1.13\right) $$
 * }
 * }

After solving the system in $$\left(4.1.13\right)$$ we obtain d:
 * {| style="width:100%" border="0"


 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |

$$\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=\left[ \begin{matrix} 4.6608 \\   -7.3471  \\   1.6735  \\ \end{matrix} \right]$$
 * $$$$
 * style= |
 * }

Construct $$U^{h}$$ and plot $$U^{h}$$ and $$u$$
Replacing $$\displaystyle d_{i}$$ in $$\left(4.1.5\right)$$ we obtain the approximated solution:


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

$$ \displaystyle U^{h}=4.6608-7.3471x-1.6735x^{2}$$


 * $$\left(4.1.14\right)$$
 * }
 * }

On the other hand, the exact solution of the PDE $$\left(4.1.1\right)$$ with the given boundary condition is:


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

$$ u\left ( x \right )=-\frac{3}{4}x^{2}-4x+\frac{19}{4} $$
 * $$\left(4.1.15\right)$$
 * }
 * }

We observe that the approximated and the exact solution are not the same. Perhaps they would be closer if we increased the number of terms. Figure 1.1 we show the exact and the approximated solution.



Repeat 1.1 to 1.6) for i = 4 and i = 6
Now we repeat the procedure using n = 4 and n=6 in Eq.1.4

As the procedure is the same shown before, we develop a routine using Matlab to perform the products, the integration and to solve the system of equations. The routine is shown below in the appendix. Figures 1.2 and 1.3 show the results of n=4 and n=6, respectively. Figure 1.4 shows the error versus the number of terms.



We see that with increase in terms, the error increases as well.

Appendix
 Matlab Code: 

Contributing Members
Solved and posted by James Roark 22:14, 28 February 2011 (UTC)

=Problem 4.2: Weak form for the equation of heat conduction= Fish and Belystschko page73, Problem 3.5.

Obtain the weak from for the equation of heat conduction with the boundary conditions $$ T\left(0\right) = 100$$ and $$ q\left(10\right) = h\bar{T} $$. The condition on the right is a convection condition.

Multiplying by the weight function
Multiplying the first two equations by the weight funtion we can write:

Applying integration by parts
Applying integration by parts to $$\left(4.2.2a\right)$$ we ca write:

Therefore we can write that:

Using the natural boundary condition we can write:

Which yields

Weak form for 1D heat conduction problems
Therefore the weak form for 1D heat conduction problems can be written as follows:

Contributing Members
Solved and posted by Cedric Adam, 1 March 2011

=Problem 4.3: Weak form for heat conduction problem in a circular plate=

Given
Fish and Belystschko page73, Problem 3.6.

Given the strong form for the heat conduction problem in a circular plate:

Natural boundary condition :

Essential boundary condition :

where $$ \displaystyle R $$ is the total radius of the plate, $$\displaystyle s $$ is the heat source per unit length along the plate radius, $$\displaystyle  T $$ is the temperature and $$\displaystyle  k $$ is the conductivity.

Find
Assume $$ k, ~s~ $$ and $$ ~R~ $$ are given:

1. Construct the weak form for the above strong form.

2. 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.

3. Solve the differential equation with the boundary conditions and show that the temperature distribution along the radius is given by

Construction of the weak form
First we multiply the equation in the strong form by the weight function and integrate over the domain.

The middle term cancels out when evaluated from $$\displaystyle 0 to  R $$ which yields

Quadratic trial (candidate) solutions
The trial solution $$ T\left ( r \right )$$ and the weight function $$ w\left ( r \right )$$ can be written as follows:

Using the boundary conditions, the coefficients of the trial solutions can be solved. First the you take the derivative of the trial solution and set it to the natural boundary condition.

Therefore $$\displaystyle \alpha_1 = 0 $$

Now use the essential boundary

Therefore

which solves for our two temperature functions

Now solve for the weight functions using $$ w \left( R \right) = 0 $$

Which yields

And the derivative of the weight function is

Before substitution, when the middle term is evaluated from $$ \left( 0 to R \right) $$ it goes to zero

Therefore we can cancel out the middle term and substitute in the rest

Integrating we get:

After evaluating the domain, the equation simplifies to:

Factoring out $$\displaystyle \beta_1 $$ and $$\displaystyle \beta_2 $$ we get:

Assuming that $$\displaystyle \beta_1 $$ and $$\displaystyle \beta_2 $$ equal zero, $$\displaystyle \alpha_2 $$ can be solved:

Now $$\displaystyle \alpha_2 $$ can be substituted into $$ T \left( r \right) $$ to obtain the solution of the weak form

Solve the differential equation
Using the strong form the differential equation can be solved by taking the integral:

Using the natural boundary condition $$ \frac{dT}{dr} \left(r = 0 \right) = 0 \rightarrowc_1 = 0 $$

Integrate once again we get:

Using the essential boundary condition $$ T \left( r = R \right) = 0 $$ solves for:

substituting in $$\displaystyle c_2 $$ yields the temperature distribution

Contributing Members
Solved and posted by Harrison Sheffield 12:19, 1 March 2011 (UTC)

=Problem 4.4: Weak form for circular bar in torsion=

Given
Fish and Belystschko page73, Problem 3.7.



The strong form for the circular bar in torsion (Figure 4.1) is the following:


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


 * style="width:95%" |
 * style="width:95%" |

$$\frac{d}{dx}\left( JG\frac{d\phi }{dx} \right)+m=0,\quad \quad \quad 0\le x\le 1,$$


 *  $$\displaystyle (4.4.1) $$
 * }
 * }

where $$\displaystyle JG=(2+3x),\quad m(x)=5x .$$

Natural boundary condition:


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


 * style="width:95%" |
 * style="width:95%" |

$$\frac{d\phi }{dx}(x=1)=-6,$$


 * <p style="text-align:right;"> $$\displaystyle (4.4.2) $$
 * }
 * }

Essential boundary condition:


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


 * style="width:95%" |
 * style="width:95%" |

$$\phi \left( x=0 \right)=\bar{\phi }=4,$$


 * <p style="text-align:right;"> $$\displaystyle (4.4.3) $$
 * }
 * }

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

Find
1. Construct the weak form for the circular bar in torsion.

2. Integrate the differential equation given above. Find the integration constants using boundary conditions.

3. Plot the result of $$\phi(x)$$.

Construct the weak form for the circular bar in torsion
As the weight function must vanish on the essential boundaries, we consider all smooth weight functions w(x) such that w(0)=0. The trial solutions must satisfy the essential boundary condition $$\phi \left( x=0 \right)=4,$$.

We start by multiplying the governing equation and the natural boundary condition over the domains where they hold by an arbitrary weight function:
 * {| style="width:100%" border="0"



$$ \int_{0}^{1} w [\frac{\mathrm{d} }{\mathrm{d} x}( (2+3x)\frac{\mathrm{d} \phi}{\mathrm{d} x})+5x]dx = 0 \forall w(x), $$


 * <p style="text-align:right;">$$\displaystyle (4.4.4) $$
 * }
 * }


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



$$ [w(\frac{\mathrm{d} \phi}{\mathrm{d} x}+6 )]_{x=1} = 0   \forall w(1). $$


 * <p style="text-align:right;">$$\displaystyle (4.4.5) $$
 * }
 * }

Next we integrate (4.4.4) by parts, then we get the following result:


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



$$ \int_{0}^{1} w [\frac{\mathrm{d} }{\mathrm{d} x}((2+3x)\frac{\mathrm{d} \phi}{\mathrm{d} x})]dx = \left. \left( w(2+3x)\frac{\text{d}\phi }{\text{d}x} \right) \right|_{x=0}^{x=1} - \int_{0}^{1} \frac{\mathrm{d} w }{\mathrm{d} x}(2+3x)\frac{\mathrm{d} \phi}{\mathrm{d} x}dx $$


 * <p style="text-align:right;">$$\displaystyle (4.4.6) $$
 * }
 * }

We have constructed the weight functions so that w(0)=0; therefore, the first term on the RHS of (4.4.6) vanishes at x=0. Substituting (4.4.6) into (4.4.4) gives:


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



$${{\left. \left( w(2+3x)\frac{\text{d}\phi }{\text{d}x} \right) \right|}_{x=1}}-\int_{0}^{1}{\frac{\text{d}w}{\text{d}x}}(2+3x)\frac{\text{d}\phi }{\text{d}x}dx+5\int_{0}^{1}{wxdx}=0 \forall w(x)\quad with\quad w(0)=0.$$


 * <p style="text-align:right;">$$\displaystyle (4.4.7) $$
 * }
 * }

Substituting (4.4.5) into the first term of (4.4.7) gives:

Integrate the differential equation and find the integration constants using boundary conditions
(Eq.4.1) yields:

Apply the boundary condition, $$\frac{d\phi }{dx}(x=1)=-6$$ gives:

So,

Apply the boundary condition $$\phi \left( 0 \right)=4  $$ gives:

Then we have the following result:

MATLAB Source for Plot
The following MATLAB code was used to plot the result.

Contributing Members
Solved and posted by Yu Hao 1 March 2011

=Problem 4.5: Various Basis Functions satisfying the Constraint Breaking Solution=

Given

 * a) Constraint Breaking Solution which needs to be satisfied is $$b_{0}(\beta)\neq 0$$ and $$b_{i}(\beta)=0$$, i=1,2,3,....,n.
 * b) Consider $$\Omega=[-2,4]$$

Find
For $$\Gamma _{g}=\left \{ \beta\right \}$$ and each of the following Families of Basis Functions Find the Corresponding Basis Functions Satisfying the Given Constraint Breaking Solution. Plot original and required basis function for $$j=1,2,3$$
 * 1. $$F_{p}:=\left \{ x^{j}, j=0,1,2,...\right \}$$.
 * 2. $$F_{c}:=\left \{ \cos(jx), j=0,1,2,... \right \}$$.
 * 3.$$F_{s}:=\left \{ 1,\sin(jx), j=1,2,... \right \}$$.
 * 4. $$F_{f}:= \left \{\cos(jx),j=0,1,2,...\sin(jx),j=1,2,...\right\}$$
 * 5. $$F_{e}:=\left \{ e^{jx},j=0,1,2,.. \right \}$$.
 * 6. Show that $$\left \{ e^{jx},j=0,1,2,.. \right \}$$ is Linear Independent on $$\Omega$$.

1) Find the Corresponding Polynomial Basis Function $$F_{p}:=\left \{ x^{j}, j=0,1,2,...\right \}$$ Satisfying the Given Constraint Breaking Solution
If we choose $$b_{i}\left ( x \right )=\left(x-\beta\right)^{j}$$, it satisfies the given Constraint Breaking Solution as below,

and

So the Polynomial Basis Function satisfying the given constraint breaking solution is,


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

$$b_{j}\left ( x \right )=\left \{\left(x-\beta\right)^{j}, j=0,1,2,...\right \}$$ Plot Between the original and the required basis function is plotted in the MATLAB using the following code.
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * <p style="text-align:right;">$$\displaystyle (4.5.3) $$
 * style= |
 * }

Appendix
 Matlab Code: 



2) Find the Corresponding Cosine Basis Function $$F_{c}:=\left \{ \cos(jx), j=0,1,2,... \right \}$$ Satisfying the Given Constraint Breaking Solution
If we select $$b_{j}\left ( x \right )=\cos(j(x-\beta)), j=0$$

It will yield, $$b_{0}\left ( \beta\right)=\cos(0)=1\neq0$$

and further

$$b_{j}\left ( x \right )=\cos(\frac{\pi}{2}(2j+1)\frac{x}{\beta}), j=1,2,...$$

it will yield,

$$b_{j}\left ( \beta \right )=\cos(n\frac{\pi}{2})=0, n = 1, 3, 5, ...$$

So the Cosine Basis Function satisfying the given constraint breaking solution is,


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

$$b_{j}\left ( x \right )=\left \{\cos(j(x-\beta)), j=0\right \}$$
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |

&

$$b_{j}\left ( x \right )=\left \{\cos(\frac{\pi}{2}(2j+1)\frac{x}{\beta}), j=1,2,...\right \}$$ Plot of the original and the required Basis Function is plotted in the MATLAB using the following code
 * <p style="text-align:right;">$$\displaystyle (4.5.4) $$
 * style= |
 * }

Appendix
 Matlab Code: 



3) Find the Corresponding Sine Basis Function $$F_{s}:=\left \{1,\sin(jx), j=1,2,... \right \}$$ Satisfying the Given Constraint Breaking Solution
If we select,

$$b_{j}\left ( x \right )=\sin(j\pi\frac{x}{\beta}), i=1,2,...$$

It yields,

$$b_{j}\left ( \beta \right )=\sin(n\pi)=0, n = 1,2,3,...$$

So the Sine Basis Function satisfying the given constraint breaking solution is,


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

$$b_{j}\left ( x \right )=\left \{1,\sin(j\pi\frac{x}{\beta}), j=1,2,...\right \}$$
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * <p style="text-align:right;">$$\displaystyle (4.5.5) $$
 * style= |
 * }

Appendix
 Matlab Code: 



===4) Find the Corresponding Fourier Basis Function $$F_{f}:=\left \{\cos(jx),j=0,1,2,... \sin(jx),j=1,2,...\right\}$$ Satisfying the Given Constraint Breaking Solution===

If we select, According to $$\displaystyle (4.5.4)$$

$$b_{j}\left ( x \right )=\cos(j(x-\beta)), j=0$$

It will yield, $$b_{0}\left ( \beta\right)=\cos(0)=1\neq0$$

and further

$$b_{j}\left ( x \right )=\cos(\frac{\pi}{2}(2j+1)\frac{x}{\beta}), j=1,2,...$$

it will yield,

$$b_{j}\left ( \beta \right )=\cos(n\frac{\pi}{2})=0, n=1,3,5,...$$

Also,

If we select, According to $$\displaystyle (4.5.5)$$,

$$b_{j}\left ( x \right )=\sin(j\pi\frac{x}{\beta}), i=1,2,...$$

It yields,

$$b_{j}\left ( \beta \right )=\sin(n\pi)=0,n=1,2,3,...$$

Eventually,

The Fourier Basis Function Satisfying the given Constraint Breaking Solution is,


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

$$b_{j}\left ( x \right )=\left \{1,\cos(\frac{\pi}{2}(2j+1)\frac{x}{\beta}), j=1,2,...,
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |

\sin(j\pi\frac{x}{\beta}), j=1,2,...\right \}$$ Plot of the original and the required Basis Function is plotted in the MATLAB using the following code.
 * <p style="text-align:right;">$$\displaystyle (4.5.6) $$
 * style= |
 * }

Appendix
 Matlab Code: 







===5) Find the Corresponding Exponential Basis Function $$F_{e}:=\left \{ e^{jx},j=0,1,2,... \right\}$$ Satisfying the Given Constraint Breaking Solution===

If we select,

$$b_{j}\left ( x \right)=\left \{e^{jx\beta}, i=0\right \}$$

and further,

$$b_{j}\left ( x \right)=\left \{e^{j(x-\beta)}-1, i=1,2,...\right \}$$

Eventually, The Fourier Basis Function Satisfying the Given Constraint Breaking Solution is,


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

$$b_{j}\left ( x \right )=\left \{1,e^{j(x-\beta)}-1, i=1,2,...\right \}$$ Plot of the original and the required Basis Function is plotted in the MATLAB using the following code.
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * <p style="text-align:right;">$$\displaystyle (4.5.7) $$
 * style= |
 * }

Appendix
 Matlab Code: 



6) Show that $$\left \{ e^{jx},j=0,1,2,.. \right \}$$ is Linear Independent on $$\Omega$$
We have,

Calculating $$\Gamma (F)$$ for the given interval $$\Omega [-2,4]$$,we get

Therefore,

To check the Linear Independent Behavior we need to check the Determinant of this matrix. So calculating the determinant of this matrix in the MATLAB Using the Following code,

Appendix
 Matlab Code: 

The Determinant arrives as $$1.1145\times (10)^{9}$$ which is a Non-zero number.

So we can say that, $$\left \{ e^{jx},j=0,1,2,.. \right \}$$ is Linear Independent on $$\Omega$$.
=Problem 4.6: Strong and weak form for elastic bar with variable distributed spring=

Given
Fish and Belystschko page74, Problem 3.9.

Consider an elastic bar with a variable distributed spring $$\displaystyle p\left( x \right)$$ along its length as shown in Figure 6.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 $$\displaystyle A\left( x \right)$$, Young’s modulus $$\displaystyle  E\left( x \right)$$ with body force $$\displaystyle  b\left( x \right)$$ and boundary conditions as shown in Figure 6.1

Find

 * 1. Construct the strong form.
 * 2. Construct the weak form.

Construct the strong form
From figure 6.1, construct a free body diagram Balance the forces in figure 6.2 to obtain:


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

$$\left[ p\left( x+\frac{dx}{2} \right)+b\left( x+\frac{dx}{2} \right) \right]dx+F\left( x+dx \right)-F\left( x \right)=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.1) $$
 * }
 * }

Rearranging:


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

$$\left[ p\left( x+\frac{dx}{2} \right)+b\left( x+\frac{dx}{2} \right) \right]+\frac{F\left( x+dx \right)-F\left( x \right)}{dx}=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.2) $$
 * }
 * }

Let $$dx\to 0$$, therefore


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

$$\left[ p\left( x \right)+b\left( x \right) \right]+\frac{dF}{dx}=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.3) $$
 * }
 * }

Note that $$F=\sigma A\left( x \right)=E\left( x \right)\varepsilon A\left( x \right)=E\left( x \right)A\left( x \right)\frac{du}{dx}$$


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

$$\left[ p\left( x \right)+b\left( x \right) \right]+\frac{d}{dx}\left[ E\left( x \right)A\left( x \right)\frac{du}{dx} \right]=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.4) $$
 * }
 * }

Apply boundary conditions, let $$u\left( l \right)=\bar{u}$$and if $$\bar{t}>0$$ at $$x=l$$, then $$\bar{t}=\left. E\left( x \right)\frac{du}{dx} \right|_{x=0}=0$$. We arrive at the strong form


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

$$\left[ p\left( x \right)+b\left( x \right) \right]+\frac{d}{dx}\left[ A\left( x \right)E\left( x \right)\frac{du}{dx} \right]=0$$
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * <p style="text-align:right;">$$\displaystyle (4.6.5a) $$
 * style= |
 * }
 * {| style="width:100%" border="0"

$$A\left( x \right)\left( E\left( x \right)\frac{du}{dx}+\bar{t} \right)=0$$
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * <p style="text-align:right;">$$\displaystyle (4.6.5b) $$
 * style= |
 * }

Construct the weak form
Let $$p\left( x \right)+b\left( x \right)=\hat{b}\left( x \right)$$in equation 6.5a.


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

$$\frac{d}{dx}\left[ A\left( x \right)E\left( x \right)\frac{du}{dx} \right]+\hat{b}\left( x \right)=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.6) $$
 * }
 * }

Following the procedure outlined in Fish & Belytschko section 3.2, page 47, multiply equations 6.5a and 6.5b by an arbitrary function $$w\left( x \right)$$such that $$w\left( l \right)=0$$and integrate to obtain


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

$$\int\limits_{0}^{l}{w\left( x \right)\left[ \frac{d}{dx}\left[ A\left( x \right)E\left( x \right)\frac{du}{dx} \right]+\hat{b}\left( x \right) \right]}=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.7a) $$
 * }
 * }


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

$$\left. w\left( x \right)A\left( x \right)\left( E\left( x \right)\frac{du}{dx}+\bar{t} \right) \right|_{x=0}=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.7b) $$
 * }
 * }

Rewrite equation 4.6.7a


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

$$\int\limits_{0}^{l}{w\left( x \right)\frac{d}{dx}\left[ A\left( x \right)E\left( x \right)\frac{du}{dx} \right]}dx+\int\limits_{0}^{l}{w\left( x \right)\hat{b}\left( x \right)}dx=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.8) $$
 * }
 * }

Integrate the first term of equation 4.6.8 by parts


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

$$\left. \left[ w\left( x \right)\underbrace{A\left( x \right)E\left( x \right)}_{\sigma } \right] \right|_{0}^{l}-\int\limits_{0}^{l}{\frac{dw}{dx}A\left( x \right)E\left( x \right)\frac{du}{dx}}dx+\int\limits_{0}^{l}{w\left( x \right)\hat{b}\left( x \right)}dx=0$$
 * <p style="text-align:right;">$$\displaystyle (4.6.9) $$
 * }
 * }

Since $$w\left( l \right)=0$$, this term drops out as zero. Noting that $$\left. {\bar{t}} \right|_{x=0}=\left. \sigma \right|_{x=0}$$


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

$$\int\limits_{0}^{l}{\frac{dw}{dx}A\left( x \right)E\left( x \right)\frac{du}{dx}}dx=\left. \left[ w\left( x \right)A\left( x \right)\bar{t} \right] \right|_{0}+\int\limits_{0}^{l}{w\left( x \right)\hat{b}\left( x \right)}dx$$
 * <p style="text-align:right;">$$\displaystyle (4.6.10) $$
 * }
 * }

Replacing $$\hat{b}$$ we arrive at the weak form


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

$$\int\limits_{0}^{l}{\frac{dw}{dx}A\left( x \right)E\left( x \right)\frac{du}{dx}}dx=\left. \left[ w\left( x \right)A\left( x \right)\bar{t} \right] \right|_{0}+\int\limits_{0}^{l}{w\left( x \right)\left[ b\left( x \right)+p\left( x \right) \right]}dx$$
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * style="width:2%; padding:1px; border:2px solid #FF0000" |
 * <p style="text-align:right;">$$\displaystyle (4.6.11) $$
 * style= |
 * }

Contributing Members
Solved and posted by James Roark 02:51, 2 March 2011 (UTC)

=Problem 4.7: working with calculix. From lecture 23 =

Given
Calculix is a nonlinear Finite Element code, open source, Abaqus-like

Find

 * 1. From the page http://www.dhondt.de/ download and install Calculix
 * 2. Read the manual, sign up with user group to as questions if any. Also access archive.
 * 3. Reproduce basic examples: disk, cylinder, sphere, sphere-volume, airfril.
 * 4. write a report for dummies explaining how to install and all the commands used in basic examples.

Download and instal Calculix
The program was succesfully downloaded and instaled

Read the manual, sign pu with user group t ask questions if any
The manual was downloaded but it was not succesfully opened from the page given in the problem, but the manual was found in the page http://web.mit.edu/calculix_v2.0/CalculiX/cgx_2.0/doc/cgx/index.html where it could be easily opened. From there the manual was read.

Reproduce basic examples: disk, cylinder, sphere, sphere-volume, airfril
The basics examples were reproduced: disc, cylinder, sphere, sphere_seg, and sphere_vol they are shown in above in the attached document.

Write report for "dummies" explining how to install and run cgx and all the commands used in basic examples
A document with a simple explanation about the installation and the construction of the simple examples listed above is shown in the pdf document

http://upload.wikimedia.org/wikiversity/en/3/35/Fe1.s11.team1.hw4.Basic_Calculix_Manual.pdf

It was found difficult to work with this interface because not icons or graphic interfaces are available, instead you have to remember all the commands and its sintaxis. to counter the beforementioned problem there are some better options like using a program where the pre and post process is made in a very easy way by using a completely graphic interphase and you can write your own conde in the program that you prefer. One of those options is the program GID which can be download free going to the page:

http://www.gid-usa.com/index.html

With the free version of GiD you can work with a limited number of nodes, but it is enough to solve simple problems that are usefull for your training.

Contributing Members
Solved and posted by Fernando Casanova 18:50, 28 February 2011 (UTC)

=Problem 4.8: Based off Problem 3.4=

Given
Based off Problem 3.4

Given the strong form

and

with

Find
Find $$ \tilde{M} $$

Assuming

$$ A = 1, E = 2, \bar{m} = 3 $$

Find $$ F \left(t \right) $$

Assuming

$$ n = 3 $$

$$ u^{h} \left ( p,t \right ) = g\left ( t \right ) = sin2t $$

From the strong form
From the strong form we can easily obtain the weak form. Find $$\displaystyle u(x) $$ smooth enough and satisfy the essential boundary condition:

such that:

since

Thus

From discrete weak form
From discrete weak form we get:

=Contributing Members=