NUMERICAL ANALYSIS of PROCESSES

1
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
2
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
3
ODE heat exchangers
NAP5

4
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)
5
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).
6
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
7
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
8
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).
9
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
10
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).
11
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
12
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.
13
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)

14
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 ...).
15
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
16
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
17
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)
18
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).
19
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.
20
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
21
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.

22
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
23
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
24
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).
25
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).
26
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
27
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
28
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.
29
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)
30
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.
31
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
32
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.
33
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
34
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?)
35
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
36
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
37
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).
38
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
39
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.
40
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
41
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.
42
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!
43
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
44
Beams (pipes)
NAP5
Results for element of length L
Nodal parameters (transversal displacement and
rotation)
Stiffness matrix
Geometry stiffness
Mass matrix
45
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.
46
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
47
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
48
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.
49
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)
50
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
51
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
52
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
53
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 ...
54
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 ...
55
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
56
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