DC resistivity modeling

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi Marios,

I think you want to say : V=I/(2*pi*conductivity*r)

The analytical expression that you use to find the electrical potential works for a point source.
In your 2D model a point source is equivalent to a line source in 3D. (It's not the same analytical expression)
You must use a 3D model to use the analytical expression.

To use a 2D model you must use a fourier transformation (see Loke (Res2dInv) >>> 2.5 D).

Best regards,

Yannick Fargier

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Yes, you are right.

I actually have a c++ code for a 2.5d modeling, and this is how I compare it.

I want to use comsol to calculate the potential for each differential wave number.

If I use the exactly same formation , in 3D modeling, assuming an infinite z direction, would I have the same results?

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

I'm not sure, but I hope.
I Think that assuming a finite z direction (5 times the maximum electrode spacing), is a good approximation.

I work a lot on ERT method with comsol (especially on sensitivity and 3D effects), so if you have questions do not hesitate to contact me.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Thanks...
For starters, can I use the model i attached in the first model?

Or, can you send me a model file from your work?

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

I think that you can't use your first model.

I took few minutes to create this model.
This model represent a wenner quadripole.
Postprocessing shows a fast representation of the "sensitivity" of a wenner quadripole.

I use to create a new application for each electrodes (here 4 electrodes = 4 applications)

When you have a lot of electrodes it implies a lot of measurements (N elec = N*N-1*N-2*N-3/8 measurments)
So, I think it's better to make combinations of pole results.

If you compute the apparent resistivity with an analytical function you will found a good approximation of the resistivity of the model (here I take 1 S/m)

I see that we work on an identical method (DC-electrical method). I'm making a Phd for the LCPC Nantes (In france). I work on the use of ERT method for the detection of infiltration in dikes. I would know what is the goal or application of your research. And I would know, without indiscretion, where are you making your studies.

I hope the model will help you.

Best regards,

Yannick

Attachments:

model_sensitivity.mph

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Thanks for the model. I will look at it and get back to you.

I have a Phd from Aristotle University of Thessaloniki.

Now I have just started a post-doc at Colorado School of Mines.

My final goal is to inject AC current so I can model Spectral IP data.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

By the way you don't use sink.
So if you have let's say 10 electrodes, (no mater what the array), you inject current to first electrode and calculate the potential to the other 9. Then inject current to the second and calculate the potential to all other.

This way you have to solve the problem only as many times as electrodes (using pole-pole combination)

Finally you have a square matrix, which each line representing different point source, and potential to all other electrodes.

So if I want to calculate for a specific array (let's say 4 electrode array) all I have to find is this

AM=current at A --> calc V at M.
AN=current at A --> calc V at N.
BM=current at B --> calc V at M.
BN=current at B --< calc V at N.

Finally array_data=k*(AM-BM-AN+BN);
where k is the geometrical factor.

The only think I do is to move A,B,M,N electrode positions accordingly to the array.

Right?

And a second question, because this is a 3D array, you don't have to solve for different wave numbers.
You can use the same model, for also 3D arrays?

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

You must have work with M. Tsourlos and M. Tsokas. I read some of their papers.

I agree with all that you say in your last post. I use to employ Comsol with Matlab to automate this.

Y.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Yes, you are right. My code is based on Tsourlos one.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

I have created the same geometry in 3D, but I can't upload it because it is 6mb. Still the solution is not the same as comsol.
Do you want me to send it via email?

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi

An annex comment, to upload, start to zip the file, and in the worst cas do a "reset model" first, but then give some clues how to mesh it in the model properties comments

Have fun Comsoling

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

sorry

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

The 3d Model

(with the correct license)

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

I don't understand why you don't find a good result ? I try to find the resistivity of your model and I find it.

Do you noticed, that your electrodes are not on the surface.
If you make this choice of electrode layout you must use this equation to find rho_a.

rho_a = dV / dV_0.

Or, you must take a more complex geometrical factor

dV_0 represents the electrical potential between two electrodes in a model of resistivity = 1 ohm.m
dV represents the electrical potential between two electrodes in a model of resistivity = 10 ohm.m (your model)

for more information see Marescot (2006) or kunetz (1966)

Y.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

I shifted the electrodes by 1 meter....

Now, inject current in first electrode, and calculate potential to the other.

Comsol finds V=0.662761
and by analytical solution is
V=1/(2*pi*conductivity*r)=0.7958.

Electrode spacing is 2 meters.

There is difference.

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

I agree with you, there is difference with the analytical solution.

I think that you must enlarge your model (1/2 space theory). ReMesh electrodes.
In this case, I find an error less than 1%. It still hight.

There is a problem of stability when you use just one electrode in ERT method. You must also create a sink of current.

One sink electrode is generally dedicated. We spell this electrode "pivot" in french.

If you use this two current electrodes you will find no error.

Y.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Basically, after some testing, you are right. I need to remesh a little bit, and the solution becomes stable.

1) Anyway, do you know a way to calculate the jacobian using comsol? (the derivatives of the measurements with respect the conductivities)?

2) Also, as an effort to reduce the size of the mesh, I want to use not homogeneous dirichlet conditions (not ground), but specify a value on the edges.
Usually, you do it using the equation
V=1/(2*pi*conductivity)*B(kr),
where B is a modified Bessel function of zero order (and I think first kind), k is the wave number and r is the distance from the source.

