Unable to get accurate solution of Laplace Equation in the presence of curved Dirichlet boundary.

[raghwendra-78]

Hi,

I am trying to calculate electric field enhancement factor at the apex of a hemi-ellipsoidal protrusion on a cathode plate in a parallel-plate diode where two conducting plate are kept at a distance 'd'. Cathode is assumed to be at zero potential and anode is at some fixed potential 'V'. Side walls are homogeneous 'Neumann' boundary.

For this purpose I am using libmesh-1.6.1 (compiled with petsc-3.15.1). Geometry and Mesh of the diode was created using Salome Platform and mesh file was exported in unv format . I adapted introduction_ex4.C file given example/introduction

Program file: introduction_ex4 -C.txt

Mesh file for just two plates:
cube-unv.txt

Mesh file for diode with hemi ellipsoid: At present I am unable to upload the file, may be due to large size of the file. Later I will post a link to the file in the thread.

In case of simple diode with two plates; cathode at zero and anode at 1000 volt, d=200 mm ( two Dirichlet boundaries and four homogeneous Neumann boundaries (side walls of the diode)), I get desired result, electric field is 5 V/mm.

When I solve the problem with hemi-ellipsoidal protrusion of height 40 mm on the cathode plate, it is expected that electric field at a point situated far-away from the protrusion should be same as that in the absence of protrusion. However I get 5.158 V/mm instead of 5.0 V/mm. I compared the field enhancement factor at the apex of the protrusion calculated using libMesh and AC/DC module of Comsol Multiphysics software as well. libMesh example program underpredicts the enhancement factor.

I shall highly appreciate any suggestion to improve the simulation result. (output log file and slice of solution of protrusion problem with mesh are attached below)

simulation-hemi-ellipsoid
log.txt

Thanks.

Raghwendra Kumar

[jwpeterson]

If I understand correctly, your quantities of interest are gradients of the solution at fixed points, i.e.

g1=meshfn->gradient(Pt1,t);
g2=meshfn->gradient(Pt2,t);
g3=meshfn->gradient(Pt3,t);
g4=meshfn->gradient(Pt4,t);

which I would expect to converge at first-order "O(h)" if you are using linear Lagrange elements (and assuming the true solution is not linear, in which case you should likely recover the true solution using this FEM discretization). However, you only mention a single value, 5.158, which does not match the value of 5 that you expect. So my questions are:

1. What happens when you refine the mesh? (Do this in both libmesh and Comsol if possible.)
2. How are gradients computed in Comsol? In libmesh, the MeshFunction approach finds the element in which the point of interest is located, and differentiates the FE solution to get the gradient, but it may be possible to get a better approximation for the gradient, especially when solving Laplace/Poisson equation.

[raghwendra-78]

I am using Comsol graphical user interface of AC/DC module for building the model and meshing. For meshing, there are two options, (i)user controlled meshing and (ii) Physics controlled meshing. I used first option, it has further options for mesh size (a) predefined and (b) custom. I used option (b) where I could give Max and Min. Mesh parameters mentioned in the table in my post. custom option has three more parameters: maximum element growth rate, resolution of curvature and resolution of narrow regions. Values for these parameters were taken to be 1.3, 0.2,1 respectively.

Salome has option of using Netgen package. It has following options apart from Max and Min. size:
second order meshing, fineness, limit size by surface curvature, optimize. I have mentioned fineness parameter in the table. Other two options, limit size by surface curvature and optimize were enabled.

I chose first order meshing, as second order option in salome with second order finite element option in libMesh with Lagrange element gives error.

So, I have kept Max. and Min size parameters same in both simulations. I am not aware about other details/algorithm used in Comsol.

[jwpeterson]

Based on your comments and my limited understanding of Comsol, it sounds like there is not a straightforward way to use the same mesh in both Comsol and libmesh. What we would need in Comsol is a "user-controlled meshing" where you simply specify the mesh to be used rather than parameters. Either that, or a way of exporting the mesh from Comsol so that you could then use it directly in libmesh. If we aren't using the same meshes, that makes it difficult to draw any conclusions from your tabulated values... I had assumed those results were using exactly the same mesh.

I chose first order meshing, as second order option in salome with second order finite element option in libMesh with Lagrange element gives error.

Would you mind opening a Github Issue with the details of this error? If it's something relatively simple to fix, I agree that using a higher-order element to represent the curved boundary is likely to give better results.

[raghwendra-78]

Thank you.

“Based on your comments and my limited understanding of Comsol, it sounds like there is not a straightforward way to use the same mesh in both Comsol and libmesh. What we would need in Comsol is a "user-controlled meshing" where you simply specify the mesh to be used rather than parameters. Either that, or a way of exporting the mesh from Comsol so that you could then use it directly in libmesh. If we aren't using the same meshes, that makes it difficult to draw any conclusions from your tabulated values... I had assumed those results were using exactly the same mesh.”

I will try to export mesh from Comsol and repeat the calculation with libmesh.

“Would you mind opening a Github Issue with the details of this error? If it's something relatively simple to fix, I agree that using a higher-order element to represent the curved boundary is likely to give better results.”

I will open a Github issue regarding second order mesh and element.

[raghwendra-78]

Updates: As suggested by Dr. Peterson, I tried using same mesh in Comsol and libMesh but could not suceed. Comsol exports its mesh in .mphbin and .mphtxt format which I can not read in libMesh. After this approach failed, I tried to import mesh exported from Salome to Comsol. I could import the mesh but could not mark faces to impose boundary conditions.

When I created second order mesh from Salome, libMesh gives following error:
ERROR: UNV-element type 22 not supported.

Do I need to open a Github issue for this error?

[jwpeterson]

When I created second order mesh from Salome, libMesh gives following error: ERROR: UNV-element type 22 not supported.

Do I need to open a Github issue for this error?

Yes, please do, and attach or provide a link to the mesh which cannot be read in. If we have a libmesh analogue for UNV element type 22, we may be able to add support for it without too much trouble.

[jwpeterson]

OK, so if I understand your table correctly, for the value of ~5.16 at Point 2, we have good agreement between libmesh and comsol, but at Point 1, which is close to the apex, we don't get good agreement. In addition, the value at Point 1 seems to vary quite a bit with the mesh in the libmesh case, but in the Comsol case the value seems almost independent of the mesh (~28). The Comsol behavior suggests the solution is already "converged" since the quantity of interest doesn't change much under the different mesh discretizations.

Is Point 1 actually a vertex of the mesh in each case? One issue I see is that, because the gradient is actually discontinuous at the mesh vertices, the MeshFunction-based approach to evaluating it can pick any element which touches that vertex to evaluate the gradient in, and you would then get different values depending on which element is picked. Comsol may have a smarter approach for recovering the gradient at vertices, which is why the value doesn't seem to depend on the mesh density.

[raghwendra-78]

Kindly suggest how to check whether point1 is vertex of mesh or not?. Apart from Linear Lagrange, what other type of elements, I can use? Please suggest, better ways to calculate gradient. I will do the post processing if required.

[jwpeterson]

Kindly suggest how to check whether point1 is vertex of mesh or not?

For debugging purposes, you could just loop over all the mesh nodes and see if any of them are "close" to Point1. It seems unlikely to me that Point1 would be a vertex of the mesh if it was chosen more or less arbitrarily, but still good to confirm one way or the other.

Please suggest, better ways to calculate gradient.

I would investigate so-called "gradient recovery" techniques for this. They are frequently used for error estimation, and @variscarey can probably suggest a good approach. In the case where Point 1 is a vertex of the mesh, simply averaging the element gradients from all the neighboring elements would likely give you a better result than what you are doing currently.