Hello everyone,
I would like to impose a Dirichlet boundary condition over a surface in 3D.
My geometry is simply a cylinder and the surface on which I would like to impose the Dirichlet condition is a disk. I used a 2D mesh (Th2D in the script below) to define the surface….but it is clearly the bad way to do it. I believe I should find the points of the 3D mesh included in my surface and affect them individually the Dirichlet condition.
Here is a part of the script… The problematic condition is the last line of the script.
Thank you in advance for your help.
int Cyz=99;
real R1 = 0.02; // Rayon du duramen
real R2 = 0.04; // Rayon de l aubier
real R4 = 0.05; // Rayon total
real H = 0.2; // Hauteur du cylindre
real kd = 5, ka = 5, ke = 20; // Conductivites thermiques
real T0 = 50.; // Temp initiale
real T = 0.002, dt = 0.0002 ;
int nt = floor(T/dt);
real[int] temps(nt);
int nx=30, ny = 30, nz = 30;
real [int] y2d(ny), z2d(nz);
real[int] xnz*nt);
real[int] uu(nz*nt);
real[int,int] u2dyz(ny,nz);
// Noms des régions
int regDuramen = 1;
int regAubier = 2;
int regEcorce = 3;
// contours
border b0(t=0, 2*pi){x=R1*cos(t); y=R1*sin(t); z = 0;}
border b1(t=0, 2*pi){x=R2*cos(t); y=R2*sin(t);z = 0;}
border b2(t=0, 2*pi){x=R4*cos(t); y=R4*sin(t);z = 0;}
// 2D mesh x-y
mesh Th2D = buildmesh(b0(nx) + b1(nx) + b2(nx));
plot(Th2D, wait=true, cmm=“Maillage 2D”);
// 3D mesh
mesh3 Th3D = buildlayers(Th2D, ny, zbound = [-H/2, H/2]);
plot(Th3D, wait=true, cmm=“Maillage 3D”);
real Duramen = Th3D(R1/2,0,0).region;
real Aubier = Th3D(R1+R2/2,0,0).region;
real Ecorce = Th3D(R1+R2+R4/2,0,0).region;
fespace Vh(Th3D, P1);
Vh u=T0, v, uold;
Vh k = kd*(region==Duramen)+ka*(region==Aubier)+ke*(region==Ecorce);
fespace Vh2(Th2D, P2);
Vh2 u2, u2e;
problem heat(u, v) =
int3d(Th3D)(u*v/dt + k*(du)*dv) + dy(u)*dy(v) + dz(u)*dz(v)))
- int3d(Th3D)(uold*v/dt)
+ on(Ecorce, u=10);
// + on(Th2D, u=10); // CL appliquee a la surface du disque (?)
