|
Dear Andy,
I just add some commentaries on your code.
Yves.
Le 19/06/2017 à 21:27, Yu (Andy) Huang a écrit :
Dear getFEM users,
Thanks for you previous help on my silly questions! Now I
manage to simulate the electric field on a toy sphere with two
layers, each layer having a different electrical conductivity.
I'm just not sure if I did it properly, because when I compare
the results to those I got from Abaqus (a commercial FEM
solver), I see some difference that I don't understand (see
attached screenshots). I used the same mesh, with the
same boundary condition and same conductivity values,
but the distribution and absolute values of voltage and field
are all different between getFEM and Abaqus. My major concern
is the way I coded the boundary conditions. The problem I'm
solving is a Laplacian equation of electric potential u, with
the following BC:
1) injecting 1 A/m^2 current density at the north pole of
the sphere: -n.J = 1
2) ground at the south pole: V = 0
3) insulation at all outer boundary: n.J = 0
4) continuity for inner boundary: n.(J1 - J2)
= 0
I only coded explicitly (1) and (2) so not sure if it's
good enough. I put part of my code below (Matlab code). Any
advice on the code is much appreciated!
% ===================================================
mesh = gfMesh('import','gmsh', 'toySphere.msh');
mfu = gf_mesh_fem(mesh, 1); % scalar-field (electric
potential)
mfE = gf_mesh_fem(mesh, 3); % 3d vector-field (electric
field)
gf_mesh_fem_set(mfu, 'fem', gf_fem('FEM_PK(3,1)')); %
P1 Lagrange
gf_mesh_fem_set(mfE, 'fem', gf_fem('FEM_PK(3,1)')); %
P1 Lagrange
mim = gf_mesh_im(mesh, gf_integ('IM_TETRAHEDRON(1)'));
% integration method
The integration method is used for both the volumic term and the
boundary conditions, so it has to be of order two
-> mim = gf_mesh_im(mesh, gf_integ('IM_TETRAHEDRON(2)'));
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mfu);
rid = gf_mesh_get(mesh,'regions');
reg1 = gf_mesh_get(mesh, 'region', rid(1));
reg2 = gf_mesh_get(mesh, 'region', rid(2));
region1 = 1;
region2 = 2;
gf_mesh_set(mesh, 'region', region1, reg1);
gf_mesh_set(mesh, 'region', region2, reg2);
Did you check that your regions were ok (for instance with
gf_plot_mesh ) ?
% governing equation and conductivities
sigma = [0.276;0.126]; % conductivity values for the two
layers in the sphere
gf_model_set(md, 'add linear generic assembly brick',
mim, [num2str(sigma(1)) '*(Grad_u.Grad_Test_u)'],region1);
gf_model_set(md, 'add linear generic assembly brick',
mim, [num2str(sigma(2)) '*(Grad_u.Grad_Test_u)'],region2);
fb1 = gf_mesh_get(mesh, 'outer faces with direction', [0
0 1], 0.1, cvid(indAnode));
fb2 = gf_mesh_get(mesh, 'outer faces with direction', [0
0 -1], 0.1, cvid(indCathode));
% the boundary condition is injecting 1 A/m^2 current
density at the north pole of the sphere, with the south
pole as ground
% here indAnode and indCathode is the index of the
element corresponding to the north and south pole
anode_area = 3;
cathode_area = 4;
gf_mesh_set(mesh, 'region', anode_area, fb1);
gf_mesh_set(mesh, 'region', cathode_area, fb2);
Did you check that your boundaries were ok (may be this is your
first graph ) ?
gf_model_set(md, 'add Dirichlet condition with
multipliers', mim, 'u', mfu, cathode_area);
gf_model_set(md, 'add initialized data','Jn',
ones(gf_mesh_get(mesh, 'nbpts'),1));
% here is the line that I suspect
mostly. I have to pass 'Jn' a vector, otherwise it won't
solve.
gf_model_set(md, 'add source term brick', mim, 'u',
[num2str(sigma(1)) '*(-Grad_u.Normal)'], anode_area, 'Jn');
The source term isthe right hand side of -J.n = f, so in your case
just 1. So you do not need any "gf_model_set(md, 'add initialized
data','Jn', ones(gf_mesh_get(mesh, 'nbpts'),1));" and you can simply
add
gf_model_set(md, 'add source term brick', mim, 'u', '1',
anode_area);
% Neumann BC of electric potential
% solve
gf_model_get(md, 'solve');
% extracted solution
u = gf_model_get(md, 'variable', 'u');
E = gf_model_get(md, 'interpolation', '-Grad_u', mfu); %
electric field
% display
figure; gf_plot(mfu, u, 'mesh','on','cvlst', get(mesh,
'outer faces')); colormap(jet); colorbar
figure; gf_plot(mfE, E, 'mesh','on','norm','on','cvlst',
get(mesh, 'outer faces')); colormap(jet); colorbar
%==============================================================================
--
Yves Renard (address@hidden) tel : (33) 04.72.43.87.08
Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------
|