User:Eml5526.s11.team4/HW5

= Problem 5.1 - G1DM1.0/D1b using WF =

Refer to lecture slide [[media:fe1.s11.mtg26.djvu|mtg-26]] for the problem statement.

Problem Statement
Solve G1DM1.0/D1b using WF with appropriate basis function (Poly, Fourier, Exp), until convergence of $$u^h(0.5)$$ to $$O(10^{-6})$$.

Boundary Conditions
The boundary conditions are

Solution
On solving the given PDE with the boundary conditions we get ,the exact equation for u(x) is as follows

Polynomial
The polynomial basis function is given by $$\mathfrak{F}_{p}=\left\{x^{j},j=0,1,...,n\right\}$$. For $$\left\{j=0,1,2,3\right\}$$ the basis functions $$\left\{b_{j}\right\}$$ are as follows: The constraint breaking solution$$\left\{\bar{b}_j\right\}$$ for the polynomial basis function is shown below: In the specific case where $$\left\{\beta=0\right\}$$ the CBS becomes:

The basis functions$$\left\{\bar{b}_j\right\}$$ satisfy the CBS requirements since$$\left\{\bar{b}_0(\beta)=\text{constant}\neq0\right\}$$ and $$\left\{\bar{b}_j(\beta)=0, for j=1,2,...,n\right\}$$.

Matlab code for polynomial function
The following figure shows the original polynomial basis functions and the CBS basis functions for $$\left\{j=1,2,3\right\}$$,We can see that convergence occurs at an n value of 8

The following figure shows the error vs n

Exponential
The exponential basis function is given by So, the basis functions $$\left\{b_j\right\}$$ are as follows: The corresponding basis functions $$\left\{\bar{b}_j\right\}$$ satisfying CBS are shown below: The basis functions$$\left\{\bar{b}_j\right\}$$ satisfy the CBS requirements


 * $$\left\{\bar{b}_0(\beta)=1\neq0\right\}$$


 * $$\left\{\bar{b}_{j}(\beta)=0, for j=1,2,...\right\}$$.

In the specific case where $$\left\{\beta=0\right\}$$ the CBS becomes:

Matlab code for exponential function
The following figure shows the original polynomial basis functions and the CBS basis functions for $$\left\{j=1,2,3\right\}$$,We can see that convergence occurs at an n value of 12

The following figure shows the error vs n

Fourier basis function

 * $$\mathfrak{F}_{f}=\left\{1, cos\frac{(2j-1)\pi x}{2}, sin(k\pi x) ,j,k=1,2,...\right\}$$

Likewise, the following figure shows the error vs n.

Plot


When n=9, the error is $$0.000004240676477$$.

The following figure is $$ u^h(x) $$

Exponential basis function

 * $$\mathfrak{F}_{e}=\left\{1, e^{i(x-1)}-1 ,i=1,2,...\right\}$$

When n = 12, $$ Error \approx 0.63\times 10^{-6} $$ and that is,


 * $$ Error\approx \begin{bmatrix}

-0.297066776970831\\ -0.114733443435346\\ 0.0165040777664232\\ 0.012641225491135\\ 0.000082171723276\\ -0.001243314977745\\ -0.000190989291357\\ 0.000100445865361\\ 0.000035683793041\\ -0.000005235050710\\ -0.000004723596919\\ -0.000000063082541 \end{bmatrix}$$

Plot
The following figure shows the error vs n. The following figure shows $$ u^h(x)$$

I will show 4 figures when n=1,2,3,4 to see how $$ u^h(x)$$ changes.

Matlab Code
= Problem 5.3 - LIBF =

Refer to lecture slide [[media:fe1.s11.mtg29.djvu|mtg-29]] for the problem statement.

Problem Statement
Repeat homework problem 5.1 by solving G1DM1.0/D1b, but use LIBF (Lagrange Interpolation Basis Functions) with uniform discretization (equidistant nodes, m = 4,6,8...). 1.) Explain how the LIBF are used as CBS (Constraint Breaking Solutions). 2.) Plot all LIBF used. 3.) Use MATLAB to integrate. 4.) Plot $$u_{m}^{h}$$ versus $$u$$ and $$u_{m}^{h}(.5)-u(.5)$$ versus $$m$$ 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$$

