# 6.3 Beam Tutorial for Calculix by Ty Beede

## Given: Tutorial by Beede on mechanicalhacks blog

A step by step tutorial explaining how to analyze a beam using Calculix can be found at [1]

## Find: Reproduce Beam Tutorial

Reproduce all steps in the tutorial by Ty Beede

## Solution

Calculix is an open source Finite Element code (with ABAQUS like input) available at http://www.dhondt.de/

More details on installation can be found at Calculix Installation Procedure

Analysis using Calculix is divided into 4 steps:

1) Buliding Geometry and Meshing
2) Exporting Mesh, Loads, and Boundary Conditions
3) Writing an Input File for the CCX Solver
4) Post Processing Results in CGX

In this tutorial we analyze a cantilever beam subjected to pressure load at one end and fixed at the other end.

### Step 1: Buliding Geometry and Meshing

Create a new file 'Beam.fbd' in the SciTE editor. To construct the 3D geometry of the beam we do the following:

(Note: code used to generate the geometry and commands required to display the geometry at each step are embedded at the bottom of each screen shot)

a) Define points and join them to create lines

b) Sweep the lines to create a surface

c) Sweep the surface to create a volume

d) Assign Element Types and mesh the geometry

### Step 2: Exporting Mesh, Loads, and Boundary Conditions

In the earlier section we created the geometry and generated a mesh.

The 'send beam abq' command used at the end of the 'beam.fbd' file, generates the node & element information upon preprocessing. This data is saved in 'beam.msh' file in the same directory.

The mesh data can be found here:beam.msh

After the mesh has been created we need to Apply the Boundary conditions and Loads:

#### Applying the Boundary Conditions

-> Use the command 'plot n all' to display all the nodes

-> Zoom in to the left end of the beam where Boundary condition is to be applied.

-> Type 'qadd fixed' and hit enter. The qadd command provides a graphical selection method that can be used interactively with the mouse. The qadd command will stay active until the 'q' button is pressed on the keyboard.

At the tip of the mouse pointer is a little rectangle outlined in black. With the 'qadd' command active the mouse is used to select items like nodes, elements, faces, points, lines, and etc. The items desired for selection must fall within the bounds of the little black rectangle shown above. When the rectangle is hovering over the desired object press the keyboard button corresponding to the item of interest in order to add it to the set. In this case the 'n' button will be used to select nodes within the bounds of the rectangle.

For adding items in bulk it will be necessary to change the size of the selection rectangle.The 'r' key on the keyboard is used to define two points of a rectangle which will become the new selection area for qadd. Simply press the 'r' button and the current cursor position becomes the new location for a corner of the selection area. Move the mouse in the horizontal and vertical directions such that it is offset from the previous point. Press the 'r' button again and a new selection rectangle has been defined.

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.

-> after mode:a has been selected hit 'n' to select all the nodes inside the rectangle to apply the fixed boundary condition.

The output window looks like this upon selecting the nodes to be fixed:

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, Change the color of the nodes and toggle the background color

Background can be toggled to black using the menu system command 'Viewing->Toggle Background Color'

Color of the node list can be changed by using the following command 'plus n fixed g'

The visualization screens appears as follows:

-> The boundary condition is saved to a separate file(.bou) using the command 'send fixed abq spc 123'

The generated file containing Boundary conditions, the command used to generate the file and the message displayed on the output window are shown:

The boundary conditions data can be found at: fixed_123.bou

In this example we apply a distributed pressure load on an element towards the right end in the Y-direction, to create a moment along the Z-direction.

-> Plot the faces of the beam using 'plot f beam' command, then use 'view elem' command.

-> Zoom in to the right end and type 'qadd load' and hit enter. Then move the mouse over the element where load needs to be applied and hit 'f' to select that face.

-> Then type 'send load abq pres 10' and hit enter, this command applies a load of 10 on the selected face and saves the load data file.

### Step 3: Writing an Input File for the CCX Solver

Open SciTE and create a new file named 'beam.inp'. This document should be located with the files from the previous articles. Add the following code:

In this example the input file contains the files names of the node, element data, boundary conditions, loads. and other information like material properties are defined in this file. (Note: It is not necessary that node data, boundary conditions, load data be declared in separate files. If the data is simple & small it can directly be typed in the input file)

After the input file is ready, we can submit the job for solving ('Ctrl+F10' or 'Tools>Solve').

a screen shot of the input code for ccx and the output window after solving the code are shown:

### Step 4: Post Processing Results in CGX

->After the solving is done, post-processing can be done by clicking on 'Tools>Post Process' or 'Shift + F10'.

->In this simulation there are two fields available for visualization. They are displacement and stress. Clicking on one of the menu items at this level will select the results field of interest.Then the particular information from the field will need to be specified. This is done by choosing an option from within the -Entity- submenu will result in a colored plot of the selected field information as shown below.

-> The deformed shape can be visualized by selecting 'Viewing>ToggleAdd-Displacement'

To scale up the displacement use the keyboard command scal d 10000.

->Another option to visualize the shape of the deformed structure is through animation. This feature will animate the model’s deformation with user specified scale. Use the menu system and select Animate->Start to view the animation. The image below shows the beam as it oscillates between -100% and 100% of the deformation amplitude.

->The cgx post-processing mode allows the user to specify a cut plane. The plane will cut through the model and show the field results for the intersection of the model and plane. The cut plane is defined by three nodes. The nodes can be specified from the menu system using the Cut menu.

The following image shows the result of cutting the beam top to bottom. The field shown is the normal stress in the x direction. The results are as we would expect, the top of the beam is in tension and the bottom is in compression.

->In addition values from the different result fields can be plotted to a graph for interpretation. The graph command is accessible from the menu system under Graph. The graph is also accessible through the command line.

The nodes are selected three at a time, then the 'tra l 10' command is used to pan the model to the left by translation of 10 units. At which point another few nodes are selected. This process is repeated as the translation command moves towards the free end. Then the last nodes are selected and the right mouse button is used to finish the selection process. Finally, the following graph is promptly displayed.

->While viewing the colored plot of a results field it is also possible to add the nodal values. The following images shows top side of the beam’s free end. The nodal values are printed in yellow using the commmand 'plus nv all y'. These values agree with the results of the graph.

# 5.5 Calculix CCX examples

## Given: open source Finite Element code (Calculix)

 Calculix is an open source Finite Element code (with ABAQUS like input) available at http://www.dhondt.de/

## Find: installation and run Calculix Crunchix(CCX)

 1) For the Disk problem, find

 a)Node info

 b)Element info

 2) Generate 3 meshes of same disk with triangular elements. (increase the number of elements)

 3) Install CCX, run examples and write a report explaining the commands used.

## Solution

### 1) Node, Element Info

 The Node Info and Element Info can be extracted by including a command SEND ALL ABQ, after the mesh has been created.

 The node info, element info is shown below:

### 2) Triangular mesh

 Triangular mesh can be generated if we set the element type to tr3 or tr6.

 Example: 'ELTY all tr3' or ' ELTY all tr6'

 We have used element type 'tr3' to generate the triangular mesh.

 The mesh size, i.e. the number of elements can be increased by increasing the number of divisions on the lines used to generate the geometry.

 The code used to generate the meshes, the meshed geometries are shown below:

### 3) CCX Installation, Basic Examples and Report

#### CCX Examples

 Basic examples demonstrating various features in CCX can be downloaded from:[CCX Examples]

 We were able to successfully run the below examples:

 After the code is run, we need to post process the beam for displacements, forces. To post process the results, click on 'Tools>PostProcess' or hit 'SHIFT+F10'. The below screen shots demonstrate how the displacements and forces can be plotted.

 Displacement Plot:

 Force Plot Plot:

 Plate Mesh

 Stress Plot

 Displacement Plot

 Meshed Geometry

 Stress Plot

 Displacement Plot

#### CCX Report

 Commands Used to generate the .inp file are explained below:

 Note: Few commands used for preprocessing were explained earlier at: [Preprocessing commands explained]
##### NSET

Keyword type: model definition

This option is used to assign nodes to a node set. The parameter NSET containing the name of the set is required (maximum 80 characters), whereas the parameter GENERATE (without value) is optional. If present, nodal ranges can be expressed by their initial value, their final value, and an increment. If a set with the same name already exists, it is reopened and complemented. The name of a set is case insensitive. Internally, it is modified into upper case and a 'N' is appended to denote it as node set.

Example:

• NSET,NSET=N1

1,8,831,208

• NSET,NSET=N2

100,N1

assigns the nodes with number 1, 8, 831 and 208 to (node) set N1 and the nodes with numbers 1, 8, 831, 208 (= set N1) and 100 to set N2

