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)) 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=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() ``` On Mar 4, 2021, at 8:02 PM, VICTOR PEREZ <viti...@gmail.com<mailto:viti...@gmail.com>> wrote: 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() -- To unsubscribe from this group, send email to fipy+unsubscr...@list.nist.gov<mailto:fipy+unsubscr...@list.nist.gov> View this message at https://list.nist.gov/fipy <Selection_004.png> -- To unsubscribe from this group, send email to fipy+unsubscr...@list.nist.gov View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to fipy+unsubscr...@list.nist.gov.
