User:Eml5526.s11.team6/hwk6

=Problem 6.1: Solving Data Set G1DM1.0/D1b Using Quadratic Lagrangian Element Basis Functions and Uniform Discretization=

Boundary Conditions
The Essential Boundary condition as applied on the boundary $$\Gamma_g = [0]$$, is

The Natural Boundary condition as applied on the boundary $$\Gamma_h = [1]$$, is

Descritization
Here we want to use Uniform Descritization such that Element nodes are equidistant with number of elements as $$ nel = 2, 4, 6, 8, ..$$

Find
1.) For $$n_{el} = 2$$, compute $$\tilde{\underline{K}}=\sum_{e=1}^2\tilde{\underline{K}}^e$$, with $$\tilde{\underline{K}}^e$$ by [[media: Fe1.s11.mtg30.djvu| (1)p.30-5]], display ,$$\tilde{\underline{K}}^e,\;e=1,2.$$

2.) Comp. $$\underline{k}^e,\underline{L}^e$$, for $$\displaystyle e = 1,2.$$

3.) Comp. $$\tilde{\underline{K}}^e={\underline{L}^e}^T\underline{k}^e\underline{L}^e$$, for $$e = 1,2.$$ Compare to 1).

4.) Plot all QLEBF for $$n_{el}=3.$$

5.) Plot $$u_{\tilde{n}}^h$$ vs. $$\displaystyle u$$, $$[u_{\tilde{n}}^h(0.5)-u(0.5)]$$ vs. $$\tilde{n}$$

Solution
When using quadratic Lagrange element basis functions for this case we will employ 1-D quadratic elements, i.e. elements consisting of 3 nodes each. Figure 1 below shows the element structure for the case when nel=2.



Element Basis Functions and Their Derivatives
The quadratic Lagrange element basis functions are

Note that these functions take these values only within their corresponding element and are zero everywhere else.

Since the discretization used is uniform these basis functions may be simplified knowing the locations of the nodes. The simplified basis functions are

Then it follows that the derivatives of these functions are

Constructing the Global Stiffness Matrix
From equation [[media:Fe1.s11.mtg30.djvu|(1)p.30-5]]

Note: This matrix denotes the contribution due to element 1 to the global stiffness matrix. Since the basis functions $$b_4(x)$$ and $$b_5(x)$$ are zero within this element, all the terms in the fourth and fifth rows and columns are equivalently zero.

Similarly for the contribution from element 2:

Note: This matrix denotes the contribution due to element 2 to the global stiffness matrix. Since the basis functions $$b_1(x)$$ and $$b_2(x)$$ are zero within this element, all the terms in the first and second rows and columns are equivalently zero.

Now, constructing the global stiffness matrix using equation [[media:Fe1.s11.mtg30.djvu|(2)p.30-4]]

Computing the Element Stiffness Matrices
For the purpose of comparison we will now compute the global stiffness matrix by computing the element stiffness matrices and then using scatter matrices to scatter the coefficients into the proper place in the global stiffness matrix.

From equation [[media:Fe1.s11.mtg31.djvu|(1)p.31-3]]

A scatter matrix is a matrix that relates the local element degrees of freedom to the global degrees of freedom as seen in equation [[media:Fe1.s11.mtg31.djvu|(3)p.31-1]]

Due to the similarity of the structure of the element stiffness matrix and the element stiffness contribution matrix previously calculated, it can be seen that the element stiffness matrices are equivalent to the non-zero sub-matrices of the stiffness contribution matrices. i.e.

And by inspection the scatter matrices are found to be

Alternate Computation of the Global Stiffness Matrix
The element contribution to the global stiffness matrix is then computed via equation [[media:Fe1.s11.mtg31.djvu|(2)p.31-1]]

Now, as before,

Note:    $$\displaystyle (Eq. 6.1.21) = (Eq. 6.1.11)$$

QLEBF Used for nel=3
The following MATLAB code was used to generate the QLEBF used for nel=3, shown in figure 2 below



Numerical Solution
This problem was then solved using the QLEBF and increasing the number of elements, $$ nel = 2,4,6...$$ until the error converged below 10-6. The following MATLAB code was used to develop figures 3 and 4 below, showing a comparison between the numerical and exact solution and the error convergence, respectively.





