Compare commits

...

2 Commits

Author SHA1 Message Date
Moonwalker777 150d42b0d8 Ignore 2022-05-19 14:53:08 +02:00
Moonwalker777 4a0865c53d Exercise 4 Assignment Data
Handed out Information and Chapter 5 of mentioned book.
2022-05-19 14:52:41 +02:00
7 changed files with 190 additions and 0 deletions

BIN
.DS_Store vendored Normal file

Binary file not shown.

2
.gitignore vendored
View File

@ -17,3 +17,5 @@ Exercise_02/**.txt
*.log
*.out
Exercise_03/.DS_Store
.DS_Store
.DS_Store

107
Exercise_03/NSSC_1.m Executable file
View File

@ -0,0 +1,107 @@
%% Euler forward Dirichlet
L = 1; % domain size (in the lecture notes this is denote h)
T = 2; % time limit (max time)
f = @(x,t) 0*x.*t; % rhs of the more general equation `u_t - d u_xx = f`
c1 = @(t) 1+0*t; % _right_ boundary condition
c2 = @(t) 0*t; % _left_ boundary condition
u0 = @(x) 0*x; % initial values
D = 0.5; % diffusion parameter `d` in `u_t - d u_xx = f`
%uex = @(x,t) cos(x).*exp(t);
N = 10; % nr. of _space_ discretization points
K = 200; % nr. of _time_ discretization points
[x, t, u] = Dirichlet_EA(L, N, T, K, c1, c2, f, u0, D);
% Report stability condition `D Delta T / (Delta x)^2 > 0.5`
Delta_T = T / K;
Delta_x = L / N;
d = D * Delta_T / Delta_x^2;
fprintf("Stability Condition: 0.5 >= D * Delta_T / Delta_x^2 = %f\n", d)
if d > 0.5
fprintf("-> NOT Stable\n")
else
fprintf("-> Stable\n")
end
figure(1)
for ii = 1:K+1 % iterates time
hold on
plot(x, u(:, ii)');
xlim([0 L])
pause(0.05);
hold off
end
% 3D plot of space solution over time
space = linspace(0,L,101);
time = linspace(0,T,201);
[xx,yy] = meshgrid(time,space);
%exsol = uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)
%% Eulero forward Mixed BC
L=2*pi;
T=5;
f=@(x,t) 0*x.*t;
c1=@(t) 1+0*t;
c2=0;
u0=@(x) 0*x;
D=1.1;
%uex=@(x,t) cos(x).*exp(t);
N=25;
K=200;
[x,t,u]=Mixed_EA(L,N,T,K,c1,c2,f,u0,D);
figure(1)
for ii=1:K+1
plot(x,u(:,ii)');
xlim([0 L])
ylim([0 1.5])
pause(0.02);
end
space=linspace(0,L,101);
time=linspace(0,T,201);
[xx,yy]=meshgrid(time,space);
%exsol=uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)
%% Eulero Backward Mixed BC
L=2*pi;
T=5;
f=@(x,t) 0*x.*t;
c1=@(t) 1+0*t;
c2=0;
u0=@(x) 0*x;
D=1;
%uex=@(x,t) cos(x).*exp(t);
N=25;
K=200;
[x,t,u]=Mixed_EI(L,N,T,K,c1,c2,f,u0,D);
figure(1)
for ii=1:K+1
plot(x,u(:,ii)');
xlim([0 L])
ylim([0 1.5])
pause(0.02);
end
space=linspace(0,L,101);
time=linspace(0,T,201);
[xx,yy]=meshgrid(time,space);
%exsol=uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)

BIN
Exercise_04/Zienkiewicz_Chapter_5.pdf (Stored with Git LFS) Normal file

Binary file not shown.

View File

@ -0,0 +1,39 @@
"""
NSSCII - FEM.
Input to be parsed through.
SI-units to be used:
+ T in K
+ L in m
+ k in W/(mK)
+ q in W/m^2 - ad Neumann
+ P in W - ad nodal forces
"""
# Group number.
groupnr = 1
# Length in x- and y-direction.
L = 0.01
# Thickness (z-direction).
hz = 0.0005
# Thermal conductivity (k=k_xx=k_yy, k_xy = 0.).
k = 429.
# Factor c for modifying thermal conductivity k for
# elements in elements_to_be_modified.
c = 10.
# Elements to be modified.
elements_to_be_modified = [
41-47,
59-63,
77-79,
95
]
# Boundary conditions.
q(y=0) = 2000000.
T(y=L) = 293.

BIN
Exercise_04/nsscii-exc4-fem.pdf (Stored with Git LFS) Normal file

Binary file not shown.

36
Exercise_04/print_HTP.py Normal file
View File

@ -0,0 +1,36 @@
def print_HTP(H, T, P, filename="output.txt"):
"""
Print matrices to .txt-file (name of file = filename).
H... overall assembled stiffness matrix
T... nodal temperature vector
P... nodal force vector
Make sure, that your system of equations is sorted by
ascending node numbers, i.e., N1 N2 ... N100.
"""
F = open(filename, 'w')
F.write("Stiffness matrix H: \n")
for row in H:
for col in row:
outline = "{0:+8.4e},".format(col)
F.write("{0:11s}".format(str(outline)))
F.write("\n")
F.write("Temperature T: \n")
for row in T:
for col in row:
outline = "{0:+8.4e},".format(col)
F.write("{0:11s} \n".format(str(outline)))
F.write("Force vector P: \n")
for row in P:
for col in row:
outline = "{0:+8.4e},".format(col)
F.write("{0:11s} \n".format(str(outline)))
F.close()
return None