Skip to content

Commit 51bde06

Browse files
committed
nafems gnl example
1 parent f412e89 commit 51bde06

133 files changed

Lines changed: 2527 additions & 2070 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

doc/CODE_OF_CONDUCT.markdown

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ chaptersDepth: 1
2020
codeBlockCaptions: false
2121
cref: false
2222
crossrefYaml: pandoc-crossref.yaml
23-
date: 2025-08-13
23+
date: 2025-08-18
2424
eqLabels: arabic
2525
eqnBlockInlineMath: false
2626
eqnBlockTemplate: |

doc/FAQ.markdown

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ chaptersDepth: 1
2020
codeBlockCaptions: false
2121
cref: false
2222
crossrefYaml: pandoc-crossref.yaml
23-
date: 2025-08-13
23+
date: 2025-08-18
2424
eqLabels: arabic
2525
eqnBlockInlineMath: false
2626
eqnBlockTemplate: |

doc/feenox.1

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
.\" Automatically generated by Pandoc 3.7.0.2
22
.\"
3-
.TH "FEENOX" "1" "2025\-08\-12" "FeenoX" "FeenoX User Manual"
3+
.TH "FEENOX" "1" "2025\-08\-15" "FeenoX" "FeenoX User Manual"
44
.SH NAME
55
FeenoX \- a cloud\-first free no\-X uniX\-like finite\-element(ish)
66
computational engineering tool

examples/.gitignore

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -125,3 +125,9 @@ steel-alum.msh
125125
steel-alum-alum.vtu
126126
steel-alum-steel.vtu
127127

128+
nafems-gnl-cantilever*.vtu
129+
nafems-gnl-cantilever*.msh
130+
nafems-gnl-cantilever*.pdf
131+
nafems-gnl-cantilever*.dat
132+
133+
steel-*.pvsm

examples/.md

Lines changed: 0 additions & 12 deletions
This file was deleted.