=Problem 6.2: Solving Data Set G1DM1.0/D1 Using Quadratic Lagrangian Element Basis Functions and Uniform Discretization=

Boundary Conditions
The Essential Boundary condition as applied on the boundary $$\Gamma_g = [1]$$, is

The Natural Boundary condition as applied on the boundary $$\Gamma_h = [0]$$, is

Descritization
Here we want to use Uniform Descritization such that Element nodes are equidistant with number of elements as $$ nel = 2, 4, 6, 8, ..$$

Find
1.) For $$n_{el} = 2$$, compute $$\tilde{\underline{K}}=\sum_{e=1}^2\tilde{\underline{K}}^e$$, with $$\tilde{\underline{K}}^e$$ by [[media: Fe1.s11.mtg30.djvu| (1)p.30-5]], display ,$$\tilde{\underline{K}}^e,\;e=1,2.$$

2.) Comp. $$\underline{k}^e,\underline{L}^e$$, for $$\displaystyle e = 1,2.$$

3.) Comp. $$\tilde{\underline{K}}^e={\underline{L}^e}^T\underline{k}^e\underline{L}^e$$, for $$e = 1,2.$$ Compare to 1).

4.) Plot all QLEBF for $$n_{el}=3.$$

5.) Plot $$u_{\tilde{n}}^h$$ vs. $$\displaystyle u$$, $$[u_{\tilde{n}}^h(0.5)-u(0.5)]$$ vs. $$\tilde{n}$$

Solution
When using quadratic Lagrange element basis functions for this case we will employ 1-D quadratic elements, i.e. elements consisting of 3 nodes each. Figure 1 below shows the element structure for the case when nel=2.



Element Basis Functions and Their Derivatives
The quadratic Lagrange element basis functions are

Note that these functions take these values only within their corresponding element and are zero everywhere else.

Since the discretization used is uniform these basis functions may be simplified knowing the locations of the nodes. The simplified basis functions are

Then it follows that the derivatives of these functions are

Constructing the Global Stiffness Matrix
From equation [[media:Fe1.s11.mtg30.djvu|(1)p.30-5]]

Note: This matrix denotes the contribution due to element 1 to the global stiffness matrix. Since the basis functions $$b_4(x)$$ and $$b_5(x)$$ are zero within this element, all the terms in the fourth and fifth rows and columns are equivalently zero.

Similarly for the contribution from element 2:

Note: This matrix denotes the contribution due to element 2 to the global stiffness matrix. Since the basis functions $$b_1(x)$$ and $$b_2(x)$$ are zero within this element, all the terms in the first and second rows and columns are equivalently zero.

Now, constructing the global stiffness matrix using equation [[media:Fe1.s11.mtg30.djvu|(2)p.30-4]]

Computing the Element Stiffness Matrices
For the purpose of comparison we will now compute the global stiffness matrix by computing the element stiffness matrices and then using scatter matrices to scatter the coefficients into the proper place in the global stiffness matrix.

From equation [[media:Fe1.s11.mtg31.djvu|(1)p.31-3]]

A scatter matrix is a matrix that relates the local element degrees of freedom to the global degrees of freedom as seen in equation [[media:Fe1.s11.mtg31.djvu|(3)p.31-1]]

Due to the similarity of the structure of the element stiffness matrix and the element stiffness contribution matrix previously calculated, it can be seen that the element stiffness matrices are equivalent to the non-zero sub-matrices of the stiffness contribution matrices. i.e.

And by inspection the scatter matrices are found to be

Alternate Computation of the Global Stiffness Matrix
The element contribution to the global stiffness matrix is then computed via equation [[media:Fe1.s11.mtg31.djvu|(2)p.31-1]]

Now, as before,

Note:    $$\displaystyle (Eq. 6.2.21) = (Eq. 6.2.11)$$

QLEBF Used for nel=3
The following MATLAB code was used to generate the QLEBF used for nel=3, shown in figure 2 below



Numerical Solution
This problem was then solved using the QLEBF and increasing the number of elements, $$ nel = 2,4,6...$$ until the error converged below 10-6. The following MATLAB code was used to develop figures 3 and 4 below, showing a comparison between the numerical and exact solution and the error convergence, respectively. The error converged quickly, with only 8 elements being needed to reach the desired convergence.





