Temperature dependent material data from file

1

Hi,

for my time dependent heat convection-diffusion equation I have some material data from files: density, conductivity and specific heat. All depend on the current local temperature, i.e. local solution of the PDE. I tried the following, that works well with constant material data (as Constant() or Expression()):

...

# material data from file, int(temperature) == index
cpX    = np.genfromtxt('analyse/cpX.ana')[:,1]
rhoX   = np.genfromtxt('analyse/densityX.ana')[:,1]
lmbdaX = np.genfromtxt('analyse/conductX.ana')[:,1]

...

class userMaterial(UserExpression):

    def __init__(self, temp, valArray, **kwargs):
        super().__init__(kwargs)
        self.temp     = temp
        self.valArray = valArray

    def eval(self, value, x):
        intT     = max(1,min(int(self.temp(x)),1600))
        value[0] = self.valArray[intT]

    def value_shape(self):
        return ()

...

U = Function(V)
U_data = np.array(U.vector())*0+100
U.vector().set_local(U_data)

lambdaU = userMaterial(temp=U,valArray=lmbdaX,degree=0)
rhoU    = userMaterial(temp=U,valArray=rhoX,degree=0)
cpU     = userMaterial(temp=U,valArray=cpX,degree=0)

...

E_thd = lambdaU/(rhoU*cpU) # thermal diffusivity
E_rcp = 1.0/(rhoU*cpU)     # 1 / rho * cp

...

for i in bcond:
    if 'Dirichlet' in bcond[i]:
        bc = DirichletBC(V,bcond[i]['Dirichlet'],boundary_markers,i)
        bcs.append(bc)

    elif 'Neumann' in bcond[i]:
        if bcond[i]['Neumann'] != 0:
            g = bcond[i]['Neumann']
            integrals_N.append(E_dt*E_rcp*g*v*ds(i))

    elif 'Robin' in bcond[i]:
        r,s = bcond[i]['Robin']
        integrals_R_a.append(E_dt*E_rcp*r*u*v*ds(i))
        integrals_R_L.append(E_dt*E_rcp*r*s*v*ds(i))

...

#===========================================================
# PDE
#
#  rho*cp*grad(T,t)+rho*cp*w*grad(T,x)=div(lambda*grad(T,x),x)
#
#===========================================================

# weak problem formulation

F    =   u*v*dx + E_dt*E_thd*dot(grad(u),grad(v))*dx - u_n*v*dx \
       + E_dt*inner(E_w,grad(u))*v*dx \
       + sum(integrals_R_a) + sum(integrals_N) - sum(integrals_R_L)

a, L = lhs(F), rhs(F)

...

t = 0
while t < tFinal:

    # Update current time
    t += dt
    print(t,"s /",tFinal,'s')

    # Compute solution
    solve(a == L, U, bcs)

    # Save to file and plot solution
    vtkfile << (U, t)
    # plot(u)

    # Update previous solution
    u_n.assign(U)

    # update material data
    lambdaU.temp = U
    rhoU.temp    = U
    cpU.temp     = U

...

Unfortunately, this doesn’t work and is also extremely slow.

How can measured, temperature dependend table data be applied to the PDE?

2

Write it into a function once, save it as a checkpoint using XDMFFile.write_checkpoint such that it can be reloaded in parallel with basically no cost using XDMFFile.read_checkpoint for later experiments.

2 Likes
3

Wow, this sounds as you know what you do. But I don’t understand much, since I’m new to FEniCS :frowning:

Could you add some pseudo code?

What I thought meanwhile is to use a MeshFunction and update it only on the vertices to avoid the expencive call to temp(x). Would this be possible?

4

See for instance: Read mesh from XDMF file (write_checkpoint) - #3 by dokken

5

Hmmm … since the material data change after each time step, I don’t think that storing the data once is sufficient. I just learn to understand how Function and UserExpression operates. I now switched over to:

...

lambdaU = Function(V)
rhoU    = Function(V)
cpiU    = Function(V)

lambdaU.vector()[:]  = updateMatdata(lambdaX,U)
rhoU.vector()[:]        = updateMatdata(rhoX,U)
cpiU.vector()[:]        = updateMatdata(cpiX,U)

...

which avoids calling a function at every time step and I think this is at the moment ok for me. updateMatdata ist more or less comparable to class userMaterial(UserExpression). I hope I use the Function statement the right way.