Thanks, it worked for me. El viernes, 5 de marzo de 2021 a las 11:19:08 UTC-5, [email protected] escribió:
> 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 <[email protected]> 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 [email protected] > > View this message at https://list.nist.gov/fipy > <Selection_004.png> > > > -- To unsubscribe from this group, send email to [email protected] View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to [email protected].