=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

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)







$$ 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).



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.



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



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





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.



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



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.



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.



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) -



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

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



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



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 -



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



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

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



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.



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.



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.

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.





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)



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.





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.



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.



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.



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

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 -

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



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.



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.



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



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.

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.



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.

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.



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.



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.





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.



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.



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.



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-



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



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)

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



Options within the animation menu are as follows -



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.



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



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.



Finally, the following graph is promptly displayed for 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.



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.



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

=Problem 6.4 Verification of Divergence Theorem=

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

$$\displaystyle (Eq. 6.4.1-1) $$ Also the given domain is a square as shown in below Figure 6.4.1-1



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 (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 (Eq. 6.4.1-4) $$ From $$Eq. 6.4.1-3 $$ and $$ Eq. 6.4.1.-4 $$, we can say that -

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

Given
The vector field is given as follows -

$$\displaystyle (Eq. 6.4.2-1) $$ on the domain shown in figure--

Verify the divergence theorem.The curved boundary of the domain is a parabola.



To Find
To verify the Divergence theorem for the given vector field and the given domain as per Figure 6.4.2-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.2-2) $$

Solution
We have, $${{q}_{x}}= \displaystyle 3x^2y+y^3 $$  and secondly  $${{q}_{y}}= \displaystyle 3x+y^3$$ as the given vector field over the domain given by Figure 6.4.2-1 Now we can write -

Integrating over the given domain, we get LHS for $$ Eq.6.4.2-2 $$ as -

Verifying with| Wolfram Alpha

Thus :

$$\displaystyle (Eq. 6.4.2-5) $$

To find RHS:

I)The computed boundary integral on the straight line AB is

where $$ \overline n^{(1)}=-\overline j, d\Gamma=dx, y=0   on   AB $$

II)The integral is also calculated on the curved portion which is parabola with the same boundaries AB.The equation of given parabola is

$$\displaystyle \int\limits_{AB }{\overline{q }.}\overline{n}^{(2)}d\Gamma= \int\limits_{AB } \left ( \right ) \left[ \begin{array}{c} {n_x \ } \\ {n_y \ } \end{array} \right] ds $$

where $$\displaystyle d \Gamma =ds  $$



Since normal vector is on the curved surface, we consider a small portion ds of the curve and then find outward normal to that portion of the curve. From the figure, it can be seen that outward normal is given as

Substituting the value of normal vector in the equation,we get,

From the equation of parabola, $$ \displaystyle dy=-2xdx \quad and \quad y=4-x^2 $$

Verifying with| Wolfram Alpha

$$\displaystyle (Eq. 6.4.2-11) $$

Thus divergence theorem is verified.

Given/Find
Refer to lecture slide [[media:Fe1.s11.mtg35.djvu|35-2]] for problem assignment. This slide refers to Fish and Belytschko Problem 6.3 on page 149.

Using the divergence theorem prove

Solution
First, we begin with the divergence theorem shown in Equation 6.6 of Page 134 of the Fish and Belytschko text, i.e.,

For this derivation, we shall consider the simple case of two dimensions in cartiesian coordinates. As a result, the del operator and vectors can be expressed in the x (or $$ \hat{i} $$ ) and y (or $$ \hat{j} $$ ) directions. Comparing (Eq. 6.4.3.1) with the RHS of (Eq. 6.4.3.2) yields the following relation for the vector q for this particular case,

Now, substitute the new relation for the vector q into (Eq. 6.4.3.2)

The derivative of a constant is zero, and integrating over the body still yields zero, therefore,

=Problem 6.5 Solve G2DM1.0/D1 using 2D LIBF(Lagrange Interpolation Basis Function)=

Given
The G2Dm1.0/D1 PDE

The G2Dm1.0/D1 PDE can be reduced to where $$\left[ \right]$$ is an identity matrix.

The essential boundary condition is $$g = 2\ \ \ \ on\ \ \ \ \partial \Omega $$.

Find
1) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) and verify results with any FE codes with detailed documentation.

2) Plot resulting solution in 3d with contours.

Solution
1) For Eq.6.5.1, we can write the weak form

