Loading...

PPT – NUMERICAL ANALYSIS of PROCESSES PowerPoint presentation | free to download - id: 649152-ODNhY

The Adobe Flash plugin is needed to view this content

NUMERICAL ANALYSIS of PROCESSES

NAP5

Models described by ordinary differential

equations boundary problems Method of weighted

residuals and finite elements Heat exchangers

(making use analytical solution) Pipelines

(pressure and flowrates using finite element

methods) Trusses, beams and rotationally

symmetric shells

Rudolf Žitný, Ústav procesní a zpracovatelské

techniky CVUT FS 2010

Models by ODE boundary problem

NAP5

Some systems are described by the same ODE as

before (ie as systems of equations with first

derivatives), except that the initial state is

not completely defined because not all conditions

are specified for the same value of the

independent variable (we speak not about initial,

but boundary conditions). While the initial

problem generally refers to non-stationary

processes (the independent variable is time)

finding a steady state is typical for boundary

problems when independent variable is a spatial

coordinate. Examples

Temperature profiles in sheltube and plate heat

exchangers

Pressures and flowrates in branched pipelines

(steady state)

Static deformation of loaded beam structures or

rotational thin wall shells

ODE heat exchangers

NAP5

ODE heat exchangers

NAP5

The simplest possible case of hat exchanger with

two parallel streams (coordinate x is measured

from one end of HE). This is example solved in

the course of Heat processes. Enthalpy balance

for CO-CURRENT

trhis is INITIAL problem, can be solved by

previous methods

T1

x

T2

L

Enthalpy balance for COUNTER-CURRENT

BOUNDARY problem, missing temperature at x0

T1

x

T2

(k is heat transfer coefficient related to unit

length of channel, Wi heat capacity of streams)

ODE heat exchangers

NAP5

Equations for temperature profiles in a matrix

form are as follows

special case for co-current flow

The matrix A is constant (for constant k),

the ODE system is therefore linear and can be

transformed from a coupled to a separated system

by linear transformation

where U is the eigenvector matrix of A

satisfying

a matrix multiplied by its eigenvector gives the

same eigenvector (scaled by eigenvalue)

is a diagonal matrix of eigenvalues ?1 ?2 (how

many equations so many eigenvalues and

eigenvectors). Calculation of ? and U is

accomplished in MATLAB by single standard

function U,Leig(A).

ODE heat exchangers

NAP5

Substituting transformation to the ODE system

Square DIAGONAL matrix (as soon as the U are

eigenvectors)

results to uncoupled system of ODE for

transformed temperatures (vector Z )

ODE for Z are independent, because ? is

diagonal

Solution of ODE

Backtransformation

coefficients dj are determined by boundary

conditions

ODE heat exchangers

NAP5

Specific case of 2channels and co-current

arrangement of HE

Eigenproblem of A2 x 2 can be solved

analytically

?

Boundary conditions for x0

and solution

.

the result is the same, only the order of vectors

is exchanged, and eigenvectors are normalised to

the unit Eucleid norm. It has no effect upon the

results.

verify eigenproblem solution in MATLABu (k1,

w11, w22) a-1 10.5 -0.5 u,leig(a)

u -0.8944 -0.7071 0.4472 -0.7071 l

-1.5000 0 0 0

ODE heat exchangers

NAP5

Generalisation for more complicated HE and HEN

(Heat Exchanger Networks)

Heat exchangers have often more than 2 streams

(generally N-streams) and each stream (channel)

is divided into M sub-channels (eg because the

inlet and outlet ports of various streams are in

different locations or heat transfer coefficient

changes along the flow direction). Let us assume

that all N streams is divided into M sub-channels

(Mgt N), the position of the input subchannels is

given by the vector x ' position and outputs

the vector x''. Subchannel i (i 1,2, ..., M)

corresponds to the temperature Ti(x), heat

capacity Wi (constant), and the heat transfer

coefficient kij with all other sub-channels (also

constant). This data is sufficient to write an

enthalpy balance of subchannels M and to

formulate system of M-differential equations

This system is solved by previous procedure

(finding the eigenvalues and eigenvectors of

A and solving separated equations). The

result is a vector of temperature profiles in M

subchannels, expressed by M-coefficients of yet

unknown vector d

From this equation we can determine the input and

output temperatures of M subchannels

(substituting the coordinates x ', x''). Some

of the inputs/outputs of subchannels are

input/output channels (ie the whole heat

exchanger), some express only and internal

linking. These links of subchannels are

represented by binary matrix GMxM whose rows

correspond to the inputs and outputs of the

column (Gij1 means that the input of i-th

sub-channel is the output of j-th sub-channel),

by the matrix of inputs G 'MxN (G'ij1 means

that the input of i-th subchannel is the input

of j-th stream) and by the matrix of outputs

