ff3d-users
[Top][All Lists]

## [ff3d-users] Problem with elliptic hollow cylinder section

 From: Tobias Nähring Subject: [ff3d-users] Problem with elliptic hollow cylinder section Date: Thu, 22 Apr 2004 21:44:13 +0200 User-agent: Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.6) Gecko/20040401 Debian/1.6-4

```Note:
A short summary of the problem with the corresponding pictures and
files is available at: http://www.tn-home.de/Soft/Fem3d/index.html.

I want to solve the Laplacian equation in the domain

G := {(x,y,z)|   1 < (x/rx)^2 + (y/ry)^2 < a^2;   0 < z < L }

which resembles a section of a hollow cylinder with elliptical
cross-section. The constants are: rx=1, ry=0.5, a=2, L=5.

The boundary condition is

u(x,y,z) = (|x,y|-1)*(a-|x,y|)*(L-z)/L

with

|x,y|:=sqrt( (x/rx)^2 + (y/ry)^2 ).

That means effectively that u(x,y,z)=0 for z=L, and for (x,y,z) on the
cylinder mantles. Only at z=0 one has some nice non-zero boundary
condition.

Because of the symmetry of the boundary condition the
computational domain can be confined to the quarter

Q:={(x,y,z) in G| x>0, y>0}

of G with natural boundary conditions for u on the boundaries at x=0
and y=0.

I tried to solve this problem with FreeFEM3D without much luck.

The results of my attempts are shown at
http://www.tn-home.de/Tobias/Soft/Fem3d/index.html

ff3d accepts the input file and the corresponding pov file (included below)
and the solver options are set to the save values

cg(epsilon=1e-8) and cg(maxiter=2000)

but the output is rather erroneous. I have played
already with other solver settings, such as bicg(...) and
krylov(type=bicg) without any luck. The grid size is as large as
possible on my computer with 512MB RAM. (If I enlarge the grid size
then ff3d freezes during "Marching cube ...".

With the current grid size I get the error message:

.
warning : meshing object plane {
(0, 1, 0),2
}
a problem occured ...
Check that:
- the object is not to small compare to the background mesh
- the object intersects the background mesh

I don't know how to handle that.

The ff3d computed solution of the problem are graphically represented
in the figures:

http://www.tn-home.de/Tobias/Soft/Fem3d/

Is there any chance to get a better solution?

One more remark: There would be a way to specify the mantle surface of
a cylinder section: the `open' directive of the povray object
cylinder (as described in the povray manual).

The input file for ff3d:
//////////////////////////////////////////////////////////////////////
// Model of a cable with ellipsoidal cross-section.

vector Grid=(70,60,30);
double rx = 1.0; // radius of the inner conductor in x-direction
double ry = 0.5; // radius of the inner conductor in y-direction
double a=2.0; // The cross-section boundary of the inner conductor
// scaled by a gives the cross-section boundary of the
// shielding.
double CableW = a*rx;
double CableH = a*ry;

double CableL = 5.0; // length of the cable
double Border = 0.1; // Border added around the modelled piece of the
// cable to build a fictitious domain

scene myScene = pov("cable-ellipse-annulus-num.pov");

```
mesh myMesh = structured( Grid, (-Border,-Border,-2.0*Border), (a*rx+Border,a*ry+Border,CableL+2.0*Border)); // Comment: in z-direcion we provided a bit more space for the fictitious domain (just guessing).
```
domain myDomain = domain(myScene, inside(<1,0,0>));
mesh myCylMesh = surface(myDomain,myMesh);

function norm2d=sqrt((x/rx)^2+(y/ry)^2);

// Realizing the boundary conditions with method-type penalty doesnot
// work (no convergence).
solve(u) in myDomain by myMesh
cg(epsilon=1e-8),
cg(maxiter=2000)
{
u=(norm2d-1)*(a-norm2d)*(CableL-z)/CableL on myCylMesh;
}

//////////////////////////////////////////////////////////////////////
// Results output for postprocessing via gnuplot:

double myx; double myy; double myz;
double mydx = CableW/x(Grid);
double mydy = CableH/y(Grid);
double mydz = CableL/z(Grid);
double myxmax=CableW+0.5*mydx;
double myymax=CableH+0.5*mydy;
double myzmax=CableL+0.5*mydz;

for( myz=0.0; myz <= myzmax; myz = myz + mydz ) {
for( myx=0.0; myx <= myxmax; myx = myx + mydx ) {
for( myy=0.0; myy <= myymax; myy = myy + mydy ) {
```
cout << "@" << myx << ", " << myy << ", " << myz << ", " << u(myx,myy,myz) << "\n";
```   }
cout << "@\n";
}
cout << "@\n";
}
//////////////////////////////////////////////////////////////////////

The pov file (named cable-ellipse-annulus-num.pov):
//////////////////////////////////////////////////////////////////////
// Model of a cable with ellipsoidal cross-section.
// rx=1.0 radius of the inner conductor in x-direction
// ry=0.5 radius of the inner conductor in y-direction
// a=2.0 the inner conductor scaled by a gives the shielding
// CableL=5 length of the cable.
// Since we impose symmetric boundary conditions we need only
// to model a quarter of the cable. (0 < x < a*rx, 0 < y < a*ry)
intersection {
difference {
cylinder {
<0,0,0>,
<0,0,1>,
2.0 // a==2.0
}
cylinder{
<0,0,-0.01>,
<0,0,1.01>,
1.0
}
plane{
<0,1,0>, 0
}
plane{
<1,0,0>, 0
}
}
box{
<0,0,0>,
<2,2,1>
}
scale <1.0,0.5,5.0> // <rx,ry,CableL>
pigment { color rgb <1,0,0> }
}

```