And then we choose Lagrange Interpolation Basis Function for each element with n=m=2


 * $$N_{I}(x,y)=L_{i,m}(x)\cdot L_{j,n}(y),\;\; where \;\;I=i+(j-1)m$$

After substituting Eq.6.5.4 into Eq.6.5.3 in which Lagrange Interpolation Basis Function are both trial function and weight function, we integrate over the whole domain to get stiffness matrix and achieve final results using the following Matlab codes

The resulting stiffness and force matrices for the n=m=2 case is as follows:

$$K_{EF}=\begin{bmatrix} 0.6667 & -0.1667 & 0 & -0.1667 & 0 & 0 & 0 & 0 & -0.3333\\ -0.1667 & 1.3333 & -0.1667 & -0.3333 & -0.3333 & 0 & 0 & 0 & -0.3333\\ 0 & -0.1667 & 0.6667 & 0 & -0.1667 & 0 & 0 & 0 & -0.3333\\ -0.1667 & -0.3333 & 0 & 1.3333 & 0 & -0.1667 & -0.3333 & 0 & -0.3333\\ 0 & -0.3333 & -0.1667 & 0 & 1.3333 & 0 & -0.3333 & -0.16667 & -0.3333\\ 0 & 0 & 0 & -0.1667 & 0 & 0.6667 & -0.16667 & 0 & -0.3333\\ 0& 0 & 0 & -0.3333 & -0.3333 & -0.1667 & 1.3333 & -0.16667 & -0.3333\\ 0 & 0 & 0 & 0 & -0.1667 & 0 & -0.1667 & 0.6667 & -0.3333\\ -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & -0.3333 & 2.6667 \end{bmatrix}$$

$$F=\begin{bmatrix} 0.1110x10^{-15}& -0.2220x10^{-15}& 0.1110x10^{-15}& 0x10^{-15}& -0.1110x10^{-15}& 0.1110x10^{-15}& 0.1110x10^{-15}& 0x10^{-15}& \end{bmatrix}^{T}$$

Knowing that Kd=F and solving for the d matrix yielded:


 * $$d=

\begin{bmatrix} 2& 2& 2& 2& 2& 2& 2& 2& 2& \end{bmatrix}^{T}$$

The result from the MATLAB code is as follows:

The following is the matlab code used to generate the following plots as well as calculate the above matrices.

2)



Based on the above results we see that regardless of how many elements we use, we get a uniform temperature distribution. The exact solution is easily obtained, because this is static problem with uniform boundary conditions. The 3D contour plots are flat surfaces with no contour due to the uniform temperature distribution all the way around the perimeter of the plate. If the temperature varied along one of the edges or at a node, then the plots would have some contour and an uneven temperature distribution.

=Problem 6.6 Solve G2DM1.0/D1 using 2D LLEBF(Linear Lagrange Element Basis Function)=

Given
For G2Dm1.0/D1, the PDE can be reduced to where $$\left[ \right]$$ is an identify matrix.

The essential boundary condition is $$g = 2\ \ \ \ on\ \ \ \ \partial \Omega $$.

Find
1) Solve G2DM1.0/D1 using 2D LLEBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) and verify results with any FE codes with detailed documentation.

2) Use distorted quadrilateral meshes, and compare results with uniform meshes.

Solution
1) For Eq.6.6.1, we can write the weak form

And then we choose Linear Lagrange basis function for each element

After substituting Eq.6.6.3 into Eq.6.6.2 in which Linear Lagrange basis function are both trial function and weight function, we integrate over the whole domain to get stiffness matrix and achieve final results using the following Matlab codes









From the above results we can find out that whatever how many elements we use, we can get a uniform temperature distribution finally. The reason for we can get the exact solution easily is that this static problem has quite easy uniform boundary conditions and therefore we can foresee the final results that the whole area will finally get same temperature as the environment.

The exactly same results can be obtained using codes from the book of Fish and Belytschko with some modify.

2)

We also can modify the mesh2d.m to use distorted elements and the same result in 1) can be achieved.





What is the difference between our own codes and F&B codes?
a) We use direct integrate over the whole domain instead of using Gauss integration. The basis function is given in Eq.6.6.2 and it's much easier to find out the methodology behind codes and the link between weak form (Eq.6.6.1)and the codes. Therefore our codes can illustrate how to use LLEBF and why LLEBF can be used as CBS.