Solution
Given the above information, the differential equation can be written as: $$\frac{\partial}{\partial x}\left[(2+3x)\frac{\partial u}{\partial x}\right]+5x=0$$ The general Lagrange Interpolation basis function is given by: $$L_{i,m}(x)=\prod^{m}_{k=1}\frac{x-x_k}{x_i-x_k}$$ where $$k\neq i$$ For the case where $$ m = 4 $$ the LIBF are as follows: $$L_{1,4}(x)=\frac{(x-x_{2})(x-x_{3})(x-x_{4})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})}$$ $$L_{2,4}(x)=\frac{(x-x_{1})(x-x_{3})(x-x_{4})}{(x_{2}-x_{1})(x_{2}-x_{3})(x_{2}-x_{4})}$$ $$L_{3,4}(x)=\frac{(x-x_{1})(x-x_{2})(x-x_{4})}{(x_{3}-x_{1})(x_{3}-x_{2})(x_{3}-x_{4})}$$ $$L_{4,4}(x)=\frac{(x-x_{1})(x-x_{2})(x-x_{3})}{(x_{4}-x_{1})(x_{4}-x_{2})(x_{4}-x_{3})}$$ Since the nodes are equidistant, $$x_1=0,\;x_2=\frac{1}{3},\;x_3=\frac{2}{3},\;x_4=1$$ Substituting these values into the LIBF equations and reducing produces: $$L_{1,4}(x)=\frac{(9x^{3}-18x^{2}+11x-2)}{-2}$$ $$L_{2,4}(x)=\frac{27x^{3}}{2}-\frac{45x^{2}}{2}+9x$$ $$L_{3,4}(x)=\frac{-27x^{3}}{2}+18x^{2}-\frac{9x}{2}$$ $$L_{4,1}(x)=\frac{9x^{3}}{2}-\frac{9x^{2}}{2}+x$$ This procedure was repeated for the cases of m = 6 and m = 8 using equidistant nodes. 1.) These LIBF all satisfy the CBS requirements such that, $$\mathbf{L_{1,4}(x=1)=L_{2,4}(x=1)=L_{3,4}(x=1)=0\;\; AND\;\; L_{4,4}(x=1)=1}$$ 2.) The LIBF for m = 4 are plotted below: The LIBF for m = 6 are plotted below: The LIBF for m = 8 are plotted below: 3.) The exact solution can be directly calculated by integration and applying the boundary conditions. Recall: $$\frac{\partial}{\partial x}\left[(2+3x)\frac{\partial u}{\partial x}\right]+5x=0$$ Integrating once W.R.T. 'x' produces: $$\left[(2+3x)\frac{\partial u}{\partial x}\right]+\frac{5x^{2}}{2}+C_1=0$$ Recall the natural boundary condition: $$\Gamma _{h}=\left\{1\right\}, h = 12$$ In other words: $$n(x=1)a_{2}(x+1)\frac{\partial u}{\partial x}(x=1)=12$$ Which can be applied with the results: $$\mathbf{C_1 = -14.5}$$ Which produces: $$\left[(2+3x)\frac{\partial u}{\partial x}\right]+\frac{5x^{2}}{2}-14.5=0$$ Solving the equation for 'du', integrating both sides, and applying the essential boundary condition yields: $$\mathbf{C_2 = .90651}$$ So the exact solution (strong form) can be represented by: $$u= 4.46296ln(3x+2)-\frac{5x(3x-4))}{36}+.90651$$ m=4, the exact solution and approximate solution using LIBF are compared in the following figure: The according matlab code is: m=6, the exact solution and approximate solution using LIBF are compared in the following figure: The according matlab code is:

m=8, the exact solution and approximate solution using LIBF are compared in the following figure: The according matlab code is: The error is decreased by increasing the nodes of basis function(LIBF), which can be clearly in the following figure: The according matlab code is:

Problem Statement
Continue to learn Calculix by:

a) Extract information from the “disk” example, including node information: node numbers with coordinates and element information: element number with element nodes.

b) Generate 3 meshes for the disk with triangular element with number of elements increasing.

c) Install ccx module and write a tutorial report.

a) extract information
Firstly run the Calculix command window and type the following command to open .fbd file in “build” mode.

cd disc cgx –b disc.fbd

Now we can only see the geometry, then generate mesh by the following command,

elty all qu4 mesh all plot me all

Then select: viewing->Toggle element edges you can see the shape of your elements (Fig. a-1).



Type,

plot n all

you can see the distribution of the nodes (Fig. a-2)



To see the information for a particular node or element, you may use the “ prnt ” command, such as,