##### BOUNDARY

Keyword type: step or model definition

This option is used to prescribe boundary conditions. This includes:

temperature, displacements and rotations for structures total temperature, mass flow and total pressure for gas networks temperature, mass flow and static pressure for liquid networks temperature, mass flow and fluid depth for channels static temperature, velocity and static pressure for 3D-fluids. For liquids and structures the total and static temperature virtually coincide, therefore both are represented by the term temperature.

The following degrees of freedom are being used:

for structures:

1: translation in the local x-direction

2: translation in the local y-direction

3: translation in the local z-direction

4: rotation about the local x-axis

5: rotation about the local y-axis

6: rotation about the local z-axis

11: temperature

for gas networks:

1: mass flow

2: total pressure

11: total temperature

for liquid networks:

1: mass flow

2: static pressure

11: temperature

for liquid channels:

1: mass flow

2: fluid depth

11: temperature

for 3D-fluids:

1: velocity in the local x-direction

2: velocity in the local y-direction

3: velocity in the local z-direction

8: static pressure

11: static temperature

##### ELEMENT

Keyword type: model definition

With this option elements are defined. There is one required parameter, TYPE and one optional parameter, ELSET. The parameter TYPE defines the kind of element which is being defined. The following types can be selected:

General 3D solids

C3D4 (4-node linear tetrahedral element)

C3D6 (6-node linear triangular prism element)

C3D8 (3D 8-node linear isoparametric element)

C3D8R (the C3D8 element with reduced integration)

C3D15 (15-node quadratic triangular prism element)

C3D20 (3D 20-node quadratic isoparametric element)

C3D20R (the C3D20 element with reduced integration)

ABAQUS 3D solids for heat transfer (names are provided for compatibility)

DC3D4: identical to C3D4

DC3D6: identical to C3D6

DC3D8: identical to C3D8

DC3D10: identical to C3D10

DC3D15: identical to C3D15

DC3D20: identical to C3D20

Shell elements

S6 (6-node triangular shell element)

S8R (the S8 element with reduced integration)

Plane stress elements

CPS6 (6-node triangular plane stress element)

CPS8 (8-node quadratic plane stress element)

CPS8R (the CPS8 element with reduced integration)

Plane strain elements

CPE6 (6-node triangular plane strain element)

CPE8 (8-node quadratic plane strain element)

CPE8R (the CPS8 element with reduced integration)

Axisymmetric elements

CAX6 (6-node triangular axisymmetric element)

CAX8R (the CAX8 element with reduced integration)

Beam elements

B32 (3-node beam element)

B32R (the B32 element with reduced integration)

Special elements

D (3-node fluid element)

GAPUNI (2-node unidirectional gap element)

Example:

• ELEMENT,ELSET=Eall,TYPE=C3D20R

1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,

16,17,18,19,20

defines one 20-node element with reduced integration and stores it in set Eall.

##### EQUATION

Keyword type: model definition

With this option, a linear equation constraint between arbitrary displacement components at any nodes where these components are active can be imposed. The equation is assumed to be homogeneous, and all variables are to be written on the left hand side of the equation. The first variable is considered to be the dependent one, and is subsequently eliminated from the equation, i.e. the corresponding degree of freedom does not show up in the stiffness matrix. This reduces the size of the matrix. A node can only be used once as the dependent node in an equation or in a SPC.

First line:

• EQUATION

Following lines, in a set: First line of set:

Number of terms in the equation. Following lines of set (maximum 12 entries per line):

Node number of the first variable.

Degree of freedom at above node for the first variable.

Value of the coefficient of the first variable.

Node number of the second variable.

Degree of freedom at above node for the second variable.

Value of the coefficient of the second variable.

Etc..

Continue, giving node number, degree of freedom, value of the coefficient, etc. Repeat the above line as often as needed if there are more than four terms in the *EQUATION. Specify exactly four terms per line for each constraint, except for the last line which may have less than four terms.

Example:

• EQUATION

3 3,2,2.3,28,1,4.05,17,1,-8.22

defines an equation of the form , where u, v and w are the displacement for degree of freedom one, two and three, respectively

##### MATERIAL

Keyword type: model definition

This option is used to indicate the start of a material definition. A material data block is defined by the options between a *MATERIAL line and either another *MATERIAL line or a keyword line that does not define material properties. All material options within a data block will be assumed to define the same material. If a property is defined more than once for a material, the last definition is used. There is one required parameter, NAME, defining the name of the material with which it can be referenced in element property options (e.g. *SOLID SECTION). The name can contain up to 80 characters.

Material data requests outside the defined ranges are extrapolated in a constant way and a warning is generated. Be aware that this occasionally occurs due to rounding errors.

First line:

• MATERIAL

Enter the NAME parameter and its value.

Example:

• MATERIAL,NAME=EL

starts a material block with name EL.

##### ELASTIC

Keyword type: model definition, material

This option is used to define the elastic properties of a material. There is one optional parameter TYPE. Default is TYPE=ISO, other values are TYPE=ORTHO and TYPE=ENGINEERING CONSTANTS for orthotropic materials and TYPE=ANISO for anisotropic materials. All constants may be temperature dependent. For orthotropic and fully anisotropic materials, the coefficients satisfy the equation:

$S_{IJ}=D_{IJKL}E_{KL} (Eq 149)$

where is the second Piola-Kirchhoff stress and is the Lagrange deformation tensor. For linear calculations, these reduce to the generic stress and strain tensors.

First line:

• ELASTIC

Enter the TYPE parameter and its value, if needed

Following line for TYPE=ISO:

Young's modulus.

Poisson's ratio.

Temperature.

Repeat this line if needed to define complete temperature dependence.

Example:

• ELASTIC,TYPE=ORTHO

500000.,157200.,400000.,157200.,157200.,300000.,126200.,126200., 126200.,294.

defines an orthotropic material for temperature T=294. Since the definition includes values for only one temperature, they are valid for all temperatures.

##### SOLID SECTION

Keyword type: model definition

This option is used to assign material properties to 3D, plane stress, plane strain and axisymmetric element sets. The parameters ELSET and MATERIAL are required, the parameter ORIENTATION is optional. The parameter ELSET defines the element set to which the material specified by the parameter MATERIAL applies. The parameter ORIENTATION allows to assign local axes to the element set. If activated, the material properties are applied to the local axis. This is only relevant for non isotropic material behavior. For plane stress and plane strain elements the thickness can be specified on the second line. Default is 1.

First line:

• SOLID SECTION

Enter any needed parameters.

Second line (only relevant for plane stress, plane strain and axisymmetric elements; can be omitted for 3D elements):

thickness (plane stress and plane strain elements)

Example:

• SOLID SECTION,MATERIAL=EL,ELSET=Eall,ORIENTATION=OR1

assigns material EL with orientation OR1 to all elements in (element) set Eall.

##### ELSET

Keyword type: model definition

This option is used to assign elements to an element set. The parameter ELSET containing the name of the set is required (maximum 80 characters), whereas the parameter GENERATE (without value) is optional. If present, element ranges can be expressed by their initial value, their final value, and an increment. If a set with the same name already exists, it is reopened and complemented. The name of a set is case insensitive. Internally, it is modified into upper case and a 'E' is appended to denote it as element set.

First line:

• ELSET

Enter any needed parameters and their values.

Following line if the GENERATE parameter is omitted:

List of elements and/or sets of elements previously defined to be assigned to this element set (maximum 16 entries per line).

Example:

• ELSET,ELSET=E1,GENERATE

20,25

• ELSET,ELSET=E2

E1,50,51

assigns the elements with numbers 20, 21, 22, 23, 24 and 25 to element set E1 and the elements with numbers 20, 21, 22, 23, 24, 25 (= set E1), 50 and 51 to element set E2.

##### SURFACE

Keyword type: model definition

This option is used to define surfaces made up of nodes or surfaces made up of element faces. A mixture of nodes and element faces belonging to one and the same surface is not possible. There are two parameters: NAME and TYPE. The parameter NAME containing the name of the surface is required. The TYPE parameter takes the value NODE for nodal surfaces and ELEMENT for element face surfaces. Default is TYPE=ELEMENT

Example:

• SURFACE,NAME=left,TYPE=NODE

part,

1,

8

assigns the nodes with number 1, and 8 and the nodes belonging to node set part to a surface with name left.

Example:

• SURFACE,NAME=new

38,S6

assigns the face 6 of element 38 to a surface with name new

##### CONTACT PAIR

Keyword type: model definition

This option is used to express that two surfaces can make contact. There is one required parameter: INTERACTION, and two optional parameters: SMALL SLIDING and ADJUST. The dependent surface is called the slave surface, the independent surface is the master surface. Surfaces are defined using the *SURFACE. The dependent surface can be defined as a nodal surface (option TYPE=NODE on the *SURFACE keyword) or as an element face surface (default for the *SURFACE card), whereas the independent surface has to be defined as an element face surface.

The INTERACTION parameter takes the name of the surface interaction (keyword *SURFACE INTERACTION) which applies to the contact pair. The surface interaction defines the nature of the contact (hard versus soft contact..)

First line:

• CONTACT PAIR

enter the required parameter INTERACTION and its parameter.

Following line:

Name of the slave surface (can be nodal or element face based).

Name of the master surface (must be based on element faces).

Example:

dep,ind

defines a contact pair consisting of the surface dep as dependent surface and the element face surface ind als independent surface. The name of the surface interaction is IN1. All slave nodes for which the clearance is smaller than or equal to 0.01 will moved onto the master surface.

##### SURFACE INTERACTION

Keyword type: model definition

This option is used to indicate the start of a surface interaction definition. A surface interaction data block is defined by the options between a *SURFACE INTERACTION line and either another *SURFACE INTERACTION line or a keyword line that does not define surface interaction properties. All surface interaction options within a data block will be assumed to define the same surface interaction. If a property is defined more than once for a surface interaction, the last definition is used. There is one required parameter, NAME, defining the name of the surface interaction with which it can be referenced in surface interaction property options (e.g. *SURFACE BEHAVIOR). The name can contain up to 80 characters.

Surface interaction data requests outside the defined ranges are extrapolated in a constant way. Be aware that this occasionally occurs due to rounding errors.

First line:

• SURFACE INTERACTION

Enter the NAME parameter and its value.

Example:

• SURFACE INTERACTION,NAME=SI1

starts a material block with name SI1.

##### BEAM SECTION

Keyword type: model definition

This option is used to assign material properties to beam element sets. The parameters ELSET, MATERIAL and SECTION are required, the parameters ORIENTATION, OFFSET1 and OFFSET2 are optional. The parameter ELSET defines the shell element set to which the material specified by the parameter MATERIAL applies. The parameter ORIENTATION allows to assign local axes to the element set. If activated, the material properties are applied to the local axis. This is only relevant for non isotropic material behavior.

The parameter SECTION defines the cross section of the beam and can take the value RECT for a rectangular cross section and CIRC for an elliptical cross section. A rectangular cross section is defined by its thickness in two perpendicular directions, an elliptical cross section is defined by the length of its principal axes. These directions are defined by specifying direction 1 on the third line of the present keyword card

First line:

• BEAM SECTION

Enter any needed parameters.

Second line:

thickness in 1-direction

thickness in 2-direction

Third line:

global x-coordinate of a unit vector in 1-direction (default:0)

global y-coordinate of a unit vector in 1-direction (default:0)

global z-coordinate of a unit vector in 1-direction (default:-1)

Example:

• BEAM SECTION,MATERIAL=EL,ELSET=Eall,OFFSET1=-0.5,SECTION=RECT

3.,1.

1.,0.,0.

assigns material EL to all elements in (element) set Eall. The reference line is in 1-direction on the back surface, in 2-direction on the central surface. The thickness in 1-direction is 3 unit lengths, in 2-direction 1 unit length. The 1-direction is the global x-axis.

##### STEP

Keyword type: step

This card describes the start of a new STEP. PERTURBATION, NLGEOM, INC, INCF and TURBULENCE MODEL are the optional parameters.

First and only line:

• STEP

Enter any needed parameters and their values

Example:

• STEP,INC=1000,INCF=20000,TURBULENCE MODEL=SST

starts a step and increases the maximum number of thermomechanical increments to complete the step to 1000. The maximum number of 3D fluid increments is set to 20000 and for the turbulence model the SST model was chosen.

##### STATIC

Keyword type: step

This procedure is used to perform a static analysis. The load consists of the sum of the load of the last *STATIC step and the load specified in the present step with replacement of redefined loads.

First line:

• STATIC

Enter any needed parameters and their values. Second line (only relevant for nonlinear analyses; for linear analyses, the step length is always 1) Initial time increment. This value will be modified due to automatic incrementation, unless the parameter DIRECT was specified (default 1.). Time period of the step (default 1.). Minimum time increment allowed. Only active if DIRECT is not specified. Default is the initial time increment or 1.e-5 times the time period of the step, whichever is smaller. Maximum time increment allowed. Only active if DIRECT is not specified. Default is 1.e+30.

Example:

• STATIC,DIRECT

.1,1.

defines a static step and selects the SPOOLES solver as linear equation solver in the step (default). If the step is a linear one, the other parameters are of no importance. If the step is nonlinear, the second line indicates that the initial time increment is .1 and the total step time is 1. Furthermore, the parameter DIRECT leads to a fixed time increment. Thus, if successful, the calculation consists of 10 increments of length 0.1.

Keyword type: step

This option allows concentrated forces to be applied to any node in the model which is not fixed by a single or multiple point constraint. Optional parameters are OP, AMPLITUDE, TIME DELAY, LOAD CASE and SECTOR. OP can take the value NEW or MOD. OP=MOD is default and implies that the concentrated loads applied to different nodes are kept over all steps starting from the last perturbation step. Specifying a force in a node for which a force was defined in a previous step replaces this value. OP=NEW implies that all previously applied concentrated loads are removed. If multiple *CLOAD cards are present in a step this parameter takes effect for the first *CLOAD card only.

First line:

Enter any needed parameters and their value.

Following line:

Node number or node set label.

Degree of freedom.

Repeat this line if needed.

Example:

1000,3,10.3

removes all previous point load forces and applies a force with magnitude 10.3 and amplitude A1 (shifted in positive time direction by 20 time units) for degree of freedom three (global if no transformation was defined for node 1000, else local) of node 1000.

##### NODE PRINT

Keyword type: step

This option is used to print selected nodal variables in file jobname.dat. The following variables can be selected:

Displacements (key=U)

Structural temperatures, total temperatures in networks and static temperatures in 3D fluids (key=NT)

Pressures in networks (key=PN). These are the total pressures for gases, static pressures for liquids and liquid depth for channels. The fluid section types dictate the kind of network.

Static pressures in 3D fluids (key=PS)

Velocities in 3D fluids (key=V)

Mass flows in networks (key=MF)

External forces (key=RF)

External concentrated heat sources (key=RFL)

First line:

• NODE PRINT

Enter the parameter NSET and its value.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:

• NODE PRINT,NSET=N1

RF

requests the storage of the reaction forces in the nodes belonging to (node) set N1 in the .dat file.

##### EL PRINT

Keyword type: step

This option is used to print selected element variables in an ASCII file with the name jobname.dat. Some of the element variables are printed in the integration points, some are whole element variables. The following variables can be selected:

Integration point variables

true (Cauchy) stress (key=S)

strain (key=E). This is the Lagrangian strain for (hyper)elastic materials and incremental plasticity and the Eulerian strain for deformation plasticity.

equivalent plastic strain (key=PEEQ)

equivalent creep strain (key=CEEQ; is converted internally into PEEQ since the viscoplastic theory does not distinguish between the two; consequently, the user will find PEEQ in the dat file, not CEEQ)

the energy density (key=ENER)

the internal state variables (key=SDV)

heat flux (key=HFL)

Whole element variables

the internal energy (key=ELSE)

the kinetic energy (key=ELKE)

the volume (key=EVOL)

First line:

• EL PRINT

Enter the parameter ELSET and its value.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:

• EL PRINT,ELSET=Copper

E

requests to store the strains at the integration points in the elements of set Copper in the .dat file.

##### NODE FILE

Keyword type: step

This option is used to print selected nodal variables in file jobname.frd for subsequent viewing by CalculiX GraphiX. The following variables can be selected:

Displacements (key=U)

Displacements: magnitude and phase (key=PU, only for *STEADY STATE DYNAMICS calculations) and *FREQUENCY calculations with cyclic symmetry).

Maximum displacements orthogonal to a given vector at all times for *FREQUENCY calculations with cyclic symmetry. The components of the vector are the coordinates of a node stored in a node set with the name RAY. This node and node set must have been defined by the user (key=MAXU).

Temperatures (key=NT). This includes both structural temperatures and total fluid temperatures in a network.

Temperatures: magnitude and phase (key=PNT, only for *STEADY STATE DYNAMICS calculations).

External forces (key=RF)

External concentrated heat sources (key=RFL)