b) We also use some codes from F&B codes but our codes are more general which can increase used elements just by change the n value in input_file_16ele.m. We also change the post-process that can plot 3D results with contour which is the requirement in the lecture.

Why the codes can have exact solution even with four elements?
We think it's a static problem and the boundary condition is so easy that we can even foresee the results before we calculate. It is obvious that we can have uniform temperature distribution which is identical with the environment. The essence of the codes is to solve Laplace equations with boundary conditions and the LLEBF can guarantee exact results on every nodes.

= Problem 6.7 Rules of Matrix algebra=

Solution
Hence, Therefore, Now, Hence, from eq. (6.7.6) and eq. (6.7.11)

Now, Using eq. (6.7.8), we can write Verified using | Wolfram Alpha

Hence, from eq.(6.7.14) and eq.(6.7.15) we can say that

Syntax for Wolfram Alpha:

Multiplication of two matrices:

{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}*{{1, 3, 5}, {1, -4, 1}, {2, 5, 8}}

Transpose of those multiplied matrices:

transpose[{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}*{{1, 3, 5}, {1, -4, 1}, {2, 5, 8}}]

Inverse of the above:

invert{transpose[{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}*{{1, 3, 5}, {1, -4, 1}, {2, 5, 8}}]}

Take inverse of A and then transpose it:

transpose[invert{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}]

Take transpose of A and then invert it:

invert{transpose[{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}]}

=Problem 6.8=

Solution
Nodes and weights for the 5-point Gauss–Legendre formula according to NIST Handbook are as follows:

For the table given in the problem statement, the calculated value of each term is as follows

In Fish and Belytschko's book on page 89 The position of Gauss points and corresponding weights are as follows

Hence, from the above tables we can conclude that the table given in problem statement (from [[media:fe1.s11.mtg36.djvu|Mtg 36 (a,x)]]) is reasonably accurate when compared to the above tables from NIST handbook and Fish and Belytschko's book on page 89.

2.Use Gauss quadrature to obtain exact values of the following integrals. Verify by analytical integration:
(c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB code.

Solution
a) We have $$2n_{gp}-1\ge 2 \Rightarrow \quad take \quad n_{gp}=2$$

Hence, $$(W_1,{\xi}_1) ,\quad (W_2,{\xi}_2)$$ can be calculated as follows:

Now, a=0, b=4

Hence,

Hence, the integral becomes,

Analytically, the integration can be calculated as follows:

We get the value of Gauss integration from MATLAB code as 25.33333 which is same as we got in eq. (6.8.7) by Gauss quadrature and also same as eq.(6.8.8) by Analytical integration.

b) The given integral is

We have

$$2n_{gp}-1\ge 4 \Rightarrow \quad take \quad n_{gp}=3$$

Hence, from the table of $$ {W_i,x_i}$$ in the first part of this problem we get,

Hence integral becomes,

Analytically, the integration can be calculated as follows:

We get the value of Gauss integration from MATLAB code as 1.7333 which is same as we got in eq. (6.8.11) by Gauss quadrature and also same as eq.(6.8.12) by Analytical integration.

3.Use three point Gauss Quadrature to evaluate the following integrals. Compare to the analytical itegral.
(c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB code.

Solution
We have

$$ \quad n_{gp}=3$$

Hence, from the table of $$ {W_i,x_i}$$ in the first part of this problem we get,

Hence, the integral becomes,

a)

We get the same results by integrating this function with Wolfram Alpha. Hence it is same as the analytical solution of the integral.

We get the value of Gauss integration from MATLAB code as 0 which is same as we got in eq. (6.8.16) by Gauss quadrature and also same as by Analytical integration in Wolfram Alpha.

b)

We have

$$ \quad n_{gp}=3$$

Hence, from the table of $$ {W_i,x_i}$$ in the first part of this problem we get,

Hence, the integral becomes,

We get the result by integrating this function with Wolfram Alpha as I=1 which is different than which is calculated with three point qudrature as I=1.5299 (eq.6.8.18).

We get the value of Gauss integration from MATLAB code as 1.53 which is same as we got in eq. (6.8.18) by Gauss quadrature but it is different than Analytical integration in Wolfram Alpha.

