User:EML4500.f08.JAMAMA/FE/mcdaniel

= 7.3 =

Given/Find: Static Solution of the Circle

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

Let Omega a circle with a one-unit radius. Find the static solution such that: $$T=4$$ $$f(\underline x)=1$$ in Omega $$g=0$$ on $$\Gamma _g=\partial \Omega $$
 * style="width:95%" |
 * style="width:95%" |

all to 10^-6 accuracy at the center. Quadrilateral and Triangular elements should be used but for only 1/4 of the circle. Plot the deformed shape in 3D. $$
 *  $$ \displaystyle
 * }

Solution: Static Solution of the Circle Using Quadrilateral and Triangular Elements to an Accuracy of 10^-6

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

The Partial Differential Equation for a vibrating membrane is the following via transient analysis: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }


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

T(\textrm{div \ grad \ u})+\textrm{f}=\rho \frac{\partial ^2u}{\partial t^2} $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *  $$ \displaystyle
 * }


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

Where: $$ u(\underline x,t)=$$ Transverse Displacement $$T=\quad$$ Tension (Force/Length) = Constant $$\textrm{f}(\underline x,t)=$$ Distributed Transverse Force $$\rho (\underline x)=$$ Mass Density (Mass/Area) & $$\textrm{div \ grad\ u}=\frac{\partial }{\partial x_i}\frac{\partial u}{\partial x_j}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }


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

The Weighted Residual Form is as follows: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }


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

\int_{\Omega } w\left [ T\frac{\partial ^2u}{\partial x_i\partial x_j}+\textrm{f}-\rho \frac{\partial ^2u}{\partial t^2} \right ]d\Omega =0 $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *  $$ \displaystyle
 * }


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

Which can thereby be broken into three separate terms: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }


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

\int_{\Omega } wT\frac{\partial ^2u}{\partial x_i\partial x_j}d\Omega =\mathbb{K} $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *  $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\int_{\Omega } w\textrm{f}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\int_{\Omega } w\rho \frac{\partial ^2u}{\partial t^2}\ d\Omega=\tilde{m}(w,u) $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

This final equation will be used later in the Weak Form. For now, the first term above (set equal to $$\mathbb{K}$$), is equivalent to: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\mathbb{K}=\int_{d\Omega } \frac{\partial }{\partial x_i}(wT\frac{\partial u}{\partial x_j})\ d\Omega-\int_{\Omega } \frac{\partial w}{\partial x_i}T\frac{\partial u}{\partial x_j}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

And by using divergence theorem, the first term of the equation above may be converted to yield: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\mathbb{K}=\int_{\partial \Omega } n_iwT\frac{\partial u}{\partial x_j}d(\partial \Omega)-\int_{\Omega } \frac{\partial w}{\partial x_i}T\frac{\partial u}{\partial x_j}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Now since $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\partial \Omega =\Gamma _g\cup \Gamma _h\ \textrm{but }\ w=0\ \ \textrm{on}\ \ \Gamma _g $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

and the following are true: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\frac{\partial w}{\partial x_i}=\textrm{grad}\ w=\bigtriangledown w $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\frac{\partial u}{\partial x_j}=\textrm{grad}\ u=\bigtriangledown u $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Then the K equation becomes: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\mathbb{K}=\int_{\Gamma _h}(n_i\ wT\frac{\partial u}{\partial x_j}\ \partial \Gamma _h)-\int_{\Omega }\bigtriangledown wT\bigtriangledown u\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Otherwise written as: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\mathbb{K}=\int_{\Gamma _h}w\ \underline nT\cdot \textrm{grad \ u}\ \partial \Gamma _h-T\int_{\Omega }\bigtriangledown w\bigtriangledown u\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

where the last term is conveniently used to write ultimately the weak form of the weighted residual form which is the following: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\tilde{\textrm{f}}(w)=\tilde{m}(w,u)+\tilde{k}(w,u)\ \ \ \forall w \ \ \ s.t.\ \ \ w=0\ \ \ on\ \ \ \Gamma _g $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\tilde{m}(w,u)=\int_{\Omega } w\rho \frac{\partial ^2u}{\partial t^2}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\tilde{k}(w,u)=T\int_{\Omega }\bigtriangledown w\bigtriangledown u \ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

\tilde{\textrm{f}}(w)=\int_{\Gamma _h}w \underline n \ T\ \textrm{grad u}\ \partial \Gamma _h+\int_{\Omega }w\textrm{f}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

From the Discrete Weighted Form, $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\underline M_{FF}\underline d_F+\underline K_{FF}\underline d_F=\underline F_F-(\underline M_{FE}g+\underline K_{FE}g)$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

But because g=0 from the boundary conditions, as well as M due to the static requirement, the displacements may be derived: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\underline K \underline d=\underline F$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\therefore \underline d=\underline K^{-1}\underline F$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Now to solve for the displacements: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$N_I(x,y)=L_{i,m}(x)\cdot L_{j,n}(y) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Where $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$I=i+(j-1)m $$ and $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$L_{i,m}(x)=\prod_{K=1, K\neq i}^{m}\left ( \frac{x-x_k}{x_i-x_k} \right ) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

And setting the Shape Function to the following: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$N=\begin{bmatrix} \frac{1}{4}(1-\xi _1)(1-\xi _2)\\ \frac{1}{4}(1+\xi _1)(1-\xi _2)\\ \frac{1}{4}(1-\xi _1)(1+\xi _2)\\ \frac{1}{4}(1-\xi _1)(1+\xi _2) \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The deformed shape coordinates may further be derived from these following two equations, where the summations are taken across all coordinate points taken from Abaqus output files: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$X_{1}^{e}=\varphi_{1}^{e}(\xi )=\sum_{I=1}^{nel}N_I(\xi )x_{1,I}^{e} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$(y=)X_{2}^{e}=\varphi_{2}^{e}(\xi )=\sum_{I=1}^{nel}N_I(\xi )x_{2,I}^{e} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Using the following Jacobian matrix: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$J=\left [ \frac{\partial x_i}{\partial \xi _j} \right ]=\begin{bmatrix} \frac{\partial x_1}{\partial \xi _1} & \frac{\partial x_1}{\partial \xi _2}\\ \frac{\partial x_2}{\partial \xi _1} & \frac{\partial x_2}{\partial \xi _2}\\ \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The following equations may now be utilized: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\underline B_I(\xi )=\bigtriangledown _xN_{I}^{e}(\xi )=(\underline J^e)^{-T}(\xi )\bigtriangledown_\xi N_{I}^{e}(\xi ) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$k_{I3}^{e}=\int_{\omega =el}N_I(\xi )\cdot T\cdot \underline B_J (\xi )J(\xi )d(el) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$F_i=\int_{\omega =el}N_I(\xi )f\ det(J(\xi ))d(el) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Implementing MATLAB at this point, in conjunction with output data files from Abaqus which contain element coordinate points for three separate quadrilateral meshes and three separate triangular meshes (included), the final displacements may be derived, to an accuracy of 10^-6: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The above MATLAB script for the 6 different meshes generate the following results: $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }