User:Eml5526.s11.team4/HW6

= Problem 6.1 - QLEBF with G1DM1.0/D1b =

Problem Statement
Solve data set G1DM1.0/D1b using Quadrilateral Lagrangian Element Basis Functions(QLEBF)with uniform discretization (equidistant element nodes) nel (number of elements) = 2,4,6,8,... This problem is similar to Homework problems 5.1, 5.3, and 5.7. Recall: The data set G1DM1.0/D1b is defined by an elastic bar with variable cross-section: (Image from Charles Cook, HW1.1) $$\Omega=\left\{0,1\right\}$$ $$\Gamma_{g}=\left\{0\right\}, g = 4$$ $$\Gamma _{h}=\left\{1\right\}, h = 12$$ $$a_{2}\left\{x\right\}=2+3x$$ $$f\left\{x,t\right\}=5x$$ $$\frac{\partial^{2} u}{\partial t^{2}}=0$$ 1.) For nel = 2, compute $$\mathbf{\tilde{K}}=\sum_{e=1}^{2}\mathbf{\tilde{K}^{e}}$$ with $$\mathbf{\tilde{K}^{e}}$$ from Equation (1) from lecture notes page 30-5. Display $$\mathbf{\tilde{K}^{e}}$$, e = 1,2. 2.) Compute $$\mathbf{k}^{e}, \mathbf{L}^{e}$$ for e = 1,2. 3.) Compute $$\mathbf{\tilde{K}^{e}}=(\mathbf{L}^{e})^{T}\mathbf{k}^{e}\mathbf{L}^{e}$$ for e = 1,2. Compare to the results from sections 1. 4.) Plot all QLEBF for nel = 3. 5.) Plot $$u_{\tilde{n}}^{h}$$ versus u,$$ \left [u_{\tilde{n}}^{h}(.5)-u(.5) \right ]$$ versus $$\tilde{n}$$.

solution
1.For nel=2, comp $$\tilde{K}=\sum\limits_{e=1}^{2}$$ ,with by (1)p30.5,display $${{\tilde{K}}^{e}}$$, e=1,2 e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix} { 0} & {0} & {0} & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Then,sum on e=1,2:

$$\tilde{K}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {32.666666666666664}  &{-20.666666666666668} & {2.833333333333334}\\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

2.comp $${{k}^{e}}$$ , $${{L}^{e}}$$, for e=1,2 e=1:

$${{k}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668}  \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  \\ \end{matrix} \right]$$

$${{L}^{1}}=\left[ \begin{matrix} 1 & 0 & 0 & 0 & 0 \\   0 & 1 & 0 & 0 & 0  \\   0 & 0 & 1 & 0 & 0  \\ \end{matrix} \right]$$

e=2:

$${{k}^{2}}=\left[ \begin{matrix} {17.833333333333332} & {-20.666666666666668} & {2.833333333333334} \\   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

$${{L}^{2}}=\left[ \begin{matrix} 0 & 0 & 1 & 0 & 0 \\   0 & 0 & 0 & 1 & 0  \\   0 & 0 & 0 & 0 & 1  \\ \end{matrix} \right]$$

3.comp $${{{\tilde{K}}}^{e}}={{L}^{e}}^{T}{{k}^{e}}{{L}^{e}}$$, for e=1,2. compare to 1). e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix}  {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix}   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Compare with result on 1), we can see that the result is the same by using these two methods. 4.plot all QLEBF for nel=3 5.plot $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$,when number of element is 2,3,4,6,respectively. plot $$\displaystyle \left[ u_{nel}^{h}(0.5)-u(0.5) \right]$$ vs number of element

The error at x=0.5 can be listed in the following table: From the table, it can be found that the error will be lower than 10^(-6) when the number of elements is greater than n=14.

Comments:
It is observed that the QLEBF method converges to the accuracy of 10^(-6) with much less element number than LLEBF method in homework 5.7. This problem was solved by shengfeng yang

= Problem 6.2 - QLEBF with G1DM1.0/D1 =

Problem Statement
Solve data set G1DM1.0/D1 using Quadrilateral Lagrangian Element Basis Functions(QLEBF)with uniform discretization (equidistant element nodes) nel (number of elements) = 2,4,6,8,... This problem is similar to Homework problems 5.2, 5.4, and 5.6. Recall: The data set G1DM1.0/D1 is defined by an elastic bar with variable cross-section: (Image from Charles Cook, HW1.1) $$\Omega=\left\{0,1\right\}$$ $$\Gamma_{g}=\left\{1\right\}, g = 4$$ $$\Gamma _{h}=\left\{0\right\}, h = 12$$ $$a_{2}\left\{x\right\}=2+3x$$ $$f\left\{x,t\right\}=5x$$ $$\frac{\partial^{2} u}{\partial t^{2}}=0$$ 1.) For nel = 2, compute $$\mathbf{\tilde{K}}=\sum_{e=1}^{2}\mathbf{\tilde{K}^{e}}$$ with $$\mathbf{\tilde{K}^{e}}$$ from Equation (1) from lecture notes page 30-5. Display $$\mathbf{\tilde{K}^{e}}$$, e = 1,2. 2.) Compute $$\mathbf{k}^{e}, \mathbf{L}^{e}$$ for e = 1,2. 3.) Compute $$\mathbf{\tilde{K}^{e}}=(\mathbf{L}^{e})^{T}\mathbf{k}^{e}\mathbf{L}^{e}$$ for e = 1,2. Compare to the results from sections 1. 4.) Plot all QLEBF for nel = 3. 5.) Plot $$u_{\tilde{n}}^{h}$$ versus u,$$ \left [u_{\tilde{n}}^{h}(.5)-u(.5) \right ]$$ versus $$\tilde{n}$$.

solution
1.For nel=2, comp $$\tilde{K}=\sum\limits_{e=1}^{2}$$ ,with by (1)p30.5,display $${{\tilde{K}}^{e}}$$, e=1,2 e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix} { 0} & {0} & {0} & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Then,sum on e=1,2:

$$\tilde{K}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\  {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {32.666666666666664}  &{-20.666666666666668} & {2.833333333333334}\\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

2.comp $${{k}^{e}}$$ , $${{L}^{e}}$$, for e=1,2 e=1:

$${{k}^{1}}=\left[ \begin{matrix} {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668}  \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  \\ \end{matrix} \right]$$

$${{L}^{1}}=\left[ \begin{matrix} 1 & 0 & 0 & 0 & 0 \\   0 & 1 & 0 & 0 & 0  \\   0 & 0 & 1 & 0 & 0  \\ \end{matrix} \right]$$

e=2:

$${{k}^{2}}=\left[ \begin{matrix} {17.833333333333332} & {-20.666666666666668} & {2.833333333333334} \\   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$

$${{L}^{2}}=\left[ \begin{matrix} 0 & 0 & 1 & 0 & 0 \\   0 & 0 & 0 & 1 & 0  \\   0 & 0 & 0 & 0 & 1  \\ \end{matrix} \right]$$

3.comp $${{{\tilde{K}}}^{e}}={{L}^{e}}^{T}{{k}^{e}}{{L}^{e}}$$, for e=1,2. compare to 1). e=1: $${{{\tilde{K}}}^{1}}=\left[ \begin{matrix}  {10.833333333333334} & {-12.666666666666666} & {1.833333333333333} & {0} & {0} \\   {-12.666666666666666} & {29.333333333333332} & {-16.666666666666668} & {0} & {0} \\   { 1.833333333333333} & {-16.666666666666668} & {14.833333333333334}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\ \end{matrix} \right]$$ e=2: $${{{\tilde{K}}}^{2}}=\left[ \begin{matrix}   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} & {0}  & {0} & {0}\\   { 0} & {0} &   {17.833333333333332} & {-20.666666666666668} & {2.833333333333334}  \\   { 0} & {0} &   {-20.666666666666668} & { 45.333333333333336} & {-24.666666666666668}  \\   { 0} & {0} &   { 2.833333333333334} & {-24.666666666666668} & {21.833333333333332}  \\ \end{matrix} \right]$$ Compare with result on 1), we can see that the result is the same by using these two methods. 4.plot all QLEBF for nel=3 5.plot $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$,when number of element is 2,3,4,6,10,respectively.

plot $$\displaystyle \left[ u_{nel}^{h}(0.5)-u(0.5) \right]$$ vs number of element

The error at x=0.5 can be listed in the following table: From the table, it can be found that the error will be lower than 10^(-6) when the number of elements is greater than n=8.

Comments:
It is observed that the QLEBF method converges to the accuracy of 10^(-6) with much less element number than LLEBF method in homework 5.8. This problem was solved by shengfeng yang

=Problem 6.3 Reproduce a tutorial=

Assigned in the lecture slide Mtg33