4. Effect on accuracy of integration due to change in number of points in Gauss Quadrature
The integral $$\displaystyle \int\limits_{-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

a. one-point quadrature is employed;

b. three-point quadrature is employed.

Check your calculations against MATLAB code.

Solution
a.

We have

$$ \quad n_{gp}=1$$

Hence, from the table of $$ {W_i,x_i}$$ in the first part of this problem we get,

Hence, the integral becomes

b.

We have

$$ \quad n_{gp}=3$$

Hence, from the table of $$ {W_i,x_i}$$ in the first part of this problem we get,

Hence, the integral becomes,

We get the value of Gauss integration from MATLAB code as 4 for both n=1 and n=3 which is same as we got in eq. (6.8.20)of part (a) and (eq.6.8.22) of part (b) by Gauss quadrature and also it is same as Analytical integration in Wolfram Alpha.

Hence, the accuracy of the given integral is not affected for one point Gauss Quadrature and for three point Gauss Quadrature.

=Problem 6.9: Reproduce F&B Examples 8.3, 8.4 and Solve F&B Problem 8.6=

Objectives
We were asked to use a finite element analysis software or routine to:

1.) Reproduce example 8.3, FB, p. 200

2.) Reproduce example 8.4, FB, p. 205

3.) Use the same codes to solve problem 8.6, FB, p. 212

Solution
The MATLAB finite element codes developed by F&B, corresponding to chapter 8, were used to reproduce and solve the following problems, these codes are available at | the FB companion site and are also shown below.

The basic procedure for using these codes is as follows:


 * define the problem parameters, such as material properties, boundary conditions, mesh and element structure in an appropriate input file (i.e. input_file_64ele.m)
 * specify desired plots created in this input file
 * define the problem geometry, nodal coordinates, and element connectivity matrix either in this input file or in the mesh2d.m subroutine
 * specify in the preprocessor.m file which input file to load
 * if any new variables are introduced, add them to the include_flags.m file to establish them as global variables
 * run the heat2d.m file to solve the problem and produce the plots specified in the input file

The function assembly.m assembles passed element stiffness and force matrices into the current global stiffness and force matrices. Called iteratively for each element.

BmatHeat2D.m computes the basis function derivative matrix.

gauss.m gets Gauss points and weight factors.

get_flux.m calculates the flux vectors.

heat2d.m calls the necessary functions and .m files to carry out the computations.

heat2Delem.m calculates the element stiffness and force matrices as well as the elemental source term.

include_flags.m initializes global variables

input_file_1ele.m is the input file configured to compute example 8.3 with 1 element.

input_file_16ele.m is the input file containing parameters to solve example 8.3 with 16 elements.

input_file_64ele.m is the input file containing parameters to solve example 8.3 with 64 elements.

mesh2d.m constructs the nodes and mesh for example 8.3.

NmatHeat2D.m constructs the matrix of basis functions.

plotmesh.m is used to plot the mesh/element structure.

postprocess.m displays the calculated results.

preprocessor.m specifies which input file to load.

setup_ID_LM.m sorts the nodes which are constrained on the essential boundary condition from those which are not.

solvedr.m partitions the K and F matrices and solves for the degrees of freedom.

src_and_flux.m computes the nodal flux vectors.

Example 8.3
Consider a plate with geometry shown in figure 1 below. The conductivity is isotropic with $$D = k\bigl[\begin{smallmatrix} 1&0\\0&1\end{smallmatrix}\bigr]$$ with k = 5 W°C-1. The temperature is prescribed at T=0 at the bottom and right edges. Heat fluxes $$ \bar{q} = 0$$ and $$ \bar{q} = 20 Wm^{-1}$$ on the right and top boundaries, respectively. A constant heat source of $$\displaystyle s = 6 Wm^{-2}$$ is applied over the plate. Solve the problem using 16 quadrilateral finite elements.



The accompaniment codes come already configured to solve this problem, one must only simply specify in the preprocessor.m file that the input_file_16ele.m is the input file to use (comment out the others and save the file). Running heat2d.m displays these results:





Example 8.4
Consider a plate with geometry shown in figure 4 below, with prescribed temperature at x = ± b and prescribed flux at y = ± b.



The essential boundary conditions are given by

