Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Setting diffusion coefficient using GMSH face definition #844

Closed
Noxical opened this issue Mar 29, 2022 · 4 comments
Closed

Setting diffusion coefficient using GMSH face definition #844

Noxical opened this issue Mar 29, 2022 · 4 comments

Comments

@Noxical
Copy link

Noxical commented Mar 29, 2022

I want to set the diffusion coefficient / thermal conductivity in a wall based on the specific wall material isolation property.

I cant figure out how to do this using GMSH. My code is attached below. It works only by setting the D value using the X, Y faceCenters. However, I would prefer to use GMSH to specifiy the correct material within the line loop.

from fipy import *
import numpy as np

mesh = Gmsh2D('''
SetFactory("OpenCASCADE");

ns = 0.01;

x0 = 0;
x1 = 1;
x2 = 0.9;

x3 = 0.8;
x4 = 0.7;

y0 = 0;
y1 = 1;
y2 = 0.1;

y3 = 0.2;
y4 = 0.3;

Point(1) = {x0, y0, 0, ns};
Point(2) = {x1, y0, 0, ns};
Point(3) = {x1, y1, 0, ns};
Point(4) = {x2, y1, 0, ns};
Point(5) = {x2, y2, 0, ns};
Point(6) = {x0, y2, 0, ns};

Point(8) = {x3, y1, 0, ns};
Point(9) = {x3, y3, 0, ns};
Point(10) = {x0, y3, 0, ns};

Point(11) = {x4, y1, 0, ns};
Point(12) = {x4, y4, 0, ns};
Point(13) = {x0, y4, 0, ns};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 1};

Line(7) = {4, 8};
Line(8) = {8, 9};
Line(9) = {9, 10};
Line(10) = {10, 6};

Line(11) = {8, 11};
Line(12) = {11, 12};
Line(13) = {12, 13};
Line(14) = {13, 10};

Line Loop(1) = {1,2,3,7,11,12,13,14,10,6};
Plane Surface(1) = {1};
Physical Surface("wall") = {1};

Physical Line("shell") = {1,2,3,4,5,6};

Physical Line("outer") = {1,2};
Physical Line("inner") = {12,13};

'''% locals())

X,Y = mesh.faceCenters
T = CellVariable(name="isotherm", mesh=mesh, value=1.)

D=FaceVariable(mesh=mesh, value=1.)

mask = mesh.physicalFaces["shell"]

D.setValue(0.1, where=mask) # this doesn work

D.setValue(1, where= (X>0.9) + (Y<0.1))
D.setValue(2, where= (X<0.9) + (Y>0.1))

T.constrain(1400, mesh.physicalFaces['inner'] )
T.constrain(40, mesh.physicalFaces['outer'] )

eq = TransientTerm(var=T) == DiffusionTerm(D, var=T)

steps = 2 #time steps
dt = 600

for i in range(steps):

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

vi = Viewer(T)

@guyer
Copy link
Member

guyer commented Mar 29, 2022

Please see the use of Physical Surface in the documentation for Gmsh2D

@Noxical
Copy link
Author

Noxical commented Mar 29, 2022

I have looked at the documentation, but this only mentions the ability to change the cellVariables using gmsh physical surface definition, this works.
However there is no clear example of how to do this to the diffusion coefficients, which are faceVariables. I would expect i will work the same, but it doens't seem to be straight forward.

@guyer
Copy link
Member

guyer commented Mar 29, 2022

Gmsh defines elements and does not provide any ready mechanism to label the internal faces within a domain.

The simplest thing is to define your diffusion coefficient as a CellVariable and let FiPy interpolate it. By default, it will calculate the arithmetic mean at the faces, e.g.,

D=CellVariable(mesh=mesh)
D.setValue(1., where=mesh. physicalCells["zone1"])
D.setValue(2., where=mesh. physicalCells["zone1"])
eq = TransientTerm(var=T) == DiffusionTerm(D, var=T)

is completely equivalent to

eq = TransientTerm(var=T) == DiffusionTerm(D.arithmeticFaceValue, var=T)

You could specify a different interpolation scheme, if you wished:

eq = TransientTerm(var=T) == DiffusionTerm(D.harmonicFaceValue, var=T)

A more physically valid thing to do for diffusion coefficients might be

Ea = CellVariable(mesh=mesh)
Ea.setValue(Ea1, where=mesh. physicalCells["zone1"])
Ea.setValue(Ea2, where=mesh. physicalCells["zone2"])
D = D0 * numerix.exp(-(Ea / (k * T)).arithmeticFaceValue))
eq = TransientTerm(var=T) == DiffusionTerm(D, var=T)

@Noxical
Copy link
Author

Noxical commented Apr 1, 2022

Thanks. This solves the problem.

@Noxical Noxical closed this as completed Apr 1, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants