User:Eml5526.s11.team6.tupsakhare/hwk3

HOMEWORK 3 -Problem 3.3

Given

Referred from - Pg 37, Reference (1)

Consider the truss shown in the following FIGURE 3.3.1(PLEASE INSERT PROBLEM STATEMENT FIG HERE)

-Nodes A and B are fixed.

-A force F = 10 N acts in the positive x- direction at node C. Coordinates of the joints are given in meters.

-Young's modulus E = 10¹¹ Pa

-Cross sectional area for all bars , A = 2. 10^-2 m²

Find

a. Number the elements and nodes

b. Assemble the global stiffness and force matrix

c. Partition the system and solve for the nodal displacements

d. Compute the stresses & reactions

Solution

Number the Element and Nodes

Node and element numbering is done in the following fashion. Refer to the FIGURE 3.3.2(PLEASE INSERT NUMBERERD FIGURE HERE) as given below-

Assemble the Global Stiffness and Force Matrix

We use the Following formula Referred from - Pg 30, Reference (1) for Elemental Stiffness Matrix of each element,

$\displaystyle {\begin{aligned}K^{e}={k^{e}}\displaystyle \left[{\begin{array}{ccccccccccccccc}(cos\phi ^{e})^{2}&cos\phi ^{e}.sin\phi ^{e}&-(cos\phi ^{e})^{2}&-cos\phi ^{e}.sin\phi ^{e}\\cos\phi ^{e}.sin\phi ^{e}&(sin\phi ^{e})^{2}&-cos\phi ^{e}.sin\phi ^{e}&-(sin\phi ^{e})^{2}\\-(cos\phi ^{e})^{2}&-cos\phi ^{e}.sin\phi ^{e}&(cos\phi ^{e})^{2}&cos\phi ^{e}.sin\phi ^{e}\\-cos\phi ^{e}.sin\phi ^{e}&-(sin\phi ^{e})^{2}&cos\phi ^{e}.sin\phi ^{e}&(sin\phi ^{e})^{2}\end{array}}\right]\end{aligned}}$

$\displaystyle (Eq.3.1)$

Where, ${k^{e}}={\frac {{A^{e}}{E^{e}}}{l^{e}}}$ and the superscript ^e denotes the element number.
Now in this case as seen from the given data, A and E for all the 4 elements have same but distinct values. So we can eliminate the superscript ^e for A and E of each element.

Element 1

• Element 1 goes from Node 1 to Node 3

• Orientation : 90 degrees to the Positive X-direction, $\phi ^{(1)}=90^{o}$ . Hence , $cos(90^{o})=0,sin(90^{o})=1$ and ${{l}^{(1)}}=1$

Using Eq. 3.3.1 for Element 1,

$\displaystyle {\begin{aligned}K^{(1)}={\frac {AE}{{l}^{(1)}}}\displaystyle \left[{\begin{array}{ccccccccccccccc}0&0&0&0\\0&1&0&-1\\0&0&0&0\\0&-1&0&1\end{array}}\right]\end{aligned}}$

Element 2

• Element 2 goes from Node 2 to Node 4

• Orientation : 90 degrees to the Positive X-direction, $\phi ^{(2)}=90^{o}$ . Hence , $cos(90^{o})=0,sin(90^{o})=1$ and ${{l}^{(2)}}=1$

Using Eq. 3.3.1 for Element 2,

$\displaystyle {\begin{aligned}K^{(2)}={\frac {AE}{{l}^{(2)}}}\displaystyle \left[{\begin{array}{ccccccccccccccc}0&0&0&0\\0&1&0&-1\\0&0&0&0\\0&-1&0&1\end{array}}\right]\end{aligned}}$

Element 3

• Element 3 goes from Node 2 to Node 3

• Orientation : 135 degrees to the Positive X-direction, $\phi ^{(3)}=135^{o}$ . Hence , $cos(135^{o})={\frac {-1}{\sqrt {2}}},sin(135^{o})={\frac {1}{\sqrt {2}}}$ and $l^{(3)}={\sqrt {2}}$

Using Eq. 3.3.1 for Element 3,

$\displaystyle {\begin{aligned}K^{(3)}={\frac {AE}{l^{(3)}}}\displaystyle \left[{\begin{array}{ccccccccccccccc}{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}\\-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}\\-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}\\{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}\end{array}}\right]\end{aligned}}$

Element 4

• Element 4 goes from Node 3 to Node 4

• Orientation : 0 degrees to the Positive X-direction,$\phi ^{(4)}=0^{o}$, $cos(0^{o})=1,sin(0^{o})=0$, $l^{(4)}=1$,

Using Eq. 3.3.1 for Element 4,

$\displaystyle {\begin{aligned}K^{(4)}={\frac {AE}{l^{(4)}}}\displaystyle \left[{\begin{array}{ccccccccccccccc}1&0&-1&0\\0&0&0&0\\-1&0&1&0\\0&0&0&0\end{array}}\right]\end{aligned}}$

GLOBAL STIFFNESS MATRIX

Using above Elemental Matrices, we construct the Global Stiffness Matrix for the given Truss System as follows,

$\displaystyle {\begin{aligned}K=AE\displaystyle \left[{\begin{array}{ccccccccccccccc}0&0&0&0&0&0&0&0\\0&1&0&0&0&-1&0&0\\0&0&{\frac {1}{2{\sqrt {2}}}}&{\frac {-1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&0&0\\0&0&-{\frac {1}{2{\sqrt {2}}}}&1+{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&0&-1\\0&0&-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&1+{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&-1&0\\0&-1&{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&1+{\frac {1}{2{\sqrt {2}}}}&0&0\\0&0&0&0&-1&0&1&0\\0&0&0&-1&0&0&0&1\end{array}}\right]\end{aligned}}$

FORCE MATRIX

The force matrix contains reactions and externally applied force components at various nodes.
Let R_x = Reaction at node in x -direction, R_y = Reaction at node in y -direction, subscripts denote the node number
Also node 3 is free to move in both x- and y- direction and is free of any external forces.
And node 4 is free to move in both x- and y- direction and has an external force of 10 N applied in the x- direction as per the given data in the problem statement.

The Force Matrix is assembled as follows -

$\displaystyle {\begin{aligned}F=\displaystyle \left[{\begin{array}{ccccccccccccccc}r_{1x}\\r_{1y}\\r_{2x}\\r_{2y}\\0\\0\\10\\0\end{array}}\right]\end{aligned}}$

NODAL DISPLACEMENT MATRIX

The noddal displacement matrix contains displacement components at various nodes.
Let u_x = Displacement at node in x -direction, u_y = Displacement at node in y -direction, subscripts denote the node number
Nodes 1 and 2 are fixed so they have all the components of displacement as zeros (0).
Also node 3 and 4 are free to move in both x- and y- direction.

The Nodal Displacement Matrix is assembled as follows -

$\displaystyle {\begin{aligned}d=\displaystyle \left[{\begin{array}{ccccccccccccccc}0\\0\\0\\0\\u_{3x}\\u_{3y}\\u_{4x}\\u_{4y}\end{array}}\right]\end{aligned}}$

Partition the system and solve for the Nodal Displacements.

Now, using the Global stiffness matrix (K), Force matrix (F) and the Nodal displacement matrix (d), we can write the Global System as follows,

$\displaystyle {\begin{aligned}AE{\begin{bmatrix}0&0&0&0&0&0&0&0\\0&1&0&0&0&-1&0&0\\0&0&{\frac {1}{2{\sqrt {2}}}}&{\frac {-1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&0&0\\0&0&-{\frac {1}{2{\sqrt {2}}}}&1+{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&0&-1\\0&0&-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&1+{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&-1&0\\0&-1&{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&1+{\frac {1}{2{\sqrt {2}}}}&0&0\\0&0&0&0&-1&0&1&0\\0&0&0&-1&0&0&0&1\end{bmatrix}}.{\begin{bmatrix}0\\0\\0\\0\\u_{3x}\\u_{3y}\\u_{4x}\\u_{4y}\end{bmatrix}}={\begin{bmatrix}r_{1x}\\r_{1y}\\r_{2x}\\r_{2y}\\0\\0\\10\\0\end{bmatrix}}\end{aligned}}$

$\displaystyle (Eq.3.3.2)$

As the displacements for first two nodes are prescribed, we partition the systems after four rows and columns.

System Partition

Also we use the following equation from - Pg 23, 'A first course in finite elements' - Fish & Belytschko -2007 and compare our global system with it to get the nodal reactions and elemental stresses using $\displaystyle (Eq.3.3.2)$ & $\displaystyle (Eq.3.3.3)$ -

$\displaystyle {\begin{aligned}\displaystyle \left[{\begin{bmatrix}K_{E}&K_{EF}\\K_{EF}^{T}&K_{F}\end{bmatrix}}{\begin{bmatrix}{\bar {d}}_{E}\\d_{F}\end{bmatrix}}={\begin{bmatrix}r_{E}\\f_{F}\end{bmatrix}}\right]\end{aligned}}$

$\displaystyle (Eq.3.3.3)$

On comparison then we get the following partitioned matrices -

$\displaystyle {\begin{aligned}\displaystyle \left[{\bar {d_{E}}}={\begin{bmatrix}0\\0\\0\\0\end{bmatrix}},\ d_{F}={\begin{bmatrix}u_{3x}\\u_{3y}\\u_{4x}\\u_{4y}\end{bmatrix}},\ r_{E}={\begin{bmatrix}r_{1x}\\r_{1y}\\r_{2x}\\r_{2y}\end{bmatrix}},\ f_{F}={\begin{bmatrix}0\\0\\10\\0\end{bmatrix}},\ K_{F}={\begin{bmatrix}1+{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&-1&0\\-{\frac {1}{2{\sqrt {2}}}}&1+{\frac {1}{2{\sqrt {2}}}}&0&0\\-1&0&1&0\\0&0&0&1\end{bmatrix}},\ K_{EF}={\begin{bmatrix}0&0&0&0\\0&-1&0&0\\-{\frac {1}{2{\sqrt {2}}}}&{\frac {1}{2{\sqrt {2}}}}&0&0\\{\frac {1}{2{\sqrt {2}}}}&-{\frac {1}{2{\sqrt {2}}}}&0&-1\end{bmatrix}}\right]\end{aligned}}$

Nodal Displacements

We solve the following system of equations to get the nodal displacement matrix,

${\textbf {d}}_{F}={{\textbf {K}}_{F}}^{-1}({\textbf {f}}_{F}-{{\textbf {K}}_{EF}}^{T}{\textbf {d}}_{E})$

$\displaystyle (Eq.3.3.3)$ -

But ${\bar {d}}_{E}$ is a Zero Matrix, which will yield, $d_{F}=K_{F}^{-1}.f_{F}$ Solving using the above equation, the unknown nodal displacements are

$\displaystyle {\begin{aligned}\displaystyle \left[{\textbf {d}}_{F}=1/AE{\begin{bmatrix}38.2885\\10.0000\\48.2885\\0\\\end{bmatrix}}=1x10^{-7}{\begin{bmatrix}0.1914\\0.05\\0.2414\\0\\\end{bmatrix}}\right]\end{aligned}}$

$\displaystyle (Eq.3.3.4)$ -

Compute the nodal reactions and Elemental stresses

Nodal reactions

${\textbf {r}}_{E}={{\textbf {K}}_{E}}{\textbf {d}}_{E}+{\textbf {K}}_{EF}{\textbf {d}}_{F}$

$\displaystyle (Eq.3.3.5)$ -

Using the above equation the unknown reaction values are

$\displaystyle {\begin{aligned}\displaystyle \left[{\begin{matrix}{\textbf {r}}_{E}&=&{\begin{bmatrix}0\\-10\\-10\\10\end{bmatrix}}\end{matrix}}\right]\end{aligned}}$

$\displaystyle (Eq.3.3.6)$

Elemental Stresses

The stresses in the elements are given by the following expression from page 34 of the reference (1)

${\begin{aligned}\sigma ^{e}={\frac {E^{e}}{l^{e}}}{\begin{bmatrix}-cos\phi ^{e}&-sin\phi ^{e}&cos\phi ^{e}&sin\phi ^{e}\end{bmatrix}}\mathbf {d} ^{e}\end{aligned}}$

(3.21)

For element 1:

$cos90^{\circ }=0,\quad sin90^{\circ }=1,$

$\mathbf {d} ^{(3)}=l^{(3)}/AE{\begin{bmatrix}u_{1x}\\u_{1y}\\u_{2x}\\u_{2y}\end{bmatrix}},\qquad \sigma ^{(3)}=1/A{\begin{bmatrix}0&-1&0&1\end{bmatrix}}{\begin{bmatrix}0\\0\\38.2843\\10\end{bmatrix}}=500N/m^{2}$

For element 2:

$cos135^{\circ }=-0.707,\quad sin135^{\circ }=0.707,$

$\mathbf {d} ^{(2)}=l^{(2)}/AE{\begin{bmatrix}u_{2x}\\u_{2y}\\u_{3x}\\u_{3y}\end{bmatrix}},\qquad \sigma ^{(2)}=1/A{\begin{bmatrix}0&-1&-0&1\end{bmatrix}}{\begin{bmatrix}0\\0\\48.2843\\0\end{bmatrix}}=0N/m^{2}$

For element 3:

$cos90^{\circ }=0,\quad sin90^{\circ }=1,$

$\mathbf {d} ^{(1)}=l^{(1)}/AE{\begin{bmatrix}u_{2x}\\u_{2y}\\u_{4x}\\u_{4y}\end{bmatrix}},\qquad \sigma ^{(1)}=1/A{\begin{bmatrix}0.5&-0.5&-0.5&0.5\end{bmatrix}}{\begin{bmatrix}0\\0\\38.2843\\10\end{bmatrix}}=-1000N/m^{2}$

For element 4:

$cos0^{\circ }=1,\quad sin0^{\circ }=0,$

$\mathbf {d} ^{(4)}=l^{(4)}/AE{\begin{bmatrix}u_{3x}\\u_{3y}\\u_{4x}\\u_{4y}\end{bmatrix}},\qquad \sigma ^{(4)}=1/A{\begin{bmatrix}-1&0&1&0\end{bmatrix}}{\begin{bmatrix}38.2843\\1\\48.2843\\0\end{bmatrix}}=10/(2*10^{-2})=500N/m^{2}$

Problem 6.3

Given

Calculix Tutorials by Beede in four parts posted on Mechanical Hacks Blog as per the links provided below :
-Calculix FEA Beam Part 1: Buliding Geometry and Meshing
-Calculix FEA Beam Part 2: Exporting Mesh, Loads, and Boundary Conditions
-Calculix FEA Beam Part 3: Writing an Input File for the CCX Solver
-Calculix FEA Beam Part 4: Post Processing Results in CGX

Find

To reproduce the steps for modeling and analysis of a beam problem from above mentioned Tutorials through CGX and CCX Finite Element software

Solution

Calculix FEA Beam Part 1: Building Geometry and Meshing

Start up data

$1.$ open the SciTE text editor, then “save as” a blank file with name beam.fbd. Firing up cgx with the current file as input is available from the Tools->Pre Process menu item or the F10 key.
Details are shown as per Figure 6.3.1(1), Figure 6.3.2(1) and Figure 6.3.3 (1)

Figure 6.3.1(1): Saving a newly opened file

Figure 6.3.2(1): Creation of batch file name : beam.fbd

Figure 6.3.3(1): Pre-processing through Tools menu

$2.$It is easier to define geometry with a batch file. Interactive mode is a bit more useful when it comes to working with the mesh like assigning boundary conditions, loads, and constraints to a set of specific nodes.

Add the following text to the beam.fbd file as shown in the below Figure 6.3.4 (1).

Figure 6.3.4(1): The beam.fbd file

The first two lines of the batch file define points in cgx model space. The first point is named p1 located at (x,y,z) = 0.
The third line of code defines a cgx line with name l1 that starts at point p1 and ends a point p2 and has a division of 25 .

Hit the F10 Key or choose Tools->Pre Process from the menu.

Figure 6.3.5(1): Pre-processing

The following result is displayed in the below Figure 6.3.6 (1).

Figure 6.3.6(1): Pre-processed line between points P1 and P2

In addition the window is divided into two distinct areas. On the right side top corner is the graphics display area bordered by a thin black line. Inside the graphics display area clicking and holding the mouse has the following functionality :

  1. Left mouse button = rotate

  2. Middle mouse button = zoom

  3. Right mouse button = pan

Moving the mouse out of the area surrounded by the thin black border the mouse has different functionality. In this area the left mouse button is used to bring up a menu which gives the user a graphical front end to some of the cgx commands which can alternatively be typed on the keyboard. Only a small subset of the commands are available through the pop up menu.

A simple 100 x 10 x 1 rectangular beam

We now proceed to generate a simple 100 x 10 x 1 rectangular beam meshed with 1x1x1 cubes.
We revise the generated batch file $beam.fbd$ with more points and lines as shown in the following Figure 6.3.7 (1) and Figure 6.3.8 (1) .

Figure 6.3.7(1): Line segments with multiple node points

Figure 6.3.8(1): Line segments with multiple node points

These extra segments will allow for meshing at the desired resolution of 1x1x1 size elements. The reason for this is that each line segment is internally limited to 99 divisions.
We also note that quadratic elements mesh at a resolution of 2 line divisions equals 1 mesh element side length. Linear elements mesh at a rate of 1 division equals 1 mesh element length.

Figure 6.3.9(1): Keyboard commands

The output window is capturing all the cgx program output and keyboard input. While the cgx window is the active window any keyboard input will be shown on the bottom line of the output window as it can be seen that two line input is provided in above figure 6.3.9 (1) -
-plot pa all tells cgx to display all points p with their names a. More specifically it displays all the points contained within the cgx set named all. Users can define their own sets with arbitrary members to manage the problem definition and processing.
-The second command plus la all is used to add the lines with their names to the existing plot. If instead plot la all had been typed after the first command the points and their names would have been eliminated from the new plot view.

Building a surface from the line using a sweep command

We now proceed to build a surface from the line using a sweep command.
We revise the generated batch file $beam.fbd$ as shown in the following Figure 6.3.10 (1) .

Figure 6.3.10(1): surface from the line

The text file is updated to create a set containing the four lines and then to sweep the set in the y direction for a length of 10 units.

The command seta is used to create a new set named lines. The l flag indicates that the items to add are lines. Several seta commands can be used to incrementally add elements, nodes, or geometry to the same set.

The swep command sweeps the set named lines via translation in (0,10,0) with 10 divisions and stores the newly created lines and points into a set labeled sweplines.

The plus sa all command shows all surfaces with their labels.

Use of plot ld all command

We use the plot ld all command to show only the lines with their divisions. The divisions in Figure 6.3.11(1) define the mesh resolution of 100 x 10 in the x and y- directions respectively.

Figure 6.3.11(1): Use of plot ld all

The following operation as shown in Figure 6.3.12-1(1) adds another swep to turn the 2D surface into a 3D body. The .fbd file and results are shown .

[[File:S.11.Team6.tupsakhare-Figure- 6.3.12-1-(1).png S.11.Team6.tupsakhare-Figure- 6.3.12-1-(1).png|thumb|centre|600px|Figure 6.3.12-1(1): Additional swep ]]

In figure 6.3.12-2 (1) the menu system toggles the background color as explained in the middle of above Figure 6.3.12-1(1) -

Figure 6.3.12-2(1): Toggling background color

Element types available to mesh with using cgx

The following element types are available to mesh with using cgx -

  1. be2: 2 node 1D linear beam element

  2. be3: 3 node 1D quadratic beam element

  3. tr3: 3 node 2D linear triangular shell element

  4. tr6: 6 node 2D quadratic triangular shell element

  5. qu4: 4 node 2D linear quadrangle element

  6. qu8: 8 node 2D quadratic quadrangle element

  7. he8: 8 node 3D linear brick element

  8. he20: 20 node 3D quadratic brick element

  9. pe6: 6 node 3D linear penta element

 10. pe15: 15 node 3D quadratic penta element

Meshing using interactive mode of CGX

We load the beam defined in the previous step in cgx and keep the mouse must be within the cgx window. Then provide with following commands

elty all he8 mesh all plot m all

The elty all he8 command indicates that we will be meshing the set all with he8 elements.

Figure 6.3.13(1): Meshing the beam using 8 node 3D linear brick element

Various views of the meshed model

Different graphical views of the model can be generated as shown in below Figure 6.3.14(1) :

Figure 6.3.14(1): View form change commands

1. Viewing -> Toggle Model Edges in the area outside of the graphics display part of cgx : as shown in Figure 6.3.14(1) below -

Figure 6.3.15(1): Toggle model edges

2. Viewing->DOTS and Viewing->Toggle Element Edges :as shown in Figure 6.3.15(1) below -

Figure 6.3.16(1): Toggle element edges

In this section , the mesh which was created based on the defined geometry and selected element type.

Calculix FEA Beam Part 2: Exporting Mesh, Loads, and Boundary Conditions

Creation of Mesh-File

Now we define into the batch file beam.fbd, a set named beam, mesh it, and save the nodes and elements in Abaqus format for input into ccx.
The following Figure 6.3.1 (2) shows the beam.fbd file after adding the new commands

Figure 6.3.1(2): Defining a set named beam

On line 1 the command seto beam has been used. This command tells cgx to open a set named beam. While beam is the open set geometry created using cgx commands will be added to this set.
Finally the set is closed on line 16 using the setc beam command. Lines 16 and 17 define the element type and build the mesh.

Now that the mesh has been created the command send beam abq is used to export the mesh created from set beam to the abq (Abaqus) format. The result of this command will reside in the same directory as the open .fbd file that created it, the working directory. Its name will be beam.msh.

Next hit the F10 key or Tools->Pre Process to process the batch file. The following window as given in Figure 6.3.2(2) appears and the mesh file is written to disk.

Figure 6.3.2(2): The mesh file written to disk.

Contents of the Mesh-file

The most interesting parts are the start of the node and element definitions respectively as shown in below figure 6.3.3(2) and Figure 6.3.4(2) -
These commands capture the shape and number ordering of the mesh.

Figure 6.3.3(2): The mesh file part-1

Line number 1 above starts with the command *NODE, NSET=Nbeam which declares the start of a node list. It also assigns a name to this set of nodes, Nbeam. Grouping together nodes and elements into logical units is a convenience of management.

Figure 6.3.4(2): The mesh file part-2

In the second image the declaration for the elements is observed on line number 2224. This command *ELEMENT, TYPE=C3D8, ELSET=Ebeam defines linear 3D bricks with 8 nodes and groups them in the Ebeam set. In cgx he8 elements turn into C3D8 elements within ccx.

Note - To define nodes and elements for a low resolution mesh the same commands shown in the images above can be written by hand in an .inp file.

To export some loads and boundary conditions

This will require selecting some nodes within interactive mode to create a set. Then the set will be exported using the send command.

With the beam.fbd file open in cgx issue the following commands to view only the nodes of the mesh. The commands are shown at the bottom of the Figure 6.3.4 (2) above.

The command plot n all is used to display the nodes, n, from the set all. The results should look like the image below. The rot -z command rotates the view to the -z direction.
In this case the set all contains everything defined in the cgx file, the same result could also be obtained by using the command plot n beam because beam is the only set meshed thus far.

Figure 6.3.5(2): Mesh with nodes displayed

Figure 6.3.6(2): Zoomed in view of the above mesh

Graphical selection

We use the command qadd fixed and press the enter button so that the mouse cursor has changed to look like Figure 6.3.7(2) below -

Variable size Selection of the graphical area can be done using command - qadd and r-button on the keyboard .
This is shown in below Figure 6.3.8(2)

Figure 6.3.8(2): Graphical area selection

In order to select all the nodes of the left hand edge the selection area will be made large enough to encompass the entire edge. In addition there are two selection modes used by the qadd command. The default selection mode adds one item each time the corresponding keyboard button is pressed. This works good for selecting one of a particular item with a small selection rectangle. If the desire is to bulk add all the items located within the selection rectangle the mode:a should be used. While the qadd command is active press the a button on the keyboard to active bulk add. The following Figure 6.3.9(2) and Figure 6.3.10(2) show the larger selection rectangle and the results of pressing the a button to active mode:a. The results of this selection will be shown in the output window as a numbered list of nodes which were added to the set. This is shown at the bottom of the figures below.

Figure 6.3.9(2): Graphical area selection

Figure 6.3.10(2): Graphical area selection

Now that the nodes have been added to the set fixed the qadd command can be quit using the q button. In order to visualize the nodes contained in the set fixed the view can be manipulated to look like the following Figure 6.3.11(2) and Figure 6.3.11a(2). Changing the color by adding a color character at the end of the plus or plot command is helpful to confirm the set contains the desired items.

Figure 6.3.11(2): Nodes visible through background color toggling

The command view edge off is equivalent to Viewing->Toggle Model Edges.
Next the plus l beam w command adds the lines from the set beam to the plot in white.
The last command issued adds the nodes contained in the new set fixed. These are plotted in green using the plus n fixed g command.

Boundary Conditions

Now that the nodes for the fixed end of the beam have been grouped together they can be exported with a boundary condition. In this case each node will be constrained such that displacements (x,y,z) = 0, that is they are fixed. This is accomplished using the send command shown in the Figure 6.3.12(2) below.

Figure 6.3.12(2): Boundary Conditions ans Constraints

The command send fixed abq spc 123 is used to write the boundary condition to a file. In this case fixed is the name of the set of nodes, abq defines the Abaqus format, spc stands for single point constraint, and 123 indicates the degrees of freedom to fix.

We notice the output below the send command. The boundary conditions are written to a file named fixed_123.bou in the working directory. Part of it is shown in figure 6.3.13(2) below. Also we notice the output statement ready which indicates processing of the file has finished.

Figure 6.3.13(2): Boundary Conditions ans Constraints

Loads

We define a distributed pressure load on an element face. Use the plot f beam command to view the element faces as shown in Figure 6.3.14(2).

Figure 6.3.14(2):Pressure load

Figure 6.3.14(2):Commands -Pressure load

Next the send command is used to write a pressure load to the load.dlo file. The send command is used to output a pressure, pres, with magnitude equal to 10 as follows -

send load abq pres 10 please wait for 'ready' write file: load.dlo ready ready

The resulting load.dlo file has the following contents within -

Figure 6.3.15(2):Pressure load

Line number 2 defines a load of magnitude 10 applied to the 6th face, P6, of element number 751. To confirm the correct element has been selected examine the model using the plot ea beam command. Also, note the Viewing->Toggle Culling Back/Front mode is used to see inside of the elements. This is necessary because the element labels are hidden behind the outside faces.

Figure 6.3.16(2):Element labels

In order to confirm the face selected in the previous step is indeed number 6 it will be necessary to turn on node number labels and zoom in on the area of element 751.

Figure 6.3.17(2):Element 751

Following Figure 6.3.18(2) shows the order of the nodes for element 751.

Figure 6.3.18(2):Order of the nodes of the Element 751

prnt e 751 elem: 751 1673 1674 1675 1676 1677 1678 1679 1680

We compare this to the element numbering guide shown in the Calculix user manual. Node numbers 1673, 1674, 1675, 1676 are located on the right hand face at x=100. These node numbers correspond to 1, 2, 3, 4 in the figure below. This figure has been copied directly from the Calculix manual. It also contains details for different element shapes. Thus it is clear to see face 6 consists of nodes 1673, 1676, 1677, 1680 and is located exactly where it should be.

The other loading method mentioned above occurs at a specific node, point loading. Loading at a point contrasts with loading a face. When a point is loaded a force vector consisting of three components for each of the x,y,z directions is applied to it. When a pressure is defined it is considered to act normal to the element face and is distributed across the surface area. If the orientation of the face or its area change the corresponding direction and magnitude of the resulting force will update as well. In the case of point loading the force will be constant in magnitude and direction as the object changes shape. Use the following commands to define a new set with the two nodes at the upper right corner of the beam where x=100, y=10 and z=0 or z=1. These nodes are a member of element 751 shown above.

seta loadnodes n 1676 1673 send loadnodes abq force 0 -5 0 Please wait for 'ready' write file:loadnodes.frc ready ready

The first command adds the nodes to a new set. This is followed by a send command to output the force in the Abaqus format with the vector components (0,-5,0). The resulting loadnodes.frc file contents are shown in the following image. Each node of the set is listed with the direction of the applied force and followed by the magnitude.

Figure 6.3.19(2):Forces based on loadnodes

Now the mesh, boundary condition, and loading has been defined and exported to files and it will be necessary to write a Calculix solver (ccx) input file (.inp) which includes these exported files to define the problem.

Calculix FEA Beam Part 3: Writing an Input File for the CCX Solver

Input file

Here we start by writing an input file and solving for the results using the finite element method.
The input file describes the finite element problem as an ordered set of commands with parameters.
We Open SciTE and create a new file named beam.inp. This file is located with the files from the previous articles. Then we add the following code from the figure 6.3.1(3) below to this file.

Figure 6.3.1(3): Meshing the beam using 8 node 3D linear brick element

The first two lines of the input file declare a *HEADING where the user can write comments and notes on the following line.
The next command, *INCLUDE, inserts the contents of the file named beam.msh into the input file.
The next command, *INCLUDE, inserts the contents of the file named beam.msh into the input file.
On lines 6 thru 9 the material model is defined.

Solid Section Definition

This assigns the elastic material EL to the element set Ebeam.

Figure 6.3.2(3): Reference for consistent set of units

Loading Step

In this step the loading is applied using the *STATIC method with only one step used.
This indicates the force will be slowly applied in a quasi-static sense. This will neglect the mass inertia of the structure.
a distributed load is applied to the tip of the beam. *NODE PRINT and *EL PRINT commands are used to print the results to a .dat text file.

Output

With the beam.inp file ,we open in SciTE choose Tools->Solve or we use the keyboard shortcut Ctrl+F10 to process this file with ccx. This will solve the finite element model as defined and output the requested files.

The output of the ccx program appears in the SciTE output window. It provides detailed diagnostics information. At the bottom it finishes the job and exits with code 0. This successful indication means a solution was found and the results have been written to disk. The ccx output should look like the following Figure 6.3.3(3) and Figure 6.3.4(3). Progress can be monitored by watching the output as ccx attempts to find a convergent solution or steps through a number of loading conditions.

Figure 6.3.3(3): SciTE Output

Figure 6.3.4(3): SciTE Output

Von-Mises stress distribution

Using the cgx menu system click on Datasets -> 2 STRESS 1.000000 , we use the menu system again to click on Datasets -> Entity -> Mises. This displays the Von Mises stress distribution for the solution at step 1 in figure 6.3.5(3) below .

Figure 6.3.5(3): Von Mises stress distribution at step 1

Calculix FEA Beam Part 4: Post Processing Results in CGX

This article will explore some of that functionality by examining the results of the cantilever beam problem from the previous article.

Open the input file beam.inp in SciTE.
Choose Tools -> Post Process from the menu system.
This launches cgx with the beam.frd results file.

Figure 6.3.1(4): SciTE Output

Note-Menu option is only available when working with an .frb file containing some simulation results e.g. Stress, Displacement etc, as seen from above figure 6.3.1(4).

plotting for various entities

The stresses and the displacement are shown in contour or graph format as below in figure 6.3.2 (4) and figure 6.3.3 (4)
An option from within the -Entity- sub-menu results in a colored plot of the selected field information.

Figure 6.3.2(4): Plotting

Plot magnification

To scale up the displacement use the keyboard command scal d 10000. This will scale the deformed shape by a factor of 10000.
This is visible in Figure 6.3.3(4) below-

Figure 6.3.3(4): Magnified displacement - deformed shape

We toggle the displacement to visualize the deformed shape. The resulting displacement from this simulation is quite small. It is not noticeably visible.

Figure 6.3.4(3): Toggle displacement

The viewing mode can be returned to the default gray beam by choosing Viewing->Show All Elements With Light as shown in below figure 6.3.4(4)

Figure 6.3.5(4): Elements With Light

Animation

Another option to visualize the shape of the deformed structure is through animation as shown in below Figure 6.3.6(4)

Figure 6.3.6(4): Animated displacement

Options within the animation menu are as follows -

Figure 6.3.7(4): Options within the animation sub-menu

To specify a cut plane

The cut plane is defined by three nodes. The nodes can be specified from the menu system using the Cut menu.

Figure 6.3.8(4): Cut plane

Result of cutting the beam top to bottom is getting the top of the beam is in tension and the bottom in compression as seen from the below figure 6.3.9(4).

Figure 6.3.9(4): Result of introducing Cut plane

Total displacement at nodes

As shown in Figure 6.3.10(4), we plot the total displacement of nodes along the top of the beam at y=10 and z=0. They span along the x-axis starting at the fixed end and progressing towards the free end. The desired field for plotting must be the currently active field selected from Datasets.

Figure 6.3.10(4): Total displacement at various nodes

Finally, the following graph is promptly displayed for displacement values at various nodes.

Figure 6.3.11(4): Displacement values at various nodes

Nodal values of entities

The nodal values are printed in yellow using the command plus nv all y. The following Figure 6.3.12(4) shows top side of the beam’s free end, near x=100, y=0, z=0.

Figure 6.3.12(4): Nodal values of stresses

Reduced set of nodes

Viewing the entire beam is difficult when several thousand nodal values are print on the screen. Limiting the nodal values to an interesting set is likely the best use of this tool.

Figure 6.3.13(4): Reduced set of nodes

This completes the 4 steps of beam analysis using CGX and CCX of calculix software

Problem 6.4-1

Given

Referring to page 2 of [Mtg 35]for problem assignment. This slide refers to Fish and Belytschko Problem 6.1 on page 148. The vector field is given as follows -

$\displaystyle {\begin{aligned}&{{q}_{x}}=-{{y}^{2}}\\&{{q}_{y}}=-2xy\\\end{aligned}}$

$\displaystyle (Eq.6.4.1-1)$

Also the given domain is a square as shown in below Figure 6.4.1-1

Figure 6.4.1-1: Domain for illustration of Divergence Theorem

Find

To verify the Divergence theorem for the given vector field and the given domain as per Figure 6.4.1-1
We need to prove that : $\int \limits _{\Omega }{{\overline {\nabla }}.}{\overline {q}}d\Omega =\oint \limits _{\Gamma }{{\overline {q}}.{\overline {n}}}d\Gamma$

$\displaystyle (Eq.6.4.1-2)$

Solution

We have, ${{q}_{x}}=-{{y}^{2}}$ and secondly ${{q}_{y}}=-2xy$ as the given vector field over the domain given by Figure 6.4.1-1
Now we can write - ${\overline {\nabla }}.{\overline {q}}={\frac {\partial {{q}_{x}}}{\partial x}}+{\frac {\partial {{q}_{y}}}{\partial y}}=0-2x=-2x$
Integrating over the given domain, we get LHS for $Eq.6.4.1-2$ as -

$\int \limits _{\Omega }{{\overline {\nabla }}.}{\overline {q}}d\Omega =$ $\int \limits _{-1}^{1}{\int \limits _{-1}^{1}{2x.dydx=}}\int \limits _{-1}^{1}{\left(2x.\left[y\right]_{-1}^{1}\right)}dx=-2\left[{{x}^{2}}\right]_{-1}^{1}=0$

Thus :

$\displaystyle \int \limits _{\Omega }{{\overline {\nabla }}.}{\overline {q}}d\Omega =0$

$\displaystyle (Eq.6.4.1-3)$

Also we have for RHS of$Eq.6.4.1-2$ as -
$\oint \limits _{\Gamma }{{\overline {q}}.{\overline {n}}}d\Gamma$ = $\int \limits _{AB}^{}{2xy.}d\Gamma +\int \limits _{BC}^{}{(-{{y}^{2}}).}d\Gamma +\int \limits _{CD}^{}{(-2xy).}d\Gamma +\int \limits _{DA}^{}{{{y}^{2}}.}d\Gamma =\int \limits _{-1}^{1}{2xy.}dx-\int \limits _{-1}^{1}{({{y}^{2}}).}dy+\int \limits _{1}^{-1}{2xy.}(-dx)+\int \limits _{1}^{-1}{{{y}^{2}}.}(-dy)$

Hence, RHS = $\oint \limits _{\Gamma }{{\overline {q}}.{\overline {n}}}d\Gamma$ = $\int \limits _{-1}^{1}{2xy.}dx-\int \limits _{-1}^{1}{({{y}^{2}}).}dy-\int \limits _{-1}^{1}{2xy.}dx+\int \limits _{-1}^{1}{{{y}^{2}}.}dy$ = $\left(2y\int \limits _{-1}^{1}{x.}dx-\int \limits _{-1}^{1}{({{y}^{2}}).}dy-2y\int \limits _{-1}^{1}{x.}dx+\int \limits _{-1}^{1}{{{y}^{2}}.}dy\right)=0$

Thus :

$\displaystyle \oint \limits _{\Gamma }{{\overline {q}}.{\overline {n}}}d\Gamma =0$

$\displaystyle (Eq.6.4.1-4)$

From $Eq.6.4.1-3$ and $Eq.6.4.1.-4$, we can say that -

$\displaystyle \int \limits _{\Omega }{{\overline {\nabla }}.}{\overline {q}}d\Omega =\oint \limits _{\Gamma }{{\overline {q}}.{\overline {n}}}d\Gamma =0$

This proves that the given vector field as per $Eq.6.4.1-1$in the stated domain of Figure 6.4.1 satisfies the Divergence Theorem

Problem 6.10

Given

2D space with boundary convection

Consider a heat conduction problem in 2-dimensional space with boundary convection as shown in the below Figure 6.10.1.

Figure 6.10.1 :Problem Domain and Boundary Conditions for heat transfer with boundary heat convection

Boundary Conditions

Above 2D space has temperature $T(x,y)$ as space unknown to solve for.
We can state it as a boundary value problem with following boundary conditions specified at 3 different kind of boundaries as : ${\Gamma _{g}}$, ${\Gamma _{h}}$ and ${\Gamma _{H}}$

Essential BC on ${\Gamma _{g}}$

$\displaystyle T(x,y)={T_{0}}$

$\displaystyle (Eq.6.10.1)$

Note : $T(x,y)$ can also can be defined as a varying quantity along ${\Gamma _{g}}$

Natural BC on ${\Gamma _{h}}$

Refer to the Figure 6.10.3 below. We specify the heat flux on boundary ${\Gamma _{h}}$ as ${q_{n}}={q_{0}}$.

Figure 6.10.3 :Figure- 6.10.3- heat flux natural BC on boundary gamma-h

This is as we have : ${q_{n}}={\overline {q}}.{\overline {n}}$ with ${\overline {n}}$ as the unit normal vector on ${\Gamma _{h}}$
Hence we can also write : ${q_{n}}=(-k.{\overline {\nabla }}T).{\overline {n}}=-k({\overline {\nabla }}T.{\overline {n}})$
OR we then have :

$\displaystyle {q_{0}}=k({\overline {\nabla }}T.{\overline {n}})$

$\displaystyle (Eq.6.10.2)$

Note- $\nabla T$ is positive as explained in the - Assumption- subsection below

We are to consider only the normal component of ${\overline {q}}$, i.e. ${q_{n}}$, for heat transactions as that is the only component causing heat flow in or out of the element. The tangential component of ${\overline {q}}$, i.e. ${q_{t}}$, is just circulating over the boundary.

Assumption

${\overline {n}}$ is unit vector pointing in the outward direction. This implies that $\nabla T$ is positive and temperature increases as we go out away from the surface towards inside of the body.
Also we note that ${\overline {n}}$ is a vector quantity normal to the $T(x,y)$ contours in 2D space.
As $T(x,y)$ increases from inside to outside, heat flux ${q_{n}}$ goes into the volume from outside of the boundary ${\Gamma _{h}}$
This implies that when ${q_{0}}$ is positive, heat comes into the system

Convection BC on ${\Gamma _{H}}$

We have as the convection boundary where heat transfer occurs via heat convection
Hence we can say that the normal component ${q_{n}}$ of the heat flux ${\overline {q}}$ equals the heat transfer due to heat convection at ${\Gamma _{H}}$. Hence we write-

${q_{n}}=h(T-Ta)$ where , we have -
$T$ : Temperature at the surface of ${\Gamma _{H}}$
${T_{a}}={T_{\infty }}=$ Say, Ambient temperature
$h$ :Coefficient of heat convection at ${\Gamma _{H}}$
Thus :

$\displaystyle -k({\overline {\nabla }}T.{\overline {n}})=h(T-Ta)$

$\displaystyle (Eq.6.10.3)$

Solution

Derivation of the strong form- 2D heat transfer

Deriving the Equilibrium Equation for 2D- heat transfer

We first derive the equilibrium equation using elemental plane surface as shown in Figure- 6.10.3 below. The elemental surface has length as 'dx' and height as 'dy'. The inflow and outflow of heat is shown symbolically with arrows at various sections.

Note-
Heat flux ' ${\underline {q}}$ ' is a vector quantity in real 2D space defined per unit surface area and also
$f=$ Heat generation per unit volume for the given elemental surface
We assume unit dimension along the direction perpendicular to the screen or paper

Figure- 6.10.3- 2D Heat Transfer Equilibrium Element

Now the equilibrium of heat can be written in following manner -

$({({q_{x}})_{x}}-{({q_{x}})_{x+dx}})(dy.1)+({({q_{y}})_{y}}-{({q_{y}})_{y+dy}})(dx.1)+f(dx.dy.1)=\rho c{\frac {\partial T}{\partial t}}$

$\displaystyle (Eq.6.10.4)$

Now we use Taylor Series Expansion technique and can then write -

${({q_{x}})_{x+dx}}={({q_{x}})_{x}}+{\left({\frac {\partial {q_{x}}}{\partial x}}\right)_{x}}dx+higher\_order\_terms$
OR we can say with the higher order terms neglected :

$({({q_{x}})_{x}}-{({q_{x}})_{x+dx}})=-{\left({\frac {\partial {q_{x}}}{\partial x}}\right)_{x}}dx=-\left({\frac {\partial {q_{x}}}{\partial x}}\right)dx$

$\displaystyle (Eq.6.10.5)$

Similarly we can write for y- direction -

$({({q_{y}})_{y}}-{({q_{y}})_{y+dy}})=-{\left({\frac {\partial {q_{y}}}{\partial y}}\right)_{y}}dy=-\left({\frac {\partial {q_{y}}}{\partial y}}\right)dy$

$\displaystyle (Eq.6.10.6)$

And now $Eq.6.10.4$ can be modified as below -

$-\left({\frac {\partial {q_{x}}}{\partial x}}\right)dxdy-\left({\frac {\partial {q_{y}}}{\partial y}}\right)dxdy+f.dxdy=\rho c{\frac {\partial T}{\partial t}}$

$\displaystyle (Eq.6.10.7)$

Canceling out $dxdy$ throughout as we have used it for elemental surface and it has to be true over the entire 2D space,

$\left({\frac {\partial {q_{x}}}{\partial x}}\right)+\left({\frac {\partial {q_{y}}}{\partial y}}\right)=f-\rho c{\frac {\partial T}{\partial t}}$

$\displaystyle (Eq.6.10.8)$

$Eq.6.10.8$ can be written with different notation in the form of Divergence of ${\overline {q}}$ as -

${\overline {\nabla }}.{\overline {q}}=f-\rho c.\nabla T$

$\displaystyle (Eq.6.10.9)$

Note that $Eq.6.10.9$ is valid for 2D as well as 3D heat transfer with ${\overline {q}}$ as 2D or 3D vector respectively in real 2D or 3D space.
Here we are solving for temperature $T(x,y)$ so we use 2D space with $T(x,y)$

Constitutive Equation for 2D- heat transfer

Here we use the Fourier Law for 2D space which is stated as  : ${\overline {q}}=-\left[k\right]{\overline {\nabla }}T$
It can be written as -

$\left({\begin{array}{ccccccccccccccc}{q_{x}}&{}\\{q_{y}}&{}\end{array}}\right)=\left({\begin{array}{ccccccccccccccc}{k_{xx}}&{k_{xy}}\\{k_{yx}}&{k_{yy}}\end{array}}\right)\left({\begin{array}{ccccccccccccccc}{\frac {\partial T}{\partial x}}&{}\\{\frac {\partial T}{\partial y}}&{}\end{array}}\right)$

$\displaystyle (Eq.6.10.10)$

Note- Usually, ${k_{xy}}=0$ and ${k_{yx}}=0$ for most of the materials as we do not see flow of heat perpendicular to the direction of the temperature gradient line.
Also , we have  : ${k_{xx}}={k_{yy}}$ for isotropic materials and ${k_{xx}}\neq {k_{yy}}$ for anisotropic materials
We consider isotropic materials here hence we can assume : ${k_{xx}}={k_{yy}}$ . Hence $Eq.6.10.10$ now becomes -

$\left({\begin{array}{ccccccccccccccc}{q_{x}}&{}\\{q_{y}}&{}\end{array}}\right)=k\left({\begin{array}{ccccccccccccccc}1&0\\0&1\end{array}}\right)\left({\begin{array}{ccccccccccccccc}{\frac {\partial T}{\partial x}}&{}\\{\frac {\partial T}{\partial y}}&{}\end{array}}\right)$

Hence we write then :

${\overline {q}}=-k{\overline {\nabla }}T$

$\displaystyle (Eq.6.10.11)$

Governing Equation for 2D- heat transfer

Using $Eq.6.10.9$ and $Eq.6.9.10$ we can say that :

${\overline {\nabla }}.(-k{\overline {\nabla }}T)=f-\rho c.\nabla T$ OR

$-k({\overline {\nabla }}.{\overline {\nabla }}T)-f=-\rho c.\nabla T$

Thus the strong form for 2D - Heat Transfer can be written as follows-

$\displaystyle {\frac {\partial }{\partial {x}_{i}}}\left({{K}_{ij}}{\frac {\partial T}{\partial {x}_{j}}}\right)+f\left(x,t\right)=\rho c{\frac {\partial T}{\partial t}}$ for $\forall (x,y)\in V(Volume\_Field)$

$\displaystyle (Eq.6.10.11)$

Deriving the Weak Form for 2D- heat transfer with boundary convection

Let us select the weighting function $w(x,y)$ such that $w(x,y)=0$ at boundary ${{\Gamma }_{g}}$ where the essential boundary condition is applied. This is actually obvious to say as the essential boundary condition is specified at ${{\Gamma }_{g}}$, the value of $T(x,y)$ will be constant there and will no longer require any weighting function multiplication to approximate it in the given domain. Also this will eliminate the unknown flux at the boundary ${{\Gamma }_{g}}$

Integrating over volume (for general 3D - case) after multiplying with the weighting function w(x) , in the form of weighted residual, we get -

$\int {w\left({\frac {\partial }{\partial {{x}_{i}}}}\left({{K}_{ij}}{\frac {\partial T}{\partial {{x}_{j}}}}\right)+f\left(x,t\right)-\rho c{\frac {\partial T}{\partial t}}\right)}dv{=0}$

$\displaystyle (Eq.6.10.12)$

Note- For 2D- case : Put $dz$ = $1$ which will give us the Area Integral
Now we use Green's Theorem to perform the integration by parts.
Our aim to get rid of the derivative on ${\overline {\nabla }}.({{K}_{ij}}{\overline {\nabla }}T)$ and a derivative should appear on $w(x)$
We state the general Green's Theorem for scalar fields $\varphi$ , T - is as follows -

$\int \limits _{S}{\varphi }{{\nabla }^{2}}.Tds=\int \limits _{V}{}\varphi ({\overline {\nabla }}T.{\overline {n}})ds-\int \limits _{V}{}\nabla \varphi .\nabla Tdv$

$\displaystyle (Eq.6.10.13)$

Here, ${{\nabla }^{2}}$ serves as Laplacian or Divergence of the gradient.
Also , ${\overline {n}}$ is the unit normal at boundary and is variable at each point on the boundary.
Thus ,we have a volume integral converted to a volume plus a surface integral
Hence, $Eq.6.10.12$ can now be modified as follows -

$\int \limits _{V}{\left[\left({\overline {\nabla }}.({{K}_{ij}}{\overline {\nabla }}T)\right)w(x)\right]}dv+\int \limits _{V}{f.w(x)dv}=\int \limits _{V}{w(x).\rho c{\frac {\partial T}{\partial t}}}dv$

$\int \limits _{S}{({{K}_{ij}}.{\overline {\nabla }}T.{\overline {n}}).w(x)ds}-\int \limits _{V}{\left[\left({\overline {\nabla }}.({{K}_{ij}}{\overline {\nabla }}T)\right)w(x)\right]}dv+\int \limits _{V}{f.w(x)dv}=\int \limits _{V}{w(x).\rho c{\frac {\partial T}{\partial t}}}dv$

$\displaystyle (Eq.6.10.14)$

This is general equation for 3D as well as 2D case and can be converted to specific case depending upon whether the gradient is 2D or 3D
Now again we note that in the weak form the natural boundary conditions get pulled in automatically, for example- the surface integral in above equation will be split into integral over the 3 boundaries - ${{\Gamma }_{g}}$ , ${{\Gamma }_{H}}$ , and${{\Gamma }_{h}}$ , as follows :

$\int \limits _{S}{({{K}_{ij}}.{\overline {\nabla }}T.{\overline {n}}).w(x)ds}=\int \limits _{{\Gamma }_{g}}{({{K}_{ij}}.{\overline {\nabla }}T.{\overline {n}}).w(x)ds}+\int \limits _{{\Gamma }_{h}}{{{q}_{0}}.w(x)ds}+\int \limits _{{\Gamma }_{H}}{h({{T}_{a}}-T).w(x)ds}$

$\displaystyle (Eq.6.10.15)$

Note- In the above Equation, the first term on the RHS will vanish as the weighting function $w(x)$ vanishes on ${{\Gamma }_{g}}$ as explained. Also here we flipped the signs for convection boundary condition on boundary on ${{\Gamma }_{g}}$
Now, we put $Eq.6.10.15$ into $eq.6.10.14$ and keep the terms required to get the stiffness matrix $\left[K\right]$ along with the terms involving temperature on LHS and all the terms with known values on the RHS.

$\displaystyle \int \limits _{V}{\underbrace {w(x).\rho c{\frac {\partial T}{\partial t}}dv} _{{\tilde {m}}(w,T)}+\underbrace {\left[\int \limits _{V}{\overline {\nabla }}T.{{K}_{ij}}{\overline {\nabla }}w(x)dv+\int \limits _{{\Gamma }_{H}}{hTw(x)ds}\right]} _{{\tilde {k}}(w,T)}}=\underbrace {\int \limits _{{\Gamma }_{h}}{{q}_{0}}w(x)ds+\int \limits _{{\Gamma }_{H}}{h{{T}_{a}}}w(x)ds+\int \limits _{V}{f}w(x)dv} _{{\tilde {f}}(w)}$

$\displaystyle (Eq.6.10.16)$

$Eq.6.10.16$ is the weak form for the given 2D Heat Transfer domain with convection boundary specified as ${{\Gamma }_{H}}$

This equation is of the following form -

$\displaystyle {\tilde {m}}(w,u)+{\tilde {k}}(w,u)={\tilde {f}}(w)\quad \forall w$ such that $\displaystyle w=0\quad on\quad {{\Gamma }_{g}}$

It is the continuous form of the weak equation and we need to convert it into a discrete form using the method of Desensitization in FEM .
Also note that -
1. In $Eq.6.10.16$ , temperature is a continuous function of field parameters (x,y).
2. This equation is good for 2D as well as 3D domains.
3. We can handle anisotropic materials by converting conductivity 'k' into matrix form as explained .

Descritization of the domain

System of equations in descritized form

Now the temperature distribution in 2D space can be written as follows with subscript 'e' denoting to an element i.e. temperature at a point (x,y) in the 2D space is a sum over that of the elemental nodes. Also here $N_{i}$ denotes shape function and $T_{i}$ denotes nodal temperature values.

$T(x,y)=\sum \limits _{e}{{{N}_{i}}{{T}_{i}}}$ $\displaystyle (Eq.6.10.17)$

This can be written in matrix form as follows where the nodal values matrix is shown by $\left\{{{X}_{e}}\right\}$

$\left\{T\right\}={{\left\{N\right\}}^{T}}\left\{{{X}_{e}}\right\}=\left\{{{N}_{1}}\right.{{N}_{2}}{{N}_{3}}\left.{{N}_{4}}\right\}\left({\begin{matrix}{{T}_{1}}&{}\\{\begin{aligned}&{{T}_{2}}\\&{{T}_{3}}\\&{{T}_{4}}\\\end{aligned}}&{}\\\end{matrix}}\right)$ $\displaystyle (Eq.6.10.18)$

Also applying the same for the weighting function, we can write :
This can be written in matrix form as follows -

$\left\{w\right\}={{\left\{N\right\}}^{T}}\left\{\partial {{X}_{e}}\right\}=\left\{{{N}_{1}}\right.{{N}_{2}}{{N}_{3}}\left.{{N}_{4}}\right\}\left({\begin{matrix}\partial {{T}_{1}}&{}\\{\begin{aligned}&\partial {{T}_{2}}\\&\partial {{T}_{3}}\\&\partial {{T}_{4}}\\\end{aligned}}&{}\\\end{matrix}}\right)$ $\displaystyle (Eq.6.10.19)$

We note that we have got gradient terms in the weak form as seen in $Eq.6.10.16$. The ${\overline {\nabla }}T$ can be written as :

${\overline {\nabla }}T=\left({\begin{matrix}{\frac {\partial T}{\partial x}}&{}\\{\frac {\partial T}{\partial y}}&{}\\\end{matrix}}\right)=\left[B\right]\left\{{{X}_{e}}\right\}$ $\displaystyle (Eq.6.10.20)$

Here $\left[B\right]$ denotes the matrix of derivatives of the shape functions as shown in below $Eq.6.10.22$ .

${\overline {\nabla }}T=\left({\begin{matrix}{\frac {\partial {{N}_{1}}}{\partial x}}&{\frac {\partial {{N}_{2}}}{\partial x}}\\{\frac {\partial {{N}_{1}}}{\partial y}}&{\frac {\partial {{N}_{2}}}{\partial y}}\\\end{matrix}}{\begin{matrix}{\frac {\partial {{N}_{3}}}{\partial x}}&{\frac {\partial {{N}_{4}}}{\partial x}}\\{\frac {\partial {{N}_{3}}}{\partial y}}&{\frac {\partial {{N}_{4}}}{\partial y}}\\\end{matrix}}\right)\left\{{\begin{aligned}&{{T}_{1}}\\&{{T}_{2}}\\&{{T}_{3}}\\&{{T}_{4}}\\\end{aligned}}\right\}$ $\displaystyle (Eq.6.10.21)$

Similar equation can be written for the weighting function as follows -

${{\left[{\overline {\nabla }}w\right]}^{T}}={{\left[B\right]}^{T}}{{\left\{\partial {{X}_{e}}\right\}}^{T}}$ $\displaystyle (Eq.6.10.22)$

Now the LHS of the weak form in $Eq.6.10.16$ can be written as for a single element -

$\int \limits _{V}{\left[{\overline {\nabla }}T.k{\overline {\nabla }}w(x)\right]}dv=\sum \limits _{o}^{{n}_{e}}{\int \limits _{V}{{{\left\{\partial {{X}_{e}}\right\}}^{T}}{{\left[B\right]}^{T}}\left[K\right]}\left[B\right]\left\{{{X}_{e}}\right\}dv}$ $\displaystyle (Eq.6.10.23)$

Hence, LHS becomes for $Eq.6.10.16$ as :

$\int \limits _{V}{\left[{\overline {\nabla }}T.k{\overline {\nabla }}w(x)\right]}dv={{\left\{\partial {{X}_{e}}\right\}}^{T}}\left[{{K}_{e1}}\right]\left\{{{X}_{e}}\right\}$ $\displaystyle (Eq.6.10.24)$

$\left[{{K}_{e1}}\right]=\int \limits _{V}{{{\left[B\right]}^{T}}\left[K\right]}\left[B\right]dv$ $\displaystyle (Eq.6.10.25)$

Depending upon various boundary conditions the RHS of $Eq.6.10.16$ can be written as $\left[{{F}_{e}}\right]{{\left\{\partial {{X}_{e}}\right\}}^{T}}$
Hence the discrete form Equation for first term of the stiffness matrix becomes -

$\displaystyle {{\left\{\partial {{X}_{e}}\right\}}^{T}}\left[{{K}_{e1}}\right]\left\{{{X}_{e}}\right\}=\left[{{F}_{e}}\right]{{\left\{\partial {{X}_{e}}\right\}}^{T}}$ $\displaystyle (Eq.6.10.26)$

where, $\left[{{K}_{e1}}\right]=\int \limits _{V}{{{\left[B\right]}^{T}}\left[K\right]}\left[B\right]dv$

The second term of the stiffness can be written in discrete form as follows-

$\displaystyle \int \limits _{{\Gamma }_{H}}{hTw(x)ds}=\int \limits _{}{h\left\{\delta {{X}_{e}}\right\}}\underbrace {{\left\{N\right\}}^{T}} _{2\times 1}\left\{N\right\}\left\{{{X}_{e}}\right\}ds$ $\displaystyle (Eq.6.10.27)$

Now the mass matrix can be written in discretized form from the first term of the continuous weak form of $Eq.6.10.16$ as follows -

$\displaystyle \int \limits _{V}{\underbrace {w(x).\rho c{\frac {\partial T}{\partial t}}dv} _{{\tilde {m}}(w,T)}}=\int \limits _{V}{\rho c}\underbrace {{\left\{N\right\}}^{T}} _{2\times 1}\left\{N\right\}{\frac {\partial T}{\partial t}}dv$ $\displaystyle (Eq.6.10.28)$

Now the first term on RHS can be discretized as -

$\displaystyle \int \limits _{{\Gamma }_{h}}{{q}_{0}}w(x)ds={{\left\{\delta {{X}_{e}}\right\}}^{T}}\left[\int \limits _{{\Gamma }_{h}}{{{q}_{0}}\underbrace {\left\{N\right\}} _{2\times 1}ds}\right]$ $\displaystyle (Eq.6.10.29)$

Now the second term on RHS can be discretized as -

$\displaystyle \int \limits _{{\Gamma }_{H}}{h{{T}_{a}}}w(x)ds=\int \limits _{{\Gamma }_{H}}{h{{T}_{a}}}{{\left\{\delta {{X}_{e}}\right\}}^{T}}\left\{N\right\}{\frac {Element\_length}{2}}ds$ $\displaystyle (Eq.6.10.30)$

And the last term on RHS of the weak form can be discretized as -

$\displaystyle \int \limits _{V}{f}w(x)dv=f{{\left\{\delta {{X}_{e}}\right\}}^{T}}\left\{N\right\}\det[J]dv$ $\displaystyle (Eq.6.10.31)$

Capacitance matrix

${\tilde {m}}(w,u)=\int \limits _{\Omega }{\rho c{{[N]}^{T}}[N]{\frac {\partial u}{\partial t}}d\Omega }$

Conductivity matrix

${\tilde {k}}(w,u)=\int \limits _{\Omega }{{{[B]}^{T}}[B]d\Omega }+\int \limits _{{\Gamma }_{h}}{H{{[N]}^{T}}[N]d{{\Gamma }_{h}}}$

Heat source matrix

${\tilde {f}}(w)=-\int \limits _{{\Gamma }_{h}}{h{{[N]}^{T}}d{{\Gamma }_{h}}}+\int \limits _{{\Gamma }_{h}}{H{{[N]}^{T}}{{u}_{\infty }}d{{\Gamma }_{h}}}+\int \limits _{\Omega }{{{[N]}^{T}}f(x,t)d\Omega }$

Contour and 3 D plots

$\displaystyle u_{p}^{h}(x)={{d}_{1}}*N_{1}+{{d}_{2}}N_{2}+{{d}_{3}}N_{3}+\cdots +{{d}_{9}}N_{9},_{}^{}j=0,1,2,\cdots ,9$ This is the trial solution with 9 node LIBF element
We then derive the trial solution using LIBF obtained is as follows-
${\begin{aligned}&U_{F}^{h}=(x)0.52*x*y*(x-1)*(y-1)-1*x*(1.0*x+1.0)*(1.0*y+1.0)*(y-1)\\&\quad +0.68*y*(1.0*x+1.0)*(1.0*y+1.0)*(x-1)-0.17*x*(1.0*y+1.0)*(x-1)*(y-1)\\&\quad -1.000*y*(1.0*x+1.0)*(x-1)*(y-1)+0.85*(1.0*x+1.0)*(1.0*y+1.0)*(x-1)*(y-1)\\&\quad +0.55*x*y*(1.0*x+1.0)*(1.0*y+1.0)+0.55*x*y*(1.0*x+1.0)*(y-1)-0.37*x*y*(1.0*y+1.0)*(x-1)\\\end{aligned}}$

Matlab Code

% Heat 2D Diffusion in File heat2d.m clear all; L = 1.0; W = L; Tend = 80; maxk = 300; dt = Tend/maxk; n = 20; u(1:n+1,1:n+1,1:maxk+1) = 70; dx = L/n; dy = W/n; h = dx; b = dt/(h*h); cond = 0.0; spheat = 0.0; rho = 1; a = cond/(spheat*rho); alpha = a*b; for i = 1:n+1 x(i) =(i-1)*h; y(i) =(i-1)*h; end for k=1:maxk+1 time(k) = (k-1)*dt; for j=1:n+1 u(1,j,k) =300*(k<120)+ 70; end end for k=1:maxk for j = 2:n for i = 2:n u(i,j,k+1) =0.*dt/(spheat*rho)+(1-4*alpha)*u(i,j,k) + alpha*(u(i-1,j,k)+u(i+1,j,k)+u(i,j-1,k)+u(i,j+1,k)); end end end mesh(x,y,u(:,:,maxk)') % Generates Sequence of 3D plots heat2d; lim =[0 1 0 1 0 400]; for k=1:5:200 mesh(x,y,u(:,:,k)') title ('heat versus space at different times' ) axis(lim); k = waitforbuttonpress; end

The 3d and contour plots are as follows-

3D plot 1

3d plot 2

contour plot 1

Contributing Members

Team- Contribution Table
Problem	Mtg	Assigned To	Solved By	Typed By	Proofread By
6.1	31-4	Shawn, Chris	Shawn	Chris	Chris, Nakul
6.2	31-4	Nakul, William	Nakul	Nakul	William, Shreyas
6.3	33-1	Anand	Anand	Anand	Shreyas, Shawn
6.4	35-2	Anand,Shreyas,Chris	All of 3	All of 3	Shawn, Li
6.5	35-4	Chris,Shawn	Chris	Shawn	William, Li
6.6	35-4	Shreyas, Li	Both	Both	Nakul, Anand
6.7	36-1	William,Nakul	William	William	Anand, Li
6.8	36-5	William,Nakul	Nakul	Nakul	Shreyas, Shawn
6.9	37-1	All except Anand	All	All	Anand
6.10	37-2	Anand	Anand	Anand	Shreyas, Shawn