How do I get internal boundary conditions? letting one square remain within another.

[VICTOR PEREZ]

for the moment I have the conditions of the edge of the square. but the intense conditions could not. In the image I show how the internal border conditions should be.

#!/usr/bin/env python3

'''

include lib

'''

from fipy import *

import time

nx = ny = 30.

nz = 6.

dx = dy = dz = 1.

L1 = dx*nx

L2 = dy*ny

L3 = dz*nz

mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)

x, y, z = mesh.cellCenters

#mask = mesh.cellCenters

#mask = mesh.facesBack == 0.5

#print("mask: ", mask)

#mask0 = mesh.cellCenters.value[0] > True

#mask1 = mesh.cellCenters.value[1]

#mask2 = mesh.cellCenters.value[2]

#print("mask0: ", mask0)

#print("mask1: ", mask1)

#print("mask2: ", mask2)

T = CellVariable(mesh=mesh, name = r"TEMPERATURE", value=23.0)

#T[mask] = 30.

'''

value_constranin = [24, 24]

value_constranin_C1 = value_constranin[1]

print(value_constranin_C1)

'''

value_C1, value_C2 = 23, 23

X, Y, Z = mesh.faceCenters

C1 = ((mesh.facesLeft & (Y < L1))

| (mesh.facesTop & (X < L1)))

C2 = ((mesh.facesRight & (Y < L2))

| (mesh.facesBottom & (X < L2)))

T.constrain(value_C1, C1)

T.constrain(value_C2, C2)

T.constrain(60, mesh.facesBack)

'''

START: Propiedades termicas del material (Ceramica)

'''

k = 1.3 # [W/m·K] [J/m.K.s] conductividad termica

rho = 2300 # [kg/m**3] desidad

cp = 840 # [J/(kg·K)] calor especifico

alpha0 = k/rho/cp#[m2/s] difusividad termica

print(alpha0)

alpha = FaceVariable(mesh=mesh, value=alpha0)

'''

END: Propiedades termicas del material (Ceramica)

'''

largeValue = 1e+7

print("Alpha: ",alpha*largeValue)

D = 0.1

eq = TransientTerm() == DiffusionTerm(coeff=largeValue*alpha)

# - ImplicitDiffusionTerm(mask*largeValue)

'''

START: Leyenda de la grafica

'''

viewer = MayaviClient(vars=T,

limits={'T min': 20, ' T max': 70},

datamin=23, datamax=70,

titel="TEMPERATURE")

'''

END: Leyenda de la grafica

'''

from builtins import range

timeStepDuration =1 #10*0.9*dx**2/(2.*D)

steps = 10

t = 0.

start_time = time.time()

#print("start: ", start_time)

print("start: t1 = ", timeStepDuration)

for step in range(steps):

eq.solve(var=T, dt=timeStepDuration)

t = timeStepDuration + t

#print("Time: ",step, time.time() - start_time)

print("Time: t ", t)

viewer.plot()

[Guyer, Jonathan E. Dr. (Fed)]

Victor -

An internal constraint, as explained at

   https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#internal-fixed-value

is applied by a combination of implicit and explicit sources, not diffusion terms, e.g.,

eq = (TransientTerm() == DiffusionTerm(coeff=alpha0)

      - ImplicitSourceTerm(coeff=mask*largeValue) 

      + mask * largeValue * 30.)

The idea is that, where `mask` is True, this equation effectively becomes

eq = (0 == - ImplicitSourceTerm(coeff=mask*largeValue) + mask * largeValue * 30.)

or

eq = (ImplicitSourceTerm(coeff=mask) == mask * 30.)

which is an implicit way to write

T = 30

when `mask` is True. 

When `mask` is False, then these sources go away and the transient diffusion equation dominates.

A few other changes are needed to get this working:

- Although `mesh.cellCenters` is a rank-1 `CellVariable` and `mesh.cellCenters[0]` is a rank-0 `CellVariable`, unfortunately, [tuple-unpacking](https://www.w3schools.com/python/python_tuples_unpack.asp) does *not* result in x, y, and z containing `CellVariable` objects. This confuses FiPY later. Instead, you can write

    x, y, z = mesh.cellCenters[0], mesh.cellCenters[1], mesh.cellCenters[2]

- You need a mask. To mask out the middle 1/3 cube of your domain, I would write

  mask = ((x > L1 / 3) & (x < 2 * L1 / 3)
          & (y > L2 / 3) & (y < 2 * L2 / 3)
          & (z > L3 / 3) & (z < 2 * L3 / 3))

- Your time step is not consistent with your diffusivity, I would try

  timeStepDuration = 10*0.9*dx**2/(2.*alpha0)

Putting it all together:

```

#!/usr/bin/env python3

'''

include lib

'''

from fipy import *

import time

nx = ny = 30.

nz = 6.

dx = dy = dz = 1.

L1 = dx*nx

L2 = dy*ny

L3 = dz*nz

mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)

x, y, z = mesh.cellCenters[0], mesh.cellCenters[1], mesh.cellCenters[2]

mask = ((x > L1 / 3) & (x < 2 * L1 / 3)

        & (y > L2 / 3) & (y < 2 * L2 / 3)

        & (z > L3 / 3) & (z < 2 * L3 / 3))

eq = (TransientTerm() == DiffusionTerm(coeff=alpha)

      - ImplicitSourceTerm(mask*largeValue)

      + mask*largeValue*30.)

'''

START: Leyenda de la grafica

'''

viewer = MayaviClient(vars=T,

limits={'T min': 20, ' T max': 70},

datamin=23, datamax=70,

title="TEMPERATURE")

'''

END: Leyenda de la grafica

'''

from builtins import range

timeStepDuration = 10*0.9*dx**2/(2.*alpha0)

steps = 10

t = 0.

start_time = time.time()

#print("start: ", start_time)

print("start: t1 = ", timeStepDuration)

for step in range(steps):

    eq.solve(var=T, dt=timeStepDuration)

    t = timeStepDuration + t

    #print("Time: ",step, time.time() - start_time)

    print("Time: t ", t)

    viewer.plot()

```

--
To unsubscribe from this group, send email to fipy+uns...@list.nist.gov
 
View this message at https://list.nist.gov/fipy
<Selection_004.png>

[VICTOR PEREZ]

Thanks, it worked for me.