Static temperatures in networks and 3D fluids (key=TS)

Total temperatures in networks and 3D fluids (key=TT)

Mass flows in networks (key=MF). The mass flow through a fluid element is stored in the middle node of the element. In the end nodes the mass flow is not unique, since more than two element can be connected to the node. For end nodes the sum of the mass flow leaving the node is stored. Notice that at nodes where mass flow is leaving the network the value will be wrong if no proper exit element (with node number 0) is attached to that node.

Velocities in 3D fluids (key=V)

Mach numbers in 3D compressible fluids (key=MACH)

Static pressures in liquid networks and 3D fluids (key=PS)

Total pressures in gas networks and 3D fluids (key=PT)

Fluid depth in channel networks (key=DEPT)

Critical depth in channel networks (key=HCRI)

Pressure coefficient in 3D compressible fluids (key=CP)

First line:

• NODE FILE

Enter any needed parameters and their values.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:

• NODE FILE,FREQUENCY=2,TIME POINTS=T1

RF,NT

requests the storage of reaction forces and temperatures in the .frd file every second increment. In addition, output will be stored for all time points defined by the T1 time points sequence.

##### EL FILE

Keyword type: step

This option is used to save selected element variables averaged at the nodal points in a frd file (extension .frd) for subsequent viewing by CalculiX GraphiX. The following element variables can be selected:

true (Cauchy) stress (key=S). For beam elements this tensor is replaced by the section forces if SECTION FORCES is selected.

stress: magnitude and phase (key=PHS, only for *STEADY STATE DYNAMICS calculations and *FREQUENCY calculations with cyclic symmetry). maximum worst principal stress at all times for *FREQUENCY calculations with cyclic symmetry. It is stored for nodes belonging to the node set with name STRESSDOMAIN. This node set must have been defined by the user with the *NSET command. The worst principal stress is the maximum of the absolute value of the principal stresses (key=MAXS).

strain (key=E). This is the Lagrangian strain for (hyper)elastic materials and incremental plasticity and the Eulerian strain for deformation plasticity.

equivalent plastic strain (key=PEEQ)

equivalent creep strain (key=CEEQ; is converted internally into PEEQ since the viscoplastic theory does not distinguish between the two; consequently, the user will find PEEQ in the frd file, not CEEQ)

the energy density (key=ENER)

the internal state variables (key=SDV)

heat flux (key=HFL)

First line:

• EL FILE

Enter any needed parameters and their values.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:

• EL FILE

S,PEEQ

requests that the (Cauchy) stresses and the equivalent plastic strain is stored in .frd format for subsequent viewing with CalculiX GraphiX.

##### END STEP

Keyword type: step

This option concludes the definition of a step.

First and only line:

• END STEP

Example:

• END STEP

concludes a step. Each *STEP card must at some point be followed by an *END STEP card.

--Eml5526.s11.team5.vijay 01:48, 24 March 2011 (UTC)

# 5.6 Quad Lagrangian Element Basis Function

## Given: Quad Lagrangian Basis Function

Lagrangian element with 3 nodes per element.

## Find

Plot the basis functions for each node

## Solution

Lagrange basis functions for an element of order m can be written as follow:

$L_{i,m}(x)=\prod_{k=1}^{m} \frac{x-x_{k}}{x_{i}-x_{k}} \forall k\ne i$

For a Quad lagrangian Basis function there are 3 nodes per element (m=3). We can write the Quad Lagrangian Basis Function as follows:

$L_{1,3}(x)= (\frac{x-x_{2}}{x_{1}-x_{2}})(\frac{x-x_{3}}{x_{1}-x_{3}})$

$L_{2,3}(x)= (\frac{x-x_{1}}{x_{2}-x_{1}})(\frac{x-x_{3}}{x_{2}-x_{3}})$

$L_{3,3}(x)= (\frac{x-x_{1}}{x_{3}-x_{1}})(\frac{x-x_{2}}{x_{3}-x_{2}})$

Let us consider the element shown below:

In the above element $X_{1}=0, X_{2}=0.5, X_{3}=1$

Substituting the above values, the basis functions at each node are as follows:

--Eml5526.s11.team5.vijay 01:49, 24 March 2011 (UTC)

# HW 1.1

Theory:

In elastodynamic case, there are 2 major contributors that affect body's motion:

• Static Forces
• Body's Inertial effects (resistance to motion caused by Static forces)

In accord to conservation laws as well as Newton's 3rd law, these 2 opposing contributors must be equal to each other.

$F_{static} = F_{dynamic} \Rightarrow \vartriangle F = ma$

where:

F = force (N)

m = mass(kg)

a = acceleration ($\frac{m}{s^2}$)

The assumed differential element shape is of trapezoid with force acting on both of its bases. Dynamic equation that resulted is shown below:

$-p(x) + f(x + \frac {\vartriangle x}{2})\vartriangle x + p(x+\vartriangle x) = ( \hat m(x) * \vartriangle x) a(t)$

where:

$\hat m$ = mass per unit length $(\frac{kg}{m})$

$\frac {p(x+\vartriangle x) - p(x)}{\vartriangle x} + f(x + \frac {\vartriangle x}{2}) = \frac {( \hat m(x) *\vartriangle x)}{\vartriangle x} a(t) = \hat m(x) a(t)$

$\Rightarrow a(t) = \frac{\vartriangle v(t)}{\vartriangle t}$

$\Rightarrow \vartriangle v(t) = \frac{u(t + \vartriangle t) - u(t)}{\vartriangle t}$

Therefore:

$\frac {p(x+\vartriangle x) - p(x)}{\vartriangle x} + f(x + \frac {\vartriangle x}{2}) = \hat m(x) \frac{1}{\vartriangle t}\frac{u(t + \vartriangle t) - u(t)}{\vartriangle t}$

$\lim_{\vartriangle x \to 0} \frac {p(x+\vartriangle x) - p(x)}{\vartriangle x} + f(x + \frac {\vartriangle x}{2}) = \lim_{\vartriangle t \to 0} \hat m(x) \frac{1}{\vartriangle t}\frac{u(t + \vartriangle t) - u(t)}{\vartriangle t}$

$\Rightarrow \lim_{\vartriangle x \to 0} \frac {p(x+\vartriangle x) - p(x)}{\vartriangle x} = \frac {\partial p}{\partial x}$

$\Rightarrow \lim_{\vartriangle t \to 0} \frac{1}{\vartriangle t} = \frac{1}{\partial t}$

$\Rightarrow \lim_{\vartriangle t \to 0} \hat m \frac{u(t + \vartriangle t) - u(t)}{\vartriangle t} = \hat m(x) \frac {\partial u}{\partial t}$

Therefore, substituting back into equation we get:

$\frac {\partial p}{\partial x} + f(x) = \hat m(x) \frac{1}{\partial t} \frac {\partial u}{\partial t} \Rightarrow \frac {\partial p}{\partial x} + f(x) = \hat m(x) \frac {\partial^2 u}{\partial t^2}$



Force = Stress * Area,

$p(x) = \sigma(x) * A(x)$

Stress = Strain * Modulus of Elasticity,

$\sigma(x) = \epsilon(x) * E(x)$

where

$\epsilon(x) = \frac {\partial u}{\partial x}$

From the above relations:

$p(x) = E(x) * A(x) * \frac {\partial u}{\partial x}$

Mass / Unit Length = Density * Area

$\hat m(x) = \rho(x) * A(x)$

Substituting the above relations we get

$\frac {\partial}{\partial x} \left[E(x)A(x)\frac{\partial u}{\partial x}\right] + f(x) = \rho(x) A(x) \frac {\partial^2 u} {\partial t^2}$



If we consider a case in which the bar has a rectangular cross section

$A(x) = h(x) * b$

where h(x) = Height of bar at x,

b = width of the bar (constant) [i.e.Independent of x,t]

Replacing A(x) with h(x)*b in our equation we get:

$\frac {\partial}{\partial x} \left[E(x)h(x)*b*\frac{\partial u}{\partial x}\right] + f(x) = \rho(x) h(x)*b* \frac {\partial^2 u} {\partial t^2}$

b is independent and can be brought outside the partial derivative of the first term

The equation is divided by b to give

$\frac {\partial}{\partial x} \left[E(x)h(x)\frac{\partial u}{\partial x}\right] + \frac{f(x)}{b} = \rho(x) h(x) \frac {\partial^2 u} {\partial t^2}$

Contribution: --Eml5526.s11.team5 20:37, 20 January 2011 (UTC) --Eml5526.s11.team5.vijay 06:07, 21 January 2011 (UTC)

|}

# HW 2.6

