179 lines
5.2 KiB
Python
179 lines
5.2 KiB
Python
#!/usr/bin/python3
|
|
|
|
from fenics import *
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from sys import argv
|
|
|
|
from scipy.constants import epsilon_0, mu_0
|
|
|
|
primeri = ('koax','izmaknjen','potlacen','kvadratni','stripline')
|
|
|
|
if argv[1] in primeri:
|
|
problem = argv[1]
|
|
else:
|
|
print('Neznan primer. Nadaljujem s koax.')
|
|
problem = 'koax'
|
|
shrani = True #spravi vse rezultate v fajle?
|
|
|
|
mreza = Mesh('geometrija/{}.xml'.format(problem))
|
|
domena = MeshFunction('size_t', mreza,
|
|
'geometrija/{}_physical_region.xml'.format(problem))
|
|
rob = MeshFunction('size_t', mreza,
|
|
'geometrija/{}_facet_region.xml'.format(problem))
|
|
|
|
dx = Measure('dx', domain=mreza, subdomain_data=domena)
|
|
ds = Measure('ds', domain=mreza, subdomain_data=rob)
|
|
|
|
V = FunctionSpace(mreza, 'P', 2)
|
|
|
|
napetost = 100
|
|
zila = DirichletBC(V, Constant(napetost), rob, 2)
|
|
oklop = DirichletBC(V, Constant(0), rob, 3)
|
|
r_pogoji = [zila, oklop]
|
|
eps_r = 2.55
|
|
|
|
u = TrialFunction(V)
|
|
v = TestFunction(V)
|
|
a = dot(grad(u), grad(v))*dx
|
|
L = Constant(0)*v*dx
|
|
u = Function(V)
|
|
solve(a == L, u, r_pogoji)
|
|
|
|
plot(mreza, linewidth=0.5)
|
|
p = plot(u)
|
|
cp = plt.colorbar(p)
|
|
cp.set_label('Elektricni potencial V (V)')
|
|
if shrani:
|
|
plt.savefig('slike/{}-V.pdf'.format(problem),bbox_inches='tight')
|
|
plt.show()
|
|
|
|
E = -grad(u)
|
|
plot(mreza, linewidth=0.5)
|
|
p = plot(E, zorder=3)
|
|
cp = plt.colorbar(p)
|
|
cp.set_label('Elektricno polje E (V/m)')
|
|
if shrani:
|
|
plt.savefig('slike/{}-E.pdf'.format(problem),bbox_inches='tight')
|
|
plt.show()
|
|
|
|
Z_w = np.sqrt(mu_0/epsilon_0/eps_r)
|
|
cross_z = lambda w: as_vector((-w[1],+w[0]))
|
|
H = cross_z(E)/Z_w
|
|
|
|
plot(mreza, linewidth=0.5)
|
|
p = plot(H, zorder=3)
|
|
cp = plt.colorbar(p)
|
|
cp.set_label('Magnetno polje H (A/m)')
|
|
if shrani:
|
|
plt.savefig('slike/{}-H.pdf'.format(problem),bbox_inches='tight')
|
|
plt.show()
|
|
|
|
nhat = FacetNormal(mreza)
|
|
cross_t = lambda w,x: w[0]*x[1]-w[1]*x[0]
|
|
k = cross_t(nhat, H)
|
|
tok = assemble(-k*ds(2))
|
|
Z0 = napetost/tok
|
|
|
|
C = epsilon_0*eps_r/napetost**2 * assemble(E**2 * dx(1))
|
|
L = mu_0/tok**2 * assemble(H**2 * dx(1))
|
|
|
|
print(tok, Z0, C, L)
|
|
if shrani:
|
|
with open('izracuni/{}-out.txt'.format(problem), 'w') as out:
|
|
out.write('$I$,{},A\n'.format(tok))
|
|
out.write('$Z_0$,{},\\ohm\n'.format(Z0))
|
|
out.write('$C$,{},pF\\per\\metre\n'.format(C*1e12))
|
|
out.write('$L$,{},nH\\per\\metre'.format(L*1e9))
|
|
|
|
# nadaljni del kode relevanten samo za koaksilani vod
|
|
# izračuna analitske vrednosti in jih za primer polj
|
|
# primerja z analitsko resitvijo
|
|
if problem == 'koax':
|
|
b = 3
|
|
a = 0.5
|
|
lg = np.log(b/a)
|
|
u_anal = lambda r: napetost/lg*np.log(r/a)
|
|
E_anal = lambda r: napetost/lg/r
|
|
H_anal = lambda r: E_anal(r)/Z_w
|
|
tok_anal = 2*np.pi*napetost/lg/Z_w
|
|
Z0_anal = 1/2/np.pi*Z_w*lg
|
|
C_anal = 2*np.pi*eps_r*epsilon_0/lg
|
|
L_anal = mu_0/2/np.pi*lg
|
|
|
|
print(tok_anal, tok)
|
|
print(Z0_anal, Z0)
|
|
print(C_anal, C)
|
|
print(L_anal, L)
|
|
if shrani:
|
|
with open('izracuni/{}-out.txt'.format(problem), 'w') as out:
|
|
out.write('velicina,FEM,analično,enote,napaka\n')
|
|
out.write('$I$,{},{},A,'.format(tok,tok_anal))
|
|
out.write('{}\n'.format((tok-tok_anal)/tok_anal*100))
|
|
out.write('$Z_0$,{},{},\\ohm,'.format(Z0,Z0_anal))
|
|
out.write('{}\n'.format((Z0-Z0_anal)/Z0_anal*100))
|
|
out.write('$C$,{},{},pF\\per\\metre,'.format(C*1e12,C_anal*1e12))
|
|
out.write('{}\n'.format((C-C_anal)/C_anal*100))
|
|
out.write('$L$,{},{},nH\\per\\metre,'.format(L*1e9,L_anal*1e9))
|
|
out.write('{}'.format((L-L_anal)/L_anal*100))
|
|
|
|
x = np.linspace(a, b, 100)
|
|
y = np.zeros(len(x))
|
|
|
|
primer_ana_e = E_anal(x)
|
|
primer_fem_e = np.array(list(map(project(E),zip(x,y))))[:,0]
|
|
plt.plot(x, primer_ana_e, label='Analitsko')
|
|
plt.plot(x, primer_fem_e, '--', label='FEniCS')
|
|
plt.xlabel('Polmer (m)')
|
|
plt.ylabel('Elektricno polje E (V/m)')
|
|
plt.grid()
|
|
plt.legend()
|
|
if shrani:
|
|
plt.savefig('slike/{}-anal-e.pdf'.format(problem),bbox_inches='tight')
|
|
plt.show()
|
|
|
|
plt.plot(x, (primer_ana_e-primer_fem_e)/primer_ana_e*100)
|
|
plt.xlabel('Polmer (m)')
|
|
plt.ylabel('Relativna napaka %')
|
|
plt.grid()
|
|
if shrani:
|
|
plt.savefig('slike/{}-napaka-e.pdf'.format(problem),bbox_inches='tight')
|
|
plt.show()
|
|
|
|
primer_ana_h = H_anal(x)
|
|
primer_fem_h = np.array(list(map(project(H),zip(x,y))))[:,1]
|
|
plt.plot(x, primer_ana_h, label='Analitsko')
|
|
plt.plot(x, primer_fem_h, '--', label='FEniCS')
|
|
plt.xlabel('Polmer (m)')
|
|
plt.ylabel('Magnetno polje H (V/m)')
|
|
plt.grid()
|
|
plt.legend()
|
|
if shrani:
|
|
plt.savefig('slike/{}-anal-h.pdf'.format(problem),bbox_inches='tight')
|
|
plt.show()
|
|
|
|
plt.plot(x, (primer_ana_h-primer_fem_h)/primer_ana_h*100)
|
|
plt.xlabel('Polmer (m)')
|
|
plt.ylabel('Relativna napaka %')
|
|
plt.grid()
|
|
if shrani:
|
|
plt.savefig('slike/{}-napaka-h.pdf'.format(problem),bbox_inches='tight')
|
|
plt.show()
|
|
|
|
if problem == 'stripline':
|
|
def We(W, b):
|
|
if W/b > 0.35:
|
|
return W
|
|
else:
|
|
We_b = W/b-(0.35-W/b)**2
|
|
return We_b*b
|
|
|
|
W = 0.8e-3
|
|
b = 1.6e-3
|
|
tmp = b/(We(W,b)+0.441*b)
|
|
Z0strip = 30*np.pi/np.sqrt(eps_r)*tmp
|
|
print(Z0strip)
|
|
print(Z0)
|
|
napaka = (Z0strip-Z0)/Z0*100
|
|
print(napaka)
|