I am using a c algorithm to calculate this?
Is any other way to calculate it with (matlab i guess?).

Thanks.

P.S. I am sending again the model I am using, just in case someone needs it!!!

Attachments:

full_3d_set_ver2.zip

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi,

2) I believe that COMSOL has two Bessel function (first and second kind) in his functions library.
Do you have an article citation for this particular point.

I don't search in this direction but COMSOL propose infinite elements to reduce the size of your model.
If you find something please share your discoveries.
To mesh your model, you must have a fine mesh near electrodes and coarse Elsewhere (mesh >> free mesh parameter >> point, and free mesh parameter >> global >> coarser) It will reduce hightly the number of mesh.

1) You can calculate Fréchet derivates (Jacobian matrix) with the adjoint state method (Park & Van, 1991), Or with the sensitivity equation method (McGillivray and Oldenburg, 1990). Comsol can make a sensitivity study. (Look in solver>>sensitivity or optimization. If you find also something interesting please share your discoveries.

Sensitivty = d V_i / d rho_j = postint ( J_emdc * J_emdc ). For a pole pole configuration (Park & Van, 1991). I use this exemple in the first model that I send to you.

P.S.: I can send you the articles if you want

for my part, I'm trying to resolve the inverse problem in 3D ERT. I want to use a quasi newton approach.

(J' Wd J + µ C) dm = J' Wd( dcalc- dpred) + ... J'=transpose of J, and J = sensitivity matrix.

I want to regularize my problem (matrix C) but I don't know how I can construct this matrix because inversion cell = tetrahedron of different size and shape.

Do you have advices for my problem ?

Thanks in advance.

Yannick

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Quasi newton is an update technique. You need to find it for the first time using e.g. the sensitivity technique, and after that in each iteration you update it like this

If num_param is the number of cells (unknows)
num_mes is the number of measurements.

dx is the update for the model in each iteration
store the d_calc for the first iteration as d_oldmes

after you calculate analytically for the first iteration (and probably and the second for accuracy), then you must update J like this

dd=0;
for i=1:num_param
dd=dd+dx(i)*dx(i);
end

dm=log10(d_calc) - log10(d_oldmes) - array_jacobian*dx;

for i=1:num_mes
for j=1:num_param
array_jacobian(i,j)=array_jacobian(i,j) + dm(i)*dx(j)/dd;
end
end

I have tested it, and the results are poor, I always prefer to calculate it in each iteration.

Regarding the regularization, (if i get it right), you should not search solution for each tetrahedron.
This is what I do.
Create several cubes side by side, which each cube has a specific conductivity. This is your unknown.
Make as many as you cover the whole 3d space. (be careful, do not create too many, usually depends on electrode separation).
Now let comsol create tetrahedron elements. Of course each element inside cube should have the same conductivity as the cube.

Now, because the cubes you created have known dimensions (create them with the same dimensions) you can create the C matrix.

Consider creating three C matrices, for each dimension (or unless you want regularization in specific dimensions)
You need Cx, Cy, and Cz.
Generally C is a matrix that has 1 in it's diagonal and -1 in another diagonal. All others are zeros.

Abderrahim Jardani

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi Marios,

algorithm of the QN what you wrote it's not correct,
for more details about the QN see Loke and Barker Paper.
abdel

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

From www.geoelectrical.com/gaussnew.pdf

dd=0;
for i=1:num_param
dd=dd+dx(i)*dx(i);
end

This is pT*p
you could write it in matlad as
dd=dx'*dx;

dm=log10(d_calc) - log10(d_oldmes) - array_jacobian*dx;

This is (Dy - B*p)

for i=1:num_mes
for j=1:num_param
array_jacobian(i,j)=array_jacobian(i,j) + dm(i)*dx(j)/dd;
end
end
This is u=(Dy - B*p)/pT*p

actually you can write it in matlab as

array_jacobian=array_jacobian + dm*dx'./dd

I don't see the error

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Another Issue.
If I consider electrodes placed in a container (box), so boundaries conditions are electric insulation everywhere, and I don't use a sink, comsol can't find a solution.
In this case, I have to use a sink. But this way I need to solve for every combination of electrodes.
Is any way to avoid this? (Without using sink, so I have to solve the problem only as times as the electrodes).

Thanks

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi marios,

I think there is any way to avoid this problem without using sink. If you don't create a sink or a special boundary condition (ground for exemple) there is no continuity, so you can't inject current.

As I said, you can create a "pivot" electrode. You define one additional electrode and you affect always a sink condition to this electrode.

Then, you solve direct problems for every electrode (additional electrode excepted ) . ex : 10 electrodes and one sink electrode = 10 direct problems.

Finally, you must combine solutions to calculate the electrical potentials. V_M = V_DP1 - V_DP2.
The electrical potential at M = the potential at M, result of the first direct problem (DP1) - the potential at M, result of the second direct problem (DP2).

V_MN = (V_DP1_M - V_DP2_M) - (V_DP1_N - V_DP2_N)

The influence of the "pivot" electrode is deleted during this operation.

Y.

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi marios,

I think there is any way to avoid this problem without using sink. If you don't create a sink or a special boundary condition (ground for exemple) there is no continuity, so you can't inject current.

As I said, you can create a "pivot" electrode. You define one additional electrode and you affect always a sink condition to this electrode.

Then, you solve direct problems for every electrode (additional electrode excepted ) . ex : 10 electrodes and one sink electrode = 10 direct problems.

Finally, you must combine solutions to calculate the electrical potentials. V_M = V_DP1 - V_DP2.
The electrical potential at M = the potential at M, result of the first direct problem (DP1) - the potential at M, result of the second direct problem (DP2).

V_MN = (V_DP1_M - V_DP2_M) - (V_DP1_N - V_DP2_N)

The influence of the "pivot" electrode is deleted during this operation.

Y.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

1) You can calculate Fréchet derivates (Jacobian matrix) with the adjoint state method (Park & Van, 1991), Or with the sensitivity equation method (McGillivray and Oldenburg, 1990). Comsol can make a sensitivity study. (Look in solver>>sensitivity or optimization. If you find also something interesting please share your discoveries.

Sensitivty = d V_i / d rho_j = postint ( J_emdc * J_emdc ). For a pole pole configuration (Park & Van, 1991). I use this exemple in the first model that I send to you.

Hi Yannic,

Did you find an another way to calculate the sensitivity? IIn your attached example, you calculate (J_source dot J_receiver), but the integration is skipped. Did you had any luck?

thanks in advance.

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hello Marios,

As I said in my previous reply, I use the "postint" function. (with comsol V3.5)
There is two way to use this function :

1) declare a volume with a string :

Si,j = postint( fem , ' (x>10)*(z>10)* (Jx_emdc1*Jx_emdc2+Jy_emdc1*Jy_emdc2+Jz_emdc1+Jz_emdc2)' , 'dl' , [1] );

2) Or create subdomains, and integrate current density over this subdomain :