Type,

prnt n 44 prnt e 30

Then you will see (Fig. a-3),



For more information about command “ prnt ” click here.

To acquire all the node and information i.e. nodal coordinates and element component, we type in the command window,

send all abq

By this action a file called “all.msh” is created and contains all nodes number and coordinates, also all element’s composition.

Opening all.msh you’ll see the following,

The node coordinates (Fig. a-4):



The element information(Fig. a-5), since we define element type “ qu4 ” i.e. “4 node shell(displayed S4 here)”, thus every element consists of 4 nodes.



b) Generate 3 meshes for the disk with triangular element with number of elements increasing
To generate new mesh, firs we need to delete the current one:

del me all

Then we create the new mesh with triangular element.

elty all tr3 mesh all plot me all

where the “ tr3 ” stands for triangular element. Alternative choices are “ tr3u ” and “ tr6 ”.

We’ll see the mesh like this (Fig b-1):



After typing “ send all abq ” and reading it, we know the number of elements of the above mesh is 384.

To generate a denser mesh i.e. with more elements, we need to change the division of lines.

Again we delete the current mesh by “ del me all ” command and then type,

plot ld all

This action will yield the following picture (Fig. b-2) telling you the division of the lines.



To get a denser mesh, we need to divide the lines into more parts. Type “ qdiv ”(with Enter) and then “ a ” (no Enter) and then type “ r ” twice (no Enter) to create a selecting window. Use the window to capture lines and change the division of them. Note that division number of 2 digits must type space in advance. Press “ q ” to quit.

Then we have change the division like (Fig b-3),



Now we can regenerate mesh,

mesh all plot me all

Remember to toggle element edge and you’ll see the denser mesh (Fig b-4), this mesh is with 864 elements.



By the same procedure we can divide the line more thoroughly and hence a very dense mesh (Fig b-5), with 1536 elements.



The preparation
For the Windows edition of Calculix, the cgx module and the ccx module are bundled together. So there’s no need to install ccx separately.

To learn the operation of ccx we can go through this tutorial to analyze a simple case: the stress/strain analysis of a disk fixed on the edge with concentrated force acting on the center.

For some reasons – more realistic, more options of element type – it’s better to replace the current 2D disk model with a 3D one. Again here we use “ del me all ” to delete mesh.

Type this command to show the label of surfaces (Fig. c-1):

plot sa all



Now we’re going to use these 4 surfaces to sweep along vertical direction to create a disk with a finite thickness. But in advance we should define them in the same set in order to manipulate them all simultaneously.

Type this command to add a new set: “set1”, with Enter button.

qadd se1

Then press “ a ” and “ r ” twice to make a selecting window (no need to press Enter). When point the window covering all 4 surfaces, press “ s ” for "surface". Finally press “ q ” to quit.

Type:

swep se1 se2 tra 0 0.25 0

This action make the original disk swep along “+y” direction making a 3D disk with 0.25 thickness (Fig. c-2).



Then we mesh de disk by,

elty all he20 mesh all plot me all

where “he20” stands for “20nodes, brick element”, and you can select your element type here.

The mesh is shown below (Fig. c-3).



To add constraint (the disk edge fixed), we need to create a set to contain nodes lying on the edge surface. First type “ plot n all ” and type “ qadd fix ” and like before press “ a ” and then “ r ” twice to select. Press “ n ” if the window is at right position. You can utilize the drop-down menu “orientation->+y direction” to help you locate. Type command “ plus n fix g ” to check whether you have chosen correctly.

When done, the cgx should display Fig. c-4.



If you mistakenly select undesired nodes you can use “ qrem ” command to remove them. The operating procedure of “ qrem ” is like “ qadd ”.

When the constraint job is over, type

send fix abq nam

to restore the selected node number into a file – “fix.nam” - for later use.

Then we can do the same procedure by adding a concentrated force in onto the center of the plate, but no need to restore them (or actually “it” since it’s a single point) into a file. We can read from the command window that it’s node #118. We can just remember this number and announce it in the ccx input file later.

Some steps before we move on to the ccx module:


 * Type “ send all abq ” creating “all.msh” file for ccx module.
 * Type “ send all fbd ” creating a .fbd file restoring geometric entities, also serve as a backup of previous work.
 * Type “ send all frd ” creating a .frd file containing node and element information.