G''NxM (G''ij1 means that the output of

i-th stream is and the output of j-th

sub-channel).

ODE heat exchangers

NAP5

N-coefficients of d is determined from the

boundary conditions (the prescribed inlet

temperature at N streams). The requirement that

the outlet temperature of subchannels are

simultaneously input temperature subchannels

downstream (it is defined by matrix G),

represents the missing M-N conditions. Inlet

temperature of all subchannels are therefore

either the heat exchanger inlet temperature or

outlet temperature of the downstream sub-channel,

ie

Previous relationship is a system of M linear

algebraic equations for the coefficients d and

for example the output temperatures of the whole

heat exchanger can be expressed as follows

Xing Luo, Meiling Li, Wilfried Roetzel A general

solution for one-dimensional multistream heat

exchangers and their networks. International

Journal of Heat and Mass Transfer 45 (2002)

26952705

ODE heat exchangers

NAP5

The previous method is appropriate when the ODE

are linear and if there exists an analytical

solution. If the coefficients of heat transfer

are not constant (but are independent of the

temperature) the exchanger channels can be

divided into smaller subchannels with constant

mean values k (see above). If the coefficient k

is temperature dependent, it would be necessary

to repeat the process iteratively (and in each

iteration to solve their own eigenproblem for the

matrix A depending on temperature). It is

alternatively possible to use numerical

integration (repeated solution of the initial

problem) Missing initial condition (x0) must be

estimated and gradually refined in repeated

integration (eg by Runge Kutta) so as to best

meet the remaining boundary conditions (for xL).

This is called the shooting method and it's

actually an optimization problem of minimization

the norm of the difference between the prescribed

boundary conditions and numerical prediction for

xL (in MATLAB the function fminsearch optimizes

the parameters of missing initial

conditions). You can try a different approach of

estimating the missing initial conditions at x0,

execution of integration from x0 to xL, and

then carry out the solution in the reverse

direction from xL to x0, but with the predicted

values (xL) replaced with prescribed boundary

terms. The whole procedure must be repeated

several times (zig zag method) and not always

converge (see next example).

ODE heat exchangers

NAP5

T190

Example of two channels HE of length L1,

cocurrent and countercurrent arrangement solved

by Runge Kutta in MATLABu (ode45)

definition of ODE for counter current (differs

only by sign)

definition ODE for cocurrent

function dtdxdts(x,t) dtdxzeros(2,1) w11

w22 k1 dtdx(1)-k/w1(t(1)-t(2)) dtdx(2)k/w

2(t(1)-t(2))

function dtdxdtp(x,t) dtdxzeros(2,1) w11

w22 k1 dtdx(1)-k/w1(t(1)-t(2)) dtdx(2)-k/w

2(t(1)-t(2))

T240

cocurrent arrangement is INITIAL problem

(therefore easy)

gtgt solode45(_at_dts,0 1,90 40) gtgt

plot(sol.x,sol.y)

T1

countrercurrent by the zig-zag method

sode45(_at_dtp,0,1,90,50) for

i110 nlength(s.x) sode45(_at_dtp,1,0,s.y(1,n)

,40) nlength(s.x) sode45(_at_dtp,0,1,90,s.y(2

,n)) end plot(s.x,s.y)

T2

ODE dryers, extractors

NAP5

Flow-through (continuous) devices such as

contactors, reactors or dryers (co-current,

counter) are solved almost as well as heat

exchangers. Only the balancing equations of heat

transfer are added to equation analogous to the

mass transfer. So, again the numerical

integration methods (Runge Kutta), method of

shooting, etc. are used Example Bruce D.M.,

Giner S.A. Mathematicall Modelling of Grain

Drying in Counter-flow Beds. J.Agric.Engng.Res.

(1993),pp.143-161 Solution by Euler, Runge Kutta,

shooting method.

Pipelines networks

NAP5

The piping network is made up of straight

sections, bends and turns. In a broader sense

also by machines like pumps, accumulators,

valves. The aim is to determine the pressure

profile along the pipe and flow distribution in

branches.

Different design procedure must be applied in the

case of

- nonstationary flow (effect of acceleration)
- compressible fluid or deformable walls (hydraulic

shock)

Pipeline network

NAP5

The basic element is a tube (of any section which

may vary in the axial direction). Flow

characteristics in a certain place of the tube

are the pressure p(x) and the characteristic

velocity u. Axial velocity profile is fairly

generally (even for compressible flow) described

by the Navier Stokes equation

The second term of the right side expresses the

frictional losses which are a function of flow

rate, viscosity and the geometry of the channel.

These factors will be included into the

coefficient k

D equivalent diameter, A cross section, m mass

flow-rate