Problem statement
Reproduce all steps in the Beede’s tutorial

Part1: Buliding Geometry and Meshing
Opening the SciTE text editor inplanted in the bConverged Calculix pacage, and paste the following content into it then save as an .fbd file, like beam.fbd

pnt p1 0 0 0 pnt p2 25 0 0 pnt p3 50 0 0 pnt p4 75 0 0 pnt p5 100 0 0 line l1 p1 p2 25 line l2 p2 p3 25 line l3 p3 p4 25 line l4 p4 p5 25 seta lines l l1 l2 l3 l4 swep lines sweplines tra 0 10 0 10 seta surfaces s A001 A002 A003 A004 swep surfaces swepsurface tra 0 0 1 1 elty all he8 mesh all plot m all

The beam created:



This beam.fbd will reproduce the same geometry property as stated in the tutorial. For detailed explanation of the lines, please refer to the the Calculix cgx online documentary or previous wikiversity tutorials by our team.

Part2: Exporting Mesh, Loads, and Boundary Conditions
Open the SciTE editor again and paste such content, then overwrite the previous beam.fbd file. The new file pretty much the same geometry and mesh as the previous but use different name for the mesh and exports the nodal information to a .msh file.

seto beam pnt p1 0 0 0 pnt p2 25 0 0 pnt p3 50 0 0 pnt p4 75 0 0 pnt p5 100 0 0 line l1 p1 p2 25 line l2 p2 p3 25 line l3 p3 p4 25 line l4 p4 p5 25 seta lines l l1 l2 l3 l4 swep lines sweplines tra 0 10 0 10 seta surfaces s A001 A002 A003 A004 swep surfaces swepsurface tra 0 0 1 1 setc beam elty beam he8 mesh beam send beam abq rot -z frame

Then followed the steps of the tutorial as it use the “interactive” style of operating Calculix.

The nodal information in beam.msh file:



Setting nodes for the fixed end:



Boundary information:



The nodal information for the load:



Part3: Writing an Input File for the CCX Solver
Use SciTE to create a new file named beam.inp and write down the following content,

*HEADING Model: Calculix Beam Input File *INCLUDE,INPUT=beam.msh *BOUNDARY *INCLUDE,INPUT=fixed_123.bou *MATERIAL,NAME=EL *ELASTIC 30000000, 0.3 *SOLID SECTION,ELSET=Ebeam,MATERIAL=EL *STEP *STATIC *DLOAD *INCLUDE,INPUT=load.dlo *NODE FILE U *EL FILE S *END STEP

Then click on the SciTE menu “Tools->Solve”, when the “Job done!” is prompt we can click on the menu Tools->Post Process thus we open the result file. Choose “Datasets->STRESS” and then “Entities->Mises” than we can see the result just same as the tutorial.



For detailed explanation on the .inp file, please refer to the Calculix ccx online documentary or our team’s previous homework. .

Part4: Post Processing Results in CGX
By manipulating the cgx module for post processing we have the following results:

The displacement distribution when the bending is seen:



The displacement distribution at a certain plain:



The Graph function of post processing:



Problem Statement
Recall matrix algebra: (1) $$(\mathbf{A}\mathbf{B})^{T}=\boldsymbol{B}^{T}\boldsymbol{A}^{T}$$ (2) $$(\boldsymbol{A}^{-1})^{T}=(\boldsymbol{A}^{T})^{-1}=\boldsymbol{A}^{-T}$$ Given the matrices: $$ \boldsymbol{A}=\begin{bmatrix} 1 & 1 & 1\\ 2 & -1 & 3\\ 3 & 2 & 6 \end{bmatrix} $$ $$ \boldsymbol{B}=\begin{bmatrix} 1 & 3 & 5\\ 1 & -4 & 1\\ 2 & 5 & 8 \end{bmatrix} $$ Verify the matrix algebra properties shown above (properties 1 and 2). Find and explain the syntax of Wolframalpha

Solution
Condition 1

(1a) is equal to (1b), this proves that $$(\mathbf{A}\mathbf{B})^{T}=\boldsymbol{B}^{T}\boldsymbol{A}^{T}$$

The wolfram alpha calculation for $$(\mathbf{A}\mathbf{B})^{T}$$ is given in the following pdf file

Click on the pdf file to see the second page.

The wolfram alpha calculation for $$\boldsymbol{B}^{T}\boldsymbol{A}^{T}$$ is given in the following pdf file

Click on the pdf file to see the second page.