After the all.msh file is generated, open it with WordPad to check if the element type is correct. This is because the cgx and ccx use different acronym for the same element type. And also some element type is compatiable with cgx but not accepted by ccx. This is particularly true in 2D cases. Make sure the ccx can recognize the element type appeared in all.msh. You may check it here.

The final look of the disk (Fig. c-4.5)



Running ccx
The ccx module requires an input file to run the analysis. This kind of input file is .inp format and is also called “input deck”. The input deck is generally consists of two parts: The model definition and some steps. See Fig. c-5.



And the model definition and step definition are all composed of some “keyword” parts, these the keyword with associated data is called a “card”. To help you understand these concepts see Fig. c-6.



The input deck is often written using Sci-TE text editor on Windows OS. Open an Sci-TE and paste such code onto it. And then save it with .inp suffix. This file is identical to what displayed in Fig. c-6

*Heading MODEL: DISC *INCLUDE, INPUT=all.msh *INCLUDE, INPUT=fix.nam *NSET, NSET=LOAD 118 *MATERIAL, Name=steel *ELASTIC 28000000, 0.3 *SOLID SECTION, Elset=Eall, Material=steel *STEP *STATIC *BOUNDARY Nfix,1,2,3 *CLOAD LOAD, 2, 5 *NODE FILE U *EL FILE S, E *END STEP

The heading has no effect on analysis. Let’s go to the model definition part first.

In the model definition we have,

*INCLUDE, INPUT=all.msh *INCLUDE, INPUT=fix.nam

This tells the ccx module the mesh information. The keyword “ INCLUDE ” tells ccx to read additional file which makes this file concise. You can also paste the content of additional file directly here.

*NSET, NSET=LOAD 118

This is the node where the concentrated force acts on and we name it LOAD for later use.

*MATERIAL, Name=steel *ELASTIC 28000000, 0.3 *SOLID SECTION, Elset=Eall, Material=steel

This is the Material description. “28000000” and “0.3” means Young’s Modulus and Poisson’s Ratio respectively. Since this is only a simple tutorial, we only need these parameter by assuming the material linearly elastic.

Then ahead is the step definition. Since we have the static case, only one step is needed for analysis.

*STEP *STATIC

This declare the end of model definition and the begin of step definition, it also tells the ccx module the analysis type is “static”.

*BOUNDARY Nfix,1,2,3 *CLOAD LOAD, 2, 5

This declare the boundary conditions: the edge is fixed and there’s concentrated force. Remember the set of nodes we defined before? Its name’s “fix” just as above, the following digits “123” is the direction of constraint. Keyword “ CLOAD ” stands for concentrated force(load), and “ 2 ” stands for the “y” direction while “ 5 ” stands for the magnitude.

*NODE FILE U *EL FILE S, E *END STEP

This asks the ccx for outputting data into a .frd file. U, S, E stands for data type as displacement, stress and strain respectively.

Then click on the Sci-TE toolbar -> tool -> solve or type in the command window:

ccx disc

Wait until the “Job finished!” appear.

Post Processing
Click on the Sci-TE toolbar -> post process or type in the command window “ cgx –v disc.frd ” to view the result.

The von Mises stress distribution (Fig. c-7),



The XX component of stress distribution (Fig. c-8),



The YY component of strain distribution (Fig. c-9),



As we can we the results generally reflect the fact that the load is concentrated. But the shapes of distribution look weird. This maybe caused my the unsufficiency of nodes and elements - in this 3D analysis we only allocated around 600 nodes and 100 elements. This may also be caused by my inappropriate choice of element type.

After Tutorial
For more information on the ccx module, please refer to the ccx online documentary.

And more examples available.

d) Plot figures showing comparison between exact and numerical solution, convergence vs number of element
From HW5.1/HW5.3 we have known the exact solution is,

The comparison of $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$:






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


We can see that the convergence becomes better when number of elements increase.

Alternative Matlab code from Fangfang
Fangfang has developed an alternative way of coding and it's with nice outcome too, see part d).

d) Plot figures showing comparison between exact and numerical solution, convergence vs number of element
From HW5.2/HW5.4 we have known the exact solution is,

The comparison of $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$:






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


We can see that the convergence becomes better when number of elements increase.

Plots by Fangfang's code
Here's 3 plots about $$\displaystyle {{u}^{h}}$$ vs $$\displaystyle u$$ when $$\displaystyle nel=4,6,8$$ generated by Fangfang's matlab code. We can it's also with nice accuracy.