dArcy Weissbach friction coefficient, in laminar

e.g. 64/Re

Friction loss

Integrating this equation in the pipe is obtained

Bernoulli equation for unsteady compressible

fluid flow (this equation will be uses later ...).

Pipeline network FEM, MWR

NAP5

Finite Element Method Method of

Weighted Residuals

Consider a special case a steady flow-rate,

constant cross section. Derive the previous

equation

There is only the pressure p(x) in this equation,

and the mass flow rate (which is constant) was

eliminated. The left side represents the friction

loss, the right side buoyancy (when the density

depends on the temperature and the temperature

along the channel change). The coefficient k is

constant in the laminar flow regime and decreases

with increasing Reynolds number in the turbulent

regime.

This ordinary differential equation will be used

for explanation of principles of the method of

weighted residuals (MWR) and its special variant

of finite element method (FEM). What is the

residue of differential equations? What remains

to the right side, when all terms will transfer

to the left

If p(x) is exact solution the residual function

res(x) is zero

Pipeline network FEM, MWR

NAP5

Finite Element Method Method of

Weighted Residuals

The method of weighted residuals collected from

the possible options p(x) such that nullifies

integrals of residuals multiplied by the selected

weight (test) function w(x)

The solution is approximated by a linear model, a

linear combination of basis functions Nj(x)

(which may even be normalised polynomials defined

in elements)

N unknown parameters pj requires at least N

equations. Therefore, it is necessary to select

at least N different weighting function wi(x).

The result is a system of N linear algebraic

equations for N unknowns

Method of Weighted Residuals

NAP5

- The weighting function wi(x,y,z) are proposed in

advance, more or less independently to solutions.

They can be of different types and correspond to

the different numerical methods - Spectral methods (analytical weight functions

wi(x)) - Finite Element Methods (continuous weight

functions) - Finite Volume Methods (discontinuous, stepwise

functions) - Collocation Methods (zero residuum in nodes.

Weights are delta functions)

(orthogonal polynomials or goniometric functions

are used as base and weight functions. Analytical

weights are used also in the Boundary Element

Methods BEM.)

(example is Galerkins method using the same base

and weight functions. FEM is preffered in

structural analysis, see software ANSYS)

(also Control Volume Method. Most frequently used

in fluid mechanics, see software FLUENT)

xi

(also Finite Difference Method. Suitable for

complicated problems in simple geometries, easy

programming)

Finite Elements and Volumes

NAP5

System (e.g.pipes) is divided into finite

elements with nodal points on the element

boundaries (element is a section of pipe between

two nodes). Each node (i) is assigned one

weighting function wi(x), one basis function

Ni(x) and the one calculated value (pressure pi).

Solution inside a finite element is defined by

interpolation of nodal values at boundary. At

finite volumes there is only one (internal) node

associated with the finite volume and therefore

the solution is approximated by a constant.

Remark A prototype of FVM are compartment

models (aim was calculation of temperatures and

concentrations the same in the whole

compartment). FEM aims to calculate values at

element boundaries, i.e. parameters of streams

connecting compartments. FVM satisfies the

balances of control volume (it follows from

principle of control volume), while FEM need not

satisfy the conservation exactly (and this is

disadvantage of FEM).

Finite Element Method

NAP5

Basis functions Ni(x) are generally constructed

so that they are equal to 1 at node i and zero at

all other nodes (1, ..., N). Basis function Ni(x)

is a non-zero only in the element that contains

the node number i It is possible to use identical

base and weight functions (for example

element-wise linear polynomial pressure and

weights). Identification of weight and base

functions (wi(x)Ni(x)) is a characteristics of

Galerkins method. Integrand of the residual

integral is the product of weight function and

the second derivative of base function. Because

the weight function in FEM is continuous its

first derivative exists and it is possible to

transform the integral by per partes method

Matrix coefficients Kij and right side vector bi

are integrals over the entire pipe network (L is

the total length of all branches). These

integrals are calculated as the sum of integrals

over individual finite elements.

Finite Element Method

NAP5

In a general finite element only two base and

weighting functions corresponding to two nodes i,

j are nonzero. The element therefore contributes

only to 4 entries of the matrix (Kii, Kjj Kij

Kji). Because the derivatives of linear Ni(x) are

constant the integral contribution of element is

calculated easily as

Assembly algorithm element matrix addition into

the global system matrix (and right side vector).

In the cycle through all elements a procedure

calculating the local matrix for a given flow

rate and geometry is called and the matrix is

added to the global array. Correspondence between

the local indices (1,2) and the global indexes

must be defined in the matrix connectivity c

whose rows correspond to the elements and in the

two columns are indices i, j nodes of each

element. This is actually defined by the network

topology.

subroutine calculating Ke