The natural boundary conditions on $$\Gamma_q$$ are $$(\bar{q}=-k\mathbf{n}^T\nabla T)$$

Lastly the source term is defined as

The MATLAB code used for example 8.3 was adapted to solve this problem by rewriting the input file. The results are shown below.









Problem 8.6
Consider a chimney constructed of two isotropic materials: dense concrete (k = 2.0 W°C-1 ) and bricks (k = 0.9 W°C-1 ). The temperature of the hot gasses on the inside surface of the chimney is 140 °C, whereas the outside is exposed to the surrounding air, which is at T = 10 °C. The dimensions of the chimney (in meters) are shown below. For the analysis, exploit the symmetry and consider 1/8 of the chimney crosssectional area. Consider a mesh of eight elements as shown below. Determine the temperature and flux in the two materials. Analyze the problem with 2 x 2, 4 x 4 and 8 x 8 quadrilateral elements for 1/8 of the problem domain. A 2 x 2 finite element mesh is shown in Figure 8.23. Symmetry implies insulated boundary conditions on edges AD and BC. Note that element boundaries have to coincide with the interface between the concrete and bricks.



The previous MATLAB code was once again adapted to solve the given problem. Below are the results for the cases of nel = 4, 16, and 64.













=Problem 6.10-Solution to a general 2-D problem (Heat Transfer)with Lagrange Interpolation Basis Functions(LIBF) =

Given
Problem statement initiated in mtgs - [[media:Fe1.s11.mtg37.djvu|37-2]] and [[media:Fe1.s11.mtg37.djvu|38-2]]

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.



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 (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.2 below. We specify the heat flux on boundary $${\Gamma _h}$$ as $${q_n} = {q_0}$$.

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 (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 (Eq. 6.10.3) $$

Given domain with 2-dimensional model data set 2
The following data set needs to be solved using 2D LIBFs -

Solution
Referenced from Section 6.3 (pg 142 )and Section 8.1 (pg 189 ) of text book - ' A first course in finite element method' by authors Fish & Belytschco

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

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

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

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

Similarly we can write for y- direction - And now $$ Eq. 6.10.4 $$ can be modified as below - Canceling out $$dxdy $$ throughout as we have used it for elemental surface and it has to be true over the entire 2D space, $$ Eq. 6.10.8 $$ can be written with different notation in the form of Divergence of $$\overline q $$ as - 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 - Note- Usually, $$=0 $$ and $$=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 : $$= $$ for isotropic materials and $$ \ne  $$ for anisotropic materials We consider isotropic materials here hence we can assume : $$= $$. Hence $$ Eq. 6.10.10 $$ now becomes - Hence we write then :

Governing Equation for 2D- heat transfer
Using $$ Eq.6.10.9 $$ and $$ Eq. 6.9.10 $$ we can say that :

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

$$\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 - 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 -

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 - 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 :

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 (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 -

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.

Developing the semi-discrete equations
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.

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

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

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 : Here $$ \left[ B \right]$$ denotes the matrix of derivatives of the shape functions as shown in below $$ Eq. 6.10.22 $$.

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

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

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

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 -

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-

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 -

Now the first term on RHS can be discretized as -

Now the second term on RHS can be discretized as -

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

Capacitance matrix

Conductivity matrix

Heat source matrix

Solving the given general 2D data set
Referring to mtg - [[media:Fe1.s11.mtg37.djvu|41-5]], the discrete form equations for the given data set as explained in the 'Given' section above, can be written as follows -

Note- In all the matrices below, we have overlap at Node 13 of the domain as shown in Figure 6.10.5 below - Now the conduction contribution from $${{\Gamma }_{H}} $$ can be written as follows -

=Contributing Members=

Eml5526.s11.team6.deshpande 17:48, 6 April 2011 (UTC)

Eml5526.s11.team6.vork 17:52, 6 April 2011 (UTC)

Eml5526.s11.team6.kurth 17:58, 6 April 2011 (UTC)

Eml5526.s11.team6.joglekar 19:12, 6 April 2011 (UTC)

Eml5526.s11.team6.gravois 19:38, 6 April 2011 (UTC)

Eml5526.s11.team6.tupsakhare 20:06, 6 April 2011 (UTC)

eml5526.s11.team6.lee 20:09, 6 April 2011 (UTC)