## Given

 $F=\left \{ 1,cos i\omega x, sin i\omega x \right \} , i = \left \{ 1, 2 \right \}$ $\displaystyle$

 on the interval, $\Omega =[0,T]$ $\displaystyle$

 i.e.,$F=\left \{ 1,cos \omega x, cos 2\omega x, sin \omega x, sin 2\omega x \right \} = \left \{b_1(x), b_2(x), b_3(x), b_4(x), b_5(x) \right \}$ $\displaystyle$

## Find

 1)Construct $\Gamma_{5x5}(F)$. Observe the properties of $\Gamma$ $\displaystyle$

 2) Find the det $\Gamma$(F) $\displaystyle$

 3) Conclude F is orthogonal Basis, i.e., $\Gamma_{ij} = \delta_{ij}$ $\displaystyle$

## Trigonometric Properties Used

 $\implies \int_{0}^{T} cos (n \omega x) = 0, \forall n = \left \{ 1, 2, 3, ....\right \}, \omega = \frac{2\pi}{T}$ $\implies \int_{0}^{T} sin (n \omega x) = 0, \forall n = \left \{ 1, 2, 3, ....\right \}, \omega = \frac{2\pi}{T}$
 $sin(x)cos(y) = \left [sin(x+y) + sin(x-y) \right ]/2,$ $cos(x)sin(y) = \left [sin(x+y) - sin(x-y) \right]/2,$ $cos(x)cos(y) = \left [cos(x-y) + cos(x+y) \right]/2,$ $sin(x)sin(y) = \left [cos(x-y) - cos(x+y)\right]/2.$

## Solution

### 1)

 The Gram Matrix is defined as follows:, $\Gamma_{ij}==\int_\Omega b_i(x)b_j(x)\,dx$ $\displaystyle$

 It Can be observed that $\Gamma_{ij}=\Gamma_{ji}$, i.e. The $\Gamma$ Matrix is symmetric $\displaystyle$

 [Since $\int_\Omega b_i(x)b_j(x)\,dx = \int_\Omega b_j(x)b_i(x)\,dx$] $\displaystyle$

 Now Let us calculate each element in the $\Gamma$ Matrix: $\displaystyle$

 $\Gamma_{1,1} = \int_{0}^{T} dx = \left [ x \right ]_{0}^{T} = T$ $\displaystyle$

 $\Gamma_{1,2} = \Gamma_{2,1} = \int_{0}^{T} cos \omega x dx = 0$ $\displaystyle$

 $\Gamma_{1,3} = \Gamma_{3,1} = \int_{0}^{T} cos 2\omega x dx = 0$ $\displaystyle$

 $\Gamma_{1,4} = \Gamma_{4,1} = \int_{0}^{T} sin \omega xdx = 0$ $\displaystyle$

 $\Gamma_{1,5} = \Gamma_{5,1} = \int_{0}^{T} sin 2\omega x dx = 0$ $\displaystyle$

 $\Gamma_{2,2} = \int_{0}^{T} cos^{2} \omega x dx = \int_{0}^{T} \frac {cos 2\omega x +1}{2} = \frac {T}{2}$ $\displaystyle$

 $\Gamma_{2,3} = \Gamma_{3,2} = \int_{0}^{T} cos \omega x cos 2\omega x dx = \frac {1}{2} \int_{0}^{T} \left [ cos \omega x + cos 3\omega x \right ] dx = 0$ $\displaystyle$

 $\Gamma_{2,4} = \Gamma_{4,2} = \int_{0}^{T} cos \omega x sin \omega x dx = \frac {1}{2} \int_{0}^{T} sin 2\omega x dx = 0$ $\displaystyle$

 $\Gamma_{2,5} = \Gamma_{5,2} = \int_{0}^{T} cos \omega x sin 2\omega x dx = \frac {1}{2} \int_{0}^{T} \left [ sin 3\omega x + sin \omega x \right ] dx = 0$ $\displaystyle$

 $\Gamma_{3,3} = \int_{0}^{T} cos^{2} \omega 2x dx = \frac {1}{2} \int_{0}^{T} \left [ cos 4\omega x + 1 \right ] dx = \frac {T}{2}$ $\displaystyle$

 $\Gamma_{3,4} = \Gamma_{4,3} = \int_{0}^{T} cos 2\omega x sin \omega x dx = \frac {1}{2} \int_{0}^{T} \left [ sin 3\omega x - sin \omega x \right ] dx = 0$ $\displaystyle$

 $\Gamma_{3,5} = \Gamma_{5,3} = \int_{0}^{T} cos 2\omega x sin 2\omega x dx = \frac {1}{2} \int_{0}^{T} sin 4\omega x dx = 0$ $\displaystyle$

 $\Gamma_{4,4} = \int_{0}^{T} sin^{2} \omega x dx = \frac {1}{2} \int_{0}^{T} \left [ 1 - cos 2\omega x \right ] dx = \frac{T}{2}$ $\displaystyle$

 $\Gamma_{4,5} = \Gamma_{5,4} = \int_{0}^{T} sin \omega x sin 2\omega x dx = \frac {1}{2} \int_{0}^{T} \left [ cos \omega x - cos 3\omega x \right ] dx = 0$ $\displaystyle$

 $\Gamma_{4,4} = \int_{0}^{T} sin^{2} 2\omega x dx = \frac {1}{2} \int_{0}^{T} \left [ 1 - cos 4\omega x \right ] dx = \frac{T}{2}$ $\displaystyle$
 $a) \Gamma_{5x5} \left (F \right) = \begin{bmatrix} T & 0 & 0 & 0 & 0\\ 0 & \frac{T}{2} & 0 & 0 & 0\\ 0 & 0 & \frac{T}{2} & 0 & 0\\ 0 & 0 & 0 & \frac{T}{2} & 0\\ 0 & 0 & 0 & 0 & \frac{T}{2} \end{bmatrix}$
 b) It is observed that the $\Gamma$ Matrix is $\rightarrow$ Diagnol $\rightarrow$ Symmetric

### 2)

 $\det \left ( \Gamma \right ) = \frac {T^{5}}{16}$

### 3)

 The $\Gamma$ Matrix is Not orthogonal

--Eml5526.s11.team5.vijay 19:01, 1 February 2011 (UTC)--

# Proof for Equations in 2.6

Consider $\int_{0}^{T}cos (2n\pi /T)x$

$= \left [ \frac{T}{2n\pi}sin (2n\pi /T)x \right ]_{0}^{T}$

$= \frac{T}{2n\pi} \left [ sin (2n\pi /T)*T - sin (2n\pi /T) *0 \right ]$

$= \frac{T}{2n\pi} \left [ sin (2n\pi) - sin 0 \right ]$

$= \frac{T}{2n\pi} \left [ 0 - 0 \right ]$

Since, $sin(2n\pi) = 0, \forall n = \left \{ 1, 2, 3, ....\right \},$ $sin 0 = 0$

$= 0$

$\implies \int_{0}^{T}cos (2n\pi /T)x = 0, \forall n = \left \{ 1, 2, 3, .... \right \}$

Let $\omega = 2 * \pi / T$,

 $\implies \int_{0}^{T} cos (n \omega x) = 0, \forall n = \left \{ 1, 2, 3, ....\right \}, \omega = \frac{2\pi}{T}$

Similarly consider

$\int_{0}^{T} sin (n \omega x) \forall n = \left \{ 1, 2, 3, ....\right \}, \omega = \frac{2\pi}{T}$

$= \int_{0}^{T}sin (2n\pi /T)x$

$= \left [ - \frac{T}{2n\pi}cos (2n\pi /T)x \right ]_{0}^{T}$

$= - \frac{T}{2n\pi} \left [ cos (2n\pi /T)*T - cos (2n\pi /T) *0 \right ]$

$= - \frac{T}{2n\pi} \left [ cos (2n\pi) - cos 0 \right ]$

$= - \frac{T}{2n\pi} \left [ 1 - 1 \right ]$

Since, $cos(2n\pi) = 1, \forall n = \left \{ 1, 2, 3, ....\right \},$ $cos 0 = 1$

$= 0$

 $\implies \int_{0}^{T} sin (n \omega x) = 0, \forall n = \left \{ 1, 2, 3, ....\right \}, \omega = \frac{2\pi}{T}$

# HW 3.3: 2D Truss Problem

## Problem Statement

a 3-Bar structure is subjected to a load P(=1000N) at C as shown in Figure 3.3.1.

The cross sectional areas are given as:

BC: 0.02

BF,BD: 0.01

All the bars are made of the same material whose Youngs modulus $E = 10^{11}$ Pa.