Si,j = postint( fem , ' (Jx_emdc1*Jx_emdc2+Jy_emdc1*Jy_emdc2+Jz_emdc1+Jz_emdc2)' , 'dl' , [...] );

3) Use the meshintegrate function :

... good luck

Create blocks (subdomains) is very time consuming but gives an accurate result, I choose this way to create my sensitivity matrix.

If you find THE BEST way to compute the sensitivity matrix with comsol, I would like to know it.

kind regards,

Yannick Fargier

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Thanks again for your reply.

I am creating subdomains, where I calculate the integral (comsol 4.2, results-->derived values --> Volume integral) of let's say a wenner array like this

(ec.Jx-ec4.Jx)*(ec2.Jx-ec3.Jx) + (ec.Jy-ec4.Jy)*(ec2.Jx-ec3.Jy) + (ec.Jz-ec4.Jz)*(ec2.Jz-ec3.Jz), in each domain.

ec -> A electrode
ec4 -> B electrode
ec2 -> M electrode
ec3 -> N electrode.

The problem is by creating those subdomain, comsol is slow but ok.
If I need to do that.] for 50 electrodes, then I need to create 50 different physics. 50 physics with blocks makes things SLOW.

If for every source location, I store V and the three components of J, would comsol be able to integrate (without storing all of the different physics)?

