3
$\begingroup$

So, I am trying to solve a simple problem using the FEM method.

The distribution of Voltage over a region,

The region is,

  << NDSolve`FEM`
bm = ToBoundaryMesh[
   "Coordinates" -> {{0., 0.}, {1., 0.}, {12/10, 0.}, {12/10, 1.}, {1., 
      1.}, {0., 1.}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        5}, {5, 6}, {6, 1}, {5, 2}}]}];
em = ToElementMesh[bm, 
   "RegionMarker" -> {{{0.5, 0.5}, 1}, {{1.2, 0.5}, 2}}];

which look like this,

     Show[{
  em["Wireframe"["MeshElementStyle" -> Directive[EdgeForm[Red], Thin],
     "MeshElementMarkerStyle" -> Blue]],
  bm["Wireframe"["MeshElementStyle" -> Black]]
  }, ImageSize -> 300]

enter image description here

I have two equations,

op1 = Laplacian[v1[x, y], {x, y}];
op2= Laplacian[v2[x, y], {x, y}];

bc1=DirichletCondition[v1[x, y] == 4, x == 12/10 && y==1];
bc2=DirichletCondition[v2[x, y] == 0, x == 0];

omg=376.991;
dd=100;
d=25*10^-6;

 eef={4.18382,4.22462,9.95342,12.8216,13.9473,14.4967,14.8028,14.9896,15.1113,15.195,15.2551,15.2998,15.3341,15.3611,15.3829,15.4007,15.4155,15.428,15.4387,15.4478,15.4556,15.4624,15.4682,15.4733,15.4777};
    ee1 = ListInterpolation[eef, {{0, Length[eef] - 1}}, InterpolationOrder -> 1]

ee1 has dependence on v1[x,y]-v2[x,y], as ee1[Abs[v1[x,y]-v2[x,y]]]

{v1fun,v2fun} = NDSolveValue[{op1==If[ElementMarker == 1, ((I omg ee1[
     Abs[v1[x, y] - v2[x, y]]])/(dd d)) (v1[x, y] - 
     v2[x, y]), 0], op2==If[ElementMarker == 1, (
   (I omg ee1[
     Abs[v1[x, y] - v2[x, y]]])/(dd d)) (v2[x, y] - 
     v1[x, y]), 0],bc1,bc2},{v1[x, y],v2{x,y}}, {x, y} ∈ em]

The output contour plot of (Abs[v1[x,y]-v2[x,y]]) should be the distribution of voltage over the entire mesh region with the highest values at the x=0 and x=1.2.

I am getting different errors in versions 11.1 and 13.1. I am trying to solve complex equations in NdSolveValue and FEM.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ I think the main problem is that it is not possible to find the derivative of the Abs of complex values. The derivative is needed for linearization. Maybe you can replace the Abs with a smooth function for which a derivative can be found. $\endgroup$
    – user21
    Feb 10 at 8:29
  • $\begingroup$ @user21 Thanks! I will follow the suggestion. However, the answer given below is also very close to solution. That could be an another way as well. $\endgroup$
    – a019
    Feb 10 at 8:33
  • $\begingroup$ Absolutely, I just wanted to express that I am not convinced that the reason is the large numbers, it's a mathematical problem, I think. $\endgroup$
    – user21
    Feb 10 at 8:34
  • 1
    $\begingroup$ Btw. I added code such that you will not see the Part message any longer (but you will see the other messages). $\endgroup$
    – user21
    Feb 10 at 8:36

1 Answer 1

Reset to default
4
$\begingroup$

Parameter k = omg/dd/d of this model is about 150796, together with function ee1 it gives about 10^7 on the right side op1 and op2. It is too high for nonlinear FEM, why NDSolve failed on this example. But the problem could be solved with linear FEM. Therefore we use the false transition method as follows

<< NDSolve`FEM`
bm = ToBoundaryMesh[
   "Coordinates" -> {{0., 0.}, {1., 0.}, {12/10, 0.}, {12/10, 
      1.}, {1., 1.}, {0., 1.}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        5}, {5, 6}, {6, 1}, {5, 2}}]}];
em = ToElementMesh[bm, 
   "RegionMarker" -> {{{0.5, 0.5}, 1}, {{1.2, 0.5}, 2}}];

op1 = Laplacian[v1[x, y], {x, y}];
op2 = Laplacian[v2[x, y], {x, y}];

bc1 = DirichletCondition[v1[x, y] == 4, x == 12/10 && y == 1];
bc2 = DirichletCondition[v2[x, y] == 0, x == 0];

omg = 376.991;
dd = 100;
d = 25*10^-6;

eef = {4.18382, 4.22462, 9.95342, 12.8216, 13.9473, 14.4967, 14.8028, 
   14.9896, 15.1113, 15.195, 15.2551, 15.2998, 15.3341, 15.3611, 
   15.3829, 15.4007, 15.4155, 15.428, 15.4387, 15.4478, 15.4556, 
   15.4624, 15.4682, 15.4733, 15.4777};
ee1 = ListInterpolation[eef, {{0, Length[eef] - 1}}, 
  InterpolationOrder -> 1];
k = omg/dd/d

V1[0][x_, y_] := 4; V2[0][x_, y_] := 0;

Do[{V1[i], V2[i]} = 
   NDSolveValue[{op1 == 
      If[ElementMarker == 
        1, ((I k ee1[Abs[V1[i - 1][x, y] - V2[i - 1][x, y]]])) (v1[x, 
           y] - v2[x, y]), 0], 
     op2 == If[
       ElementMarker == 
        1, ((I  k ee1[Abs[V1[i - 1][x, y] - V2[i - 1][x, y]]])) (v2[x,
            y] - v1[x, y]), 0], bc1, bc2}, {v1, 
     v2}, {x, y} \[Element] em];, {i, 5}]

Visualization on every step shows stable solution at i>2

Table[{DensityPlot[Abs[V1[i][x, y]], {x, y} \[Element] em, 
   ColorFunction -> "Rainbow", PlotRange -> All, 
   PlotLegends -> Automatic, PlotLabel -> "v1"], 
  DensityPlot[Abs[V2[i][x, y]], {x, y} \[Element] em, 
   ColorFunction -> "Rainbow", PlotRange -> All, 
   PlotLegends -> Automatic, PlotLabel -> "v2"]}, {i, 5}]

Figure 1

Update 1. With bc2 = DirichletCondition[v2[x, y] == 0, x == 0&& y==0] we have pictures Figure 2

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks!! In the problem, V1=4 is at the right up corner of the mesh and V2=0 is at the left down corner of the mesh. However, in density plots, both V1 and V2 at the right up corner. V2 should be at the opposite side of the V1. $\endgroup$
    – a019
    Feb 10 at 8:28
  • 1
    $\begingroup$ @a019 Do you mean that your bc2=DirichletCondition[v2[x, y] == 0, x == 0] is wrong and it should be bc2=DirichletCondition[v2[x, y] == 0, x == 0&&y==0]? $\endgroup$
    – Alex Trounev
    Feb 10 at 10:13
  • $\begingroup$ Yes, I thought x==0 will be enough, however it should be bc2=DirichletCondition[v2[x, y] == 0, x == 0&&y==0] $\endgroup$
    – a019
    Feb 10 at 10:16
  • 1
    $\begingroup$ @a019 See update to my answer. $\endgroup$
    – Alex Trounev
    Feb 10 at 10:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.