## Find

a) Construct the global stiffness matrix and Load matrix

b) Partition the matrices and solve for the x,y displacements at B, x displacement at D.

c) Find the stress in the 3 bars BC, BD, BF.

d) Find the Reactions at the nodes C, D, F.

## Solution

$\bold {a)}$

The structure can be divided into 3 elements:

Element 1 = Bar BF

Element 2 = Bar BC

Element 3 = Bar BD (Element numbers are shown in the above figure)

Before we construct the Stiffness matrix for each bar element, let us review the element stiffness matrix for a bar element having 2 degrees of freedom (X,Y displacements) and inclined at an angle $\phi$ with the X-axis.

 $K^{e}=k^{e} \begin{bmatrix} \cos^2 \phi^e & \cos \phi^e \sin \phi^e & -\cos^2 \phi^e & -\cos \phi^e \sin \phi^e \\ \cos \phi^e \sin \phi^e & \sin^2 \phi^e & -\cos^2 \phi^e &-\sin^2 \phi^e \\ -\cos^2 \phi^e & -\cos \phi^e \sin \phi^e & \cos^2 \phi^e & \cos \phi^e \sin \phi^e \\ -\cos \phi^e \sin \phi^e & -\sin^2 \phi^e & \cos^2 \phi^e &\sin^2 \phi^e \\ \end{bmatrix}$


(ref: Fish and Belytschko P30, Eq (2.47))

where

 $k^e= \frac{A^eE^e}{l^e}$ $\displaystyle$

$A^e=$ Cross sectional area of the bar element

$E^e=$ Young's Modulus of the bar element

$l^e=$ Length of the bar element

$\bold {Element 1:}$

Element 1 connects global nodes 1,4.

$\phi^{(1)}= \frac{\pi}{4}$

$\quad A^{(1)}= 0.01$

$\quad E^{(1)}=10^{11}$

$l^{(1)}=\sqrt{2}$

 $\implies K^{(1)}=\frac{10^{9}}{\sqrt{2}} \begin{matrix} \begin{bmatrix} \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2}\\ \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2}\\ -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & \frac{1}{2}\\ -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{bmatrix} \begin{matrix} \\ (1)\\ \\ (2)\\ \\ \end{matrix} \\ \begin{matrix} &(1) & & (2) & \end{matrix} \end{matrix}$


$\bold {Element 2:}$

Element 2 connects global nodes 1,3.

$\phi^{(2)}= \frac{\pi}{2}$

$\quad A^{(2)}= 0.02$

$\quad E^{(2)}=10^{11}$

$\quad l^{(2)}=1$

 $\implies K^{(2)}=2*10^{9} \begin{matrix} \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & -1\\ 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 1\\ \end{bmatrix} \begin{matrix} \\ (1)\\ \\ (3)\\ \\ \end{matrix} \\ \begin{matrix} &(1) & & (3) & \end{matrix} \end{matrix}$


$\bold {Element 3:}$

Element 3 connects global nodes 1,4.

$\phi^{(3)}= \frac{3\pi}{4}$

$\quad A^{(3)}= 0.01$

$\quad E^{(3)}=10^{11}$

$l^{(3)}=\sqrt{2}$

 $\implies K^{(3)}=\frac{10^{9}}{\sqrt{2}} \begin{matrix} \begin{bmatrix} \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2}\\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}\\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}\\ \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \end{bmatrix} \begin{matrix} \\ (1)\\ \\ (4)\\ \\ \end{matrix} \\ \begin{matrix} &(1) & & (4) & \end{matrix} \end{matrix}$


Direct Assembly Gives the global matrices:

 $\implies K={10^{9}} \begin{matrix} \begin{bmatrix} \frac{1}{\sqrt2}& 0 & -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & 0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2}\\ 0 & 2+\frac{1}{\sqrt2} & -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & -2 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} \\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 \\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -2 & 0 & 0 & 0 & 2 & 0 & 0\\ -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} \\ \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \\ \end{bmatrix} \begin{matrix} \\ (1)\\ \\ (2)\\ \\ (3)\\ \\ (4)\\ \\ \end{matrix} \\ \begin{matrix} (1)& & &(2)& & &(3)& & &(4) \end{matrix} \end{matrix}$

 $d=\begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x}\\ 0\\ 0\\ 0\\ 0\\ 0 \end{bmatrix}, f=1000\begin{bmatrix} 1\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0 \end{bmatrix}, r=\begin{bmatrix} 0\\ 0\\ 0\\ r_{2y}\\ r_{3x}\\ r_{3y}\\ r_{4x}\\ r_{4y} \end{bmatrix}$

$\bold {b)}$

$\implies 10^9\begin{bmatrix} \frac{1}{\sqrt2}& 0 & -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & 0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2}\\ 0 & 2+\frac{1}{\sqrt2} & -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & -2 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} \\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 \\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -2 & 0 & 0 & 0 & 2 & 0 & 0\\ -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} \\ \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \\ \end{bmatrix} \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x}\\ 0\\ 0\\ 0\\ 0\\ 0 \end{bmatrix} = \begin{bmatrix} 1000\\ 0\\ 0\\ r_{2y}\\ r_{3x}\\ r_{3y}\\ r_{4x}\\ r_{4y} \end{bmatrix}$

It can be observed that partitioning needs to be done after 3 rows and 3 columns.

The reduced system of equations are:

$K_{F}= 10^9\begin{bmatrix} \frac{1}{\sqrt2}& 0 & -\frac{1}{2\sqrt2}\\ 0 & 2+\frac{1}{\sqrt2} & -\frac{1}{2\sqrt2}\\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \end{bmatrix}, K_{EF}= 10^9\begin{bmatrix} -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \\ 0 & 0 & 0\\ 0 & -2 & 0\\ -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 \\ \frac{1}{2\sqrt2}& -\frac{1}{2\sqrt2} & 0 \end{bmatrix}, K_{E}= 10^9\begin{bmatrix} \frac{1}{2\sqrt2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 0 & 0\\ 0 & 0 & 0 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} \\ 0 & 0 & 0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \end{bmatrix},$

$d_{F}= \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x} \end{bmatrix}, \bar d_{E}= \begin{bmatrix} 0\\ 0\\ 0\\ 0\\ 0 \end{bmatrix} = \begin{bmatrix} \bar u_{2y}\\ \bar u_{3x}\\ \bar u_{3y}\\ \bar u_{4x}\\ \bar u_{4y} \end{bmatrix}, f_{F}= \begin{bmatrix} 1000\\ 0\\ 0 \end{bmatrix}, r_{E}= \begin{bmatrix} r_{2y}\\ r_{3x}\\ r_{3y}\\ r_{4x}\\ r_{4y} \end{bmatrix},$

 $\implies \begin{bmatrix} K_{F} & K_{EF}^{T}\\ K_{EF} & K_{E} \end{bmatrix} \begin{bmatrix} d_{F}\\ \bar d_{E} \end{bmatrix} = \begin{bmatrix} f_{F}\\ r_{E} \end{bmatrix}$


In other words:

 $K_{F} d_{F} + K_{EF}^T \bar d_{E} = f_{F}$ $\displaystyle (Eq 3.6.1)$

 $K_{EF} d_{F} + K_{E} \bar d_{E} = r_{E}$ $\displaystyle (Eq 3.6.2)$

but,

 $d_{E}=\begin{bmatrix} 0&0&0&0&0\end{bmatrix}^{T}$, $\displaystyle (Eq 3.6.3)$

Combining Eq 3.6.1, Eq 3.6.3:

  $\quad K_{F} d_{F} = f_{F}$  $\displaystyle (Eq 3.6.4)$

Combining Eq 3.6.2, Eq 3.6.3:

  $\quad K_{EF} d_{F} = r_{E}$  $\displaystyle (Eq 3.6.5)$

from Eqn 3.6.4:

$d_{F} = K_{F}^{-1} f_{F}$

$\therefore d_{F}= \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x} \end{bmatrix} = \frac {1}{10^9}\begin{bmatrix} \frac{1}{\sqrt2}& 0 & -\frac{1}{2\sqrt2}\\ 0 & 2+\frac{1}{\sqrt2} & -\frac{1}{2\sqrt2}\\ -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \end{bmatrix}^{-1} \begin{bmatrix} 1000\\ 0\\ 0 \end{bmatrix}$

 $\therefore d_{F}= \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x} \end{bmatrix} = \frac{1}{10^6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 1+2\sqrt{2} \end{bmatrix}$

 X,Y-displacement at B(in mm) $\quad = \left( u_{1x}, u_{1y} \right)=\left( 3.3284*10^{-3}, 5*10^{-4} \right)$ X-displacement at D(in mm) $\quad = \left( u_{2x} \right)=\left( 3.8284*10^{-3} \right)$

$\bold {c)}$