Yannick Fargier

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi marios,

this is my personnal email :
yannick.fargier@ifsttar.fr

can you send me a mail, and I will send you a part of my code to compute the senstivity matrix (but not this week I will not have the time).

To resume :
create one model for all forward problem
compute each forward problem with the fem solver : 50 electrodes = 50 different forwards problems
create one super FEM object with 50 application, containing 50 small fem object (function femsol).
then compute Si,j = postint( FEM , ' (Jx_emdc1*Jx_emdc2+Jy_emdc1*Jy_emdc2+Jz_emdc1+Jz_emdc2)' , 'dl' , [...] );
then, sum the pole pole sensitivity to have quadripole sensitivity.
then do inversion, or opitimization, or ...

regards,

Yannick Fargier

Yingying Zhang

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

hi,
I am wondering if the comsol can do the 3D resistivity inversion of ERT. If it can, how? Any clue or document or model or ... is welcome. Thank you in advance.

Rodolphe Duval

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi Yingying,

I'm not sure comsol is best suited to do inversion (we speak about 100 to 50 000 inversion cell).
So I don't have a model.

But, comsol can solve the forward problem and matlab (with the LiveLink) can solve the inverse problem.

There exist a comsol model example (industrial tutoriel) that show the forward problem and how to compute fréchet derivates,

best regards,

mi ahmad

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

Hi Yingying,

I am also looking for the same solution. If you have any information, kindly share with me as well.

Many thanks

Regards,

Hafisoh

Jonathan Nyquist

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

I, too, would be interested in learning what progress you make with this problem. Comsol should give it some attention.

Marios Karaoulis

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago

There are three different replies on this topic.

1) Short answer is NO.
2) Long answer is yes with the help of maltab (you can export results there) and a parametric solver in comsol. There are several steps how to do it. Briefly for every electrode, you export the results in each node of the mesh. For every pole-pole configuration, you calculate the dot product of current densities and integrate over the element. That's how you calculate you Jacobian. Comsol can also integrate over an element, but I found it extremely slow. Then you build your normal equations and invert. Find new model, throw it back to comsol, new iteration...
3) Exotic answer: YES. Search for a paper named "Efﬁcient solution of nonlinear, underdetermined inverse problems with a generalized PDE mode

EDIT
Unfortunately, when using mapped mesh, accuracy is terrible (compared with analytical solutions). A simple fem code from "The Finite Element Method in Electromagnetics, 2nd Edition", using the some mesh (export mesh to matlab, and write the code showed in chapter 4), is much more accurate. Typical trianglura mesh gives better results.

Maybe comsol can comment out.

