Meshes

Different meshes are available when solving the diffusion or any other PDE

For a 1D system, the diffusional flux is defined as

$ J = -D \frac{\partial x}{\partial z} $

Cartesian

$ \frac{\partial x}{\partial t} = -\frac{\partial}{\partial z} (J) $

Cylindrical

$ \frac{\partial x}{\partial t} = - \frac{1}{r} \frac{\partial}{\partial r} (r J) $

Spherical

$ \frac{\partial x}{\partial t} = - \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 J) $

To show the different meshes, we could create a simple heat transfer model to apply the meshes to and compare to analytical solutions

 1  import numpy as np
 2  import matplotlib.pyplot as plt
 3  
 4  from kawin.diffusion.mesh import AbstractMesh, DiffusionPair, MeshData
 5  from kawin.GenericModel import GenericModel
 6  
 7  class ThermalModel(GenericModel):
 8      def __init__(self, mesh: AbstractMesh, k, rho, cp):
 9          super().__init__()
10          self.mesh = mesh
11          self.T0 = 298.15
12          self.k = k
13          self.rho = rho
14          self.cp = cp
15          self.D = k / rho / cp
16  
17          # Assume mesh is in temperature and convert to enthalpy
18          # Since self.mesh.y could be in any shape, we want to flatten it to (N,1)
19          T = self.mesh.flattenResponse(self.mesh.y)
20          H = self.rho * self.cp * (T - self.T0)
21          self.data = MeshData(self.mesh, record=1)
22          self.data.setResponseAtN(H)
23  
24      def getCurrentX(self):
25          return [self.data.currentY]
26      
27      def getDt(self, dxdt):
28          return 0.4*self.mesh.dz**2 / self.D / self.mesh.dims
29      
30      def getdXdt(self, t, x):
31          H_d, z_d = self.mesh.getDiffusivityCoordinates(x[0])
32          ks = self.k*np.ones(H_d.shape)
33          H_r, z = self.mesh.getResponseCoordinates(x[0])
34          T_r = self._enthalpyToTemp(H_r)
35          return [self.mesh.computedXdt([DiffusionPair(diffusivity=ks, response=T_r)])]
36      
37      def postProcess(self, time, x):
38          super().postProcess(time, x)
39          self.data.record(time, x[0])
40          return self.getCurrentX(), False
41      
42      def postSolve(self):
43          self.data.finalize()
44      
45      def _enthalpyToTemp(self, H):
46          return H/self.rho/self.cp + self.T0
47      
48      @property
49      def z(self):
50          return self.mesh.z
51      
52      @property
53      def T(self):
54          H = self.data.currentY
55          return self.mesh.unflattenResponse(self._enthalpyToTemp(H))
56      

 57  from kawin.diffusion.mesh import Cartesian1D, Cylindrical1D, Spherical1D
 58  from kawin.diffusion.mesh import MixedBoundary1D, ProfileBuilder, ConstantProfile
 59  
 60  xlim = [0.1, 0.2]
 61  ylim = [300, 1000]
 62  N = 100
 63  responses = ['T']
 64  k = 13 # W/m-K
 65  rho = 7870 # kg/m3
 66  cp = 490 # J/kg-K
 67  
 68  
 69  # Dirichlet boundary conditions on ends of mesh
 70  bc = MixedBoundary1D(responses)
 71  bc.setLBC('T', 'dirichlet', ylim[0])
 72  bc.setRBC('T', 'dirichlet', ylim[1])
 73  
 74  profile = ProfileBuilder()
 75  profile.addBuildStep(ConstantProfile(ylim[0]), 'T')
 76  
 77  meshCartesian = Cartesian1D(responses, xlim, N)
 78  meshCylindrical = Cylindrical1D(responses, xlim, N)
 79  meshSpherical = Spherical1D(responses, xlim, N)
 80  meshes = [meshCartesian, meshCylindrical, meshSpherical]
 81  labels = ['Cartesian', 'Cylindrical', 'Spherical']
 82  for i, mesh in enumerate(meshes):
 83      mesh.setResponseProfile(profile, bc)
 84      model = ThermalModel(mesh, k=k, rho=rho, cp=cp)
 85      model.solve(3600)
 86      plt.plot(model.z, model.T, label=labels[i])
 87  
 88  # Solution for cartesian: x = A*z + B
 89  # Solution for cylindrical: x = A*ln(r) + B
 90  # Solution for spherical: x = A/r + B
 91  A = (ylim[1] - ylim[0]) / (xlim[1] - xlim[0])
 92  B = ylim[0] - A*xlim[0]
 93  cart_ideal = A*meshCartesian.z + B
 94  
 95  A = (ylim[1] - ylim[0]) / np.log(xlim[1] / xlim[0])
 96  B = ylim[0] - A*np.log(xlim[0])
 97  cyl_ideal = A*np.log(meshCylindrical.z) + B
 98  
 99  A = (ylim[1] - ylim[0]) / (1/xlim[1] - 1/xlim[0])
100  B = ylim[0] - A/xlim[0]
101  sph_ideal = A/meshSpherical.z + B
102  plt.plot(meshCartesian.z, cart_ideal, color='k', linestyle='--', linewidth=0.75, label='Analytical')
103  plt.plot(meshCylindrical.z, cyl_ideal, color='k', linestyle='--', linewidth=0.75, label='Analytical')
104  plt.plot(meshSpherical.z, sph_ideal, color='k', linestyle='--', linewidth=0.75, label='Analytical')
105  plt.xlim(xlim)
106  plt.ylim(ylim)
107  plt.xlabel('z or r (m)')
108  plt.ylabel('Temperature (K)')
109  plt.legend(labels=labels+['Analytical'])
110  plt.show()

Periodic boundary conditions are also available, but is only valid for cartesian meshes

111  from kawin.diffusion.mesh import PeriodicBoundary1D, BoundedRectangleProfile
112  
113  profile = ProfileBuilder()
114  profile.addBuildStep(ConstantProfile(300))
115  profile.addBuildStep(BoundedRectangleProfile(0, 0.25, 700, 0))
116  profile.addBuildStep(BoundedRectangleProfile(0.5, 0.75, 700, 0))
117  
118  pb = PeriodicBoundary1D()
119  mesh = Cartesian1D(responses=['T'], zlim=[0,1], N=100)
120  mesh.setResponseProfile(profile, pb)
121  
122  model = ThermalModel(mesh, k=k, rho=rho, cp=cp)
123  
124  fig, ax = plt.subplots(1,2, figsize=(8,3))
125  ax[0].plot(model.z, model.T)
126  ax[0].set_xlim([0,1])
127  ax[0].set_ylim([200, 1100])
128  
129  model.solve(1000)
130  ax[1].plot(model.z, model.T)
131  
132  for i in range(2):
133      ax[i].set_xlim([0,1])
134      ax[i].set_ylim([200, 1100])
135      ax[i].set_xlabel('Distance (m)')
136      ax[i].set_ylabel('Temperature (K)')
137  fig.tight_layout()
138  plt.show()
← Nucleation Rate
kawin Framework →