The stress in each element can be computed as follows:

$\sigma^e = E^e \frac{u_{2x}^{'e}-u_{1x}^{'e}}{l^{e}} = \frac{E^e}{l^e} \begin{bmatrix} -1 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} u_{1x}^{'e}\\ u_{1y}^{'e}\\ u_{2x}^{'e}\\ u_{2y}^{'e} \end{bmatrix} =\frac{E^e}{l^e} \begin{bmatrix} -\cos \phi^e & -\sin \phi^e & \cos \phi^e & \sin \phi^e \end{bmatrix} \bold d^e$


(ref: Fish and Belytschko P30, Eq (2.47))

Computing the Stresses in each bar (element) using the above formula:

$\bold {Element 1:}$

$\bold d^{(1)}= \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x}\\ u_{2y} \end{bmatrix} = \frac{1}{10^6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 1+2\sqrt{2}\\ 0 \end{bmatrix}$

$\implies \sigma^{(1)}=\frac{10^{11}}{\sqrt{2}} \begin{bmatrix} -\frac{1}{\sqrt 2} & -\frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \end{bmatrix}* \frac{1}{10^6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 1+2\sqrt{2}\\ 0 \end{bmatrix}$

 $\implies \sigma^{(1)}=0$

$\bold {Element 2:}$

$\bold d^{(2)}= \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{3x}\\ u_{3y} \end{bmatrix} = \frac{1}{10^6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 0\\ 0 \end{bmatrix}$

$\implies \sigma^{(2)}=10^{11} \begin{bmatrix} 0 & -1 & 0 & 1 \end{bmatrix}* \frac{1}{10^6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 0\\ 0 \end{bmatrix}$

 $\implies \sigma^{(2)}=-50*10^3 = -50KPa$

$\bold {Element 3:}$

$\bold d^{(3)}= \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{4x}\\ u_{4y} \end{bmatrix} = \frac{1}{10^6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 0\\ 0 \end{bmatrix}$

$\implies \sigma^{(3)}=\frac{10^{11}}{\sqrt{2}} \begin{bmatrix} \frac{1}{\sqrt 2} & -\frac{1}{\sqrt 2} & -\frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \end{bmatrix}* \frac{1}{10^6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 0\\ 0 \end{bmatrix}$

 $\implies \sigma^{(3)}=\sqrt{2} * 10^5Pa= 141.42KPa$

$\bold {d)}$

Substituting $d_{F}$ in Eqn 3.6.2, we calculate the reactions forces ($r_{E}$ matrix)

$r_{E}= \begin{bmatrix} r_{2y}\\ r_{3x}\\ r_{3y}\\ r_{4x}\\ r_{4y} \end{bmatrix} = 10^9\begin{bmatrix} -\frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \\ 0 & 0 & 0\\ 0 & -2 & 0\\ -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & 0 \\ \frac{1}{2\sqrt2}& -\frac{1}{2\sqrt2} & 0 \end{bmatrix} * \frac{1}{10^{6}} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 1+2\sqrt{2} \end{bmatrix}$

 $r_{E}= \begin{bmatrix} r_{2y}\\ r_{3x}\\ r_{3y}\\ r_{4x}\\ r_{4y} \end{bmatrix} =10^3 \begin{bmatrix} 0\\ 0\\ -1\\ -1\\ 1 \end{bmatrix}$

 The reactions at each Node (in N): C(X-component,Y-component) = (0,-1000) D(Y-component) = 0 F(X-component,Y-component) = (-1000, 1000)

## Comparing Results with ABAQUS

The Stresses and Displacements calculated have been compared with the results generated by ABAQUS and there is a very good match as seen below.

2D Truss Model in ABAQUS, after the Loads and BCs have been applied.

The above image displays the Stress in each bar.

The above image displays the X-Displacements of the structure.

The above image displays the Y-Displacements of the structure.

--Eml5526.s11.team5.vijay 19:56, 15 February 2011 (UTC)

# 4.7 Calculix Graphics Installation, Basic Examples, & Tutorials

## Given

 Calculix is an open source Finite Element code (with ABAQUS like input) available at http://www.dhondt.de/

## Find

 1) Install CGX

 3) Reproduce Basic examples

 4)Write a report, explaining to novices how to install and run cgx

## Solution

 1) Calculix Graphics Modules (CGX) was installed successfully.

 Installation procedure for a Windows OS is described below:

 Upon visiting http://www.dhondt.de/, you will be redirected to http://www.bconverged.com/download.php for a windows executable, where there is a link 'CalculiX_2_2_win_002.zip' which down loads the required ZIP files.

 Unzip the above downloaded file and double click on the executable to install calculix (provided by bconverged).

 After successful installation, You will have access to CalculiX CrunchiX (CCX 2.2), CalculiX GraphiX(CGX 2.2), Help file (which contains manuals for CCX, CGX), test cases and a custom built SciTE, a text editor which is integrated with the other tools.

 2) We have gone through the CGX User manual (available in the help file provided with bConverged or at [2] ) and signed up with the user group.

 The user group can be found at :http://groups.yahoo.com/group/calculix/

 3) We then reproduced the basic examples (disk, cylinder, sphere, sphere_volume, airfoil) with the help of example files downloaded from the following link: http://dhondt.de/cgx_2.2.exa.tar.bz2.

 The input files required to reproduce each of the following examples are shown below:

 a) Disc

 b) Cylinder

 c) Sphere

 d) Sphere_Volume

 e) Airfoil

 4) The installation procedure for Calculix Graphics (CGX) has been described in (1).

 There are two ways in which we can write CGX input files and run them:

 a) Using CalculiX Command window

 b) Using SciTE, a text editor which is integrated with other tools

 We will use SciTE to write and run input files. This process involves three steps:

 i) Open SciTE and create a New file:

 SciTE can be opened as follows: Start> All Programs > bConverged > SciTE A text editor pops up where an input file can be edited or created. A new file can be created by using the command "Ctrl+N" or by selecting "File > New"

 A text editor pops up where an input file can be edited or created. A new file can be created by using the command "Ctrl+N" or by selecting "File > New"

 ii) Type the required commands & save the file (in .fbd format): Once a new file has been created we are ready to start writing our input file. Type the commands for geometry creation and mesh generation in the text editor. Few commands required to generate basic shapes are discussed later on. After the required commands are typed in, Save the input file giving it an appropriate name and extension (.fbd) Example: "Sphere.fbd"

 Once a new file has been created we are ready to start writing our input file. Type the commands for geometry creation and mesh generation in the text editor. Few commands required to generate basic shapes are discussed later on. After the required commands are typed in, Save the input file giving it an appropriate name and extension (.fbd) Example: "Sphere.fbd"

 iii) Pre Process the input file (.fbd): The .fbd file created using the above steps can be PreProcessed by hitting the "F10" key or by selecting "Tools > PreProcess"

 The .fbd file created using the above steps can be PreProcessed by hitting the "F10" key or by selecting "Tools > PreProcess"

 Generally the Plot window which pops up upon PreProcessing hides the mesh and shows only the wire frame geometry. To View the meshed part do the following:

 Menu > Viewing > Show All Elements With Light

 Menu > Viewing > Toggle Element Edges

 (Note: Menu is displayed by clicking on the input file name below the Plot)

 CGX Commands used in creating the basic shapes:

 PNT

 This keyword is used to define or redefine a point.

 Syntax:

 'pnt' |'!' [ ]|[ ]|[ ]|[]

 Examples:

 pnt p1 11 1.2 34

 Creates point "p1" with coordinates (11,1.2,34)

 pnt ! 11 1.2 34

 The name is chosen automatically with coordinates (11,1.2,34).

 pnt ! L1 0.25 3

 To create 3 points on a line "L1" with 0.25 spacing

 pnt ! P1 P2 0.25 3

 To create 3 points between P1, P2 with spacing 0.25

 LINE

 This keyword is used to define or redefine a line. A line depends on points. A line can only be defined if the necessary points are already defined.

 Syntax:

 'line' |'!'

 Examples:

 line l1 p1 p2 4

 To create a line "l1" by the points p1 and p2 divided into 4 segments.

 line ! p1 p2 cp 4

 To create a arc line joining points "p1", "p2" with "cp" as center.

 GSUR

 This keyword is used to define or redefine a surface in the most basic way. Each surface must have three to five lines or combined lines (see lcmb) to be mesh-able. However, the recommend amount of edges is four.

 Syntax:

 'gsur' |'!' '+|-' 'BLEND|' '+|-' '+|-' -> .. (3-5 times)

 Example:

 gsur S004 + BLEND - L002 + L00E + L006 - L00C

 will create the surface S004 with a mathematically positive orientation indicated by the + sign after the surface name. The keyword BLEND indicates that the interior of the surface will be defined according to Coons or a NURBS surface is referenced. Use a + or - in front of the lines or lcmbs to indicate the orientation.

 GBOD

 This keyword is used to define or redefine a body in the most basic way. Each body must have five to seven surfaces to be mesh-able. However, the number of recommended surfaces is six. The first two surfaces should be the top and the bottom surfaces.

 Syntax:

 'gbod' |'!' 'NORM' '+|-' '+|-' ->.. ( 5-7 surfaces )

 Example:

 gbod B001 NORM - S001 + S002 - S005 - S004 - S003 - S006

 will create a body B001. The keyword NORM is a necessary placeholder for future functionallity but has no actual meaning. Next, follow the surfaces with a sign + or - in front that indicates the orientation of each surface.

 SETA

 This keyword is used to create or redefine a set. The following entities are known: Nodes n, Elements e, Faces f, Points p, Lines l, Surfaces s, Bodies b, Nurb Surfaces S, Nurb Lines L and names of other sets se

 Syntax:

 'seta' '!'|'n'|'e'|'p'|'l'|'c'|'s'|'b'|'S'|'L'|'se' | ['n'|'e' '-' ]

 Example:

 seta dummy p p1 p2

 This create a set "dummy" and will add the points p1 and p2 to the set.

 ELTY

 This keyword is used to assign a specific element type to a set of entities .

 The element name is composed of the following parts:

 I) The leading two letters define the shape (be: beam, tr: triangle, qu: quadrangle, he: hexahedra),

 II) The number of nodes

 III) At last an attribute describing the mathematical formulation or other features (u: unstructured mesh, r: reduced integration, i: incompatible modes, f: fluid element for ccx). If the element type is omitted, the assignment is deleted. If all parameters are omitted, the actual assignments are posted.

 Syntax:

 'elty' [] ['be2'|'be3'|'tr3'|'tr3u'|'tr6'|'qu4'|->'qu8'|'he8'|'he8f'||'he8i'|'he8r'|'he20'|'he20r']

 Examples:

 elty

 will print only the sets with assigned elements. Multiple definitions are possible. For example,

 elty all

 deletes all element definitions. If the geometry was already meshed, the mesh will NOT be deleted. If the mesh command is executed again after new assignments has taken place, additional elements will be created.

 elty all he20

 assigns 20 node brick-elements to all bodies in the set all.

 elty part1 he8

 redefines that definition for all bodies in the set part1.

 elty part2 qu8

 assigns 8 node shell elements to all surfaces in set part2.

 MESH

 This keyword is used to start the meshing of the model. before using the mesh command, the element types must be defined with the elty command

 Syntax:

 'mesh' ['fast'] ['block'|'lonly'|'nolength'|'noangle'|'length'|'angle']

 PLUS
 This keyword is used to display the entities of an additional set after a plot command was used.

 Syntax:

 'plus' ['n'|'e'|'f'|'p'|'l'|'s'|'b'|'S'|'L']&['a'|'d'|'p'|'q'|'v'] -> 'w'|'k'|'r'|'g'|'b'|'y'|'m'|'i'

 Now, Let us review how the above discussed commands are used to crate a Disk.

 i) Initially Points are created by typing the point name followed by their coordinates. A screen shot of the created points is shown below along with the commands

 ii) We now have all the points required to create the disk.

 Arc lines join points on circumference that are adjacent to each other (As shown below)

 Straight lines join the circumferential points to the center of the circle (as shown below)

 iii) After all the lines have been created, we create surfaces and this is what the geometry looks like before it is meshed.

 The Disk after it has been meshed is shown in (3)

