MAT-VET-F 20002
Examensarbete 15 hp June 2020
Investigating strains on the
Oseberg ship using photogrammetry and finite element modeling
Andreas Eriksson
Erik Thermaenius
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress:
Box 536 751 21 Uppsala Telefon:
018 – 471 30 03 Telefax:
018 – 471 30 00 Hemsida:
http://www.teknat.uu.se/student
Abstract
Investigating strains on the Oseberg ship using photogrammetry and finite element modeling
Andreas Eriksson, Erik Thermaenius
The Oseberg ship is known as one of the finest surviving artifacts from the Viking age, with origins dated back to the 800s. The ship has been displayed in the Viking ship museum in Oslo since 1926. The nearly 100 years on museum display along with the over 1000 years it was buried has weakened the structure of the ship. To slow down the deterioration, several research projects has been initiated, among them the project ''Saving Oseberg''. A part of ''Saving Oseberg'' is contributing to the planning of a new museum for the ship. As a basis for the planning, the ship has been monitored with photogrammetry. This is intended as a way to visualise the deformation and displacements of the ship due to seasonal changes in indoor temperature and humidity.
In this thesis the photogrammetry data from the hull of the ship was used to make a finite element model, and through this model calculate the average strain on each element. The method was based on a previous research project conducted on the Swedish warship Vasa by a research group at the Division of Applied Mechanics at Uppsala University.
The measurements of the ship was formed into a hull by Delaunay triangulation. The strain was approximated as a Green strain and evaluated using isoparametric mapping of the elements. Through the nodal displacements, the strain was evaluated by approximating the elements as tetrahedrons and calculating the average strain from these elements between the measurements.
The result showed an oscillating behavior of the displacements, proving the proposal of seasonal depending displacements. The measured principal strains also matched to the corresponding relative humidity fluctuation during the year. The strain magnitude was relatively even throughout the ship, mostly varying between ±0.4% but certain areas were more subjected than others. A few elements on the starboard side showed very large strains through most of the measurements, this seemed very unusual and was probably the result of inaccuracies or errors in the data.
Though the ship is subjected to relative small strains and permanent displacements after annual cycles, the mechano-sorptive strains may lead to accumulated
deformation and eventually failure in the weak parts of the wood or at the high stress concentraion parts. In addition, the cyclic strain even in elastic range may cause fatigue failure in any material which could pose a large threat for the future conservation of the ship.
Populärvetenskaplig sammanfattning
Osebergsskeppet är känd som en av de bäst bevarade och mest åtråvärda artefakterna som överlevt från Vikingatiden. Det stora långskeppet med en omkrets på över 40 meter hittades i en gravplats utanför staden Tønsberg i Norge år 1904 och grävdes upp följande år av arkeologerna Haakon Shetelig och Gabriel Gustafson. Till en början förvarades skeppet i en av Oslo Universitets lokaler till dess att Vikingaskeppsmuséet slog upp sina dörrar år 1926, där det befunnit sig sedan dess.
Trots att skeppet inte brytits ned under de 1000 år det varit begravet har denna process nu ökat takten. Ett flertal forskningsprojekt har initierats för att sakta ned nedbrytningen och möjliggöra skeppets förefintlighet, bland annat det pågående ”Saving Oseberg”.
Då museet som huserar Osbergsskeppet är en gammal byggnad som saknar en modern klimatanläg- gning vilken kan reglerar inomhustemperatur och luftfuktighet till ideala förhållanden så ingår det i ”Saving Oseberg” att bidra med projekteringen till det nya muséet. En forskare har som underlag till projekteringen för den nya muséebyggnaden gjort punktmätningar över skrovet på skeppet med jämna mellanrum. Målet med dessa mätningar är att se hur de uppmätta punkterna förflyttas mellan mättillfällena vilket ger möjlighet att analysera deformationen av skeppet. Detta var syftet med det presenterade kandidatarbetet, som innebar att skapa en modell av skeppet från all mätdata och utifrån förflyttningarna i mätdatan översätta det till töjningar i skrovet på skeppet i %.
Mätdatan på skeppet analyserades med programvaran Matlab, där det sammanfogades och tri- angulerades till ett tvådimensionellt skal av skrovet. Utifrån förflyttningarna av mätpunkterna, beräknades den genomsnittliga töjningen på varje triangel (även kallat element) som bygger upp skrovet för varje mätning. Töjning visualiserades med hjälp av en färgskala, där negativ och positiv töjning färgade elementen blå respektive röda. Detta gjorde det enkelt att överblicka skeppet och identifiera s.k. ”hotspots” som utsätts för större töjningar än övriga delar på skeppet. Dessutom visualiserades medelvärdet av mätpunkternas förflyttningar vid varje mättillfälle för de två sidorna på skeppet i en graf.
Resultatet visade att det fanns delar på skeppet som var mer utsatta än andra, särskilt utsatt var relingen runt båten och delar av kölen. Detta förklarades enklast av att relingen saknar samma stöd som resten av skrovet och eftersom ungefär halva tyngden av skeppet vilar på kölen medför detta en högre belastning och därmed större töjningar. Halva tyngden eftersom den andra halvan är fördelat på stöttor runt sidorna av skeppet. När medelvärdet av förflyttningarna studerades uppmärksammades ett cykliskt beteende av förflyttningarna på babordsidan vilket troligen beror på årstidsvariationer i temperatur och luftfuktighet, då förflyttningar återfanns på ungefär samma positioner ett år senare. Styrbordsidan var dock mycket mer utsatt och hade ett delvis cyklisk beteende men att utslagen var betydligt högre. För alla skeppets delar syntes även flera element vars töjning var flera storleksordningar högre än elementen runt om. Detta var väldigt orimligt, och vad det berodde på blev inte klarlagt. Troligen något fel i den ursprungliga mätdatan.
Slutsatsen blev att skeppet utsätts för töjningar, att dessa skiljde sig på olika delar av skeppet och att förflyttningarna varierade cykliskt under året. Med viss osäkerhet inkluderat så var resultatet relativt väntat och rimligt.
Acknowledgements
We, the authors of this thesis, would first of all like to express our gratitude towards our supervi- sors, professor Kristofer Gamstedt and researcher Reza Afshar from the Division of Applied Mechanics at Uppsala University. Both for introducing the project to us and for their guidance and valuable input throughout the project. Along with their previous work on a similar project which has been a helpful inspiration for the general structure of this thesis.
Secondly, we want to thank doctoral research fellow David Hauer from the Department of Col- lection Management at Oslo University who is the researcher from ”Saving Oseberg” behind this project. His help with supplying us with the data and answering all our question regarding the data and ”Saving Osberg”-project in general was much appreciated.
We would also like to thank our project mentor, researcher Johan Gising from the Division of Drug Design and Discovery and the Division of Nanotechnology and Functional Materials at Uppsala University. His effort in managing the stage-gates, and his positive feedback and support when the project was falling behind the planned time schedule was a great help in keeping us on track.
Finally, we would like to acknowledge our dear classmates who kept us company during the whole time span of the project. From the lunches in the sun, to the activities on the nights and weekends was truly delightful highlights which kept our moral high considering the unexpected pandemic which fully broke out just days prior to the start of the project.
Without their contributions, this thesis would not had been made possible.
Andreas Eriksson Erik Thermaenius Uppsala, May 2020
Contents
1 Introduction 1
1.1 Background . . . 1
1.2 Objective . . . 2
1.3 Aims . . . 2
2 Theory 3 2.1 Strain and deformation . . . 3
2.2 Finite elements . . . 7
2.3 Isoparametric mapping . . . 10
2.4 Delaunay triangulation . . . 13
2.5 Mapping between vector spaces . . . 14
3 Method 17 3.1 Calculating the strain . . . 17
3.2 Verifying the strain calculation . . . 18
3.3 Measurements of the ship . . . 20
3.4 Utilising the measurement data . . . 22
4 Results 25
5 Discussion 38
6 Conclusions 39
Appendix A Strain analysis A.1
A.1 strain.m . . . A.1 A.2 elementInOnePlane.m . . . A.2 A.3 colourFromStrain.m . . . A.2 A.4 assemble.m . . . A.4 A.5 centerDataInSphere.m . . . A.7
Glossary
aft The direction towards the stern, or rear of a ship. 26, 28, 38 bow The direction towards the stem, or front of a ship. 26, 28, 38
engineering strain Describes strain during a deformation as the ratio between the change of length from the initial dimension and the original length. 3, 6
functional A mathematical construction that maps a vector space into a field of scalars. 7, 8 gunwale The top edge of the hull of a ship. 38, 39
port The left facing side of a ship relative the direction the ship is heading. 20, 21, 25–27, 30, 32, 34, 35, 38
shear strain Describes the strain as the change of angles between two points on an element during a deformation. 6
starboard The right facing side of a ship relative the direction the ship is heading. 20, 25, 28–30, 32, 36–38
stem The forward-facing part of a ship. 1, 20, 25, 32, 33, 38 stern The back part of a ship. 1, 20, 21, 25, 30, 31, 38
strake A strip of planking running alongside the bottom and sides of a ship, making up its’ hull.
18, 22
Starboard Port
Stem
Stern
Gunwale -
Strake 6
An illustration of the nautical terms described in the glossaries. Note that the
1 Introduction
1.1 Background
The Oseberg Ship is known to be one of the finest and most well-preserved artefacts from the Viking Era. The over 20 meters long and 5 meters broad ship of the type Karve was found in 1903 at a large burial mound near the Norwegian town of Tønsberg and excavated the following years 1904-1905 by the Norwegian archaeologist Haakon Shetelig and Swedish archaeologist Gabriel Gustafson, see Figure 1.1a. The age of the ship has not been determined exactly, researchers have seen that timber found in the grave with the ship is from the year 834 AD but the ship is thought to be older.[2]
Since the excavation, the most deteriorated artefacts of the find were treated with alum as a conser- vation process. However, the oak built ship was in better condition and could be slowly air dried, it was only surface treated with a mixture of creosote and linseed oil. [3] Completely destroyed parts were reconstructed with new wood, among them were the posts on the stem and stern. The ship was then stored at the University of Oslo temporally until the first halls of the Viking Ship Museum in Oslo were finished in 1926 and the ship was moved where it has been ever since, see Figure 1.1b.[2][4] Being displayed in the museum hall in over 90 years, the ship has slowly been degrading and shown signs of deformation and fracture formations, it was not until the 2000s that the extent of the degradation became apparent.[4]
(a) (b)
Figure 1.1: (a) The Oseberg Ship being excavated outside Tønsberg.[5] (b) The ship displayed in the Viking Ship Museum in Oslo.[6]
In 2014 a research project was launched called ”Saving Oseberg” with the aim of preventing and if possible stop the deterioration of the ship.[4] One part of the ”Saving Oseberg”-project is contributing to the planning of a new museum and the future move of the ships to the new building. The group has been monitoring the ship through photogrammetric surveillance, making point-wise position measuring of the ship on a regular basis. The researchers have observed that the measured positions of the ship are not consistent throughout the year, but instead are displaced from the original position due to seasonal changes in the indoor climate. The Viking Ship Museum housing the ship is an old building, and have not been modified with a modern air conditioning system making the indoor humidity and temperature of the museum irregular. The changes in humidity makes the wood of the ship swell and shrink, possibly causing the displacements of the measured positions.
1.2 Objective
The deformation of the ship can pose a potential threat for the future conservation of it, depending on the magnitude and character of the deformations. Elastic deformations are completely reversible, meaning that the displaced points of the ship will return to their original position after a time cycle. But even if the deformations are elastic, cyclic deformation as fatigue in elastic range may cause failure in any material. Plastic deformations are irreversible, slowly deforming the ship until inevitable fracture. Concerned with the displacements over time, the research team monitoring the ship wants to determine if the ship is subjected to irreversible deformations or if the displacements are strictly composed of elastic deformations as well as if there are ”hot spots” where the strain is larger than in other places.
A research group at the Division of Applied Mechanics at Uppsala University have conducted a similar study on the Swedish warship Vasa. Using geodetic measurements of the ship over a 12- year period, they computed the average strain of different sections on the hull of the ship through the measured position data from the original measurement to the final measurement. The building housing the ship have an advanced air conditioning system installed monitoring, among other things, the humidity in the museum which is kept constant. Their result was that the Vasa ship was slowly deformed during the monitored time and subjected to a constant stress called creep deformation.
The displacements varied in different regions and directions on the hull.[7]
The objective of this thesis is to conduct a similar study as the one conducted on the Vasa Ship.
Using the same technique as the Vasa Ship project, this project aims to visualise the displacements of the Oseberg ships hull during the measured time. The data will in this project instead be measured using photogrammetry. Furthermore, using the same technique as the Vasa Ship project, calculate the average strain of the hull, determine the strain in percentages with respect to the reference position along the different directions on the hull and the different sections of the ship.
1.3 Aims
We the authors of this thesis aim to produce a finite element model of the ship using the gathered position data with Matlab and its triangulation functions to visualise the hull. This model will be easy to implement on general position data for later use on similar projects. If the ship is further studied and more position data is added it should be possible to implement the code on this data as well. The model will calculate the average strain on finite elements composed from the data. The result will then be plotted with colourised triangular elements where the colour is decided from the scale of the strain. The plots for each time point will then be used in an animation to give a better view of the deformation over time. From the result it can be discussed where the strains are most tangible and possibly determine the cause of these strains.
2 Theory
2.1 Strain and deformation
Consider a point in an arbitrary body B in R3, P : (x, y, z). Then look at the same point in the deformed body B⇤ in R3, P⇤ : (x⇤, y⇤, z⇤). Assuming the coordinate functions (x⇤, y⇤, z⇤) are continuous and differentiable in the (x, y, z) variables. This implies that
8>
<
>:
x = x(x⇤, y⇤, z⇤), y = y(x⇤, y⇤, z⇤), z = z(x⇤, y⇤, z⇤).
Note that (x⇤, y⇤, z⇤) is considered as independent variables. It can be shown that when dealing with small displacements it is not necessary to distinguish between the two points (x, y, z) and (x⇤, y⇤, z⇤).
Now lets add the point Q : (x + dx, y + dy, z + dz) to B to make a line element ds = P Q. This is deformed in the body B⇤ into ds⇤ by adding the point Q⇤ : (x⇤+ dx⇤, y⇤+ dy⇤, z⇤+ dz⇤). See Figure 2.1.
x z
y y⇤
z⇤
x⇤ B
B⇤
•
•Q : (x + dx, y + dy, z + dz) P : (x, y, z)
•
•
P⇤: (x⇤, y⇤, z⇤)
Q⇤: (x⇤+ dx⇤, y⇤+ dy⇤, z⇤+ dz⇤)
ds
ds⇤
(u, v, w)
Figure 2.1: The line segments ds between the points P , Q in the undeformed body B and ds⇤ between the points P⇤, Q⇤ in the deformed body B⇤. Note the distance between the points P and P⇤ is (u, v, w).
The deformation of ds can be described by the engineering strain "E
" = ds⇤ ds
. (2.1)
Since there are two different bodies, B and B⇤ which do not necessarily have the same coordinate axes, each line element variable can be expressed as
dx⇤ = @x⇤
@xdx +@x⇤
@y dy +@x⇤
@z dz. (2.2)
The same follows for dy⇤ and dz⇤. Noting from Figure 2.1 the displacement of the point P to P⇤ as (u, v, w), the relation between the spaces is described as
8>
<
>:
x⇤= x + u, y⇤ = y + v, z⇤= z + w.
(2.3)
The length of the line element ds can be described as
|ds|2= dx2+ dy2+ dz2,
|ds⇤|2= dx⇤2+ dy⇤2+ dz⇤2.
There are several ways to calculate strain, among them are Greens strain which is used to calculate finite strains. The expression for Greens strain "G is
"G= 1 2
ds⇤2 ds2 ds2 .
This can also be written as the magnification factor M which is just another expression for Greens strain, see [8] for the full derivation of M.
M = 1 2
✓⇣ds⇤ ds
⌘2
1
◆
= "E+1 2"2E
= l2"xx+ lm"xy + ln"xz+ ml"yx+ m2"yy+ mn"yz+ nl"zx+ nm"zy+ n2"zz
= l2"xx+ m2"yy+ n2"zz+ 2lm"xy+ 2ln"xz+ 2mn"yz
(2.4)
where the parameters l, m and n are defined as l = dx
ds, m = dy
ds, n = dz
ds (2.5)
and the strain terms as
"xx = @u
@x+ 1 2
✓@u
@x
◆2
+
✓@v
@x
◆2
+
✓@w
@x
◆2! ,
"yy = @v
@y +1 2
✓@u
@y
◆2
+
✓@v
@y
◆2
+
✓@w
@y
◆2! ,
"zz = @w
@z +1 2
✓@u
@z
◆2
+
✓@v
@z
◆2
+
✓@w
@z
◆2! ,
"xy = "yx= 1 2
✓@v
@x+@u
@y + @u
@x
@u
@y +@v
@x
@v
@y +@w
@x
@w
@y
◆ ,
"xz = "zx= 1 2
✓@w
@x +@u
@z +@u
@x
@u
@z + @v
@x
@v
@z +@w
@x
@w
@z
◆ ,
"yz = "zy = 1 2
✓@w
@y +@v
@z + @u
@y
@u
@z + @v
@y
@v
@z +@w
@y
@w
@z
◆ .
(2.6)
Note that when the line element ds is along for example the x-axis (ds = dsx) then equation 2.5 will output l = 1 and m = n = 0. In turn equation 2.4 will output
M = Mx = "xx. Similar to equation 2.5, the deformed forms can be written as
l⇤ = dx⇤
ds⇤ = dx⇤ ds
ds
ds⇤, m⇤ = dy⇤ ds⇤ = dy⇤
ds ds
ds⇤, n⇤ = dz⇤ ds⇤ = dz⇤
ds ds ds⇤. Combining equation 2.2 and 2.3, it is known that
dx⇤ ds =⇣
1 +@u
@x
⌘l +@u
@ym +@u
@zn, dy⇤
ds = @v
@xl +⇣ 1 +@v
@y
⌘m +@v
@zn, dz⇤
ds = @w
@xl +@w
@ym +⇣
1 +@w
@z
⌘n.
Rewriting equation 2.1, it is possible to express ds
ds⇤ = 1 1 + "E.
From the equations above, the expressions for the deformed line elements is given as (1 + "E)l⇤ =⇣
1 +@u
@x
⌘l +@u
@ym +@u
@zn, (1 + "E)m⇤ = @v
@xl +⇣ 1 +@v
@y
⌘
m +@v
@zn, (1 + " )n⇤ = @w
l +@w m +⇣
1 +@w⌘ n.
This will be useful in the following. The goal is to explain how a shear strain can be expressed.
Considering two line elements ds1 and ds2, the product of these vectors can be expressed as ds1· ds2 = cos (✓), ds⇤1· ds⇤2= cos (✓⇤).
From this the shear strain between the two deformed line elements can be written as
12= (1 + "E1)(1 + "E2) cos (✓⇤)
= 2l1l2"xx+ 2m1m2"yy+ 2n1n2"zz
+ 2(l1m2+ l2m1)"xy+ 2(m1n2+ m2n1)"yz+ 2(l1n2+ l2n1)"xz.
If ds1 is oriented along the x-axis and ds2 along the y-axis, then l1 = m2 = 1and m1 = n1 = l2 = n2= 0. This yields
12= xy = 2"xy which also implies that
xy = 2"xy, yz= 2"yz, xz = 2"xz.
Assuming small displacements, the quadratic terms in equation 2.6 can be discarded. This yields the following simplification
M ⇡ "E,
"xx ⇡ @u
@x, "yy ⇡ @v
@y, "zz ⇡ @w
@z, (2.7)
"xy ⇡ 1 2
⇣ @v
@x+ @u
@y
⌘, "xz ⇡ 1 2
⇣@w
@x +@u
@z
⌘, "yz⇡ 1 2
⇣@w
@y +@v
@z
⌘.
Since the engineering strain " and shear strain is of interest, this is related to the magnification factor M by equation 2.4. When calculating the strain in practice, the following strain tensor is obtained [9]
" = 2 64
"xx "xy "xz
"yx "yy "yz
"zx "zy "zz
3 75 =
2 64
"xx xy/2 xz/2
yz/2 "yy yz/2
zx/2 zy/2 "zz
3
75 . (2.8)
It should also be mentioned that one can calculate the principle strains from this tensor. The principle strains can be explained by imagining a volume in the 3D space which is deformed into a new volume. The deformation can be in any directions defined by a global coordinate system. If a local coordinate system of the body is set so the deformation is along the coordinate axes, the strain tensor " will be a diagonal matrix. Doing this it will eliminate any shear strains and only be made up of three directional strains. An easy way to implement this is by calculating the eigenvalues of the strain tensor " as
"xx 1 "xy "xz
"yx "yy 2 "yz
"zx "zy "zz 3
= det(" I) = 0 (2.9)
2.2 Finite elements
When solving equations and computational problems numerically, the equation and computational domain needs to be discretised. Once a computational domain is discretised, it is called a mesh.
There are several different methods for discretising computational domains, most common are as point wise differences or as finite elements. These types of discretisations originates from methods for solving partial differential equations numerically, the above mentioned comes from the finite difference method and finite element method.[10]
The reason for meshing a domain is because it allows easy constructions of polynomial function spaces that is otherwise difficult. Discretisation using finite elements can utilise several different types of polygons to shape the mesh, most common are lines, triangles, tetrahedrons and bricks depending on the dimension and type of the computation domain. In this thesis the finite elements are composed of tetrahedral geometries K, these are defined by a polynomial function space P . The polynomial space is defined by a set of linear functionals Lj(·) defining the degrees of freedom, where i = 1, 2, . . . , n and n = dim (P ). The function space P is spanned by a basis {Ni}ni=1 and since they form the geometry of the finite element they will be referred to as the shape functions.
The functionals Lj(·) is used to uniquely define the shape functions by requiring them to satisfy the relation [10]
Xn i,j=1
Lj(Ni) = ij =
(1, i = j,
0, i6= j. (2.10)
Combining this with the theory of strain mentioned in Section 2.1, a numerical method for computing strain over an arbitrarily body is achieved. Assuming the body B, containing the point P in Figure 2.1 is discretised with 3D polygons and is then displaced into the body B⇤ and the point P into P⇤. The polygon element containing P is called K with the nodes ni, see Figure 2.2.[9]
x z
y •
B
• P K
• •
•
•
n1
n2
n3
n4
Figure 2.2: The element K on the body B containing the point P .
The displacement components of the element (u, v, w) are continuous functions of the spatial coordi- nates (x, y, z), where each node niof the polygon have separate displacement components (ui, vi, wi), see Figure 2.3.[9]
x z
y
u1 v1
w1
u2
v2
w2
u v
w
u4
v4
w4
u3 v3
w3
•
•
•
•P 1
4
2
• 3 K
Figure 2.3: The tetrahedron element K, showing the nodal displacements and the displacement of the point P .
For the implementation done in this thesis it suffices to use continuous functionals Li(·) with dis- continuous derivatives at the element boundaries. This type of functionals are called Lagrange functionals and are defined as
Lj(v) = v(Ni), Ni = f (x(i), y(i), z(i))
where (x(i), y(i), z(i))are selected nodal points. For a tetrahedral element with linear shape functions where the polynomial space P is the space P1(K). The simplest basis for this space will be the canonical basis {1, x, y, z} such that any of the shape functions can be expressed as
N1 = a1+ a2x + a3y + a4z
where {ak}4k=1 are coefficients unique to each shape function. For a reference element with nodal points at (0, 0, 0), (1, 0, 0), (0, 1, 0) and (0, 0, 1), see Figure 2.4, the coefficients can easily be calcu- lated.
z
y
• x
•
•
•
(0, 1, 0) (0, 0, 0)
(1, 0, 0) (0, 0, 1)
K
Figure 2.4: A uniform reference tetrahedron.
From the criteria of equation 2.10 the following system is obtained
e1 = 2 66 66 4
1 0 0 0
3 77 77 5=
2 66 66 4
L1(1) L1(x) L1(y) L1(z) L2(1) L2(x) L2(y) L2(z) L3(1) L3(x) L3(y) L3(z) L4(1) L4(x) L4(y) L4(z)
3 77 77 5
2 66 66 4
a1
a2 a3 a4
3 77 77 5=
2 66 66 4
1 0 0 0 1 1 0 0 1 0 1 0 1 0 0 1
3 77 77 5
2 66 66 4
a1
a2 a3 a4
3 77 77
5= V· a
where the matrix V is called the Vandermonde matrix. The coefficients a can be calculated as a = V 1e1 which in this case gives a = h
1 1 1 1
i|
. It implies that N1 = 1 x y z. Similar calculations for N2, N3 and N4 will give the following
N1= 1 x y z, N2= x,
N3= y, N4= z.
This is only true for the element presented in Figure 2.4 because of its specific geometry. For other geometries the shape functions would be different. Although for easy calculations it is preferred that all elements is of this geometry, and it can actually be accomplished using an isoparametric mapping which will be discussed in the following section.[10]
The shape functions are related to the displacement functions as
u(x, y, z) = X4
i=1
Ni(x, y, z)ui,
v(x, y, z) = X4
i=1
Ni(x, y, z)vi,
w(x, y, z) = X4
i=1
Ni(x, y, z)wi.
(2.11)
Note that the shape functions sum up to one.[9]
2.3 Isoparametric mapping
When dealing with curved geometries and non-uniform finite elements, the computations become advanced and cumbersome when calculating the element shape functions. This can be solved using an isoparametric mapping of the coordinates which builds up each element. This will make the shape functions equal for all elements even if their geometries may vary, making it possible to reuse the calculated shape functions for all elements. This makes the calculations a lot less complicated especially if a higher order shape function is used, for example a quadratic shape function. The higher order shape functions is very common for curved geometries and since these have more terms that need calculation, isoparametric mapping is almost essential when dealing with these.[10]
The following set up is used. Consider a finite element K called the physical element, this will usually not be uniform and is described in the global coordinate system. The element uesd for easy calculations is on the other hand ¯K called the reference element or the natural element, where the nodes are given in the local coordinate system. If the element have i = 1, 2, ..., n nodal points, each coordinate can be descried in the physical nodal points from the reference nodal points as[9][10]
x(r, s, t) = Xn i=1
x(i)Ni(r, s, t),
y(r, s, t) = Xn i=1
y(i)Ni(r, s, t),
z(r, s, t) = Xn i=1
z(i)Ni(r, s, t).
(2.12)
Ni is the shape function of each coordinate-axis and (r, s, t) is the local coordinates which maps (x(i), y(i), z(i)) of ¯K to the global coordinates of K. Therefore a description of the mapping from the reference coordinates to the physical coordinates is obtained, see Figure 2.5.[10]
z
y
• x •
•
K
•
(a)
t
s
• r
•
•
•
1
1 1
1
K¯
(b)
Figure 2.5: (a) The non-uniform element K. (b) The isoparametric mapping K¯.
This can also be written in vector notation as
P (r, s, t) = N· ¯P where ¯P =h
x1 y1 z1 x2 y2 z2 . . . i|
is the vector for the nodal points and N is the matrix of the shape functions which is written as
N = 2 64
N1 0 0 N2 0 0 N3 0 0 N4 0 0
0 N1 0 0 N2 0 0 N3 0 0 N4 0
0 0 N1 0 0 N2 0 0 N3 0 0 N4
3 75
for finite elements with four nodal points.[9]
In the same way a function v of the finite element K can also be described as equation 2.12
v(r, s, t) = Xn
i=1
viNi(r, s, t). (2.13)
Since equation 2.7 have partial derivatives of first order, the chain rule is used to derive an expression
like equation 2.13
@v
@x = @v
@r
@r
@x+ @v
@s
@s
@x +@v
@t
@t
@x,
@v
@y = @v
@r
@r
@y +@v
@s
@s
@y+@v
@t
@t
@y,
@v
@z = @v
@r
@r
@z + @v
@s
@s
@z +@v
@t
@t
@z. This can in vector form be written as
2 66 66 66 4
@v/@x
@v/@y
@v/@z 3 77 77 77 5
= 2 66 66 66 4
@r/@x @s/@x @t/@x
@r/@y @s/@y @t/@y
@r/@z @s/@z @t/@z 3 77 77 77 5
· 2 66 66 66 4
@v/@r
@v/@s
@v/@t 3 77 77 77 5
= J 1· 2 66 66 66 4
@v/@r
@v/@s
@v/@t 3 77 77 77 5
where J is the Jacobian matrix[10]
J = 2 66 66 66 4
@x/@r @y/@r @z/@r
@x/@s @y/@s @z/@s
@x/@t @y/@t @z/@t 3 77 77 77 5 .
The Jacobian can also be used to describe the differential terms ds from equation 2.2 as dx = @x
@rdr + @x
@sds +@x
@tdt, dy = @y
@rdr +@y
@sds + @y
@tdt, dz = @z
@rdr + @z
@sds +@z
@tdt,
=) ds = J|·h
dr ds dti| .
The shape functions of the two coordinate systems is also related like this [9]
2 66 66 66 4
@Ni/@x
@Ni/@y
@Ni/@z 3 77 77 77 5
= J 1· 2 66 66 66 4
@Ni/@r
@Ni/@s
@Ni/@t 3 77 77 77 5
. (2.14)
For equation 2.13 the Jacobian will be decomposed as follows J11= @x
@r =X
i
@Ni
@r x(i), J12= @y
@r =X
i
@Ni
@r y(i), J13= @z
@r =X
i
@Ni
@r z(i), J21= @x
@s =X
i
@Ni
@s x(i), J22= @y
@s =X
i
@Ni
@s y(i), J23= @z
@s =X
i
@Ni
@s z(i), J31= @x
@t =X
i
@Ni
@t x(i), J32= @y
@t =X
i
@Ni
@t y(i), J33= @z
@t =X
i
@Ni
@t z(i).
(2.15)
Following the properties of the elements, polynomial spaces can be constructed on each element.
Combined with the isoparametric mapping of a element, the construction of the nodal shape func- tions becomes fairly easy as described in the previous section.
2.4 Delaunay triangulation
Mesh generation is an essential tool within many branches of mathematics and computer science, examples ranges from image analysis, computer graphics and numerical analysis. When discretising a domain with triangular polygons for finite element analysis the meshes should be suitable for interpolation. The aim is to find a set of triangles that covers the domain and satisfy certain constrains such as the angles in the triangles should not be to small or to large, and the size of triangles should not be smaller nor larger than necessary. Delaunay triangulation and refinement is a mathematical algorithm that guarantees these constraints are met when discretising domains.[10][11]
Such a mesh generator that triangulates 2D surfaces attempts to satisfy three conditions:
Firstly, the union of the triangles is the triangulation domain, meaning that the triangles shall cover the whole discretisation domain.
Secondly, the triangles should have a sufficient shape, bounding how large and small the angles of the triangles can be. A lower bound for an angle implicitly bounds the largest angle, meaning if no angle is smaller than ✓ then no angle is larger than 180 2✓. Hence, mesh generators typically approach to bound the smallest angle.
Thirdly, the generator should offer as much control as possible over the size of the triangles. This makes it possible to trade off accuracy against speed when dealing with finite element analysis.
Large, sparse triangles offers fast computations but often have low accuracy and small, densely packed triangles offers great accuracy but are computationally expensive. Given this option, the user can make this trade off depending on the given problem. A coarse mesh with few triangles should be able to be refined with more triangles, and the opposite. This is made possible with the Delaunay refinement algorithm.[11]
Both Delaunay triangulation and refinement is built in as functions in Matlab, and therefore a sufficient meshing method for finite element problems is made possible. This function makes it possible to both triangulate a given a geometry or a shape, and given a set of data points which can be formed into a surface, see Figure 2.6 and 2.7.[10]
(a) (b)
Figure 2.6: (a) A random elliptical body. (b) The Delaunay triangulation of that body.
(a) (b)
Figure 2.7: (a) A set of 30 random data points. (b) The Delaunay triangulation of those points.
2.5 Mapping between vector spaces
A point vector in the Cartesian space can be represented by P = aˆı + bˆ|+ cˆk. Where ˆı, ˆ|, ˆk are unit vectors along the coordinate axes, in this case x, y, z. This can also be written as
P = 2 66 66 4
x y z w
3 77 77 5,
a = x/w b = y/w c = z/w
be written ash
1 0 0 1 i|
or h
2 0 0 2 i|
and similarly for the other unit vectors.
In this thesis it is important to look at translation and rotation of vector spaces in case the mea- surements of the data have different origins and coordinate axes. It will be needed to transform all vector spaces into one global space. Lets look at transformation of the space H with a 4 ⇥ 4 matrix M. This matrix will be able to represent translation, rotation, scaling and shearing of the space.
Given a point P0 represented in two different spaces, the same point will be written as P1 in the second space. The transformation will then be written as
P1= M· P0.
If the matrix only need to deal with translation it will look like
M =Trans(a, b, c) = 2 66 66 4
1 0 0 a 0 1 0 b 0 0 1 c 0 0 0 1
3 77 77 5.
In turn the points will be translated as
P1 = 2 66 66 4
1 0 0 a 0 1 0 b 0 0 1 c 0 0 0 1
3 77 77 5·
2 66 66 4
x y z 1
3 77 77 5=
2 66 66 4
x + aw y + bw z + cw
w 3 77 77 5=
2 66 66 4
x/w + a y/w + b z/w + c
1 3 77 77 5.
If the matrix instead deal with rotation only it will look like the following for rotations about the x, y and z axes by an angle ✓
M =Rot(x, ✓) = 2 66 66 4
1 0 0 0
0 cos (✓) sin (✓) 0 0 sin (✓) cos (✓) 0
0 0 0 1
3 77 77 5,
M =Rot(y, ✓) = 2 66 66 4
cos (✓) 0 sin (✓) 0
0 1 0 0
sin (✓) 0 cos (✓) 0
0 0 0 1
3 77 77 5,
M =Rot(z, ✓) = 2 66 66 4
cos (✓) sin (✓) 0 0 sin (✓) cos (✓) 0 0
0 0 1 0
0 0 0 1
3 77 77 5.
Scaling and shearing are represented by exchanging the diagonal elements by the axis scale factor.
By combining these different operations the matrix M will be able to transform one vector space
into another.[12] The transformation matrix can therefore be described by a set of variables as
M = 2 66 66 4
↵11 ↵12 ↵13 a
↵21 ↵22 ↵23 b
↵31 ↵32 ↵33 c
0 0 0 1
3 77 77 5.
3 Method
3.1 Calculating the strain
The data obtained for calculating the displacement and strain of the ship is made up of a number of measurement points. These points are given as coordinates in a 3D Cartesian domain at a certain time. The measurements has been done at several time points. This data will be analysed to obtain the movement of the data points over time. The theory presented in Section 2.2 gives us a way of joining the data points to form finite elements making up the hull of the ship and determine the strain using the nodal points. The finite element discretisation is done with Matlab using its built-in Delaunay triangulation function.
We use the theory presented in Section 2.1 to describe the strain in each finite element as a strain tensor " given by equation 2.8. Along with equation 2.7 we get an expression for the tensor we want to calculate as
" = 1 2
2 66 66 66 4
@u/@x + @u/@x @v/@x + @u/@y @w/@x + @u/@z
@v/@x + @u/@y @v/@y + @v/@y @w/@y + @v/@z
@w/@x + @u/@z @w/@y + @v/@z @w/@z + @w/@z 3 77 77 77 5
= 1
2(A + A|)
where A is given by
A = 2 66 66 66 4
@u/@x @u/@y @u/@z
@v/@x @v/@y @v/@z
@w/@x @w/@y @w/@z 3 77 77 77 5
=h
u v w i
· 2 66 66 66 4
@/@x
@/@y
@/@z 3 77 77 77 5 .
Using this we find that the partial derivatives of the displacement functions (u, v, w) can be calcu- lated by equation 2.11. For an easier calculation we use an isoparametric mapping as described in Section 2.3, then the partial derivatives of the shape functions can be rewritten as equation 2.14 using the Jacobian presented. The Jacobians elements is given by equation 2.15. By these equations along with the displacement vector of the finite element and its original positions we can easily com- pute the strain tensor, see Algorithm 1. This is implemented in a simple Matlab routine presented in Appendix A.1. The algorithm takes three points in the Cartesian 3D coordinate system and adds a fourth node in the normal from the plane between the nodes making a tetrahedral as the finite element, see Figure 3.1. From this it calculates the strain tensor in the element.
In practice, strain is a 3D phenomenon, but in this thesis the measuring data only gives information for a 2D surface. This means that the in-and-out strain of the hull cannot be interpreted correctly with this information and is therefore neglected in the analysis. By doing this procedure and adding an additional node to the finite elements, creating 3D finite elements on a 2D hull essentially makes the calculated strain in-and-out of the hull equal to zero. The actual strain can be non-zero as explained in [7]. We will therefore only analyse the strain along and across the hull. Because of this
we have also made it possible to calculate the strain along the strakes of the ship by rotating each element to the same 2D plane during calculation.
Algorithm 1 Calculate the average strain of triangular finite element
Input: Vectors with nodal points of initial and deformed triangular element X0 & X1
Output: Strain tensor of the initial element "
1: Compute normal vectors to elements ni= ((Xi(2) Xi(1))kX⇥(Xi(3) Xi(1)))
ik , i = 0, 1
2: Find centre points of triangular elements ci=E(Xi), i = 0, 1
3: Place normal vectors on element centre points to form tetrahedrons Xi(4) = ci+ ni, i = 0, 1
4: Calculate the displacement between the tetrahedrons u = X1 X0
5: Find the matrix for the derivative of the shape functions with respect to the isoparametric mapping of the elements @Nrst
6: Compute the Jacobian of transformation from isoparametric to Cartesian coordinates J =
@Nrst· X0
7: Compute the matrix for the derivative of the shape functions with respect to the reference coordinates of the elements @Nxyz= J 1· @Nrst
8: Compute the derivative of the displacement with respect to the reference coordinates @u =
@Nxyz· u
9: Compute the strain for the initial element " = 12· (@u + @u|)
•
•
• (a)
•
•
•
(b)
•
•
•
•
(c)
Figure 3.1: (a) Initial triangle element. (b) Adding normal vector. (c) Obtained tetrahedron element.
3.2 Verifying the strain calculation
The implementation of the method is done via Matlab using the function in Appendix A.1. To verify it we construct a square made up of two elements. The magnitude of the strain from the ship will be around 0.1%, therefore the square spans an area of 1000⇥1000 to show how movement of 1 unit will affect the result. If one side of the square is moved by 1 unit out from the square the expected output is 0.1%. If a shear is applied, moving one side to either direction parallel to the side by 2 units the expected output should also be 0.1%. This is tested and shown in the following figures where both input and output is presented.
y z
x
P3: (0, 1000, 1000)
P1: (0, 0, 0)
P30 : (0, 1000, 1001) P4: (0, 0, 1000)
P2: (0, 1000, 1000) P40: (0, 0, 1001)
(a)
1 0
1 0 500 1,000
0 500 1,000
x-axis
y-axis
z-axis
z-strain 0.1%
0.2 0.1 0 0.1 0.2
Strain [%]
(b)
Figure 3.2: (a) The elements sent into the analysis script exposed to a strain in the z-direction. (b) The output of the script.
y z
x
P3: (0, 1000, 1000)
P1: (0, 0, 0)
P30: (0, 1001, 1000) P4: (0, 0, 1000)
P2: (0, 1000, 0) P20: (0, 1001, 0)
(a)
1 0
1 0 500 1,000
0 500 1,000
x-axis
y-axis
z-axis
y-strain 0.1%
0.2 0.1 0 0.1 0.2
Strain [%]
(b)
Figure 3.3: (a) The elements sent into the analysis script exposed to a strain in the y-direction. (b) The output of the script.
y z
x
P3: (0, 1000, 1000)
P1: (0, 0, 0)
P30: (0, 1000, 1002) P4: (0, 0, 1000)
P2: (0, 1000, 1000) P40: (0, 2, 1000)
(a)
1 0
1 0 500 1,000
0 500 1,000
x-axis
y-axis
z-axis
y-shear 0.1%
0.2 0.1 0 0.1 0.2
Strain [%]
(b)
Figure 3.4: (a) The elements sent into the analysis script exposed to a shear in the zy-direction. (b) The output of the script.
The output is one of the elements of the strain tensor ". In Figure 3.2 it is "zz which is plotted as a colour. For Figure 3.3 it is "yy and for Figure 3.4 it is the "yz element. In these cases the other elements are zero when these specified deformations are applied. This is not shown in this thesis but can easily be seen when applying the function in Appendix A.1. These results corresponds as expected and the method is therefore correct.
3.3 Measurements of the ship
The measurements of the ship are done using photogrammetry. The measuring was divided into four sections of the ship, and all sections had fixed points on the museum floor as reference for the coordinate systems. All points on the ship and surrounding environment are fixed year around and are not removed after each measurement, guaranteeing reliable comparison between each measure- ment, see Figure 3.5b. The individual measurement sections of the ship are the stem, the stern, the starboard and port sides. A total of 15 different measurements have been conducted at the time of this thesis between August 2018 and December 2019, see Table 3.1. As mentioned in the introduction, the museum building does not have the ability to control the humidity sufficiently causing the variations in temperature and humidity between the measurements.
The measurements for each section was saved as ATOS files for the software GOM Inspect, and each section was divided into further components. GOM Inspect held information on how many points each component of every section was composed of, and gave warnings if any point was not found on later measurements, see Figure 3.5a. A threshold residual around each point could be set to help the program identify points which had been lost. A maximum residual of 10 mm was set so no point would intervene with another point, as there were no points closer than 10 mm to
Table 3.1: The dates for the measurements along with the temperature and relative humidity inside the museum at the time of the measurements.
Measurement date Temperature [ C] Relative humidity [%]
24 August 2018 21.0 60.5
19 September 2018 19.2 66.1
17 October 2018 16.5 68.9
15 November 2018 17.0 61.6
12 December 2018 15.5 49.4
10 January 2019 15.1 48.3
6 February 2019 14.6 45.6
7 March 2019 15.6 48.2
3 April 2019 16.5 42.4
29 April 2019 19.3 51.5
30 May 2019 18.3 56.8
27 June 2019 22.2 66.1
22 July 2019 22.9 69.2
21 August 2019 22.0 62.1
16 December 2019 15.5 54.7
(a)
(b)
Figure 3.5: (a) Example of how the measurement points looks in the GOM software, in this case the stern of the ship. All measurement points are taken from a set of images which can also be seen in the program. (b) One of the images which the measurement points are calculated from, in this case the port side of the stern.
Table 3.2: The measurement sections with components and number of measure- ment points from the initial measurement in August 2018.
Stem
Component Measurement points
Stem keel 10
Stem keel beam 13
Stem outer section 17
Stem strake 7 10
Stem strake 10 8
Stem strake 12 20
Stem supports 9
Stern
Component Measurement points
Stern keel 9
Stern keel beam 13
Stern outer section 13
Stern strake 7 10
Stern strake 10 9
Stern strake 12 15
Stern supports 9
Starboard
Component Measurement points
Starboard keel 10
Starboard keel beam 5 Starboard strake 7 10 Starboard strake 10 6 Starboard strake 12 20
Starboard supports 9
Port
Component Measurement points
Port keel 12
Port keel beam 6
Port strake 7 10
Port strake 10 10 Port strake 12 16
Port supports 10
As there were points missing after the residual was set to 10 mm for some of the measurements, we reconstructed those missing points as the mean value of the measurement before and after.
Reconstructed points was highlighted as black dots in the results.
3.4 Utilising the measurement data
Each data file can be analysed as earlier described. This could be of interest if a certain part of the ship is more exposed, although it would be of greater interest to look at the whole ship all at once. Since the measurements have been done in separate parts for the ship it will have to be merged into one data file for each time point. This is done by transforming the coordinate systems into one common system. If we are given two sets of data which should be joined by common points, we can utilise the theory in Section 2.5. Both data sets need to contain three common points for the following to be used. The first data set described in the space U have the points u0, u1, u2. The data set V also have points v0, v1, v2 which are equal to those of U but described by different coordinates. We want to rotate and translate the set V so that the common points are matched with set U. The first step is to define a common origin as one of these points and settingh i
new origin h
0 0 0i
by adding ru =h
0 0 0i
u0 to all points in both U and V. The rotation of V can now be performed, but first the rotation matrix R have to be identified. Since the origin in both sets are u0= v0=h
0 0 0i
, a base for these two spaces can be described by three vectors pointing out from the origin. The first two basis vectors are u1, u2 and v1, v2 translated as earlier described, the third is given by the cross product of these forming an orthogonal vector to these two. Now three basis vectors are given in set U as u1, u2, u3 = u1⇥ u2 and set V as v1, v2, v3 = v1⇥ v2. The rotation matrix is then given by forming two matrices with the basis vectors as rows and calculating R from the inverse
2 64
u1 u2
u2 3 75 = R ·
2 64
v1 v2
v3 3
75 =) R = 2 64
v1 v2
v3 3 75
1
· 2 64
u1 u2
u3 3 75 =
2 64
↵11 ↵12 ↵13
↵21 ↵22 ↵23
↵31 ↵32 ↵33 3 75 .
To move the original set V to match the common points in the original U all points in V must be translated to the new origin, rotated and then translated again to match U as
u⇤= (v⇤+ rv)R + ru.
This can also be implemented as in the 4⇥4 matrices in Section 2.5. The different systems will have parallel axes and will only be mirrored or translated between each other, therefore it will not be needed to do any scaling and only translation and rotation in the transformation matrix. In practice this is done with the Matlab routine in Appendix A.4.
Once all the parts are merged into one file for the ship, the data is triangulated using the built-in Delaunay triangulation function in Matlab. The triangulated model is saved as a matrix where each row contains the three indices from the initial data defining each triangular element. The triangulation is done with respect to two of the axes, as we are interested in a 2D hull and not a 3D body. The ship forms an oval shape and to translate this into a plane for the triangulation it was first projected onto the underside of a large sphere. This was achieved through Algorithm 2 and was implemented as a Matlab routine in Appendix A.5. The triangulation could then be done on the x and y coordinates.
Algorithm 2 Project a given data set on the boundary of a sphere
Input: A data set X with rows for each data point and four columns as (index, x, y, z), the coordinates of the origin of the data set C0, the coordinates for the new origin of the data set C1
Output: The data set centred around the new origin XC, the data set projected on the surface of a sphere XS
1: Find the radius of the sphere r = C1 C0
2: Add the radius to all the coordinates of the initial data set XC = X + r
3: Define the boundaries of the sphere in all directions as the maximum and minimum values of XC and save it as a vector R
4: Project the new data set XC on the boundary of the sphere XS = XkXC·R
Ck
Once the triangulation is obtained, the strain in each element is calculated using Algorithm 1 and the strain is then translated as a corresponding RGB colour in a chosen colourmap. This is implemented
in Algorithm 3 which inputs two sets of data, the initial set and the deformed, the triangulation of one of the sets, the chosen axis for the strain calculation and spectrum of the strain in %. This was implemented as a Matlab routine presented in Appendix A.3.
Algorithm 3 Map the colour of finite elements depending on the strain of the element
Input: Row vectors of the initial and deformed data set X0 & X1, triangulation for one of the sets T, chosen axis for strain calculation (x, y, z), vector s with spectrum of the strain in %
Output: Matrix with RGB coordinates for each element from T from a chosen colourmap depend- ing on the strain
1: Create zero matrix for allocation of colour coordinates for all elements C = zeros(size(T))
2: Calculate the strain " in each element using Algorithm 1, and save the strain for the chosen axis in a vector Ei = "i(axis, axis)
3: Choose a colourmap and save the RGB coordinates of the map in a separate vector M
4: Find indices of the colourmap limited by the specified spectrum and corresponds to the strain I =j size(M)
·|s(1)|
Ps+E·100·size(M)P s
m
5: Allocate each element in the matrix C by the indices of the colourmap by the values of I, Ci = M(Ii), if Ii > Mmax then Ii = Mmax of if Ii < 1 then Ii= 1
We only choose to focus on five of the measurements during the analysis, the measurements from January, April 3rd, June, August and December 2019 with the January measurement as reference.
The reasoning was that this gives a good idea over the seasonal dependence of the strain, as only one specific year was analysed and the measurements are somewhat even distributed and includes the different seasons and where the different peaks in temperature and humidity was located, see Figure 3.6.
Feb 2019 Apr 2019 July 2019 Oct 2019 Dec 2019
40 50 60 70 80
Relativehumidity[%]
Figure 3.6: The relative humidity measured inside the museum during 2019 along with the humidity for the chosen measurements. The highlighted markers shows the chosen measurements.
The spectrum of the strain in the plots was set between ±0.4% through initial observations of the
