Comment déterminez-vous la résistance effective d'une grille finie de résistances?

14

Avertissement: Je suis géophysicien avec une formation limitée en génie électrique. Je ne sais pas si ce problème est incroyablement facile, incroyablement complexe ou complètement insensé.

Mon objectif: déterminer la résistivité globale d'un échantillon de roche à l'aide de réseaux de résistances.

L'échantillon de roche doit être modélisé en utilisant un réseau de résistances avec certaines résistances ayant une résistance élevée (représentant la roche solide) et d'autres résistances ayant une faible résistance (représentant les voies de fluide dans la roche).

Supposons que j'ai un réseau de résistances sur une grille uniforme comme indiqué ci-dessous. Dans l'exemple illustré, chaque segment de ligne a une résistance associée étiquetée de 1 à 24 sur une grille 3 par 3. Les résistances de chaque segment de ligne sont connues.

La longueur totale de la grille est et la "zone" est A (dans ce cas, il s'agit d'un exemple 2D, donc la zone n'est également qu'une longueur). La résistivité globale de l'échantillon est alors donnée par:LA

ρbulk=ReffAL

entrez la description de l'image ici

Ma question: comment déterminer la résistance effective, du réseau?Reff

J'ai regardé en ligne mais tout ce que je peux trouver, ce sont des discussions sur les réseaux infinis , les courants de source et de puits, etc. Je ne suis pas intéressé par le courant ou la tension.

Ce problème peut-il être résolu tel quel?