Kzeros(n,n) bzeros(n,1) for e1ne

Ke,belocal(e) for i12

igc(e,i) b(ig)b(ig)be(i) for

j12 jgc(e,j)

K(ig,jg)K(ig,jg)Ke(i,j) end end end

connectivity matrix

Finite Element Method

NAP5

Physical interpretation Product of local matrix

of one element Ke by the vector of calculated

pressures results to the mass flowrate through

the element (e), more specifically flowrate from

node i to the node j

Assembly of local matrices represents requirement

that the sum of oriented flowrates through any

node is zero.

- The assembled matrix K is singular in this

phase, as soon as the boundary conditions are not

specified (and distribution of pressures is not

unique). The boundary conditions have to be

prescribed at all end-nodes of mesh - End pressures (the correponding row of K is

substituted by 1 at diagonal and the right hand

side is prescribed pressure). - As soon as the mass flowrate is prescribed at an

end-node the flowrate is substituted as the right

hand side bi.

Finite Element Method

NAP5

System of linear algebraic equations can be

solved either by direct methods (characterized by

a finite number of steps, example is Gauss

elimination method, frontal method ) or iteration

methods (e.g. Gauss Seidel overrelaxation, or

method of conjugated gradients). The system

matrix K is sparsed (in a row corresponding

to a node connecting only 2 elements are only 3

nonzero entries) and therefore special storage

methods are used for sparse matrices, e.g.

Profit is not only a memory saving but first of

all speed of computation. For example the Gauss

elimination with a full and asymetric matrix N x

N requires N3/3 floating point multiplication.

Using the band matrix algorithm the number of

operations is reduced to N.Nb2/2, where Nb is the

band width.

For small problems (for example pipe networks) it

is sufficient to use a standard Gauss eliminatiom

and to write the solution of Kpb in

MATLAB as pK\b

Finite Element Method

NAP5

Matrix of equations for nodal values of pressures

is a constant only at laminar.flows. For

turbulent or non-Newtonian flows the matrix K

is a function of Reynolds number and calculated

pressure drop (k(p)). The system is non-linear

and must be solved iteratively where k is index

of iteration (and a new matrix K must be

assembled at each iteration). Systém of

nonlinear equations can be used by previously

introduced substitution method (Picard) or by

using Newton Raphson method, that is a

straightforward generalisation of Newtons

tangent-line method. The method is based upon

Taylor expansion of residuals.

residuum ?

residuum / ?p

Finite Element Method

NAP5

Historical remark Galerkin method and FEM

In russsian literature the Galerkin method is

referred to Bubnov method published in 1914

(Bubnov I.G. Stroitelnaja mechanika korablja.

II. SPB-Tipografija morskogo ministerstva

Bubnov was a designer of battleships). Galerkin

published his works later Galerkin B.G. K voprosu

ob issledovajii naprjazenij i deformacij

v uprugom izotropnom tele, Doklady Akademii Nauk

SSSR, 1930, Ser.A, No.14, s.353-358. FEM

application of Galerkin method was developed

much later in connection to analysis of

catastrophes of british jet-liners Comet. Turner

M.J. et al Stiffness and deflection analysis of

complex structures. J. of Aeronautical Sciences,

23, pp.805-823 (1956).

Pipe networks circulation model

NAP5

Fluid flow in a model of blood circulation

(physical model in the laboratory of

cardiovascular mechanics). The model demonstrates

that the FEM need not be designed by Galerkin or

even by weighted residual methods (in this case

only Bernoulli and continuity equations are used).

Pipeline network - nodes

NAP5

Nodal parameter are pressure p and the volumetric

flowrate Q. It means that nodes cannot be located

at the branching points, where the flowrate is

not defined!

Y m

Number of unknowns 2N-Nb, where N is number of

nodes and Nb is number of end-points (where

either flowrate or pressure is prescribed as a

boundary condition)

19

18

16

28

15

17

27

14

22

20

29

7

9

12

24

25

21

6

8

5

26

10

It is not possible to define pressure and

flowrate at a node simultaneously!!!

11

13

4

3

23

Prescribed pressure p, calculated Q

n

2

Calculated pressure and flowrate

m

Auxilliary node defining geometry of branching

k

1

X m

-0.4 -0.3 -0.2 -0.1 0 0.1 0.2

0.3 0.4 0.5 0.6 0.7 0.8

Pipeline network - elements

NAP5

Two different FE two-nodes element generate 2

equations (Bernoullis equation along streamline

between nodes 1,2 and continuity equation) and a

three-nodes FE generating 3 equations (Bernoulli

1,2 and 1,3 plus continuity equation). Equations

of system are not generated by assembly of

elements contributions. Equations are generated

directly at an element level. Number of unknowns

agrees with the number of generated equations