Mayank Dhull

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 10 years ago

I have got a simple problem which makes me struggle.
I have a rectangular thin conductive sheet.
I would like to calculate with COMSOL the electrical resistance between the corners of the conductive sheet and any other point on the sheet?
and i am using joule heating module.
I have applied a voltage and now how can i measure the resistance?

Onyekachukwu Nnodim

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 9 years ago

Hi, i'm new to COMSOL. Working on image reconstruction, ERT. I would like to know how to create a number of cells (side by side, creating a square of cells in both x and y direction), so as to find the sensitivity matrix and conductivity distribution.

Thanks

Yingying Zhang

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 9 years ago

In Comsol model library, there is a model named " aquifer_characterization". Square cells are made. You can refer to this model to build yours.

Onyekachukwu Nnodim

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 9 years ago

Thank you much, Yingying Zhang

In Comsol model library, there is a model named " aquifer_characterization". Square cells are made. You can refer to this model to build yours.

Onyekachukwu Nnodim

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 9 years ago

Hi, how do i calculate my sensitivity? I have my derived values of current density norm for 8 electrodes. My formula is Vo = Vin(Rfb/Rf) where Rfb = 15k. I know my sensitivity has to be the (S = Vo of water - Vo of mixture) but i have no idea how to implement these as i 'm new to ERT. Thanks.

I'm not sure, but I hope.
I Think that assuming a finite z direction (5 times the maximum electrode spacing), is a good approximation.

I work a lot on ERT method with comsol (especially on sensitivity and 3D effects), so if you have questions do not hesitate to contact me.

Ayo Onilude

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 9 years ago

Good Day,

Working on electrical resistance tomography.
a) In the attached model, I have defined 3 layers for cylinder geometry (cyl1). This automatically gives concentric cylinders which are seen to be divided into 90 degrees (giving 4 quadrants in plane xy view). How can i vary the angles such that i decide number of section generated?

b) Secondly How can I perform a material sweep such that when a single domain is of a particular material 'a', then the rest of the domains are of a different material 'b', such that every domain one at a time will will be of material 'a' while others are 'b'. This will result to a complete solution where every domain has assumed the material 'a' while others were 'b'.

Thank you, for your much anticipated help in advance.
Kind regards
Ayo

Attachments:

cyl1.mph

Fitz Gerald

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 9 years ago

Dear Marios. I'm on ERT also with focus on sensitivity in flow pipelines.

Any ideas on designing centroid of mass and reading coordinates w.r.t to centroid variables.

Cheers

I shifted the electrodes by 1 meter....

Now, inject current in first electrode, and calculate potential to the other.

Comsol finds V=0.662761
and by analytical solution is
V=1/(2*pi*conductivity*r)=0.7958.

Electrode spacing is 2 meters.

There is difference.

Ayo Onilude

Send Private Message Flag post as spam

Please login with a confirmed email address before reporting spam

Posted: 9 years ago

Dear Ayo. I'm on ERT with focus on sensitivity in flow pipelines.

Any ideas on designing centroid of mass and reading coordinates w.r.t to centroid variables?

Cheers

Good Day,

Working on electrical resistance tomography.
a) In the attached model, I have defined 3 layers for cylinder geometry (cyl1). This automatically gives concentric cylinders which are seen to be divided into 90 degrees (giving 4 quadrants in plane xy view). How can i vary the angles such that i decide number of section generated?

b) Secondly How can I perform a material sweep such that when a single domain is of a particular material 'a', then the rest of the domains are of a different material 'b', such that every domain one at a time will will be of material 'a' while others are 'b'. This will result to a complete solution where every domain has assumed the material 'a' while others were 'b'.

Thank you, for your much anticipated help in advance.
Kind regards
Ayo

Hi Fitz,

Not sure I exactly understand your question. But from the little interpretation I could get.

Once you have designed your "centroid of mass", use "select point" to select a location on your "centroid of mass" drawing. Under geometry, "measure" should give you the coordinates you are look for.

Regards,
Ayo

Discussion Forum

DC resistivity modeling

Suggested Content