forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpnpLayeredWallqboth.py
143 lines (130 loc) · 5.62 KB
/
pnpLayeredWallqboth.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
# Finite Element Modeling with Abaqus and Python for Thermal and
# Stress Analysis
# (C) 2017, Petr Krysl
#
"""
Layered wall. Temperature boundary condition on the left, and heat flux
boundary condition on the right.
"""
from numpy import array as array
from numpy import zeros as zeros
from numpy import arange as arange
from numpy import dot as dot
from numpy import linalg as linalg
from numpy import vstack as vstack
from numpy import diff as diff
import matplotlib.pyplot as plt
Dz = 1.0 # Thickness of the slice (it cancels out in the end)
kappaL = 0.05 # thermal conductivity, layer on the left
kappaR = 1.8 # thermal conductivity, layer on the right
# Coordinates of nodes
x = array([[0, 0], [0.07, 0], [0.3, 0], [0, 0.1], [0.07, 0.1], [0.3, 0.1]])
N = 6 # total number of nodes
N_f = 6 # total number of free degrees of freedom
# Mapping from nodes to degrees of freedom
node2dof=array([6, 1, 3, 5, 2, 4])
# Connectivity of the left-region triangles (polyurethane)
connL = array([[4, 1, 5], [2, 5, 1]])
# Connectivity of the right-region triangles (concrete)
connR = array([[6, 5, 2], [2, 3, 6]])
# Connectivity of the boundary L2 element on the left
connLbdry = array([[4, 1]])
qnbarL = +30
# Connectivity of the boundary L2 element on the right
connRbdry = array([[3, 6]])
qnbarR = -30
# gradients of the basis functions with respect to the param. coordinates
gradNpar = array([[-1, -1], [1, 0], [0, 1]])
Kg = zeros((N_f, N_f)) # allocate the global conductivity matrix
Lg = zeros(N_f).reshape(N_f, 1) # allocate the global heat loads vector
# Loop over the triangles in the mesh, left domain
kappa = kappaL
zconn = array(connL)-1
for j in arange(zconn.shape[0]):
J = dot(x[zconn[j, :], :].T, gradNpar) # compute the Jacobian matrix
gradN = dot(gradNpar, linalg.inv(J)) # compute the x,y grads of the b. funcs
Ke = (kappa*Dz*linalg.det(J)/2)*dot(gradN, gradN.T) # elementwise matrix
# Element degree-of-freedom array, converted to zero base
zedof = array(node2dof[zconn[j, :]])-1
# Assemble elementwise conductivity matrix
for ro in arange(len(zedof)):
for co in arange(len(zedof)):
if (zedof[ro] < N_f) and (zedof[co] < N_f):
Kg[zedof[ro], zedof[co]] = Kg[zedof[ro], zedof[co]] + Ke[ro, co]
# Compute the elementwise load vector due to prescribed temperatures
LTe = zeros(3).reshape(3, 1)
for co in arange(len(zedof)):
if zedof[co] >= N_f:
LTe= LTe - Ke[:, co].reshape(3, 1)*Tfix[zedof[co]] # the "-" sign!
# Assemble elementwise heat load vector
for ro in arange(len(zedof)):
if zedof[ro] < N_f:
Lg[zedof[ro]] = Lg[zedof[ro]] + LTe[ro]
# Loop over the triangles in the mesh, right domain
kappa = kappaR
zconn = array(connR)-1
for j in arange(zconn.shape[0]):
J = dot(x[zconn[j, :], :].T, gradNpar) # compute the Jacobian matrix
gradN = dot(gradNpar, linalg.inv(J)) # compute the x,y grads of the b. funcs
Ke = (kappa*Dz*linalg.det(J)/2)*dot(gradN, gradN.T) # elementwise matrix
# Element degree-of-freedom array, converted to zero base
zedof = array(node2dof[zconn[j, :]])-1
# Assemble elementwise conductivity matrix
for ro in arange(len(zedof)):
for co in arange(len(zedof)):
if (zedof[ro] < N_f) and (zedof[co] < N_f):
Kg[zedof[ro], zedof[co]] = Kg[zedof[ro], zedof[co]] + Ke[ro, co]
# Compute the elementwise load vector due to prescribed temperatures
LTe = zeros(3).reshape(3, 1)
for co in arange(len(zedof)):
if zedof[co] >= N_f:
LTe= LTe - Ke[:, co].reshape(3, 1)*Tfix[zedof[co]] # the "-" sign!
# Assemble elementwise heat load vector
for ro in arange(len(zedof)):
if zedof[ro] < N_f:
Lg[zedof[ro]] = Lg[zedof[ro]] + LTe[ro]
# Boundary element heat flux term, left-hand side
qnbar = qnbarL
zconn = array(connLbdry)-1
for j in arange(zconn.shape[0]):
he = linalg.norm(diff(x[zconn[j, :], :], axis=0))
# Element degree-of-freedom array, converted to zero base
zedof = array(node2dof[zconn[j, :]])-1
Leq = -qnbar*he*Dz/2*array([[1], [1]])
# Assemble elementwise heat load vector
for ro in arange(len(zedof)):
if (zedof[ro] < N_f):
Lg[zedof[ro]] = Lg[zedof[ro]] + Leq[ro]
# Boundary element heat flux term, right-hand side
qnbar=qnbarR
zconn=array(connRbdry)-1
for j in arange(zconn.shape[0]):
he = linalg.norm(diff(x[zconn[j, :], :], axis=0))
# Element degree-of-freedom array, converted to zero base
zedof = array(node2dof[zconn[j, :]])-1
Leq = -qnbar*he*Dz/2*array([[1], [1]])
# Assemble elementwise heat load vector
for ro in arange(len(zedof)):
if (zedof[ro] < N_f):
Lg[zedof[ro]] = Lg[zedof[ro]] + Leq[ro]
# Solve for the global temperatures at the free degrees of freedom
Tg = linalg.solve(Kg, Lg)
print('Tg=', Tg)
# Set up an array of the temperatures at all the degrees of freedom
T = Tfix.copy()
T[0:N_f] = Tg # insert the computed values of the free degrees of freedom
# Plotting
plt.figure()
plt.gca().set_aspect('equal')
# setup three 1-d arrays for the x-coord, the y-coord, and the z-coord
xs = x[:, 0].reshape(N,)# one value per node
ys = x[:, 1].reshape(N,)# one value per node
ix = node2dof[arange(N)]-1
zs = (T[ix]).reshape(N,)# one value per node
triangles = vstack((connL-1, connR-1))# triangles are defined by the conn arrays
plt.tricontourf(xs, ys, triangles, zs)
plt.colorbar()
plt.title('Contour plot of temperature')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.show()