2N-Nb 2Mp3Mv

Number of elements BRANCH

Number of elements PIPE

Number of nodes

Number of end-nodes

Pipeline network - Bernoulli

NAP5

Previous solution is based upon Bernoulli

equation derived by integration of momentum

balance along streamline between nodes 1 and 2

The last term ez12 , energy loss J/kg, depends

on flowrate Q, apparent viscosity, local losses,

seez E.Fried, E.I.Idelcik Flow resistance, a

design guide for engineers, E.I.Idelcik Handbook

of hydraulic resistances.

Pipeline network MATLAB

NAP5

MATLAB possible solution mainpq.m

Input data files xyz.txt x y z

(coordinates of nodes) cb.txt i pq

(boundary conditions i-node,

pressure, -flowrate as pq) cp.txt i j d

(connectivity matrix, diameter

of pipe) cv.txt i j k l d1 d2 d3 ?

(connectivity for branching, diameter of pipes,

branching angle)

Truss, beam, shells

NAP5

In this chapter you will se similarity between

calculation of pressures and flowrates in

pipelines and calculation of forces and

sisplacements in a frame structure.

Truss, beam, shells FEM

NAP5

What is the difference between trusses and beams?

Both are rods but trusses are connected by joints

(no moment can be transfered) while beams are

welded and are deformed by bending and torsional

moments.

The same differences characterise membranes and

shells. Membrane is an analogy of truss (no

bending moments, unlike shells).

Bernoulli beam deformed cross section remains

perpendicular to the beam axis Timoshenko beam

deformed cross section remains planar but not

perpendicular to axis (effect of shear) Analogy

in shells Kirchhoff plate cross section

perpendicular to the deformed mid plane 9analogy

to Bernoulli beam) Mindlin plate is analogy to

Timoshenko beam

Truss, beam, shells FEM

NAP5

Following examples share the following

features Solution is based upon deformation

method, calculating displacements and rotations

at nodes (stresses are calculated ex post).

Deformation method corresponds to Lagrangean

principle of minimum of total potentiual energy.

There exist another principles based upon

complementary energy (Castiggliano) when primary

results are forces or hybrid principles (Hu

Waschizu) calculating both kinematic and dynamic

varibles as primary unknowns.. Basic result of

each example are element matrices (local

matrices) implemented as MATLAB functions. You

can copy these functions and use them for

solution of particular problems.

Truss FEM

NAP5

Frame construction from trusses is

computationally identical with the calculation of

pressures in pipeline networks. Instead of

pressures are displacements of nodes, the

stiffness coefficient (cross-section x Young

modulus) replaces the friction factor and

buoyancy is substituted by continuous load.

Weighted residual method results to algebraic

equations for displacements

Remark FUNCTIONALISTs would derived the same

result in a simpler way by minimising potential

energy

and also the stiffness matrix is identical

Trusses FEM

NAP5

In this way it could be possible to solve only

the simplest 1D problem of serially connected

springs in the direction x. Nodal displacements

are generally vectors ux,uy (2D case), uz (3D

case). Local matrices must be transformed to a

global coordinate system shared by all finite

elements. For example the local element matrix 2

x 2 is to be transformed to 4 x 4 matrix in

planar construction (each node is characterized

by two displacements ux, uy). Transform

from local to global coordinatge system concerns

only rotations (one rotation at 2D and three

rotations at 3D frames).

Local coordinate system (axis ?) Global system

x,y (ccos?, ssin?)

Trusses FEM

NAP5

Equations in local coordinate system

After transformation

Final matrix in thee coordinate systém x,y

gtgt tc2 cscs s2 gtgt ket -t-t t

Truss

NAP5

Example

x0 1 0 1 2 y0 0 1 1 2 con1 3 2 3 2 4

3 4 3 5 4 5 a4e-4 4e-4 4e-4 4e-4 4e-4

4e-4 nulength(x) nelength(a) n2nu kzeros(n

,n)bzeros(n,1) for e1ne

kl,igkloc(e,x,y,con,a) k(ig(14),ig(14))k

(ig(14),ig(14))kl(14,14) end for i14

k(i,)0k(i,i)1 end b(n)-100 pk\b

assembly

nulové posuvy v uzlech 1,2 a zatížení -100N v

uzlu 5 (ošklivé rešení, platí jen pro danou

topologii)

Connectivity matrix

Gauss elimination. Vector p are displacements

function kl,igkloc(ie,x,y,con,a) E200e9 aea(

ie) i1con(ie,1)i2con(ie,2) ig2i1-1 2i1

2i2-1 2i2 le((x(i1)-x(i2))2(y(i1)-y(i2))2)

0.5 c(x(i2)-x(i1))/les(y(i2)-y(i1))/le tdE