Darcy
la source
2
Je le branche sur un simulateur et je laisse le simulateur le résoudre. Vous pouvez construire votre modèle comme un circuit d'épices. Ensuite, pour trouver une résistance, utilisez simplement la loi d'Ohm (V = I * R). Spice vous indiquera le courant afin que vous puissiez calculer R.
mkeith
1
Le tout peut potentiellement être automatisé à l'aide d'épices en ligne de commande, mais pour preuve de concept, vous pouvez entrer votre circuit dans une épice gratuite telle que LTSpice. Appliquez une tension et affichez le courant. LTspice peut également afficher des fonctions simples telles qu'une tension divisée par un courant (résistance).
mkeith
Darcy, il existe un certain nombre d'approches. Je voudrais poser quelques questions avant de faire part de mes réflexions. (1) Il existe un logiciel très simple à écrire. Cherchez-vous ce genre d'approche? (2) Vous pouvez résoudre ce problème en utilisant l'analyse nodale traditionnelle. Cherchez-vous ce genre d'approche? (3) Votre problème se décompose en sommets et arêtes . (Compte tenu de votre expérience de géophysicien, je m'attends à ce que vous connaissiez la signification de ces deux termes.) Comment déterminez-vous, a priori, les valeurs que vous insérez pour les bords?
jonk
@jonk Je serais principalement intéressé par l'option (1) d'écrire moi-même un court morceau de code pour ce faire. Je détermine les résistances de bord sur la base d'une géométrie de pores a priori et d'une résistivité connue d'un minéral ou d'un fluide rocheux.
Darcy
Darcy, il y a aussi des techniques qui s'appuient sur des réseaux triangulaires irréguliers qui me viennent immédiatement à l'esprit lorsque vous écrivez «voies fluides». Avez-vous lu quelque chose sur ce sujet? Je ne sais pas quels sont vos objectifs ultimes, mais vous voudrez peut-être aussi les rechercher. Ceux-ci seraient formidables à utiliser pour les gradients vous aidant à comprendre où les "courants" auraient tendance à se concentrer. Si c'est un problème.
jonk

Réponses:

11

L'idée de base est assez simple. Vous disposez d'une matrice ( ) qui représente des "nœuds" ou sommets dans votre système. Chacun de ces nœuds est associé à une "tension" à valeur scalaire qui peut être modifiée ou mise à jour au fur et à mesure de l'algorithme. Il y aura également deux nœuds dont la tension ne peut pas être modifiée. Nous allons appliquer ici une sorte de "batterie", donc ces deux nœuds représentent les deux extrémités de cette batterie.V

Séparément, deux autres matrices ( et R h ) représentent les bords du système, horizontal et vertical. Ce sont vos valeurs de résistance, je suppose. Je ne sais pas comment vous comptez les remplir. Mais c'est votre problème. Cette technique suppose que vous pouvez également remplir ces matrices.RvRh

Selon le langage informatique que vous utilisez, vous pouvez ou non être en mesure d'utiliser des indices négatifs. Peu importe. C'est juste une question de garder à l'esprit ce à quoi vous êtes confronté.

Supposons que la longueur soit divisée en N sections L et que la "longueur" A soit divisée en N sections A. Ensuite, vous devrez construire une matrice avec ( N L + 1 )( N A + 1 ) sommets pour les valeurs de tension scalaire. (ou plus.) Vous aurez également besoin de ces deux autres matrices avec N A( N L + 1 ) bords verticaux et N L( N A + 1LNLANA(NL+1)(NA+1)NA(NL+1) bords horizontaux entre ces sommets.NL(NA+1)

Maintenant. Initialise tous les sommets avec . Choisissez l'un des sommets à gauche (au milieu, de préférence) et notez-le comme un 00VValeur V qui n'est PAS autorisée à changer. Utilisez la méthode que vous souhaitez pour cela. Choisissez l'un des sommets à droite (au milieu, de préférence) et changez sa valeur en 10V , tout en prenant note à nouveau que sa valeur n'est PAS autorisée à changer. Une technique qui fonctionne ici consiste simplement à le laisser changer normalement, puis à remplacer la valeur à chaque étape. Mais peu importe comment vous y parvenez, tant que vous y parvenez.1V

(Il existe d'autres techniques pour des raisons d'efficacité. Mais cela ne vaut probablement pas la peine de les utiliser ici.)

Passons maintenant à l'algorithme, qui est parfois appelé algorithme en damier ou rouge-noir . En vous déplaçant dans la matrice de tension de votre nœud, traitez chaque nœud où la somme des deux indices, est paire, en effectuant l'affectation simple suivante:i+j

Vi,j=Rhi,j1Rhi,j(Vi1,jRvi,j+Vi+1,jRvi1,j)Rhi,j1Rhi,j(Rvi,j+Rvi1,j)+Rvi1,jRvi,j(Rhi,j+Rhi,j1)+Rvi1,jRvi,j(Vi,j1Rhi,j+Vi,j+1Rhi,j1)Rhi,j1Rhi,j(Rvi,j+Rvi1,j)+Rvi1,jRvi,j(Rhi,j+Rhi,j1)

L'équation ci-dessus n'est rien de plus que le calcul de la tension d'un nœud central ayant quatre résistances qui s'y connectent, où les tensions aux autres extrémités des quatre résistances sont connues. La tension du nœud central est ensuite calculée à partir de l'équation ci-dessus. Étant donné que le diviseur est le même pour chaque terme, vous pouvez simplement calculer la somme des numérateurs, puis diviser une fois par le dénominateur.

Cela mettra à jour tous les nœuds où la somme est paire. Maintenant, vous effectuez la même procédure pour tous les nœuds où la somme i + ji+ji+j est impaire. Une fois ces deux étapes effectuées, vous avez terminé un cycle.

Si nécessaire, réinitialisez les deux nœuds spéciaux (pour 0V et pour 1V comme expliqué précédemment.) Ou, si vous avez protégé ces deux nœuds, il n'est pas nécessaire de les réinitialiser.

Vous êtes prêt pour le prochain cycle. Effectuez ces cycles autant de fois que vous le jugez nécessaire pour que l'état général se stabilise (et ce sera le cas.)

Lorsque vous arrêtez le processus, vous pouvez facilement déterminer la résistance en choisissant de regarder les nœuds entourant votre nœud protégé gauche ou bien de regarder les nœuds entourant votre nœud protégé droit. (Il peut être judicieux d'agrandir votre matrice juste assez [par 1 dans toutes les directions] pour avoir en fait quatre nœuds entourant l'un ou l'autre choix.) La différence de tension entre les nœuds environnants et le nœud spécial, divisée par le la résistance dans les bords entre eux vous indique la sortie / entrée actuelle de votre nœud spécial. Puisqu'il s'agit d'un nœud "batterie", ce courant doit être TOUT du courant. Puisque la tension est1V , par définition, divisant 1 par la somme de ces quatre courants que vous trouvez vous indique la résistance totale.

Je regarde un code que j'ai écrit qui totalise, avec beaucoup de commentaires, seulement 67 lignes. Il n'est donc PAS difficile d'écrire.

1V


Pourquoi faut-il séparer le système en i + j = pair et i + j = impair?

V5,5=f(V4,5,V6,5,V5,4,V5,6)V5,5V5,6=f(V4,6,V6,6,V5,5,V5,7)V5,5V5,7=f(V4,7,V6,7,V5,6,V5,8)car aucune des entrées de la fonction n'est un nœud qui a été modifié au cours de cette étape. Ensuite, vous vous retournez et calculez les suppléants, en évitant le maculage, mais maintenant en mettant à jour les suppléants. Vous devez vraiment le faire de cette façon.

De plus, la formule est-elle identique pour les étapes paires et impaires?

Oui, c'est pareil.

Tout cela peut-il être résolu en une seule étape en utilisant une sorte de système linéaire Ax = b où A est un opérateur linéaire et b fournit les conditions aux limites? En le regardant, cela semble quelque peu analogue aux méthodes de différences finies pour résoudre des équations aux dérivées partielles.

Il y a une connexion. Je pense que cela s'appelle une implémentation «sans matrice».


Voici un exemple. L'ensemble de valeurs de résistance suivant a été placé dans LTSpice pour la simulation:

entrez la description de l'image ici

1V30.225mA30.224552mA

J'ai exécuté le programme VB.NET suivant:

Module GEOGRID

    Const NL As Integer = 2
    Const NA As Integer = 2
    Const INF As Double = 1.0E+32

    Sub Main()

        Static Rh As Double(,) = New Double(NL + 2, NA + 1) {
                    {INF, INF, INF, INF},
                    {INF, 5, 21, INF},
                    {INF, 76, 10, INF},
                    {INF, 32, 22, INF},
                    {INF, INF, INF, INF}}
        Static Rv As Double(,) = New Double(NA + 1, NL + 2) {
                    {INF, INF, INF, INF, INF},
                    {INF, 61, 50, 16, INF},
                    {INF, 56, 45, 18, INF},
                    {INF, INF, INF, INF, INF}}
        Dim V As Double(,) = New Double(NL + 2, NA + 2) {
                    {0, 0, 0, 0, 0},
                    {0, 0, 0, 0, 0},
                    {0, 0, 0, 1, 0},
                    {0, 0, 0, 0, 0},
                    {0, 0, 0, 0, 0}}
        Dim PDE As Func(Of Integer, Integer, Double) = Function(ByVal i As Integer, ByVal j As Integer) (
                    Rh(i, j - 1) * Rh(i, j) * (V(i - 1, j) * Rv(i, j) + V(i + 1, j) * Rv(i - 1, j)) +
                    Rv(i - 1, j) * Rv(i, j) * (V(i, j - 1) * Rh(i, j) + V(i, j + 1) * Rh(i, j - 1))
                  ) / (
                    Rh(i, j - 1) * Rh(i, j) * (Rv(i, j) + Rv(i - 1, j)) +
                    Rv(i - 1, j) * Rv(i, j) * (Rh(i, j) + Rh(i, j - 1))
                  )
        Dim IV As Func(Of Integer, Integer, Double) = Function(ByVal i As Integer, ByVal j As Integer) 0 +
                    (V(i, j) - V(i - 1, j)) / Rv(i - 1, j) + (V(i, j) - V(i + 1, j)) / Rv(i, j) +
                    (V(i, j) - V(i, j - 1)) / Rh(i, j - 1) + (V(i, j) - V(i, j + 1)) / Rh(i, j)
        Dim idx As Integer = NA \ 2 + 1
        Dim jdx1 As Integer = NL + 1
        Dim jdx2 As Integer = 1
        For x As Integer = 1 To 1000
            For k As Integer = 0 To (NA + 1) * (NL + 1) - 1 Step 2
                Dim i As Integer = k \ (NL + 1)
                Dim j As Integer = k - i * (NL + 1) + 1
                i += 1
                If Not (i = idx AndAlso (j = jdx1 OrElse j = jdx2)) Then V(i, j) = PDE(i, j)
            Next
            For k As Integer = 1 To (NA + 1) * (NL + 1) - 1 Step 2
                Dim i As Integer = k \ (NL + 1)
                Dim j As Integer = k - i * (NL + 1) + 1
                i += 1
                If Not (i = idx AndAlso (j = jdx1 OrElse j = jdx2)) Then V(i, j) = PDE(i, j)
            Next
        Next
        Console.WriteLine("R = " & (1.0 / IV(idx, jdx1)).ToString)
        Console.WriteLine("R = " & (-1.0 / IV(idx, jdx2)).ToString)
    End Sub

End Module

R=33.0856844038614Ω

Le programme ci-dessus montre un moyen de configurer les résistances, verticales et horizontales, ainsi que la matrice de tension, de sorte qu'il simplifie certains des tests pour les nœuds et / ou les valeurs de résistance inexistants. Le code est un peu plus propre de cette façon, bien qu'il nécessite quelques éléments de tableau supplémentaires. (J'ai simplement rendu les valeurs de résistance supplémentaires infinies.) Il suffit de comparer la façon dont j'ai configuré les tableaux avec la façon dont le schéma a été présenté également, et je pense que vous serez en mesure de déterminer exactement détails ici.

J'ai également piraté les résistances et les valeurs des nœuds, bien sûr, sans en faire un programme à usage général pour lire un tableau de valeurs. Mais cette généralité est assez facile à ajouter. Et ce code devrait rendre tout ce que j'ai écrit absolument sans ambiguïté.

XX

0V1V

(D'accord. Une dernière note. Ce serait beaucoup mieux ciblé sur F # ou tout compilateur décent ciblant un système informatique massivement parallèle. Chaque calcul en "rouge" ou "noir" peut être effectué en parallèle, complètement indépendamment les uns des autres. F # rend cela trivial. Donc codé en F #, vous pouvez l'exécuter sur tous vos cœurs disponibles sans rien de spécial à faire. Cela fonctionne juste. Juste une note au cas où vous collectez BEAUCOUP de données d'une manière et que vous pourriez vouloir prendre tous les avantages d'un système multicœur.)


NOTE DE FIN:

La dérivation est assez simple à partir de KCL. Placez quatre résistances dans la disposition suivante:

schématique

simuler ce circuit - Schéma créé à l'aide de CircuitLab

Appliquer KCL:

VR1+VR2+VR3+VR4=V1R1+V2R2+V3R3+V4R4V=(V1R1+V2R2+V3R3+V4R4)(R1∣∣R2∣∣R3∣∣R4)

Certains jouer avec l'algèbre obtient le résultat que j'ai utilisé dans le code.

jonk
la source
Merci pour la bonne réponse. J'ai quelques questions de clarification. 1) Pourquoi faut-il séparer le système enje+j = pair et je+j= impair? 2) Tout peut-il être résolu en une seule étape en utilisant une sorte de système linéaireUNEX=bUNE est un opérateur linéaire et bfournit les conditions aux limites? En y regardant, cela semble quelque peu analogue aux méthodes de différences finies pour résoudre des équations aux dérivées partielles ...
Darcy
De plus, la formule est-elle identique pour les étapes paires et impaires?
Darcy
2
@Darcy J'écrirai un peu plus pour aider à couvrir ces problèmes.
jonk
Merci encore pour les détails. Une dernière question (et cela pourrait peut-être devenir une question entièrement distincte, mais je vais la poser ici): si toutes les résistances du réseau ont la même résistance (disons 1 Ohm), alors s'ensuit-il que la résistance effective devrait également être 1 Ohm? Mon intuition dit que oui, mais je n'en suis pas sûr.
Darcy
1
@Darcy Votre intuition est fausse et le résultat MATLAB est correct.
jonk
1

Vous pouvez certainement adopter l'approche d'un réseau de résistances 2D pour modéliser un problème 2D mais cela peut devenir quelque peu délicat lors du passage à 3 dimensions. Vous voudrez peut-être envisager d'utiliser une approche plus traditionnelle (de nos jours) avec des conducteurs de volume définis dans vos domaines avec une conductivité appropriée attribuée à chacun. Le code freeware FEMM ( http://www.femm.info/wiki/HomePage ) est très performant et peut être utilisé pour la 2D, la symétrie axiale et la 3D. De là, vous pouvez envisager de passer à des codes beaucoup plus performants comme SCIrun ( https://www.sci.utah.edu/), qui est un code académique pour les problèmes de conducteurs de volume d'une complexité substantielle. Je l'utilise régulièrement pour des mailles de plus d'un million de tétraèdres. Même s'il a été principalement développé pour la modélisation biologique, il devrait fonctionner parfaitement pour ce que vous faites. Les exemples de problèmes avancés dans la boîte à outils avant / inverse devraient vous aider. Les problèmes inverses peuvent également être utiles pour la tomographie par impédance. J'utilise généralement la version 4 car la version 5 est toujours en cours de réalisation. Le logiciel possède également une interface avec tetgen qui est un excellent code de construction de maillage.

Enfin, si vous n'êtes pas opposé à dépenser de l'argent, il y a toujours COMSOL, qui est très facile à utiliser (et assez cher).

wilk
la source