CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   writing own laplacian scheme with boundary conditions (https://www.cfd-online.com/Forums/openfoam-programming-development/157308-writing-own-laplacian-scheme-boundary-conditions.html)

riddim July 25, 2015 19:09

writing own laplacian scheme with boundary conditions
 
Hey there,

I think I found out that I "simply" need to write my own laplacian scheme with boundary conditions ....

http://www.cfd-online.com/Forums/ope...te-scheme.html <<< this threat was very useful to me but I have to admit I'm still a little bit stuck at the start ...
...maybe you could give me some hints how to implement this ...

My aim is to simply solve

Laplacian(A) == -µ_r*µ_0*J

µ_0 is const. and µ_r can be different from cell to cell

A * n == 0
and (A x curl(A)) * n == 0

where n is the surface normal vector.

do I have to implement these boundary conditions inside the laplacian scheme or is it possible to use the standard laplacian scheme and implement these conditions inside my solver ...?
How do i access the surface normal vector?

Where do you think I should start ...?

chriss85 July 27, 2015 05:44

If you have a variable magnetic permeability this equation is incorrect because you are neglecting a term. Please see the derivation of the formula as I don't recall it from my head. Alternatively, you can use another formulation for calculating a magnetic field. See for example this paper for the T-T0-Phi formulation: http://ieeexplore.ieee.org/xpl/login...mber%3D4787389

riddim July 28, 2015 14:14

thank you very much for your quick reply and sorry for my slow one ...

I'm very aware i'm neglecting a term here...

I already found a couple of formulations e.g. in
http://ieeexplore.ieee.org/xpl/artic...eld=Search_All
or
http://ieeexplore.ieee.org/xpl/artic...eld=Search_All

my Problem is that I didn't see how I can implement two/three equations to be solved simultaniously ...

e.g. in the first paper they suggest to solve

Laplacian(A) == -µ_r*µ_0*J

with

B_t1 == dA_n/dt2 - dA_t2/dn

(t1 = tangential component of cell 1, t2 = tan comp. of cell 2 and n = normal vector of the border)


I'm now a little bit confused where it is possible to implement this relation ...
My idea is to "correct" the tangential/normal component of the B_Field of the surrounding cells temporarily before solving the Laplacian term. Is such a temporary correction possible in my solver or do i have to write a new laplacian scheme for my purposes ...? (or do i have to implement a new surface-normal-correction-scheme ...?)

riddim July 28, 2015 15:22

ok and a little bit tricky might be that i need to correct the vectors of the A field and therefore would neet a curl⁻1 operator since B = curl(A) ...

i suppose it is possible to access the neighbour cell vectors in the laplacian scheme ... (it must ... otherwise you couldn't calculate the laplacian or the divergence or whatsoever ...) is there a list of commands that are possible somewhere ...? i didn't find it yet ..

danny123 July 29, 2015 10:29

Hello,

I am looking into that laplacian too, see

http://www.cfd-online.com/Forums/ope...auf-p_rgh.html

There is already an explicit laplacian of the form Laplacian(A) implemeted, it seems to me. Look into gaussLaplacianScheme.C line 134.


The implicit form is only of the type:

Code:

fvm::laplacian(gamma, vf)
It calls fvmLaplacianUncorrected, where the matrix seems to be built. At this point the source seems to be empty (or at least should be).

Back in fvmLaplacian,
Code:

mesh.V()*fvc::div(tfaceFluxCorrection())().internalField()
to the source term.

To be honest, I did not understand the entire code, but
Code:

deltaCoeffs
seems to be the gradient matrix coefficient.

In order to build
Code:

fvm::laplacian(vf)
I think you can copy the code of
Code:

fvm::laplacian(gamma, vf)
and remove everything related to gamma. Otherwise, an easy way is just to define a surfaceScalarfield gamma, all entries being 1 and use the laplacian as it is implemented.

If you could look into my problem, while you are looking into yours, I would appreciate your input. I am kind off stuck.

Regards,

Daniel

chriss85 July 29, 2015 10:39

You can use fvc::curl(A), but that is an explicit formulation. As for solving equations simultaneously, you have the choice of using a segregated approach (that is a loop, in which both equations are solved separately until convergence) or a coupled approach, in which both equations are discretized in such a way that only one linear system remains to be solved. However, I'm not sure how the coupled approach would be implemented in detail, openfoam-extend might be needed for that.

riddim July 30, 2015 12:40

@chriss85 hmmm - okay - the coupled approach would be pretty cool of course but I would go with the sequential approach at first ... a correct working slower solver is better than a solver not working at all :)

the segregated approach is like in the "compressible/rhoSimpleFoam" where 3 different equations are defined and one after another is done .solve() ..?

e.g. the UEqn.H is pretty simple ... (OF extend 3.0)
Code:

    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevRhoReff(U)
    );

    UEqn().relax();

    eqnResidual = solve
    (
        UEqn() == -fvc::grad(p)
    ).initialResidual();

    maxResidual = max(eqnResidual, maxResidual);

is this .relax() important?
the residuals are optional and not necessary for any solver - are they ...?


Is it correct that the functions i can use in my solver are pretty much the ones defined in src/finiteVolume/finiteVolume(/fvc)
where I have Curl, Div, Laplacian, Grad, SnGrad, ...

Within the snGradSchemes I found the expression
Code:

    const unallocLabelList& owner = mesh.owner();
    const unallocLabelList& neighbour = mesh.neighbour();

    forAll(owner, faceI)
    {
        ssf[faceI] =
            deltaCoeffs[faceI]*(vf[neighbour[faceI]] - vf[owner[faceI]]);
    }

this is how I could address neigbour-field-values within my solver/Laplacian-Scheme ...?

faceI is an iteration parameter used by the solvers ...?

riddim July 30, 2015 13:21

ok sorry for not just testing it ...

the Laplacian Scheme seems to iterate with "PatchI" and mesh.owner() + mesh.neighbour() are valid expressions ... however I don't see where I can access the actual field values ... I think I have to read a little bit more and think about which functions i really do need before being able to ask qualified questions ...

danny123 July 31, 2015 05:30

No, I think you are somewhat off-track...

The iteration parameter is vf. So, the code lines you found correspond to an explicit gradient, since the multiplication is done by those code lines.

To the other questions: relax means that your equation to be solved is relaxed. You can set relax parameters in fvSolutions. If you do not, there is no relaxation.

initialResidual() is one of the solverperformance parameters. This statement iterates the equation to find a suitable U to fit UEqn = - fvc::grad(p). The RHS is explicit, so added to the source term of the linear equation system. While doing this the initial residual (residual of your guess value) is written into a scalar eqnResidual (To be defined somewhere else). In a system of several fields to be iterated separatly (segregated approach), you need a check that at start of each iteration, where your residual is at that point. Only if all initial residuals are below a tolerance value to be defined by you, the whole system of equations has converged.

Regards,

Daniel

riddim July 31, 2015 05:37

...ok yes vf is the iteration parameter that takes care of every cell ... i thought faceI would be an "inner" iteration process because every face has to be dealt with in every large iteration step ...

thx :) sorry now i will be busy for 3 days so i'll come back here after that

have a nice weekend


All times are GMT -4. The time now is 20:03.