ae/lec2 cscs s2 kltd -td-td td

Stiffness matrix of element

Beams (pipes)

NAP5

Displacement of neutral axis of a straight beam

is described by differential equation of the

fourth order

where I is inertia moment of cross section, Fx is

axial force, k is stiffness of support and fz(x)

is continuous transversal load. Corresponding

residual integral is

Per partes integration is applied twice with the

aim to decrease order of derivatives

The terms at boundary points resulting from per

partes integration are neglected. This is

consistent with the assumption that strong

boundary conditions (fixed deformations) are

applied or that there is no loading at boundary

(free ends).

Beams (pipes)

NAP5

If you are not interested in details, skip these

gray pages

The Bernoulli beam is described by ODE, relating

bending moment M and resistive moment

Where the second derivative is a measure of

curvature (1/d2u/dx2 radius of circular arc).

Bending moment can be caused by transversal

force (reaction MFz(b-x)), but also by axial

force Fx

This equation is a basis for stability analysis

of axially loaded beams (Euler stability).

Bending moment M(x) corresponding to isolated

transversal forces is only a linear function of

x-coordinate and its second derivative is zero.

Only a continuous transversal loading or a

continuous elastic support give moments with

non-zero second derivative

Derivative of integral with variable upper bound

Beams (pipes)

NAP5

Lets return back to the weighted residual method

Twice applied per partes to the first term and

once applied per partes to the second term give

New terms corresponding to boundary conditions at

ends of beam are expressed in terms of moments

and forces

In previous derivation of Galerkin approach these

terms were neglected.

Beams (pipes)

NAP5

Strong and weak boundary conditions

Degrees of freedom (DOF) of the systém taken as a

stiff body should be eliminated by STRONG

boundary conditions (nodal displacements and

rotations). Strong BC are implemented easily by

zeroing corresponding row of the system matrix

with the exception of diagonal element that is

substituted by one and by specifying prescribed

displacement as the right hand side. If strong BC

are satisfied at all approximations than the

limit of approximations satisfies the BC too.

However, some boundary conditions cannot be

specified as strong but WEAK (conditions

containing first derivatives in case of second

order differential equations). But it is

questionable whether they will comply with the

prescribed conditions and limitations of

numerical solutions for refining the network

elements. Not always, at least not when the

differential equation simply does not prescribe

any parameter as a strong boundary condition. An

example is the equation of the beam entered the

second derivative of the deflection lines

M

Cantilever beam loaded constant torque M. Linear

basis functions. Nodes have a single parameter,

the displacement uz. Zero rotation at the fixed

(x 0) or not possible to prescribe.

The same task, but using the cubic basis

functions. Nodal parameters are displacement and

rotation. Zero rotation can be prescribed, but

the solution with a large number of elements

ignores this condition

Exact solution

In the solution shown in the figures this term is

neglected, we assume that Natural boundary

conditions hold. Is it not the cause of failure?

20 elements

1 element

2 elements

3 elements

3 elementy

20 elementu

Exact solution

Beams (pipes)

NAP5

Neglecting the term

assumes zero duz/dx at the end-point L. The

equation

of weighted residuals method is really trying to

accomplish (this is evident from the previous

figures, where an approximation of the deflection

line by linear and cubic polynomials respect this

unwanted boundary condition). When we include the

boundary term into the stiffness matrix, the

result will be improved (for cubic polynomials

solutions will be even accurate), but the

boundary condition of zero rotation is not

satisfied when using linear polynomials. It is

due to the fact that for second order

differential equations only displacements can be

applied as the strong boundary condition. The

first derivative of displacement (rotation) can

be prescribed as a strong boundary condition only

for the fourth-order equations.

Rešení uvažující clen

Linear base functions

20 elements

3 elements

Boundary conditions, which a priori ensures each

approximation solutions are called strong, while

the boundary conditions transferred to the

integral formula MWR application integration by

parts are weak. For differential equations 2s-th

order weak boundary conditions in a least

s-normal derivatives, boundary conditions with

lower derivatives are strong. For second order

differential equations are strong boundary

conditions prescribed values and weak conditions

include the first derivative, while the fourth

order equation (like bending plates) are strong

boundary condition values and first derivative

(displacements and rotations, ie kinematic

conditions), poor conditions concern until the

second or third derivative (stressful conditions,

the second derivative of corresponding moments M

EId2u / dx2, the third derivative of a rolling

force F -EI d3u / dx3). Compliance with weak

boundary conditions is a matter of correct

choices integral formula for the method of

weighted residuals and can not enforce them by

fixing the nodal parameters.

Beams (pipes)

NAP5

The final formulation of the Galerkin method of

weighted residuals for Bernoullis beam on

elastic foundation

Remark. Do not forget that this is still only