The Wolfram alpha results match with the manually calculated results of condition 1 of problem statement.This validates the proof.

Condition 2

(2a) is equal to (2b).This proves that $$(\boldsymbol{A}^{-1})^{T}=(\boldsymbol{A}^{T})^{-1}=\boldsymbol{A}^{-T}$$

The wolfram alpha calculation for $$(\boldsymbol{A}^{-1})^{T}$$ is given in the following pdf file

The wolfram alpha calculation for $$(\boldsymbol{A}^{T})^{-1}$$ is given in the following pdf file

Click on the pdf file to see the second page.

The Wolfram alpha results match with the manually calculated results of condition 2 of problem statement.This validates the proof.

More details on the appropriate syntax can be found at the following link: http://reference.wolfram.com/mathematica/tutorial/BasicMatrixOperations.html

Problem Statement
Table of Gauss points and weights (used in Gauss Quadrature integration), from wikipedia

Verify this table against Fish and Belytschko's table on p. 89, and the NIST Handbook. Perform problems 4.6, 4.7, and 4.9 on p. 91.

Verification of Gauss points and weights
Evaluating the wikipedia table into decimal values:

From Fish and Belytschko, the numerical values are listed below:

From the NIST Digital Library of Mathematical Functions, the 5-point values are:

Comparing the three tables above it is clear they are in agreement. Further, the values are confirmed with the MATLAB script at the end of this problem which can calculate the points and weights for $$n$$ points.

FB Problem 4.6
Use Gauss quadrature to obtain exact values for the following integrals. Verify by analytical integration.

For each integration the number of points needed is given by

part a
Solved with Wolfram Alpha here

For the integration of this second order polynomial we need 2 points. Using the matlab script at the end of this problem's solution to get the quadrature points and weights we integrate with the following script.

The MATLAB scripts calculates the integration to be 25.333333333333329, which is within round-off error of the analytic result found (thus exact for a numerical result).

part b
$$\int_{-1}^{1}\left(\xi^4+2\xi^2\right)d\xi = 26/15$$

Solved with Wolfram Alpha here

For a fourth-order polynomial we need three integration points. The following code finds the integration result exactly.

The numerical integration provided a result of 1.733333333333333, a numerically exact result.

part c
The MATLAB function has been created and used in the above parts. The function can be seen at the end of this problem.

FB Problem 4.7
Use three-point Gauss quadrature to evaluate the following integrals. Compare to the analytical integral.

part a
NOTE: There is a typo in the text here. The kernel uses $$\xi$$ and $$x$$ inconsistently.

Evaluated with Wolfram Alpha here.

Evaluating numerical also finds a value of zero (0.0).

part b
NOTE: There is a typo in the text here. The kernel uses $$\zeta$$ and $$\xi$$ inconsistently.

Evaluated with Wolfram Alpha here.

Here, numerical integration with three points does not match the analytical result. The numerical result with three points is 1.529961646516273. This is expected as we are not integrated a polynomial but a trig function (Gauss quadrature is based upon interpolating polynomials of increasing order).

FB Problem 4.8
The integral $$\int_{-1}^{1}\left(3\xi^3+2\right)d\xi$$ can be integrated exactly using two-point Gauss quadrature. How is the accuracy affected if

part a
one-point quadrature is employed

part b
three-point quadrature

Solution
First note that the integration with two-points results in a value of 4.0. With one-point the value is also 4.0. With three-points the value is also 4.0 as expected.

This result is somewhat unexpected and is the result of coincidence and the integral mean value theorem. The integral mean value theorem states that there is a point within a domain that is valued the same as the integration value. For this integral, this point so happens to be at zero, the location where the sample point is taken for one quadrature point. Thus, for this particular equation, with this particular domain of integration, the single point predicts the exact value correctly.

Problem Statement
Problem 1: The heat conduction problem in Example 8.1 ( Fish and Belytschko ) is modeled with 16 quadrilateral finite elements.The problem is solved manually using the finite element code given in Section 12.5

Solution
1.) The problem is solved using the MATLAB code given in section 12.5 of companion site of the text book Fish and Belytschko.On running the code we receive the following post process result plots for temperature and flux figs (a),(b). Fluxes are calculated by looping over the number of elements.For the four node quadrilateral there are four Gauss points as shown in fig (c).







2.) Reproduce Example 8.4 (page 205) from Fish and Belytschko. 3.) Complete problem 8.6 (page 212) from Fish Belytschko.