# 4.6 Elastic bar with variable distributed spring

## Given

 An elastic bar is loaded with a variable distributed spring p(x) along its length as shown in the below figure. The distributed spring imposes an axial force on the bar in proportion to the displacement.

 L- Length of bar,

 A(x)- Cross sectional area at x,

 E(x)- Young's Modulus at x,

 b(x)- Body force at x

 Boundary conditions:

 i) Displacement at x(=0) = 0

 ii) Force at x(=L) = $\bold t$

## Find

 a) Construct the weak form

 b) Construct the strong form

## Solution

### a) Strong Form

 Consider the free body diagram of a differential element for the above problem.

 Balancing the forces acting on the differential element in the X-direction (along the length of the bar):

  $\quad \Sigma F = \sigma (x)A(x)n(x) + \sigma (x+dx)A(x+dx)n(x+dx) + b(x)dx - p(x)dx = 0$  $\displaystyle (Eq 4.6.1)$

 p(x)dx is the force exerted by the spring on the differential element.

 b(x)dx is the body force acting on the differential element.

 Substituting n(x)=1, n(x+dx)=-1 and expanding the terms $\quad \sigma (x+dx), A(x+dx)$ Using Taylor's Series and neglecting Higher order terms, Eq 4.6.1 transforms into:

 $\frac {\partial \sigma}{\partial x}dx + b(x)dx - p(x)dx = 0$ $\displaystyle (Eq 4.6.2)$

 We know that $\sigma(x) = E(x)A(x)\frac {\partial u}{\partial x}$

 Substituting the above relation and cancelling the dx term, Eq 4.6.2 yields the Strong Form:

 $P(u) = \frac {\partial }{\partial x} \left[ E(x)A(x) \frac {\partial u}{\partial x} \right] + b(x) - p(x) = 0$ $\displaystyle (Eq 4.6.3)$

 The Boundary Conditions for the Strong form are:

 i) Essential Boundary Condition:

 $\quad u(x=0) = 0$

 ii) Natural Boundary Condition:

 $E(x)\frac {\partial u(x=L)}{\partial x} = \bold t$

 where $\bold t$ is the Traction applied.

### b) Weak Form

 The weak form by simplifying the following equation

  $\quad \int_{\Omega} W(x)P(u)dx=0$ 

 i.e., $\int_0^L W(x) \left( \frac {\partial }{\partial x} \left[ E(x)A(x) \frac {\partial u}{\partial x} \right] + b(x) - p(x) \right)dx = 0$

  $\underbrace{ \int_0^L W(x) \left( \frac {\partial }{\partial x} \left[ E(x)A(x) \frac {\partial u}{\partial x} \right] \right)dx}_{I} + \underbrace{\int_0^L W(x) \left(b(x) - p(x) \right)dx}_{II} = 0$  $\displaystyle (Eq 4.6.4)$

 Term I can be expanded using the chain rule of integration as follows:

 $\left[ W(x)E(x)A(x) \frac {\partial u}{\partial x} \right]_0^L - \int_0^L \left( \frac {\partial W(x)}{\partial x} E(x) A(x) \frac {\partial u(x)}{\partial x} \right)dx$

 $\left[ W(x)E(x)A(x) \frac {\partial u}{\partial x} \right]_L - \left[ W(x)E(x)A(x) \frac {\partial u}{\partial x} \right]_0 - \int_0^L \left( \frac {\partial W(x)}{\partial x} E(x) A(x) \frac {\partial u(x)}{\partial x} \right)dx$

 Applying the Natural Boundary Condition $E(x)\frac {\partial u(x=L)}{\partial x} = \bold t$ and selecting W Such that W(x=0)=0, the above expression becomes:

 $\left[ W(L)A(L) \bold t \right] - \int_0^L \left( \frac {\partial W(x)}{\partial x} E(x) A(x) \frac {\partial u(x)}{\partial x} \right)dx$

 Substituting the above expression back into Eq 4.6.4 We get the weak form:

 $\left[ W(L)A(L) \bold t \right] - \int_0^L \left( \frac {\partial W(x)}{\partial x} E(x) A(x) \frac {\partial u(x)}{\partial x} \right)dx + \int_0^L W(x) \left(b(x) - p(x) \right)dx = 0$

 The Weak Form can be rewritten as follows:

 $\int_0^L \left( \frac {\partial W(x)}{\partial x} E(x) A(x) \frac {\partial u(x)}{\partial x} \right)dx = \left[ W(L)A(L) \bold t \right] + \int_0^L W(x) \left(b(x) - p(x) \right)dx, \forall W,$ Such That W(x=0)=0

 Essential Boundary Condition: u(x=0)=0,

 Natural Boundary Condition is Already absorbed into the weak form.