very simplified model of the beam, functioning

only for small deformations, slender beams (not

considering shear forces as Timošenko model) and

does not even have axial stiffness. The model

calculates only with a bending stiffness!

Beams (pipes)

NAP5

Significant difference from the problem with

trusses is that there are now the second

derivatives of basis functions in the integrand

and therefore it is not possible to use linear

polynomials - basis functions must be constructed

so as to have a continuous first derivative, for

example the Hermite cubic polynomials

Expression of deflection line

Substituting Hermite polynomials into integrand

results to a system of linear algebraic equations

whose matrix is the sum of three matrices

KM bending stiffness, KF geometrical

stiffness, M is mass matrix

Beams (pipes)

NAP5

Results for element of length L

Nodal parameters (transversal displacement and

rotation)

Stiffness matrix

Geometry stiffness

Mass matrix

Beams (pipes) MATLAB

NAP5

function kstif,kgeom,mass,igkloc(ie,x,y,con,ei,

k) EIei(ie)E200e9 i1con(ie,1)i2con(ie,2) i

g2i1-1 2i1 2i2-1 2i2 l((x(i1)-x(i2))2(y

(i1)-y(i2))2)0.5 kstifEIE/l312 6l -12

6l6l 4l2 -6l 2l2 -12 -6l 12

-6l6l 2l2 -6l 4l2 kgeom1/(30l)36

3l -36 3l3l 4l2 -3l -l2 -36

-3l 36 -3l3l -l2 -3l 4l2 masskl/4201

56 22l 54 -13l22l 4l2 13l -3l2

54 13l 156 -22l-13l -3l2 -22l 4l2

These matrices can be used to construct several

elements, but only oriented along the axis x

(mass matrix will be use when analyzing vibration

of beams, geometric stiffness matrix for the

investigation of stability). In each node are two

calculated parameters, deflection and bending

angle of deflection. We can not use the

coordinate transformation for rotating the truss

beams as well as for rods, because the cantilever

element has zero stiffness in the x direction.

The addition of stiffness matrix of the beam

stiffness of the element type of rod is shown in

the following example.

Beams frame 1/3

NAP5

Stiffness matrix of steel beam truss

function kstif,igklocp(ie,x,y,con,ei,ea) EIei(

ie)EAea(ie) E200e9 i1con(ie,1)i2con(ie,2)

ig3i1-2 3i1-1 3i1 3i2-2 3i2-1

3i2 l((x(i1)-x(i2))2(y(i1)-y(i2))2)0.5 c

(x(i2)-x(i1))/ls(y(i2)-y(i1))/l rc s 0-s c

00 0 1zzeros(3,3) Qr zz

r kmEIE/l3EAl2/EI 0 0 -EAl2/EI

0 0 0 12

6l 0 -12 6l

0 6l 4l2 0 -6l

2l2 -EAl2/EI 0 0

EAl2/EI 0 0 0

-12 -6l 0 12 -6l

0 6l 2l2 0

-6l 4l2 kstifQ'kmQ

Previous matrix KM extended by first and

fourth rows and columns (inserting stiffness

matrix of truss)

Transformation to global coordinate system

Transformation stiffness matrix from local to

global coordinates KQKeQ

Generalised displacements

Generalised loads

Beams frame 2/3

NAP5

Nodal coordinates, connectiv ity matrix,

assembly, boundary values spacification and

solution of systém is the same as with trusses

Inertia moment for pipe

Square cross section 1x1 cm. Area 1e-4 m2,

inertia moment I8.3333e-10 m4

Inertia moment for rectangular profile

x0 0 1 y0 1 1 con1 2 2

3 ei8.3333e-10 8.3333e-10ea1e-4

1e-4 nulength(x) nelength(ea) n3nu kzeros(

n,n)bzeros(n,1) for e1ne

kl,igklocp(e,x,y,con,ei,ea)

k(ig(16),ig(16))k(ig(16),ig(16))kl(16,16)

end for i13 k(i,)0k(i,i)1 end b(n-1)-

100 pk\b

Fy-100 N

Zero displacement and rotation in node 1 and load

-100N at node 3

Beams frame 3/3

NAP5

Result is vector of nodal parameters p(19),

generalised displacements

p -0.0000 ux1 0.0000 uy1

0 ?1 0.3000 ux2 -0.0000

uy2 -0.6000 ?2 0.3000 ux3

-0.8000 uy3 -0.9000 ?3

Internal forces and moments are calculated as the

product of stiffness matrix of particular element

times displacement vector (in this element). For

example element 2

fklp(ig(16)) f 0 Fx2

100.0000 Fy2 100.0000 M2 arm

of force 1m, force 100 N 0

Fx3 -100.0000 Fy3 0.0000 M3