examples/airfoil.yaml

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
---
2+
category: laplace
3+
intro: |
4+
# Potential flow around an airfoil profile
5+
6+
The Laplace equation can be used to model potential flow, as illustrated with this example from [Prof. Enzo Dari](https://www.conicet.gov.ar/new_scp/detalle.php?keywords=enzo%2Bdari&id=32035&datos_academicos=yes) for his course "Fluid Mechanics" at [Instituto Balseiro](https://www.ib.edu.ar/).
7+
8+
For the particular case of a airfoil profile, the Dirichlet condition at the wing has to satisfy the [Kutta condition](https://en.wikipedia.org/wiki/Kutta%E2%80%93Joukowski_theorem).
9+
10+
This example
11+
12+
1. Creates a symmetric airfoil-like [Joukowsky profile](https://en.wikipedia.org/wiki/Joukowsky_transform) using the the Gmsh Python API
13+
2. Solves the steady-state 2D Laplace equation with a different Dirichlet value at the airfoil until the solution $\phi$ evaluated at the continuation of the wing tip matches the boundary value $c$. It then computes the circulation integral of the velocities
14+
a. over the profile itself
15+
b. over a circle around the airfoil, computing the unitary tangential vector
16+
i. from the internal normal variables `nx` and `ny`
17+
ii. from two functions `tx` and `ty` using the circle's equation
18+
c. around the original rectangular domain
19+
20+
It also computes the drag and the lift scalars as the integral of the pressure (computed from Bernoulli's principle) times `nx` and `nz` over the circle, respectively.
21+
22+
![Full domain](airfoil-msh1.png){width=100%}
23+
24+
![Zoom over the airfoil profile](airfoil-msh2.png){width=100%}
25+
26+
```{.python include="airfoil.py"}
27+
```
28+
29+
terminal: |
30+
$ python airfoil.py
31+
[...]
32+
Info : Writing 'airfoil.msh'...
33+
Info : Done writing 'airfoil.msh'
34+
$ feenox airfoil.fee
35+
# # step_static c e
36+
# 1 0.500000 -3.4e-03
37+
# 2 0.400000 +3.7e-03
38+
# 3 0.452067 +2.0e-09
39+
circ_profile = 2.49918
40+
circ_circle_t = 2.49975
41+
circ_circle_n = 2.49907
42+
circ_box = 2.50379
43+
pressure_drag = -0.00168068
44+
pressure_lift = 2.60614
45+
$
46+
...
47+
PROBLEM laplace 2D MESH airfoil.msh
48+
static_steps = 20
49+
50+
# boundary conditions constant -> streamline
51+
BC bottom phi=0
52+
BC top phi=1
53+
54+
# initialize c = streamline constant for the wing
55+
DEFAULT_ARGUMENT_VALUE 1 0.5
56+
IF in_static_first
57+
c = $1
58+
ENDIF
59+
BC wing phi=c
60+
61+
SOLVE_PROBLEM
62+
63+
PHYSICAL_GROUP kutta DIM 0
64+
e = phi(kutta_cog[1], kutta_cog[2]) - c
65+
PRINT TEXT "\# "step_static %.6f c %+.1e e HEADER
66+
67+
# check for convergence
68+
done_static = done_static | (abs(e) < 1e-8)
69+
IF done_static
70+
PHYSICAL_GROUP left DIM 1
71+
vx(x,y) = +dphidy(x,y) * left_length
72+
vy(x,y) = -dphidx(x,y) * left_length
73+
74+
# write post-processing views
75+
WRITE_MESH $0-converged.msh phi VECTOR vx vy 0
76+
WRITE_MESH $0-converged.vtk phi VECTOR vx vy 0
77+
78+
# circulation on profile
79+
INTEGRATE vx(x,y)*(+ny)+vy(x,y)*(-nx) OVER wing RESULT circ_profile
80+
PRINTF "circ_profile = %g" -circ_profile
81+
82+
# function unit tangent
83+
PHYSICAL_GROUP circle DIM 1
84+
cx = circle_cog[1]
85+
cy = circle_cog[2]
86+
tx(x,y) = -(y-cy)/sqrt((x-cx)*(x-cx)+(y-cy)*(y-cy))
87+
ty(x,y) = +(x-cx)/sqrt((x-cx)*(x-cx)+(y-cy)*(y-cy))
88+
89+
# circulation on the circumference with tangent
90+
INTEGRATE vx(x,y)*tx(x,y)+vy(x,y)*ty(x,y) OVER circle RESULT circ_circle_t
91+
PRINTF "circ_circle_t = %g" -circ_circle_t
92+
93+
# circulation on the outer box
94+
INTEGRATE vx(x,y)*(+ny)+vy(x,y)*(-nx) OVER circle RESULT circ_circle_n
95+
PRINTF "circ_circle_n = %g" -circ_circle_n
96+
97+
INTEGRATE -vx(x,y) OVER top GAUSS RESULT inttop
98+
INTEGRATE +vx(x,y) OVER bottom GAUSS RESULT intbottom
99+
INTEGRATE -vy(x,y) OVER left GAUSS RESULT intleft
100+
INTEGRATE +vy(x,y) OVER right GAUSS RESULT intright
101+
circ_box = inttop + intleft + intbottom + intright
102+
PRINTF "circ_box = %g" -circ_box
103+
104+
# pressure forces
105+
p(x,y) = 1.0 - vx(x,y)^2 + vy(x,y)^2
106+
INTEGRATE p(x,y)*(nx) OVER circle RESULT pFx
107+
PRINTF "pressure_drag = %g" pFx
108+
INTEGRATE p(x,y)*(ny) OVER circle RESULT pFy
109+
PRINTF "pressure_lift = %g" pFy
110+
111+
ELSE
112+
113+
# update c
114+
IF in_static_first
115+
c_last = c
116+
e_last = e
117+
c = c_last - 0.1
118+
ELSE
119+
c_next = c_last - e_last * (c_last-c)/(e_last-e)
120+
c_last = c
121+
e_last = e
122+
c = c_next
123+
ENDIF
124+
ENDIF
125+
---
126+
figures: |
127+
![Potential and velocities on the whole domain](airfoil-converged1.png){width=100%}
128+
129+
![Potential and velocities zoomed over the airfoil](airfoil-converged2.png){width=100%}
130+
...

examples/asme-expansion.ppl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
set preamble "\usepackage{amsmath}"
1+
set preamble "\\usepackage{amsmath}"
22

33
set width 14*unit(cm)
44
#set size ratio 9.0/16.0

examples/asme-expansion.yaml

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
---
2+
category: basic
3+
intro: |
4+
# On the evaluation of thermal expansion coefficients
5+
6+
When solving thermal-mechanical problems it is customary to use
7+
thermal expansion coefficients in order to take into account the mechanical
8+
strains induced by changes in the material temperature with respect to
9+
a reference temperature where the deformation is identically zero.
10+
These coefficients $\alpha$ are defined as the partial derivative
11+
of the strain $\epsilon$ with respect to temperature $T$ such that
12+
differential relationships like
13+
14+
$$
15+
d\epsilon = \frac{\partial \epsilon}{\partial T} \, dT = \alpha \cdot dT
16+
$$
17+
18+
hold. This derivative $\alpha$ is called the _instantaneous_
19+
thermal expansion coefficient.
20+
For finite temperature increments, one would like to be able to write
21+
22+
$$
23+
\Delta \epsilon = \alpha \cdot \Delta T
24+
$$
25+
26+
But if the strain is not linear with respect to the temperature---which
27+
is the most common case---then $\alpha$ depends on $T$.
28+
Therefore, when dealing with finite temperature increments $\Delta T = T-T_0$
29+
where the thermal expansion coefficient $\alpha(T)$ depends on the temperature
30+
itself then _mean_ values for the thermal expansion ought to be used:
31+
32+
$$
33+
\Delta \epsilon = \int_{\epsilon_0}^{\epsilon} d\epsilon^{\prime}
34+
= \int_{T_0}^{T} \frac{\partial \epsilon}{\partial T^\prime} \, dT^\prime
35+
= \int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime
36+
$$
37+
38+
We can multiply and divide by $\Delta T$ to obtain
39+
40+
$$
41+
\int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime \cdot \frac{\Delta T}{\Delta T}
42+
= \bar{\alpha}(T) \cdot \Delta T
43+
$$
44+
45+
where the mean expansion coefficient for the temperature range $[T_0,T]$ is
46+
47+
$$
48+
\bar{\alpha}(T) = \frac{\displaystyle \int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime}{T-T_0}
49+
$$
50+
51+
This is of course the classical calculus result of the mean value of a
52+
continuous one-dimensional function in a certain range.
53+
54+
Let $\epsilon(T)$ be the linear thermal expansion of a given material in a
55+
certain direction when heating a piece of such material from an initial
56+
temperature $T_0$ up to $T$ so that $\epsilon(T_0)=0$.
57+
58+
![Table TE2 of thermal expansion properties for Aluminum alloys from ASME II Part D (figure from [this report](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/asmeansys-potential-issue-with-thermal-expansion-calculations/).](asme-expansion-table.png){#fig:table-asme-expansion}
59+
60+
From our previous analysis, we can see that in @fig:table-asme-expansion:
61+
62+
$$
63+
\begin{aligned}
64+
A(T) &= \alpha(T) = \frac{\partial \epsilon}{\partial T} \\
65+
B(T) &= \bar{\alpha}(T) = \frac{\epsilon(T)}{T-T_0} = \frac{\displaystyle \int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime}{T - T_0} \\
66+
C(T) &= \epsilon(T) = \int_{T_0}^T \alpha(T^\prime) \, dT^\prime
67+
\end{aligned}
68+
$$
69+
70+
Therefore,
71+
72+
i. $A(T)$ can be computed out of $C(T)$
73+
ii. $B(T)$ can be computed either out of $A(T)$ or $C(T)$
74+
iii. $C(T)$ can be computed out of $A(T)$
75+
terminal: |
76+
$ cat asme-expansion-table.dat
77+
# temp A B C
78+
20 21.7 21.7 0
79+
50 23.3 22.6 0.7
80+
75 23.9 23.1 1.3
81+
100 24.3 23.4 1.9
82+
125 24.7 23.7 2.5
83+
150 25.2 23.9 3.1
84+
175 25.7 24.2 3.7
85+
200 26.4 24.4 4.4
86+
225 27.0 24.7 5.1
87+
250 27.5 25.0 5.7
88+
275 27.7 25.2 6.4
89+
300 27.6 25.5 7.1
90+
325 27.1 25.6 7.8
91+
$ feenox asme-expansion.fee > asme-expansion-interpolation.dat
92+
$ pyxplot asme-expansion.ppl
93+
$
94+
...
95+
# just in case we wanted to interpolate with another method (linear, splines, etc.)
96+
DEFAULT_ARGUMENT_VALUE 1 steffen
97+
98+
# read columns from data file and interpolate
99+
# A is the instantaenous coefficient of thermal expansion x 10^-6 (mm/mm/ºC)
100+
FUNCTION A(T) FILE asme-expansion-table.dat COLUMNS 1 2 INTERPOLATION $1
101+
# B is the mean coefficient of thermal expansion x 10^-6 (mm/mm/ºC) in going
102+
# from 20ºC to indicated temperature
103+
FUNCTION B(T) FILE asme-expansion-table.dat COLUMNS 1 3 INTERPOLATION $1
104+
# C is the linear thermal expansion (mm/m) in going from 20ºC
105+
# to indicated temperature
106+
FUNCTION C(T) FILE asme-expansion-table.dat COLUMNS 1 4 INTERPOLATION $1
107+
108+
VAR T' # dummy variable for integration
109+
T0 = 20 # reference temperature
110+
T_min = vecmin(vec_A_T) # smallest argument of function A(T)
111+
T_max = vecmax(vec_A_T) # largest argument of function A(T)
112+
113+
# compute one column from another one
114+
A_fromC(T) := 1e3*derivative(C(T'), T', T)
115+
116+
B_fromA(T) := integral(A(T'), T', T0, T)/(T-T0)
117+
B_fromC(T) := 1e3*C(T)/(T-T0) # C is in mm/m, hence the 1e3
118+
119+
C_fromA(T) := 1e-3*integral(A(T'), T', T0, T)
120+
121+
# write interpolated results
122+
PRINT_FUNCTION A A_fromC B B_fromA B_fromC C C_fromA MIN T_min+1 MAX T_max-1 STEP 1
123+
---
124+
figures: |
125+
![Column $A(T)$ from $C(T)$](asme-expansion-A.svg){width=100%}
126+
![Column $B(T)$ from $A(T)$ and $C(T)$](asme-expansion-B.svg){width=100%}
127+
![Column $C(T)$ from $A(T)$](asme-expansion-C.svg){width=100%}
128+
129+
> The conclusion (see
130+
> [this](https://www.linkedin.com/pulse/accuracy-thermal-expansion-properties-asme-bpv-code-angus-ramsay/),
131+
> [this](https://www.seamplex.com/docs/SP-WA-17-TN-F38B-A.pdf) and
132+
> [this](https://www.linkedin.com/pulse/ansys-potential-issue-thermal-expansion-calculations-angus-ramsay/) reports)
133+
> is that values rounded to only one decimal value as presented in the ASME code
134+
> section II subsection D tables are not enough to satisfy the mathematical relationships
135+
> between the physical magnitudes related to thermal expansion properties of the materials listed.
136+
> Therefore, care has to be taken as which of the three columns is chosen when using the data
137+
> for actual thermo-mechanical numerical computations.
138+
> As an exercise, the reader is encouraged to try different interpolation algorithms to
139+
> see how the results change. _Spoiler alert_: they are also highly sensible to the interpolation method used
140+
> to “fill in” the gaps between the table values.
141+
...
Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,13 @@
11
# compute axis-angle from original and transformed vectors
22

3-
VECTOR orig SIZE 3 DATA 1 0 0
4-
VECTOR trans SIZE 3 DATA 1/sqrt(3) 1/sqrt(3) -1/sqrt(3)
3+
VECTOR orig[3] DATA 1 0 0
4+
VECTOR trans[3] DATA 1/sqrt(3) 1/sqrt(3) -1/sqrt(3)
55

6-
VECTOR a SIZE 3
6+
VECTOR a[3]
77

88
a[1] = orig[2]*trans[3]-orig[3]*trans[2]
99
a[2] = orig[3]*trans[1]-orig[1]*trans[3]
1010
a[3] = orig[1]*trans[2]-orig[2]*trans[1]
1111

1212
theta = acos(vecdot(orig, trans)/(vecnorm(orig)*vecnorm(trans)))
1313
PRINTF "(%g, %g, %f), %g" a(1) a(2) a(3) theta
14-
# PRINT "(" a(1) "," a(2) "," a(3) ")," theta SEP " "

examples/azmy-full.ppl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
set preamble "\usepackage{amsmath}"
1+
set preamble "\\usepackage{amsmath}"
22

33
set width 12*unit(cm)
44
set size ratio 9.0/16.0

0 commit comments

Comments
 (0)