kl local stiffness matrix ig selects 6

corresponding displacements from the vector of

global displacements p.

Rotationally symmetric shells

NAP5

In rotationally symmetric vessels loaded eg. by

internal pressure exist normal membrane stresses

(formulas of type pr / s), but at the locations,

where the curvature or thickness of wall changes,

bending and shear stresses appears (wave of

stresses in range typically ?(rs) ). The vessel

is then to be understood as a shell. FEM systems,

e.g. ANSYS, COSMOS, ... always contain spatial

shell elements (curved triangles or quadrangles),

but for application to rotationally symmetric

shells are sufficient and important only 2-nodal

elements (representing a conical ring). Perhaps

the most successful element of this type already

suggested about 30 years ago ing.Vykutil (VUT

Brno), is implemented in Femina and is easy to

write in MATLAB

See Schneider P., Vykutil J.Aplikovaná metoda

konecných prvku, skripta VÚT Brno, 1997.

ukázky grafických výstupu programu FEMINA

(Vykutilovský element)

Rotationally symmetric shells

NAP5

Element representing the conical ring has two

nodes and each node 3 parameters displacement in

the axial direction ux, uy radial direction and

rotation ?.

r

Stiffness matrix has dimension 6x6

Derivatives of base functions

Elastic constants

Rotationally symmetric shells

NAP5

MATLAB

function kstif,bp,igkloc(ie,x,y,con,he,pe)

kstif(6x6) matice tuhosti, bp(6) zatížení tlakem,

ig(6) indexy globální soustavy ie index

elementu,x(),y() globalni souradnice, con(,2)

konektivita he() tlouštky sten, pe() vnitrní

pretlaky E200e9mi0.28hhe(ie)ppe(ie) i1con

(ie,1)i2con(ie,2) ig3i1-2 3i1-1 3i1

3i2-2 3i2-1 3i2 r1y(i1)r2y(i2) l((x(i1)-

x(i2))2(r1-r2)2)0.5 r(r1r2)/2 c(x(i2)-x(i

1))/ls(r2-r1)/l b1/l-c -s 0 c s

0 0 l/(2r) 0 0 l/(2r) 0

0 0 -1 0 0 1 0

0 sl/(2r) 0 0 sl/(2r) -s -c

l/2 s -c l/2 dEh/(1-mi2)1 mi

0 0 0 mi

1 0 0 0 0

0 h2/12 mih2/12 0

0 0 mih2/12 h2/12 0

0 0 0 0

5(1-mi)/12 kstifrlb'db plpl/8 bp-pl

s(3r1r2) plc(3r1r2) 0 -pls(3r2r1)

plc(3r2r1) 0

Pressure loading

Rotationally symmetric shells

NAP5

Solution of system Kub is vector of

generalized displacements, for element with nodes

1,2

5 generalized forces N? N/m meridian membrane

force N? N/m circumferential force M? N

meridian moment M? N circumferential moment Q

N/m transversal force

FEM 2D and 3D elements

NAP5

Truss, beam, but also the rotary shell are

typical 1D finite elements corresponding to

ordinary differential equations for displacement,

deflection or rotation of a single independent

variable (length in the direction of the

coordinate axis of the element). Analogous

equations for displacements, deflections or

rotations also applies to 2D or 3D elements (eg.

triangles, tetrahedrons, hexahedron ...). Basis

weight functions and Ni are, however, a function

of x, y, z (again it may e.g. linear or Hermite

polynomials). Application of Galerkin method and

implementation of assembly, boundary conditions

and the solution of system of algebraic equations

is the same for 1D, 2D or 3D elements. The 2D

and 3D problems differ from the 1D problems more

or less formally, because corresponding

differential equations are partial differential

equations and not ordinary differential equations

(because there is more than one independent

variable, x,y,z). A solution of partial

differential equations will be a subject of next

lecture ...

What is important

NAP5

The lecture was devoted to solving boundary value

problems described by ordinary differential

equations. It is quite extensive and quite

important, it is important to remember that a

little more than usual ...

What is important

NAP5

How to use solution of eigenproblem for

transformation of differential equations to

separated equations

Weighted residual method. What is it weight

function. Type of weight functions used in

spectral, finite element, finite volumes and

collocation methods.

FEM for solution of popeline networks

(incompressible fluid, steady state).

Finite element PIPE

Assembly of elements to global matrix

for e1ne kl,igkloc(e)

k(ig(12),ig(12))k(ig(12),ig(12))kl(12,12)

end

What is important

NAP5

What is the difference between truss, beam and

membrane, shell.

FEM for trusses

Transformation from local coordinate systém of

element to global coordinate systém x,y

FEM for beams. Cubic polynomials

Bending and geometrical stiffness, mass matrix